| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * Copyright (c) 2000-2022 Inria | ||
| 3 | * All rights reserved. | ||
| 4 | * | ||
| 5 | * Redistribution and use in source and binary forms, with or without | ||
| 6 | * modification, are permitted provided that the following conditions are met: | ||
| 7 | * | ||
| 8 | * * Redistributions of source code must retain the above copyright notice, | ||
| 9 | * this list of conditions and the following disclaimer. | ||
| 10 | * * Redistributions in binary form must reproduce the above copyright notice, | ||
| 11 | * this list of conditions and the following disclaimer in the documentation | ||
| 12 | * and/or other materials provided with the distribution. | ||
| 13 | * * Neither the name of the ALICE Project-Team nor the names of its | ||
| 14 | * contributors may be used to endorse or promote products derived from this | ||
| 15 | * software without specific prior written permission. | ||
| 16 | * | ||
| 17 | * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" | ||
| 18 | * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | ||
| 19 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | ||
| 20 | * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE | ||
| 21 | * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR | ||
| 22 | * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF | ||
| 23 | * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS | ||
| 24 | * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN | ||
| 25 | * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) | ||
| 26 | * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE | ||
| 27 | * POSSIBILITY OF SUCH DAMAGE. | ||
| 28 | * | ||
| 29 | * Contact: Bruno Levy | ||
| 30 | * | ||
| 31 | * https://www.inria.fr/fr/bruno-levy | ||
| 32 | * | ||
| 33 | * Inria, | ||
| 34 | * Domaine de Voluceau, | ||
| 35 | * 78150 Le Chesnay - Rocquencourt | ||
| 36 | * FRANCE | ||
| 37 | * | ||
| 38 | */ | ||
| 39 | |||
| 40 | #include <geogram/delaunay/LFS.h> | ||
| 41 | #include <geogram/basic/geometry.h> | ||
| 42 | #include <geogram/basic/logger.h> | ||
| 43 | #include <geogram/basic/command_line.h> | ||
| 44 | |||
| 45 | namespace { | ||
| 46 | |||
| 47 | using namespace GEO; | ||
| 48 | |||
| 49 | const index_t pid_size = 6; | ||
| 50 | int pid[pid_size][4] = { | ||
| 51 | {0, 1, 2, 3}, | ||
| 52 | {0, 2, 1, 3}, | ||
| 53 | {0, 3, 1, 2}, | ||
| 54 | {1, 2, 0, 3}, | ||
| 55 | {1, 3, 2, 0}, | ||
| 56 | {2, 3, 0, 1} | ||
| 57 | }; | ||
| 58 | |||
| 59 | /** | ||
| 60 | * \brief Computes the circumcenter and squared radius of a tetrahedron. | ||
| 61 | * \param[in] p first vertex of the tetrahedron | ||
| 62 | * \param[in] q second vertex of the tetrahedron | ||
| 63 | * \param[in] r third vertex of the tetrahedron | ||
| 64 | * \param[in] s fourth vertex of the tetrahedron | ||
| 65 | * \param[out] t_circumcenter computed circumcenter of the tetrahedron | ||
| 66 | * \param[out] squared_radius computed squared radius of the tetrahedron | ||
| 67 | * \param[out] dihedral_angles if non-nullptr, | ||
| 68 | * computed dihedral angles of the tetrahedron | ||
| 69 | * \return true if successful, false otherwise | ||
| 70 | * (for instance, if the tetrahedron is flat). | ||
| 71 | */ | ||
| 72 | ✗ | bool tetra_circumcenter_squaredradius( | |
| 73 | const vec3& p, const vec3& q, | ||
| 74 | const vec3& r, const vec3& s, | ||
| 75 | vec3& t_circumcenter, double& squared_radius, | ||
| 76 | double* dihedral_angles = nullptr | ||
| 77 | ) { | ||
| 78 | |||
| 79 | ✗ | vec3 qp = q - p; | |
| 80 | ✗ | vec3 rp = r - p; | |
| 81 | ✗ | vec3 sp = s - p; | |
| 82 | |||
| 83 | ✗ | double tet_vol_x_12 = 2.0 * dot(qp, cross(rp, sp)); | |
| 84 | |||
| 85 | // it is not a safe check, one should avoid | ||
| 86 | // to input a degenerated tetra. | ||
| 87 | ✗ | if(tet_vol_x_12 == 0.0) { | |
| 88 | ✗ | return false; | |
| 89 | } | ||
| 90 | |||
| 91 | vec3 d_t_circumcenter = | ||
| 92 | ✗ | length2(qp) * cross(rp, sp) + | |
| 93 | ✗ | length2(rp) * cross(sp, qp) + | |
| 94 | ✗ | length2(sp) * cross(qp, rp); | |
| 95 | |||
| 96 | ✗ | squared_radius = length2( | |
| 97 | ✗ | d_t_circumcenter) / (tet_vol_x_12 * tet_vol_x_12 | |
| 98 | ); | ||
| 99 | ✗ | d_t_circumcenter = (1.0 / tet_vol_x_12) * d_t_circumcenter; | |
| 100 | ✗ | t_circumcenter = p + d_t_circumcenter; | |
| 101 | |||
| 102 | ✗ | if(dihedral_angles) { | |
| 103 | // compute Dihedral angle directly | ||
| 104 | // the following computation is not optimized. | ||
| 105 | ✗ | vec3 point[4]; | |
| 106 | ✗ | point[0] = p; | |
| 107 | ✗ | point[1] = q; | |
| 108 | ✗ | point[2] = r; | |
| 109 | ✗ | point[3] = s; | |
| 110 | ✗ | for(index_t i = 0; i < pid_size; i++) { | |
| 111 | ✗ | vec3 b1 = point[pid[i][0]] - point[pid[i][2]]; | |
| 112 | ✗ | vec3 b2 = point[pid[i][1]] - point[pid[i][0]]; | |
| 113 | ✗ | vec3 b3 = point[pid[i][3]] - point[pid[i][1]]; | |
| 114 | ✗ | vec3 b2b3 = cross(b2, b3); | |
| 115 | ✗ | dihedral_angles[i] = ::fabs( | |
| 116 | ✗ | ::atan2( | |
| 117 | ✗ | length(b2) * dot(b1, b2b3), | |
| 118 | ✗ | dot(cross(b1, b2), b2b3) | |
| 119 | ) | ||
| 120 | ); | ||
| 121 | } | ||
| 122 | } | ||
| 123 | ✗ | return true; | |
| 124 | } | ||
| 125 | |||
| 126 | /** | ||
| 127 | * \brief Gets a Delaunay vertex by global index. | ||
| 128 | * \param[in] delaunay the Delaunay triangulation | ||
| 129 | * \param[in] v the index of the vertex | ||
| 130 | * \return a const reference to the vertex, as a vec3 | ||
| 131 | */ | ||
| 132 | ✗ | inline const vec3& delaunay_vertex(Delaunay* delaunay, index_t v) { | |
| 133 | ✗ | return *(const vec3*) delaunay->vertex_ptr(v); | |
| 134 | } | ||
| 135 | |||
| 136 | /** | ||
| 137 | * \brief Gets a Delaunay vertex by tetrahedron index | ||
| 138 | * and local vertex index. | ||
| 139 | * \param[in] delaunay the Delaunay triangulation | ||
| 140 | * \param[in] c the index of the tetrahedron | ||
| 141 | * \param[in] lv the local index of the vertex (0,1,2 or 3) | ||
| 142 | * in tetrahedron \p c. | ||
| 143 | * \return a const reference to the vertex, as a vec3. | ||
| 144 | */ | ||
| 145 | ✗ | inline const vec3& delaunay_tet_vertex( | |
| 146 | Delaunay* delaunay, index_t c, index_t lv | ||
| 147 | ) { | ||
| 148 | ✗ | return delaunay_vertex(delaunay, index_t(delaunay->cell_vertex(c, lv))); | |
| 149 | } | ||
| 150 | |||
| 151 | index_t tet_facet_vertex[4][3] = { | ||
| 152 | {1, 2, 3}, | ||
| 153 | {0, 3, 2}, | ||
| 154 | {3, 0, 1}, | ||
| 155 | {2, 1, 0} | ||
| 156 | }; | ||
| 157 | |||
| 158 | /** | ||
| 159 | * \brief Computes the normal to a tetrahedron facet. | ||
| 160 | * \param[in] delaunay the Delaunay triangulation | ||
| 161 | * \param[in] t the index of the tetrahedron | ||
| 162 | * \param[in] f the local index (0,1,2 or 3) of | ||
| 163 | * the facet in the tetrahedron \p t | ||
| 164 | * \return the normal to the facet \p f or tetrahedron \p t | ||
| 165 | */ | ||
| 166 | ✗ | vec3 delaunay_facet_normal( | |
| 167 | Delaunay* delaunay, index_t t, index_t f | ||
| 168 | ) { | ||
| 169 | ✗ | geo_debug_assert(f < 4); | |
| 170 | ✗ | index_t v1 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][0])); | |
| 171 | ✗ | index_t v2 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][1])); | |
| 172 | ✗ | index_t v3 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][2])); | |
| 173 | ✗ | const vec3& p1 = delaunay_vertex(delaunay, v1); | |
| 174 | ✗ | const vec3& p2 = delaunay_vertex(delaunay, v2); | |
| 175 | ✗ | const vec3& p3 = delaunay_vertex(delaunay, v3); | |
| 176 | ✗ | return cross(p2 - p1, p3 - p1); | |
| 177 | } | ||
| 178 | } | ||
| 179 | |||
| 180 | /****************************************************************************/ | ||
| 181 | |||
| 182 | namespace GEO { | ||
| 183 | |||
| 184 | ✗ | void LocalFeatureSize::init(index_t nb_pts, const double* pts) { | |
| 185 | // Note: I need PDEL here instead of ANN/BNN since I need | ||
| 186 | // to access the cells (and ANN/BNN only give me the | ||
| 187 | // neighbors) | ||
| 188 | |||
| 189 | ✗ | Logger::out("LFS") << "Delaunay" << std::endl; | |
| 190 | ✗ | Delaunay_var delaunay = Delaunay::create(3, "PDEL"); | |
| 191 | ✗ | delaunay->set_vertices(nb_pts, pts); | |
| 192 | ✗ | Logger::out("LFS") << "Done Delaunay" << std::endl; | |
| 193 | |||
| 194 | ✗ | vector<vec3> circumcenter(delaunay->nb_cells()); | |
| 195 | vector<bool> voronoi_cell_is_infinite( | ||
| 196 | delaunay->nb_vertices(), false | ||
| 197 | ✗ | ); | |
| 198 | ✗ | vector<bool> is_sliver(delaunay->nb_cells()); | |
| 199 | |||
| 200 | static constexpr index_t NO_POLE = NO_INDEX; | ||
| 201 | // Index of incident tet whose circumcenter is the pole | ||
| 202 | ✗ | vector<index_t> positive_pole(delaunay->nb_vertices(),NO_POLE); | |
| 203 | ✗ | vector<index_t> negative_pole(delaunay->nb_vertices(),NO_POLE); | |
| 204 | |||
| 205 | vector<vec3> avg_infinite_dir( | ||
| 206 | ✗ | delaunay->nb_vertices(), vec3(0, 0, 0) | |
| 207 | ✗ | ); | |
| 208 | ✗ | vector<double> dist(delaunay->nb_vertices(), 0.0); | |
| 209 | |||
| 210 | ✗ | Logger::out("LFS") << "(1) Circumcenters and slivers" << std::endl; | |
| 211 | // Step 1: compute circumcenters and check for slivers | ||
| 212 | ✗ | const double sliver_quality = sliver_angle_threshold_ / 180.0 * M_PI; | |
| 213 | ✗ | for(index_t t = 0; t < delaunay->nb_cells(); t++) { | |
| 214 | ✗ | const vec3& p = delaunay_tet_vertex(delaunay, t, 0); | |
| 215 | ✗ | const vec3& q = delaunay_tet_vertex(delaunay, t, 1); | |
| 216 | ✗ | const vec3& r = delaunay_tet_vertex(delaunay, t, 2); | |
| 217 | ✗ | const vec3& s = delaunay_tet_vertex(delaunay, t, 3); | |
| 218 | |||
| 219 | double dihedral_angle[6]; | ||
| 220 | double tet_radius; | ||
| 221 | ✗ | bool ok = tetra_circumcenter_squaredradius( | |
| 222 | p, q, r, s, | ||
| 223 | ✗ | circumcenter[t], tet_radius, dihedral_angle | |
| 224 | ); | ||
| 225 | |||
| 226 | ✗ | for(index_t a = 0; a < 6; a++) { | |
| 227 | ✗ | if(!ok) { | |
| 228 | ✗ | break; | |
| 229 | } | ||
| 230 | ✗ | ok = ok && | |
| 231 | ✗ | dihedral_angle[a] >= sliver_quality && | |
| 232 | ✗ | dihedral_angle[a] <= M_PI - sliver_quality; | |
| 233 | } | ||
| 234 | ✗ | is_sliver[t] = !ok; | |
| 235 | } | ||
| 236 | |||
| 237 | ✗ | Logger::out("LFS") << "(2) Positive poles" << std::endl; | |
| 238 | // Step 2: compute positive poles | ||
| 239 | ✗ | for(index_t t = 0; t < delaunay->nb_cells(); t++) { | |
| 240 | ✗ | if(is_sliver[t]) { | |
| 241 | ✗ | continue; | |
| 242 | } | ||
| 243 | ✗ | index_t f_inf = NO_INDEX; | |
| 244 | ✗ | for(index_t f = 0; f < 4; f++) { | |
| 245 | ✗ | if(delaunay->cell_adjacent(t, f) == NO_INDEX) { | |
| 246 | ✗ | f_inf = index_t(f); | |
| 247 | ✗ | break; | |
| 248 | } | ||
| 249 | } | ||
| 250 | ✗ | if(f_inf == NO_INDEX) { | |
| 251 | // tet t does not have facet on border | ||
| 252 | ✗ | const vec3& c = circumcenter[t]; | |
| 253 | ✗ | for(index_t lv = 0; lv < 4; lv++) { | |
| 254 | ✗ | index_t iv = index_t(delaunay->cell_vertex(t, lv)); | |
| 255 | ✗ | const vec3& v = delaunay_vertex(delaunay, iv); | |
| 256 | ✗ | double d = length2(c - v); | |
| 257 | ✗ | if(d > dist[iv]) { | |
| 258 | ✗ | dist[iv] = d; | |
| 259 | ✗ | positive_pole[iv] = t; | |
| 260 | } | ||
| 261 | } | ||
| 262 | } else { | ||
| 263 | // tet t has facet f_inf on border | ||
| 264 | // positive pole = average direction of infinite Voronoi edges | ||
| 265 | ✗ | index_t v = index_t(delaunay->cell_vertex(t, index_t(f_inf))); | |
| 266 | ✗ | voronoi_cell_is_infinite[v] = true; | |
| 267 | ✗ | avg_infinite_dir[v] += | |
| 268 | ✗ | delaunay_facet_normal(delaunay, t, index_t(f_inf)); | |
| 269 | } | ||
| 270 | } | ||
| 271 | |||
| 272 | ✗ | Logger::out("LFS") << "(3) Negative poles" << std::endl; | |
| 273 | // Step 3: compute negative poles | ||
| 274 | ✗ | std::fill(dist.begin(), dist.end(), 0.0); | |
| 275 | ✗ | for(index_t t = 0; t < delaunay->nb_cells(); t++) { | |
| 276 | ✗ | if(is_sliver[t]) { | |
| 277 | ✗ | continue; | |
| 278 | } | ||
| 279 | ✗ | bool t_is_infinite = false; | |
| 280 | ✗ | for(index_t f = 0; f < 4; f++) { | |
| 281 | ✗ | if(delaunay->cell_adjacent(t, f) == NO_INDEX) { | |
| 282 | ✗ | t_is_infinite = true; | |
| 283 | ✗ | break; | |
| 284 | } | ||
| 285 | } | ||
| 286 | ✗ | if(t_is_infinite) { | |
| 287 | ✗ | continue; | |
| 288 | } | ||
| 289 | ✗ | const vec3& c = circumcenter[t]; | |
| 290 | ✗ | for(index_t lv = 0; lv < 4; lv++) { | |
| 291 | ✗ | index_t iv = index_t(delaunay->cell_vertex(t, lv)); | |
| 292 | ✗ | const vec3& v = delaunay_vertex(delaunay, iv); | |
| 293 | ✗ | vec3 N; | |
| 294 | ✗ | if(voronoi_cell_is_infinite[iv]) { | |
| 295 | ✗ | N = -avg_infinite_dir[iv]; | |
| 296 | } else { | ||
| 297 | ✗ | N = v - circumcenter[positive_pole[iv]]; | |
| 298 | } | ||
| 299 | ✗ | double d = dot(c - v, N); | |
| 300 | ✗ | if(d > dist[iv]) { | |
| 301 | ✗ | dist[iv] = d; | |
| 302 | ✗ | negative_pole[iv] = t; | |
| 303 | } | ||
| 304 | } | ||
| 305 | } | ||
| 306 | |||
| 307 | ✗ | Logger::out("LFS") << "(4) Kd-tree" << std::endl; | |
| 308 | // Step 4: create search structure | ||
| 309 | ✗ | vector<bool> is_pole(delaunay->nb_cells(), false); | |
| 310 | ✗ | for(index_t iv = 0; iv < delaunay->nb_vertices(); iv++) { | |
| 311 | ✗ | if(negative_pole[iv] != NO_POLE) { | |
| 312 | ✗ | is_pole[index_t(negative_pole[iv])] = true; | |
| 313 | } | ||
| 314 | ✗ | if(!voronoi_cell_is_infinite[iv]) { | |
| 315 | ✗ | if(positive_pole[iv] != NO_POLE) { | |
| 316 | ✗ | is_pole[index_t(positive_pole[iv])] = true; | |
| 317 | } | ||
| 318 | } | ||
| 319 | } | ||
| 320 | |||
| 321 | ✗ | index_t nb_poles = 0; | |
| 322 | ✗ | for(index_t t = 0; t < delaunay->nb_cells(); t++) { | |
| 323 | ✗ | if(is_pole[t]) { | |
| 324 | ✗ | nb_poles++; | |
| 325 | } | ||
| 326 | } | ||
| 327 | |||
| 328 | ✗ | poles_.reserve(nb_poles * 3); | |
| 329 | ✗ | for(index_t t = 0; t < delaunay->nb_cells(); t++) { | |
| 330 | ✗ | if(is_pole[t]) { | |
| 331 | ✗ | poles_.push_back(circumcenter[t].x); | |
| 332 | ✗ | poles_.push_back(circumcenter[t].y); | |
| 333 | ✗ | poles_.push_back(circumcenter[t].z); | |
| 334 | } | ||
| 335 | } | ||
| 336 | |||
| 337 | ✗ | spatial_search_ = Delaunay::create(3, "NN"); | |
| 338 | ✗ | spatial_search_->set_vertices(poles_.size() / 3, poles_.data()); | |
| 339 | |||
| 340 | ✗ | Logger::out("LFS") << "Done init." << std::endl; | |
| 341 | ✗ | } | |
| 342 | } | ||
| 343 |