| 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/voronoi/convex_cell.h> | ||
| 41 | |||
| 42 | #ifndef STANDALONE_CONVEX_CELL | ||
| 43 | #include <geogram/numerics/predicates.h> | ||
| 44 | #include <geogram/mesh/mesh.h> | ||
| 45 | #endif | ||
| 46 | |||
| 47 | #include <iostream> | ||
| 48 | #include <fstream> | ||
| 49 | #include <vector> | ||
| 50 | #include <cmath> | ||
| 51 | #include <limits> | ||
| 52 | #include <stack> | ||
| 53 | |||
| 54 | |||
| 55 | namespace { | ||
| 56 | using namespace VBW; | ||
| 57 | |||
| 58 | /** | ||
| 59 | * \brief a class for a stack of ushorts allocated on the stack. | ||
| 60 | * \details Used by clip_by_plane_fast() internally. I do not want | ||
| 61 | * clip_by_plane_fast() to do dynamic allocation. | ||
| 62 | */ | ||
| 63 | class SmallStack_ushort { | ||
| 64 | public: | ||
| 65 | |||
| 66 | /** | ||
| 67 | * \brief SmallStack constructor. | ||
| 68 | * \param[in] buffer a buffer. Callers keeps ownership | ||
| 69 | * (in most cases it will be created by alloca()). | ||
| 70 | * \param[in] capacity the number of uints that can | ||
| 71 | * be stored in the buffer. | ||
| 72 | */ | ||
| 73 | SmallStack_ushort( | ||
| 74 | ushort* buffer, int capacity | ||
| 75 | ) : buffer_(buffer), | ||
| 76 | index_(-1), | ||
| 77 | capacity_(capacity) { | ||
| 78 | } | ||
| 79 | |||
| 80 | /** | ||
| 81 | * \brief Tests whether this stack is empty. | ||
| 82 | * \retval true if this stack is empty. | ||
| 83 | * \retval false otherwise. | ||
| 84 | */ | ||
| 85 | bool empty() const { | ||
| 86 | return (index_ == -1); | ||
| 87 | } | ||
| 88 | |||
| 89 | /** | ||
| 90 | * \brief Pushes an item on the stack. | ||
| 91 | * \param[in] val the item to be pushed. | ||
| 92 | */ | ||
| 93 | void push(ushort val) { | ||
| 94 | ✗ | ++index_; | |
| 95 | vbw_assert(index_ < capacity_); | ||
| 96 | (void)capacity_; // To silence a warning. | ||
| 97 | ✗ | buffer_[index_] = val; | |
| 98 | ✗ | } | |
| 99 | |||
| 100 | /** | ||
| 101 | * \brief Pops an item from the stack. | ||
| 102 | */ | ||
| 103 | void pop() { | ||
| 104 | vbw_assert(!empty()); | ||
| 105 | ✗ | --index_; | |
| 106 | } | ||
| 107 | |||
| 108 | /** | ||
| 109 | * \brief Gets the item on the top of the stack. | ||
| 110 | * \return the item. | ||
| 111 | */ | ||
| 112 | ushort top() const { | ||
| 113 | vbw_assert(!empty()); | ||
| 114 | ✗ | return buffer_[index_]; | |
| 115 | } | ||
| 116 | private: | ||
| 117 | ushort* buffer_; | ||
| 118 | int index_; | ||
| 119 | int capacity_; | ||
| 120 | }; | ||
| 121 | } | ||
| 122 | |||
| 123 | /*************************************************************/ | ||
| 124 | |||
| 125 | namespace VBW { | ||
| 126 | |||
| 127 | 2 | ConvexCell::ConvexCell(ConvexCellFlags flags) : | |
| 128 | 2 | max_t_(64), | |
| 129 | 2 | max_v_(32), | |
| 130 | t_(max_t_), | ||
| 131 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | t_adj_(max_t_), |
| 132 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | plane_eqn_(max_v_), |
| 133 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | v2t_(max_v_), |
| 134 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | v2e_(max_v_) |
| 135 | { | ||
| 136 | #ifndef STANDALONE_CONVEX_CELL | ||
| 137 | 2 | use_exact_predicates_ = true; | |
| 138 | #endif | ||
| 139 | 2 | nb_t_ = 0; | |
| 140 | 2 | nb_v_ = 0; | |
| 141 | 2 | first_free_ = END_OF_LIST; | |
| 142 | 2 | first_valid_ = END_OF_LIST; | |
| 143 | 2 | geometry_dirty_ = true; | |
| 144 | 2 | has_vglobal_ = ((flags & WithVGlobal) != 0); | |
| 145 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(has_vglobal_) { |
| 146 | ✗ | vglobal_.assign(max_v_,index_t(-1)); | |
| 147 | } | ||
| 148 | 2 | has_tflags_ = ((flags & WithTFlags) != 0); | |
| 149 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(has_tflags_) { |
| 150 | ✗ | tflags_.assign(max_t_,0); | |
| 151 | } | ||
| 152 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | v2t_.assign(max_v_,ushort(-1)); |
| 153 |
1/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
2 | v2e_.assign(max_v_,uchar(-1)); |
| 154 | 2 | } | |
| 155 | |||
| 156 | /***********************************************************************/ | ||
| 157 | |||
| 158 | 442 | void ConvexCell::clear() { | |
| 159 | 442 | nb_t_ = 0; | |
| 160 | 442 | nb_v_ = 0; | |
| 161 | 442 | first_free_ = END_OF_LIST; | |
| 162 | 442 | first_valid_ = END_OF_LIST; | |
| 163 | 442 | geometry_dirty_ = true; | |
| 164 | #ifdef VBW_DEBUG | ||
| 165 | // Initialize all triangle flags with something | ||
| 166 | // different from VALID_TRIANGLE. | ||
| 167 | for(index_t t=0; t<max_t(); ++t) { | ||
| 168 | set_triangle_flags(t, END_OF_LIST); | ||
| 169 | } | ||
| 170 | #endif | ||
| 171 | 442 | } | |
| 172 | |||
| 173 | /***********************************************************************/ | ||
| 174 | |||
| 175 | ✗ | void ConvexCell::init_with_box( | |
| 176 | double xmin, double ymin, double zmin, | ||
| 177 | double xmax, double ymax, double zmax | ||
| 178 | ) { | ||
| 179 | ✗ | clear(); | |
| 180 | |||
| 181 | // The vertex at infinity. | ||
| 182 | ✗ | plane_eqn_[0] = make_vec4(0,0,0,0); | |
| 183 | |||
| 184 | // Offset for the 6 bounding box plane equations. | ||
| 185 | // Here they come first (offset is zero). | ||
| 186 | index_t boff = 1; | ||
| 187 | |||
| 188 | // The equations of the six faces of the bounding box. | ||
| 189 | ✗ | plane_eqn_[boff ] = make_vec4( 1.0, 0.0, 0.0, -xmin); | |
| 190 | ✗ | plane_eqn_[boff+1] = make_vec4(-1.0, 0.0, 0.0, xmax); | |
| 191 | ✗ | plane_eqn_[boff+2] = make_vec4( 0.0, 1.0, 0.0, -ymin); | |
| 192 | ✗ | plane_eqn_[boff+3] = make_vec4( 0.0,-1.0, 0.0, ymax); | |
| 193 | ✗ | plane_eqn_[boff+4] = make_vec4( 0.0, 0.0, 1.0, -zmin); | |
| 194 | ✗ | plane_eqn_[boff+5] = make_vec4( 0.0, 0.0,-1.0, zmax); | |
| 195 | |||
| 196 | // Create the 8 triangles that correspond to the | ||
| 197 | // 8 vertices of the bounding box. | ||
| 198 | // (Unused) adjacency info. ----------------. | ||
| 199 | // Triangle vertices -. | | ||
| 200 | // v v | ||
| 201 | new_triangle( boff+2,boff+5,boff+0, 1,4,2); | ||
| 202 | new_triangle( boff+5,boff+3,boff+0, 5,0,3); | ||
| 203 | new_triangle( boff+1,boff+5,boff+2, 0,6,3); | ||
| 204 | new_triangle( boff+5,boff+1,boff+3, 7,1,2); | ||
| 205 | new_triangle( boff+4,boff+2,boff+0, 0,5,6); | ||
| 206 | new_triangle( boff+4,boff+0,boff+3, 1,7,4); | ||
| 207 | new_triangle( boff+2,boff+4,boff+1, 7,2,4); | ||
| 208 | new_triangle( boff+4,boff+3,boff+1, 3,6,5); | ||
| 209 | |||
| 210 | // We already created 6 vertices (for the 6 bounding box | ||
| 211 | // plane equations) plus the vertex at infinity. | ||
| 212 | ✗ | nb_v_ = 7; | |
| 213 | |||
| 214 | ✗ | geometry_dirty_ = true; | |
| 215 | ✗ | } | |
| 216 | |||
| 217 | |||
| 218 | ✗ | void ConvexCell::init_with_tet( | |
| 219 | vec4 P0, vec4 P1, vec4 P2, vec4 P3 | ||
| 220 | ) { | ||
| 221 | ✗ | clear(); | |
| 222 | |||
| 223 | // The vertex at infinity. | ||
| 224 | ✗ | plane_eqn_[0] = make_vec4(0,0,0,0); | |
| 225 | |||
| 226 | // Offset for the 4 plane equations. | ||
| 227 | // Plane 0 is vertex at infinity. | ||
| 228 | index_t boff = 1; | ||
| 229 | |||
| 230 | ✗ | plane_eqn_[boff ] = P0; | |
| 231 | ✗ | plane_eqn_[boff+1] = P1; | |
| 232 | ✗ | plane_eqn_[boff+2] = P2; | |
| 233 | ✗ | plane_eqn_[boff+3] = P3; | |
| 234 | |||
| 235 | // Create the 4 triangles (that correspond to | ||
| 236 | // the 4 vertices of the tetrahedron) | ||
| 237 | // (Unused) adjacency info. ----------. | ||
| 238 | // Triangle vertices -. | | ||
| 239 | // v v | ||
| 240 | new_triangle(boff+3, boff+2, boff+1, 3, 2, 1); | ||
| 241 | new_triangle(boff+3, boff+0, boff+2, 3, 0, 2); | ||
| 242 | new_triangle(boff+3, boff+1, boff+0, 3, 1, 0); | ||
| 243 | new_triangle(boff+2, boff+0, boff+1, 2, 0, 1); | ||
| 244 | |||
| 245 | // We already created 4 vertices (for the 4 facets | ||
| 246 | // plane equations) plus the vertex at infinity. | ||
| 247 | ✗ | nb_v_ = 5; | |
| 248 | |||
| 249 | ✗ | geometry_dirty_ = true; | |
| 250 | ✗ | } | |
| 251 | |||
| 252 | |||
| 253 | ✗ | void ConvexCell::init_with_tet( | |
| 254 | vec4 P0, vec4 P1, vec4 P2, vec4 P3, | ||
| 255 | global_index_t P0_global_index, | ||
| 256 | global_index_t P1_global_index, | ||
| 257 | global_index_t P2_global_index, | ||
| 258 | global_index_t P3_global_index | ||
| 259 | ) { | ||
| 260 | geo_debug_assert(has_vglobal_); | ||
| 261 | ✗ | init_with_tet(P0, P1, P2, P3); | |
| 262 | |||
| 263 | // Offset for the 4 plane equations. | ||
| 264 | // Plane 0 is vertex at infinity. | ||
| 265 | index_t boff = 1; | ||
| 266 | |||
| 267 | ✗ | vglobal_[boff ] = P0_global_index; | |
| 268 | ✗ | vglobal_[boff+1] = P1_global_index; | |
| 269 | ✗ | vglobal_[boff+2] = P2_global_index; | |
| 270 | ✗ | vglobal_[boff+3] = P3_global_index; | |
| 271 | ✗ | } | |
| 272 | |||
| 273 | |||
| 274 | /***********************************************************************/ | ||
| 275 | |||
| 276 | ✗ | void ConvexCell::save(const std::string& filename, double shrink) const { | |
| 277 | ✗ | std::ofstream out(filename.c_str()); | |
| 278 | ✗ | save(out, 1, shrink); | |
| 279 | ✗ | } | |
| 280 | |||
| 281 | |||
| 282 | 442 | index_t ConvexCell::save( | |
| 283 | std::ostream& out, global_index_t v_offset, | ||
| 284 | double shrink, bool borders_only | ||
| 285 | ) const { | ||
| 286 | |||
| 287 | vec3 g = make_vec3(0.0, 0.0, 0.0); | ||
| 288 |
1/2✓ Branch 0 taken 442 times.
✗ Branch 1 not taken.
|
442 | if(shrink != 0.0) { |
| 289 | 442 | const_cast<ConvexCell*>(this)->compute_geometry(); | |
| 290 | 442 | g = barycenter(); | |
| 291 | } | ||
| 292 | |||
| 293 | 442 | vector<index_t> v2t(nb_v(),index_t(-1)); | |
| 294 |
1/4✓ Branch 1 taken 442 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
442 | vector<index_t> t_index(nb_t(),index_t(-1)); |
| 295 | index_t nt=0; | ||
| 296 | |||
| 297 | { | ||
| 298 | 442 | index_t t = first_valid_; | |
| 299 |
2/2✓ Branch 0 taken 9684 times.
✓ Branch 1 taken 442 times.
|
10126 | while(t != END_OF_LIST) { |
| 300 | TriangleWithFlags T = get_triangle_and_flags(t); | ||
| 301 | vec4 p; | ||
| 302 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9684 times.
|
9684 | if(geometry_dirty_) { |
| 303 | ✗ | p = compute_triangle_point(t); | |
| 304 | ✗ | p.x /= p.w; | |
| 305 | ✗ | p.y /= p.w; | |
| 306 | ✗ | p.z /= p.w; | |
| 307 | p.w = 1.0; | ||
| 308 | } else { | ||
| 309 | 9684 | p.x = triangle_point_[t].x; | |
| 310 | 9684 | p.y = triangle_point_[t].y; | |
| 311 | 9684 | p.z = triangle_point_[t].z; | |
| 312 | p.w = 1.0; | ||
| 313 | } | ||
| 314 | |||
| 315 |
1/2✓ Branch 0 taken 9684 times.
✗ Branch 1 not taken.
|
9684 | if(shrink != 0.0) { |
| 316 | 9684 | p.x = shrink * g.x + (1.0 - shrink) * p.x; | |
| 317 | 9684 | p.y = shrink * g.y + (1.0 - shrink) * p.y; | |
| 318 | 9684 | p.z = shrink * g.z + (1.0 - shrink) * p.z; | |
| 319 | } | ||
| 320 | out << "v " << p.x << " " << p.y << " " << p.z << std::endl; | ||
| 321 | 9684 | t_index[t] = nt; | |
| 322 | 9684 | ++nt; | |
| 323 | 9684 | v2t[T.i] = t; | |
| 324 | 9684 | v2t[T.j] = t; | |
| 325 | 9684 | v2t[T.k] = t; | |
| 326 | 9684 | t = index_t(T.flags); | |
| 327 | } | ||
| 328 | } | ||
| 329 | |||
| 330 |
2/2✓ Branch 0 taken 6912 times.
✓ Branch 1 taken 442 times.
|
7354 | for(index_t v=1; v<nb_v(); ++v) { |
| 331 | ✗ | if(borders_only && | |
| 332 | ✗ | has_vglobal() && | |
| 333 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 6912 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
6912 | v_global_index(v) != global_index_t(-1) && |
| 334 | v_global_index(v) != global_index_t(-2) | ||
| 335 | // index_t(-2) is for fluid free bndry | ||
| 336 | ) { | ||
| 337 | ✗ | continue; | |
| 338 | } | ||
| 339 |
2/2✓ Branch 0 taken 5726 times.
✓ Branch 1 taken 1186 times.
|
6912 | if(v2t[v] != index_t(-1)) { |
| 340 | index_t t = v2t[v]; | ||
| 341 |
1/2✓ Branch 1 taken 5726 times.
✗ Branch 2 not taken.
|
5726 | out << "f "; |
| 342 | do { | ||
| 343 |
2/4✓ Branch 1 taken 29052 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 29052 times.
✗ Branch 5 not taken.
|
29052 | out << (t_index[t]+v_offset) << " "; |
| 344 | index_t lv = triangle_find_vertex(t,v); | ||
| 345 |
2/2✓ Branch 0 taken 23326 times.
✓ Branch 1 taken 5726 times.
|
29052 | t = triangle_adjacent(t, (lv + 1)%3); |
| 346 |
2/2✓ Branch 0 taken 23326 times.
✓ Branch 1 taken 5726 times.
|
29052 | } while(t != v2t[v]); |
| 347 | out << std::endl; | ||
| 348 | } | ||
| 349 | } | ||
| 350 | |||
| 351 | 442 | return nt; | |
| 352 | } | ||
| 353 | |||
| 354 | ✗ | void ConvexCell::for_each_Voronoi_vertex( | |
| 355 | index_t v, | ||
| 356 | std::function<void(index_t)> vertex | ||
| 357 | ) { | ||
| 358 | geo_debug_assert(!geometry_dirty_); | ||
| 359 | ✗ | if(v2t_[v] != END_OF_LIST) { | |
| 360 | ✗ | index_t t = index_t(v2t_[v]); | |
| 361 | do { | ||
| 362 | ✗ | vertex(t); | |
| 363 | index_t lv = triangle_find_vertex(t,v); | ||
| 364 | ✗ | t = triangle_adjacent(t, (lv + 1)%3); | |
| 365 | ✗ | } while(t != v2t_[v]); | |
| 366 | } | ||
| 367 | ✗ | } | |
| 368 | |||
| 369 | #if !defined(STANDALONE_CONVEX_CELL) && !defined(GEOGRAM_PSM) | ||
| 370 | |||
| 371 | ✗ | void ConvexCell::append_to_mesh( | |
| 372 | GEO::Mesh* mesh, double shrink, bool borders_only, | ||
| 373 | GEO::Attribute<GEO::index_t>* facet_attr | ||
| 374 | ) const { | ||
| 375 | |||
| 376 | global_index_t v_offset = mesh->vertices.nb(); | ||
| 377 | |||
| 378 | vec3 g = make_vec3(0.0, 0.0, 0.0); | ||
| 379 | ✗ | if(shrink != 0.0) { | |
| 380 | ✗ | const_cast<ConvexCell*>(this)->compute_geometry(); | |
| 381 | ✗ | g = barycenter(); | |
| 382 | } | ||
| 383 | |||
| 384 | ✗ | vector<index_t> v2t(nb_v(),index_t(-1)); | |
| 385 | ✗ | vector<index_t> t_index(nb_t(),index_t(-1)); | |
| 386 | index_t nt=0; | ||
| 387 | |||
| 388 | { | ||
| 389 | ✗ | index_t t = first_valid_; | |
| 390 | ✗ | while(t != END_OF_LIST) { | |
| 391 | TriangleWithFlags T = get_triangle_and_flags(t); | ||
| 392 | ✗ | vec4 p = compute_triangle_point(t); | |
| 393 | ✗ | p.x /= p.w; | |
| 394 | ✗ | p.y /= p.w; | |
| 395 | ✗ | p.z /= p.w; | |
| 396 | ✗ | p.w = 1.0; | |
| 397 | ✗ | if(shrink != 0.0) { | |
| 398 | ✗ | p.x = shrink * g.x + (1.0 - shrink) * p.x; | |
| 399 | ✗ | p.y = shrink * g.y + (1.0 - shrink) * p.y; | |
| 400 | ✗ | p.z = shrink * g.z + (1.0 - shrink) * p.z; | |
| 401 | } | ||
| 402 | ✗ | mesh->vertices.create_vertex(p.data()); | |
| 403 | ✗ | t_index[t] = nt; | |
| 404 | ✗ | ++nt; | |
| 405 | ✗ | v2t[T.i] = t; | |
| 406 | ✗ | v2t[T.j] = t; | |
| 407 | ✗ | v2t[T.k] = t; | |
| 408 | ✗ | t = index_t(T.flags); | |
| 409 | } | ||
| 410 | } | ||
| 411 | |||
| 412 | ✗ | for(index_t v=1; v<nb_v(); ++v) { | |
| 413 | ✗ | if(borders_only && | |
| 414 | ✗ | has_vglobal() && | |
| 415 | ✗ | v_global_index(v) != global_index_t(-1) && | |
| 416 | v_global_index(v) != global_index_t(-2) // This one for | ||
| 417 | // fluid free bndry | ||
| 418 | ) { | ||
| 419 | ✗ | continue; | |
| 420 | } | ||
| 421 | std::vector<global_index_t> facet_vertices; | ||
| 422 | ✗ | if(v2t[v] != index_t(-1)) { | |
| 423 | index_t t = v2t[v]; | ||
| 424 | do { | ||
| 425 | ✗ | facet_vertices.push_back(t_index[t]+v_offset); | |
| 426 | index_t lv = triangle_find_vertex(t,v); | ||
| 427 | ✗ | t = triangle_adjacent(t, (lv + 1)%3); | |
| 428 | ✗ | } while(t != v2t[v]); | |
| 429 | } | ||
| 430 | ✗ | if(facet_vertices.size() < 3) { | |
| 431 | continue; | ||
| 432 | } | ||
| 433 | ✗ | global_index_t f = mesh->facets.create_polygon( | |
| 434 | GEO::index_t(facet_vertices.size()) | ||
| 435 | ); | ||
| 436 | ✗ | for(index_t i=0; i<facet_vertices.size(); ++i) { | |
| 437 | ✗ | mesh->facets.set_vertex(f, i, facet_vertices[i]); | |
| 438 | } | ||
| 439 | ✗ | if(facet_attr != nullptr && facet_attr->is_bound()) { | |
| 440 | ✗ | (*facet_attr)[f] = v_global_index(v); | |
| 441 | } | ||
| 442 | } | ||
| 443 | |||
| 444 | ✗ | } | |
| 445 | |||
| 446 | #endif | ||
| 447 | |||
| 448 | /***********************************************************************/ | ||
| 449 | |||
| 450 | ✗ | bool ConvexCell::has_v_global_index(global_index_t v) const { | |
| 451 | vbw_assert(has_vglobal_); | ||
| 452 | ✗ | for(index_t i=0; i<nb_v(); ++i) { | |
| 453 | ✗ | if(vglobal_[i] == v) { | |
| 454 | return true; | ||
| 455 | } | ||
| 456 | } | ||
| 457 | return false; | ||
| 458 | } | ||
| 459 | |||
| 460 | |||
| 461 | ✗ | void ConvexCell::clip_by_plane(vec4 eqn, global_index_t j) { | |
| 462 | vbw_assert(has_vglobal_); | ||
| 463 | ✗ | clip_by_plane(eqn); | |
| 464 | ✗ | vglobal_[nb_v()-1] = j; | |
| 465 | ✗ | } | |
| 466 | |||
| 467 | 1326 | void ConvexCell::clip_by_plane(vec4 eqn) { | |
| 468 | 1326 | geometry_dirty_ = true; | |
| 469 | |||
| 470 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1324 times.
|
1326 | index_t lv = nb_v_; |
| 471 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1324 times.
|
1326 | if(lv == max_v()) { |
| 472 | 2 | grow_v(); | |
| 473 | } | ||
| 474 | 1326 | plane_eqn_[lv] = eqn; | |
| 475 | vbw_assert(lv < max_v()); | ||
| 476 | 1326 | ++nb_v_; | |
| 477 | |||
| 478 | // Step 1: Find conflict zone and link conflicted triangles | ||
| 479 | // (to recycle them in free list). | ||
| 480 | |||
| 481 | index_t conflict_head = END_OF_LIST; | ||
| 482 | index_t conflict_tail = END_OF_LIST; | ||
| 483 | |||
| 484 | // Classify triangles, compute conflict list and valid list. | ||
| 485 | // Note: This could be done by climbing from a random triangle, | ||
| 486 | // but here we prefer complete linear scan for several reasons: | ||
| 487 | // - it is more robust to numerical errors (here we are not | ||
| 488 | // using exact predicates). | ||
| 489 | // - the code is simpler. | ||
| 490 | // - and more importantly, we got no more than a few tenths of | ||
| 491 | // vertices. | ||
| 492 | // The 'climbing from a random triangle' strategy is implemented | ||
| 493 | // in clip_by_plane_fast(). We keep both implementations for now, | ||
| 494 | // until we make sure than one is more efficient than the other one. | ||
| 495 | |||
| 496 | 1326 | index_t t = first_valid_; | |
| 497 | 1326 | first_valid_ = END_OF_LIST; | |
| 498 |
2/2✓ Branch 0 taken 23904 times.
✓ Branch 1 taken 1326 times.
|
25230 | while(t != END_OF_LIST) { |
| 499 | TriangleWithFlags T = get_triangle_and_flags(t); | ||
| 500 |
2/2✓ Branch 1 taken 1708 times.
✓ Branch 2 taken 22196 times.
|
23904 | if(triangle_is_in_conflict(T,eqn)) { |
| 501 | set_triangle_flags( | ||
| 502 |
2/2✓ Branch 0 taken 298 times.
✓ Branch 1 taken 1410 times.
|
1708 | t, ushort(conflict_head) | ushort(CONFLICT_MASK) |
| 503 | ); | ||
| 504 | conflict_head = t; | ||
| 505 |
2/2✓ Branch 0 taken 298 times.
✓ Branch 1 taken 1410 times.
|
1708 | if(conflict_tail == END_OF_LIST) { |
| 506 | conflict_tail = t; | ||
| 507 | } | ||
| 508 | } else { | ||
| 509 | 22196 | set_triangle_flags(t, ushort(first_valid_)); | |
| 510 | 22196 | first_valid_ = t; | |
| 511 | } | ||
| 512 | 23904 | t = index_t(T.flags); | |
| 513 | } | ||
| 514 | |||
| 515 | 1326 | triangulate_conflict_zone(lv, conflict_head, conflict_tail); | |
| 516 | 1326 | } | |
| 517 | |||
| 518 | // This version of clip_by_plane(), with a user-defined predicate, | ||
| 519 | // is duplicated from the standard version above. May be fixed by | ||
| 520 | // moving it to the header file (to have faster invokation of the | ||
| 521 | // predicate) and make the default version call it with PCK, but | ||
| 522 | // I do not want to do that because: | ||
| 523 | // - will make the header more heavy, longer compilation time | ||
| 524 | // - not sure about the impact on performance | ||
| 525 | // - the function is short, not a big drama to duplicate it... | ||
| 526 | ✗ | void ConvexCell::clip_by_plane( | |
| 527 | vec4 eqn, global_index_t global_index, | ||
| 528 | std::function<bool(ushort,ushort)> triangle_conflict_predicate | ||
| 529 | ) { | ||
| 530 | ✗ | geometry_dirty_ = true; | |
| 531 | |||
| 532 | ✗ | index_t lv = nb_v_; | |
| 533 | ✗ | if(lv == max_v()) { | |
| 534 | ✗ | grow_v(); | |
| 535 | } | ||
| 536 | ✗ | plane_eqn_[lv] = eqn; | |
| 537 | vbw_assert(lv < max_v()); | ||
| 538 | ✗ | ++nb_v_; | |
| 539 | |||
| 540 | // Note: it is unlikely that this function is used without | ||
| 541 | // global indices (because without global indices, it would | ||
| 542 | // mean we are only using geometry, then we should use the | ||
| 543 | // default predicate), so we could probably make it mandatory | ||
| 544 | // to have global indices here. | ||
| 545 | ✗ | if(has_vglobal_) { | |
| 546 | ✗ | vglobal_[nb_v()-1] = global_index; | |
| 547 | } | ||
| 548 | |||
| 549 | |||
| 550 | // Step 1: Find conflict zone and link conflicted triangles | ||
| 551 | // (to recycle them in free list). | ||
| 552 | |||
| 553 | index_t conflict_head = END_OF_LIST; | ||
| 554 | index_t conflict_tail = END_OF_LIST; | ||
| 555 | |||
| 556 | // Classify triangles, compute conflict list and valid list. | ||
| 557 | // Note: This could be done by climbing from a random triangle, | ||
| 558 | // but here we prefer complete linear scan for several reasons: | ||
| 559 | // - it is more robust to numerical errors (here we are not | ||
| 560 | // using exact predicates). | ||
| 561 | // - the code is simpler. | ||
| 562 | // - and more importantly, we got no more than a few tenths of | ||
| 563 | // vertices. | ||
| 564 | // The 'climbing from a random triangle' strategy is implemented | ||
| 565 | // in clip_by_plane_fast(). We keep both implementations for now, | ||
| 566 | // until we make sure than one is more efficient than the other one | ||
| 567 | // (and clip_by_plane_fast() does not have a version with the user- | ||
| 568 | // defined predicate for now). | ||
| 569 | |||
| 570 | ✗ | index_t t = first_valid_; | |
| 571 | ✗ | first_valid_ = END_OF_LIST; | |
| 572 | ✗ | while(t != END_OF_LIST) { | |
| 573 | TriangleWithFlags T = get_triangle_and_flags(t); | ||
| 574 | ✗ | if(triangle_conflict_predicate(ushort(t), ushort(nb_v()-1))) { | |
| 575 | set_triangle_flags( | ||
| 576 | ✗ | t, ushort(conflict_head) | ushort(CONFLICT_MASK) | |
| 577 | ); | ||
| 578 | conflict_head = t; | ||
| 579 | ✗ | if(conflict_tail == END_OF_LIST) { | |
| 580 | conflict_tail = t; | ||
| 581 | } | ||
| 582 | } else { | ||
| 583 | ✗ | set_triangle_flags(t, ushort(first_valid_)); | |
| 584 | ✗ | first_valid_ = t; | |
| 585 | } | ||
| 586 | ✗ | t = index_t(T.flags); | |
| 587 | } | ||
| 588 | |||
| 589 | ✗ | triangulate_conflict_zone(lv, conflict_head, conflict_tail); | |
| 590 | ✗ | } | |
| 591 | |||
| 592 | |||
| 593 | ✗ | void ConvexCell::clip_by_plane_fast(vec4 P, global_index_t j) { | |
| 594 | vbw_assert(has_vglobal_); | ||
| 595 | ✗ | clip_by_plane_fast(P); | |
| 596 | ✗ | vglobal_[nb_v()-1] = j; | |
| 597 | ✗ | } | |
| 598 | |||
| 599 | ✗ | void ConvexCell::clip_by_plane_fast(vec4 P) { | |
| 600 | ✗ | geometry_dirty_ = true; | |
| 601 | ✗ | index_t lv = nb_v_; | |
| 602 | ✗ | if(lv == max_v()) { | |
| 603 | ✗ | grow_v(); | |
| 604 | } | ||
| 605 | ✗ | plane_eqn_[lv] = P; | |
| 606 | vbw_assert(lv < max_v()); | ||
| 607 | ✗ | ++nb_v_; | |
| 608 | |||
| 609 | // Step 1: Find a good seed triangle (likely to | ||
| 610 | // be in conflict). If it is not in conflict, | ||
| 611 | // it is not a big problem (it will be just a | ||
| 612 | // bit slower), so we use inexact predicates | ||
| 613 | // here. | ||
| 614 | |||
| 615 | index_t t_init = END_OF_LIST; | ||
| 616 | |||
| 617 | { | ||
| 618 | index_t t_pred = END_OF_LIST; | ||
| 619 | ✗ | index_t t = first_valid_; | |
| 620 | ✗ | if(t == END_OF_LIST) { | |
| 621 | ✗ | return; | |
| 622 | } | ||
| 623 | |||
| 624 | ✗ | auto triangle_distance = [this](index_t t_in, vec4 P_in) { | |
| 625 | ✗ | Triangle T = get_triangle(t_in); | |
| 626 | vbw_assert(T.i != VERTEX_AT_INFINITY); | ||
| 627 | vbw_assert(T.j != VERTEX_AT_INFINITY); | ||
| 628 | vbw_assert(T.k != VERTEX_AT_INFINITY); | ||
| 629 | vec4 p1 = vertex_plane(T.i); | ||
| 630 | vec4 p2 = vertex_plane(T.j); | ||
| 631 | vec4 p3 = vertex_plane(T.k); | ||
| 632 | ✗ | return det4x4( | |
| 633 | p1.x, p2.x, p3.x, P_in.x, | ||
| 634 | p1.y, p2.y, p3.y, P_in.y, | ||
| 635 | p1.z, p2.z, p3.z, P_in.z, | ||
| 636 | p1.w, p2.w, p3.w, P_in.w | ||
| 637 | ✗ | ); | |
| 638 | ✗ | }; | |
| 639 | |||
| 640 | index_t count = 100; | ||
| 641 | ✗ | double t_dist = triangle_distance(t,P); | |
| 642 | |||
| 643 | ✗ | still_walking: | |
| 644 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 645 | index_t t_next = triangle_adjacent(t, le); | ||
| 646 | ✗ | if(t_next == t_pred) { | |
| 647 | ✗ | continue; | |
| 648 | } | ||
| 649 | ✗ | double t_next_dist = triangle_distance(t_next,P); | |
| 650 | ✗ | if(t_next_dist < t_dist) { | |
| 651 | ✗ | continue; | |
| 652 | } | ||
| 653 | ✗ | --count; | |
| 654 | t_pred = t; | ||
| 655 | t = t_next; | ||
| 656 | t_dist = t_next_dist; | ||
| 657 | ✗ | if(count > 0 && t_dist < 0.0) { | |
| 658 | ✗ | goto still_walking; | |
| 659 | } | ||
| 660 | } | ||
| 661 | t_init = t; | ||
| 662 | } | ||
| 663 | |||
| 664 | // note: t_init is now one triangle, probably in conflict | ||
| 665 | // (but not always: there can be degenerate cases where count | ||
| 666 | // reach zero). When t_init is in conflict, it is in general | ||
| 667 | // not the "highest" triangle (as one expect in a Delaunay | ||
| 668 | // triangulation code, but here it is different). | ||
| 669 | |||
| 670 | |||
| 671 | // Step 2: mark all the triangles that have the same | ||
| 672 | // conflict status as t_init with CONFLICT_MASK | ||
| 673 | // (even if t_init was not in conflict, then we will | ||
| 674 | // swap confict mask) | ||
| 675 | |||
| 676 | ✗ | bool t_init_is_in_conflict = triangle_is_in_conflict( | |
| 677 | get_triangle_and_flags(t_init), P | ||
| 678 | ); | ||
| 679 | |||
| 680 | { | ||
| 681 | ✗ | ushort* buff = (ushort*)alloca(max_t_*sizeof(ushort)); | |
| 682 | SmallStack_ushort S(buff, int(max_t_)); | ||
| 683 | ✗ | S.push(ushort(t_init)); | |
| 684 | set_triangle_flags( | ||
| 685 | t_init, | ||
| 686 | get_triangle_flags(t_init) | | ||
| 687 | ✗ | ushort(CONFLICT_MASK) | | |
| 688 | ushort(MARKED_MASK) | ||
| 689 | ); | ||
| 690 | ✗ | while(!S.empty()) { | |
| 691 | index_t t = index_t(S.top()); | ||
| 692 | S.pop(); | ||
| 693 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 694 | index_t t_neigh = triangle_adjacent(t,le); | ||
| 695 | ✗ | if( | |
| 696 | (get_triangle_flags(t_neigh) & ushort(MARKED_MASK)) | ||
| 697 | == 0 | ||
| 698 | ) { | ||
| 699 | ✗ | if( triangle_is_in_conflict( | |
| 700 | get_triangle_and_flags(t_neigh), P | ||
| 701 | ) == t_init_is_in_conflict | ||
| 702 | ) { | ||
| 703 | set_triangle_flags( | ||
| 704 | t_neigh, | ||
| 705 | get_triangle_flags(t_neigh) | | ||
| 706 | ✗ | ushort(CONFLICT_MASK) | | |
| 707 | ushort(MARKED_MASK) | ||
| 708 | ); | ||
| 709 | S.push(ushort(t_neigh)); | ||
| 710 | } else { | ||
| 711 | set_triangle_flags( | ||
| 712 | t_neigh, | ||
| 713 | ✗ | get_triangle_flags(t_neigh) | | |
| 714 | ushort(MARKED_MASK) | ||
| 715 | ); | ||
| 716 | } | ||
| 717 | } | ||
| 718 | } | ||
| 719 | } | ||
| 720 | } | ||
| 721 | |||
| 722 | // Step 3: update conflict list and active triangles list | ||
| 723 | index_t conflict_head = END_OF_LIST; | ||
| 724 | index_t conflict_tail = END_OF_LIST; | ||
| 725 | { | ||
| 726 | index_t new_first_valid = END_OF_LIST; | ||
| 727 | ✗ | index_t t = first_valid_; | |
| 728 | ✗ | while(t != END_OF_LIST) { | |
| 729 | bool t_is_in_conflict = triangle_is_marked_as_conflict(t); | ||
| 730 | index_t t_next = | ||
| 731 | ✗ | index_t( | |
| 732 | get_triangle_flags(t) & | ||
| 733 | ~(CONFLICT_MASK | MARKED_MASK) | ||
| 734 | ✗ | ); | |
| 735 | // Flip conflict flag if t_init was not in conflict. | ||
| 736 | ✗ | if(!t_init_is_in_conflict) { | |
| 737 | ✗ | t_is_in_conflict = !t_is_in_conflict; | |
| 738 | } | ||
| 739 | ✗ | if(t_is_in_conflict) { | |
| 740 | set_triangle_flags( | ||
| 741 | ✗ | t, ushort(conflict_head) | ushort(CONFLICT_MASK) | |
| 742 | ); | ||
| 743 | conflict_head = t; | ||
| 744 | ✗ | if(conflict_tail == END_OF_LIST) { | |
| 745 | conflict_tail = t; | ||
| 746 | } | ||
| 747 | } else { | ||
| 748 | ✗ | set_triangle_flags(t, ushort(new_first_valid)); | |
| 749 | new_first_valid = t; | ||
| 750 | } | ||
| 751 | t = t_next; | ||
| 752 | } | ||
| 753 | ✗ | first_valid_ = new_first_valid; | |
| 754 | } | ||
| 755 | |||
| 756 | ✗ | triangulate_conflict_zone(lv, conflict_head, conflict_tail); | |
| 757 | } | ||
| 758 | |||
| 759 | /***********************************************************************/ | ||
| 760 | |||
| 761 | 1326 | void ConvexCell::triangulate_conflict_zone( | |
| 762 | index_t lv, index_t conflict_head, index_t conflict_tail | ||
| 763 | ) { | ||
| 764 | // Special case: no triangle in conflict. | ||
| 765 |
2/2✓ Branch 0 taken 298 times.
✓ Branch 1 taken 1028 times.
|
1326 | if(conflict_head == END_OF_LIST) { |
| 766 | return; | ||
| 767 | } | ||
| 768 | |||
| 769 | |||
| 770 | // Link the vertices on the border of the conflict zone. | ||
| 771 | // Consider the edge e of triangle t such that: | ||
| 772 | // t is not marked as conflict | ||
| 773 | // triangle_adjacent(t,e) is marked as conflict | ||
| 774 | // Let v1 = triangle_vertex(t, (e+1)%3) | ||
| 775 | // v2 = triangle_vertex(t, (e+2)%3) | ||
| 776 | // We set: | ||
| 777 | // v2t_[v1] = t | ||
| 778 | // v2e_[v1] = e | ||
| 779 | // Once done for all the triangles, one can traverse the | ||
| 780 | // border of the conflict zone with: | ||
| 781 | // | ||
| 782 | // v = first_v_on_border] | ||
| 783 | // do { | ||
| 784 | // t = v2t_[v]; | ||
| 785 | // e = v2t_[e]; | ||
| 786 | // do something with t,e | ||
| 787 | // v = triangle_vertex(t, (e+2)%3); | ||
| 788 | // } while(v != first_v_on_border]); | ||
| 789 | |||
| 790 | geo_debug(index_t nb = 0); // for sanity check, nbr of vertices on border | ||
| 791 | // of conflict zone. | ||
| 792 | |||
| 793 | VBW::index_t first_v_on_border = END_OF_LIST; | ||
| 794 | 298 | for( | |
| 795 | VBW::ushort t = first_triangle(); | ||
| 796 |
2/2✓ Branch 0 taken 3570 times.
✓ Branch 1 taken 298 times.
|
3868 | t != END_OF_LIST; |
| 797 | t = next_triangle(t) | ||
| 798 | ) { | ||
| 799 | vbw_assert(!triangle_is_marked_as_conflict(t)); | ||
| 800 | |||
| 801 |
2/2✓ Branch 0 taken 443 times.
✓ Branch 1 taken 3127 times.
|
3570 | if(triangle_is_marked_as_conflict(triangle_adjacent(t,0))) { |
| 802 | first_v_on_border = triangle_vertex(t,1); | ||
| 803 | 443 | v2t_[first_v_on_border] = t; | |
| 804 | 443 | v2e_[first_v_on_border] = 0; | |
| 805 | geo_debug(++nb); | ||
| 806 | } | ||
| 807 |
2/2✓ Branch 0 taken 599 times.
✓ Branch 1 taken 2971 times.
|
3570 | if(triangle_is_marked_as_conflict(triangle_adjacent(t,1))) { |
| 808 | first_v_on_border = triangle_vertex(t,2); | ||
| 809 | 599 | v2t_[first_v_on_border] = t; | |
| 810 | 599 | v2e_[first_v_on_border] = 1; | |
| 811 | geo_debug(++nb); | ||
| 812 | } | ||
| 813 |
2/2✓ Branch 0 taken 588 times.
✓ Branch 1 taken 2982 times.
|
3570 | if(triangle_is_marked_as_conflict(triangle_adjacent(t,2))) { |
| 814 | first_v_on_border = triangle_vertex(t,0); | ||
| 815 | 588 | v2t_[first_v_on_border] = t; | |
| 816 | 588 | v2e_[first_v_on_border] = 2; | |
| 817 | geo_debug(++nb); | ||
| 818 | } | ||
| 819 | } | ||
| 820 | |||
| 821 | geo_debug(index_t nb2 = 0); // for sanity check, number of vertices on border | ||
| 822 | // of conflict zone (should match nb). | ||
| 823 | |||
| 824 | // Traverse the list of edges on the border of the conflict zone | ||
| 825 | // (see previous comment block for explanations). For each edge | ||
| 826 | // on the border of the conflict zone, generate a new triangle. | ||
| 827 | // - Connect it to the triangle on the border of the conflict zone | ||
| 828 | // - Connect it to the previous new triangle | ||
| 829 | // - Special case: in the end, connect the first new triangle with | ||
| 830 | // the last one. | ||
| 831 | |||
| 832 | // check we are not in the special case with all triangles in conflict | ||
| 833 |
1/2✓ Branch 0 taken 298 times.
✗ Branch 1 not taken.
|
298 | if(first_v_on_border != END_OF_LIST) { |
| 834 | VBW::index_t v = first_v_on_border; | ||
| 835 | VBW::ushort prev_new_t = VBW::ushort(-1); | ||
| 836 | VBW::ushort first_new_t = VBW::ushort(-1); | ||
| 837 | do { | ||
| 838 | geo_debug(++nb2); | ||
| 839 | 1630 | index_t t = v2t_[v]; | |
| 840 | 1630 | index_t e = v2e_[v]; | |
| 841 | 1630 | index_t v1 = triangle_vertex(t, (e+1)%3); | |
| 842 | 1630 | index_t v2 = triangle_vertex(t, (e+2)%3); | |
| 843 | vbw_assert(v1 == v); | ||
| 844 |
2/2✓ Branch 1 taken 1332 times.
✓ Branch 2 taken 298 times.
|
1630 | VBW::ushort new_t = VBW::ushort(new_triangle(lv, v2, v1)); |
| 845 | set_triangle_adjacent(new_t, 0, t); | ||
| 846 | set_triangle_adjacent(t, e, new_t); | ||
| 847 |
2/2✓ Branch 0 taken 1332 times.
✓ Branch 1 taken 298 times.
|
1630 | if(prev_new_t == VBW::ushort(-1)) { |
| 848 | first_new_t = new_t; | ||
| 849 | } else { | ||
| 850 | set_triangle_adjacent(prev_new_t, 2, new_t); | ||
| 851 | set_triangle_adjacent(new_t, 1, prev_new_t); | ||
| 852 | } | ||
| 853 | prev_new_t = new_t; | ||
| 854 | v = v2; | ||
| 855 |
2/2✓ Branch 0 taken 1332 times.
✓ Branch 1 taken 298 times.
|
1630 | } while (v != first_v_on_border); |
| 856 | set_triangle_adjacent(prev_new_t, 2, first_new_t); | ||
| 857 | set_triangle_adjacent(first_new_t, 1, prev_new_t); | ||
| 858 | } | ||
| 859 | |||
| 860 | vbw_assert(nb2 == nb); | ||
| 861 | |||
| 862 | // Recycle triangles in conflict zone | ||
| 863 | 298 | set_triangle_flags(conflict_tail, ushort(first_free_)); | |
| 864 | 298 | first_free_ = conflict_head; | |
| 865 | } | ||
| 866 | |||
| 867 | |||
| 868 | /***********************************************************************/ | ||
| 869 | |||
| 870 | 23904 | bool ConvexCell::triangle_is_in_conflict( | |
| 871 | TriangleWithFlags T, const vec4& eqn | ||
| 872 | ) const { | ||
| 873 | #ifndef STANDALONE_CONVEX_CELL | ||
| 874 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23904 times.
|
23904 | if(use_exact_predicates_) { |
| 875 | ✗ | if(T.i == VERTEX_AT_INFINITY) { | |
| 876 | ✗ | vec3 E = make_vec3(eqn.x, eqn.y, eqn.z); | |
| 877 | vec3 n2 = vertex_plane_normal(T.j); | ||
| 878 | vec3 n3 = vertex_plane_normal(T.k); | ||
| 879 | ✗ | return(GEO::PCK::det_3d(E.data(), n2.data(), n3.data()) <= 0); | |
| 880 | } | ||
| 881 | |||
| 882 | ✗ | if(T.j == VERTEX_AT_INFINITY) { | |
| 883 | ✗ | vec3 E = make_vec3(eqn.x, eqn.y, eqn.z); | |
| 884 | vec3 n3 = vertex_plane_normal(T.k); | ||
| 885 | vec3 n1 = vertex_plane_normal(T.i); | ||
| 886 | ✗ | return(GEO::PCK::det_3d(n1.data(), E.data(), n3.data()) <= 0); | |
| 887 | } | ||
| 888 | |||
| 889 | ✗ | if(T.k == VERTEX_AT_INFINITY) { | |
| 890 | ✗ | vec3 E = make_vec3(eqn.x, eqn.y, eqn.z); | |
| 891 | vec3 n1 = vertex_plane_normal(T.i); | ||
| 892 | vec3 n2 = vertex_plane_normal(T.j); | ||
| 893 | ✗ | return(GEO::PCK::det_3d(n1.data(), n2.data(), E.data()) <= 0); | |
| 894 | } | ||
| 895 | |||
| 896 | // The triangle is in conflict with eqn if the | ||
| 897 | // result of compute_triangle_point(t) injected in eqn | ||
| 898 | // has a negative sign. | ||
| 899 | // Examining the formula in compute_triangle_point(), | ||
| 900 | // this corresponds to (minus) the 4x4 determinant of the | ||
| 901 | // 4 plane equations developed w.r.t. the 4th column. | ||
| 902 | // (see Edelsbrunner - Simulation of Simplicity for similar examples | ||
| 903 | // of computations). | ||
| 904 | |||
| 905 | vec4 p1 = vertex_plane(T.i); | ||
| 906 | vec4 p2 = vertex_plane(T.j); | ||
| 907 | vec4 p3 = vertex_plane(T.k); | ||
| 908 | |||
| 909 | ✗ | return(GEO::PCK::det_4d( | |
| 910 | p1.data(), p2.data(), p3.data(), eqn.data()) >= 0 | ||
| 911 | ✗ | ); | |
| 912 | } | ||
| 913 | #endif | ||
| 914 | |||
| 915 | double det = 0.0; | ||
| 916 | |||
| 917 | // If one of the vertices of the triangle is the | ||
| 918 | // vertex at infinity, then the triangle is in conflict | ||
| 919 | // with eqn if the oriented director vector of the intersection | ||
| 920 | // between the two planes of the two other vertices | ||
| 921 | // is opposite to the normal vector to eqn. | ||
| 922 | |||
| 923 |
2/2✓ Branch 0 taken 816 times.
✓ Branch 1 taken 23088 times.
|
23904 | if(T.i == VERTEX_AT_INFINITY) { |
| 924 | vec3 n2 = vertex_plane_normal(T.j); | ||
| 925 | vec3 n3 = vertex_plane_normal(T.k); | ||
| 926 | 816 | det = -det3x3( | |
| 927 | 816 | eqn.x, n2.x, n3.x, | |
| 928 | 816 | eqn.y, n2.y, n3.y, | |
| 929 | 816 | eqn.z, n2.z, n3.z | |
| 930 | ); | ||
| 931 |
2/2✓ Branch 0 taken 1485 times.
✓ Branch 1 taken 21603 times.
|
23088 | } else if(T.j == VERTEX_AT_INFINITY) { |
| 932 | vec3 n3 = vertex_plane_normal(T.k); | ||
| 933 | vec3 n1 = vertex_plane_normal(T.i); | ||
| 934 | 1485 | det = -det3x3( | |
| 935 | 1485 | n1.x, eqn.x, n3.x, | |
| 936 | 1485 | n1.y, eqn.y, n3.y, | |
| 937 | 1485 | n1.z, eqn.z, n3.z | |
| 938 | ); | ||
| 939 |
2/2✓ Branch 0 taken 876 times.
✓ Branch 1 taken 20727 times.
|
21603 | } else if(T.k == VERTEX_AT_INFINITY) { |
| 940 | vec3 n1 = vertex_plane_normal(T.i); | ||
| 941 | vec3 n2 = vertex_plane_normal(T.j); | ||
| 942 | 876 | det = -det3x3( | |
| 943 | 876 | n1.x, n2.x, eqn.x, | |
| 944 | 876 | n1.y, n2.y, eqn.y, | |
| 945 | 876 | n1.z, n2.z, eqn.z | |
| 946 | ); | ||
| 947 | } else { | ||
| 948 | |||
| 949 | // The triangle is in conflict with eqn if the | ||
| 950 | // result of compute_triangle_point(t) injected in eqn | ||
| 951 | // has a negative sign. | ||
| 952 | // Examining the formula in compute_triangle_point(), | ||
| 953 | // this corresponds to (minus) the 4x4 determinant of the 4 plane | ||
| 954 | // equations developed w.r.t. the 4th column. | ||
| 955 | // (see Edelsbrunner - Simulation of Simplicity for similar examples | ||
| 956 | // of computations). | ||
| 957 | |||
| 958 | vec4 p1 = vertex_plane(T.i); | ||
| 959 | vec4 p2 = vertex_plane(T.j); | ||
| 960 | vec4 p3 = vertex_plane(T.k); | ||
| 961 | 20727 | det = det4x4( | |
| 962 | 20727 | p1.x, p2.x, p3.x, eqn.x, | |
| 963 | 20727 | p1.y, p2.y, p3.y, eqn.y, | |
| 964 | 20727 | p1.z, p2.z, p3.z, eqn.z, | |
| 965 | 20727 | p1.w, p2.w, p3.w, eqn.w | |
| 966 | ); | ||
| 967 | } | ||
| 968 | 23904 | return (det > 0.0); | |
| 969 | } | ||
| 970 | |||
| 971 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9684 times.
|
9684 | vec4 ConvexCell::compute_triangle_point(index_t t) const { |
| 972 | |||
| 973 | double infinite_len = 16.0; | ||
| 974 | TriangleWithFlags T = get_triangle_and_flags(t); | ||
| 975 | |||
| 976 | // Special cases with one of the three vertices at infinity. | ||
| 977 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9684 times.
|
9684 | if(T.i == VERTEX_AT_INFINITY) { |
| 978 | vec4 Pj = vertex_plane(T.j); | ||
| 979 | vec4 Pk = vertex_plane(T.k); | ||
| 980 | ✗ | vec3 Njk = normalize(cross( | |
| 981 | make_vec3(Pj.x, Pj.y, Pj.z), | ||
| 982 | make_vec3(Pk.x, Pk.y, Pk.z) | ||
| 983 | )); | ||
| 984 | ✗ | index_t t_adj = t_adj_[t].i; // vv2t(T.k, T.j); | |
| 985 | vbw_assert(!triangle_is_infinite(t_adj)); | ||
| 986 | ✗ | vec4 result = compute_triangle_point(t_adj); | |
| 987 | ✗ | result.x += result.w * Njk.x * infinite_len; | |
| 988 | ✗ | result.y += result.w * Njk.y * infinite_len; | |
| 989 | ✗ | result.z += result.w * Njk.z * infinite_len; | |
| 990 | return result; | ||
| 991 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9684 times.
|
9684 | } else if(T.j == VERTEX_AT_INFINITY) { |
| 992 | vec4 Pk = vertex_plane(T.k); | ||
| 993 | vec4 Pi = vertex_plane(T.i); | ||
| 994 | ✗ | vec3 Nki = normalize(cross( | |
| 995 | make_vec3(Pk.x, Pk.y, Pk.z), | ||
| 996 | make_vec3(Pi.x, Pi.y, Pi.z) | ||
| 997 | )); | ||
| 998 | ✗ | index_t t_adj = t_adj_[t].j; // vv2t(T.i, T.k); | |
| 999 | vbw_assert(!triangle_is_infinite(t_adj)); | ||
| 1000 | ✗ | vec4 result = compute_triangle_point(t_adj); | |
| 1001 | ✗ | result.x += result.w * Nki.x * infinite_len; | |
| 1002 | ✗ | result.y += result.w * Nki.y * infinite_len; | |
| 1003 | ✗ | result.z += result.w * Nki.z * infinite_len; | |
| 1004 | ✗ | return result; | |
| 1005 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9684 times.
|
9684 | } else if(T.k == VERTEX_AT_INFINITY) { |
| 1006 | vec4 Pi = vertex_plane(T.i); | ||
| 1007 | vec4 Pj = vertex_plane(T.j); | ||
| 1008 | ✗ | vec3 Nij = normalize(cross( | |
| 1009 | make_vec3(Pi.x, Pi.y, Pi.z), | ||
| 1010 | make_vec3(Pj.x, Pj.y, Pj.z) | ||
| 1011 | )); | ||
| 1012 | ✗ | index_t t_adj = t_adj_[t].k; // vv2t(T.j, T.i); | |
| 1013 | vbw_assert(!triangle_is_infinite(t_adj)); | ||
| 1014 | ✗ | vec4 result = compute_triangle_point(t_adj); | |
| 1015 | ✗ | result.x += result.w * Nij.x * infinite_len; | |
| 1016 | ✗ | result.y += result.w * Nij.y * infinite_len; | |
| 1017 | ✗ | result.z += result.w * Nij.z * infinite_len; | |
| 1018 | ✗ | return result; | |
| 1019 | } | ||
| 1020 | |||
| 1021 | // Get the plane equations associated with each vertex of t | ||
| 1022 | |||
| 1023 | vec4 pi1 = vertex_plane(T.i); | ||
| 1024 | vec4 pi2 = vertex_plane(T.j); | ||
| 1025 | vec4 pi3 = vertex_plane(T.k); | ||
| 1026 | |||
| 1027 | // Find the intersection of the three planes using Cramer's formula. | ||
| 1028 | // (see Edelsbrunner - Simulation of Simplicity for other examples). | ||
| 1029 | // | ||
| 1030 | // Cramer's formula: each component of the solution is obtained as | ||
| 1031 | // the ratio of two determinants: | ||
| 1032 | // - the determinant of the system where the ith column is replaced | ||
| 1033 | // with the rhs | ||
| 1034 | // divided by: | ||
| 1035 | // - the determinant of the system. | ||
| 1036 | // | ||
| 1037 | // System of equations to be solved: | ||
| 1038 | // pi1.x * x + pi1.y * y + pi1.z * z = -pi1.w | ||
| 1039 | // pi2.x * x + pi2.y * y + pi2.z * z = -pi2.w | ||
| 1040 | // pi3.x * x + pi3.y * y + pi3.z * z = -pi3.w | ||
| 1041 | // | ||
| 1042 | // Expression of the solution given by Cramer's formula: | ||
| 1043 | // | -pi1.w p1.y pi1.z | | pi1.x pi1.y pi1.z | | ||
| 1044 | // x = | -pi2.w p2.y pi2.z | / | pi2.x pi2.y pi2.z | | ||
| 1045 | // | -pi3.w p3.y pi3.z | | pi3.x pi3.y pi3.z | | ||
| 1046 | // | ||
| 1047 | // | pi1.x -p1.w pi1.z | | pi1.x pi1.y pi1.z | | ||
| 1048 | // y = | pi2.x -p2.w pi2.z | / | pi2.x pi2.y pi2.z | | ||
| 1049 | // | pi3.x -p3.w pi3.z | | pi3.x pi3.y pi3.z | | ||
| 1050 | // | ||
| 1051 | // | pi1.x p1.y -pi1.w | | pi1.x pi1.y pi1.z | | ||
| 1052 | // z = | pi2.x p2.y -pi2.w | / | pi2.x pi2.y pi2.z | | ||
| 1053 | // | pi3.x p3.y -pi3.w | | pi3.x pi3.y pi3.z | | ||
| 1054 | |||
| 1055 | vec4 result; | ||
| 1056 | |||
| 1057 | 9684 | result.x = -det3x3( | |
| 1058 | pi1.w, pi1.y, pi1.z, | ||
| 1059 | pi2.w, pi2.y, pi2.z, | ||
| 1060 | pi3.w, pi3.y, pi3.z | ||
| 1061 | ); | ||
| 1062 | |||
| 1063 | 9684 | result.y = -det3x3( | |
| 1064 | pi1.x, pi1.w, pi1.z, | ||
| 1065 | pi2.x, pi2.w, pi2.z, | ||
| 1066 | pi3.x, pi3.w, pi3.z | ||
| 1067 | ); | ||
| 1068 | |||
| 1069 | 9684 | result.z = -det3x3( | |
| 1070 | pi1.x, pi1.y, pi1.w, | ||
| 1071 | pi2.x, pi2.y, pi2.w, | ||
| 1072 | pi3.x, pi3.y, pi3.w | ||
| 1073 | ); | ||
| 1074 | |||
| 1075 | result.w = det3x3( | ||
| 1076 | pi1.x, pi1.y, pi1.z, | ||
| 1077 | pi2.x, pi2.y, pi2.z, | ||
| 1078 | pi3.x, pi3.y, pi3.z | ||
| 1079 | ); | ||
| 1080 | |||
| 1081 | 9684 | return result; | |
| 1082 | } | ||
| 1083 | |||
| 1084 | /***********************************************************************/ | ||
| 1085 | |||
| 1086 | 4 | void ConvexCell::grow_v() { | |
| 1087 | 4 | max_v_ *= 2; | |
| 1088 | 4 | plane_eqn_.resize(max_v_); | |
| 1089 | 4 | vglobal_.resize(max_v_, global_index_t(-1)); | |
| 1090 | 4 | v2t_.resize(max_v_, ushort(-1)); | |
| 1091 | 4 | v2e_.resize(max_v_, uchar(-1)); | |
| 1092 | 4 | } | |
| 1093 | |||
| 1094 | 4 | void ConvexCell::grow_t() { | |
| 1095 | 4 | max_t_ *= 2; | |
| 1096 | 4 | t_.resize(max_t_); | |
| 1097 | 4 | t_adj_.resize(max_t_); | |
| 1098 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if(has_tflags_) { |
| 1099 | ✗ | tflags_.resize(max_t_,0); | |
| 1100 | } | ||
| 1101 | 4 | } | |
| 1102 | |||
| 1103 | /***********************************************************************/ | ||
| 1104 | |||
| 1105 | ✗ | void ConvexCell::kill_vertex(index_t v) { | |
| 1106 | ✗ | for(index_t t=0; t<nb_t(); ++t) { | |
| 1107 | Triangle T = get_triangle(t); | ||
| 1108 | ✗ | if(T.i == v) { | |
| 1109 | T.i = VERTEX_AT_INFINITY; | ||
| 1110 | } | ||
| 1111 | ✗ | if(T.j == v) { | |
| 1112 | T.j = VERTEX_AT_INFINITY; | ||
| 1113 | } | ||
| 1114 | ✗ | if(T.k == v) { | |
| 1115 | T.k = VERTEX_AT_INFINITY; | ||
| 1116 | } | ||
| 1117 | ✗ | t_[t].i = T.i; | |
| 1118 | ✗ | t_[t].j = T.j; | |
| 1119 | ✗ | t_[t].k = T.k; | |
| 1120 | } | ||
| 1121 | ✗ | } | |
| 1122 | |||
| 1123 | /***********************************************************************/ | ||
| 1124 | |||
| 1125 | 442 | void ConvexCell::compute_geometry() { | |
| 1126 |
1/2✓ Branch 0 taken 442 times.
✗ Branch 1 not taken.
|
442 | if(!geometry_dirty_) { |
| 1127 | return; | ||
| 1128 | } | ||
| 1129 | |||
| 1130 | 442 | triangle_point_.resize(nb_t()); | |
| 1131 | 442 | v2t_.assign(max_v(),END_OF_LIST); | |
| 1132 | |||
| 1133 | 442 | index_t t = first_valid_; | |
| 1134 |
2/2✓ Branch 0 taken 9684 times.
✓ Branch 1 taken 442 times.
|
10126 | while(t != END_OF_LIST) { |
| 1135 | TriangleWithFlags T = get_triangle_and_flags(t); | ||
| 1136 | 9684 | vec4 p = compute_triangle_point(t); | |
| 1137 | 9684 | triangle_point_[t] = make_vec3(p.x/p.w, p.y/p.w, p.z/p.w); | |
| 1138 | 9684 | v2t_[T.i] = ushort(t); | |
| 1139 | 9684 | v2t_[T.j] = ushort(t); | |
| 1140 | 9684 | v2t_[T.k] = ushort(t); | |
| 1141 | 9684 | t = index_t(T.flags); | |
| 1142 | } | ||
| 1143 | |||
| 1144 | 442 | geometry_dirty_ = false; | |
| 1145 | } | ||
| 1146 | |||
| 1147 | ✗ | inline double triangle_area(vec3 p1, vec3 p2, vec3 p3) { | |
| 1148 | ✗ | double Ux = p2.x - p1.x; | |
| 1149 | ✗ | double Uy = p2.y - p1.y; | |
| 1150 | ✗ | double Uz = p2.z - p1.z; | |
| 1151 | ✗ | double Vx = p3.x - p1.x; | |
| 1152 | ✗ | double Vy = p3.y - p1.y; | |
| 1153 | ✗ | double Vz = p3.z - p1.z; | |
| 1154 | ✗ | double Wx = Uy * Vz - Uz * Vy; | |
| 1155 | ✗ | double Wy = Uz * Vx - Ux * Vz; | |
| 1156 | ✗ | double Wz = Ux * Vy - Uy * Vx; | |
| 1157 | ✗ | return 0.5 * ::sqrt(Wx*Wx + Wy*Wy + Wz*Wz); | |
| 1158 | } | ||
| 1159 | |||
| 1160 | ✗ | double ConvexCell::facet_area(index_t v) const { | |
| 1161 | vbw_assert(v < nb_v()); | ||
| 1162 | vbw_assert(!geometry_dirty_); | ||
| 1163 | |||
| 1164 | ushort t1t2[2]; | ||
| 1165 | index_t cur=0; | ||
| 1166 | double result = 0.0; | ||
| 1167 | |||
| 1168 | ✗ | if(v2t_[v] != END_OF_LIST) { | |
| 1169 | ✗ | index_t t = v2t_[v]; | |
| 1170 | index_t count = 0; | ||
| 1171 | do { | ||
| 1172 | ✗ | if(cur < 2) { | |
| 1173 | ✗ | t1t2[cur] = ushort(t); | |
| 1174 | } else { | ||
| 1175 | ✗ | result += triangle_area( | |
| 1176 | ✗ | triangle_point_[t1t2[0]], | |
| 1177 | ✗ | triangle_point_[t1t2[1]], | |
| 1178 | triangle_point_[t] | ||
| 1179 | ); | ||
| 1180 | ✗ | t1t2[1] = ushort(t); | |
| 1181 | } | ||
| 1182 | ✗ | ++cur; | |
| 1183 | index_t lv = triangle_find_vertex(t,v); | ||
| 1184 | ✗ | t = triangle_adjacent(t, (lv + 1)%3); | |
| 1185 | ++count; | ||
| 1186 | ✗ | geo_assert(count < 100000); | |
| 1187 | ✗ | } while(t != v2t_[v]); | |
| 1188 | } | ||
| 1189 | |||
| 1190 | ✗ | return result; | |
| 1191 | } | ||
| 1192 | |||
| 1193 | 16204 | inline double tet_volume(vec3 p1, vec3 p2, vec3 p3, vec3 p4) { | |
| 1194 | 16204 | double Ux = p2.x - p1.x; | |
| 1195 | 16204 | double Uy = p2.y - p1.y; | |
| 1196 | 16204 | double Uz = p2.z - p1.z; | |
| 1197 | |||
| 1198 | 16204 | double Vx = p3.x - p1.x; | |
| 1199 | 16204 | double Vy = p3.y - p1.y; | |
| 1200 | 16204 | double Vz = p3.z - p1.z; | |
| 1201 | |||
| 1202 | 16204 | double Wx = p4.x - p1.x; | |
| 1203 | 16204 | double Wy = p4.y - p1.y; | |
| 1204 | 16204 | double Wz = p4.z - p1.z; | |
| 1205 | |||
| 1206 | 16204 | double UVx = Uy * Vz - Uz * Vy; | |
| 1207 | 16204 | double UVy = Uz * Vx - Ux * Vz; | |
| 1208 | 16204 | double UVz = Ux * Vy - Uy * Vx; | |
| 1209 | |||
| 1210 | 16204 | return ::fabs( | |
| 1211 | 16204 | UVx * Wx + UVy * Wy + UVz * Wz | |
| 1212 | 16204 | ) / 6.0; | |
| 1213 | } | ||
| 1214 | |||
| 1215 | ✗ | double ConvexCell::volume() const { | |
| 1216 | vbw_assert(!geometry_dirty_); | ||
| 1217 | double result = 0.0; | ||
| 1218 | |||
| 1219 | ushort t_origin = END_OF_LIST; | ||
| 1220 | ✗ | for(index_t v=0; v<nb_v_; ++v) { | |
| 1221 | ✗ | if(v2t_[v] == END_OF_LIST) { | |
| 1222 | ✗ | continue; | |
| 1223 | } | ||
| 1224 | ✗ | if(t_origin == END_OF_LIST) { | |
| 1225 | t_origin = v2t_[v]; | ||
| 1226 | ✗ | continue; | |
| 1227 | } | ||
| 1228 | ushort t1t2[2]; | ||
| 1229 | index_t cur=0; | ||
| 1230 | ✗ | index_t t = v2t_[v]; | |
| 1231 | |||
| 1232 | index_t count = 0; | ||
| 1233 | do { | ||
| 1234 | ✗ | if(cur < 2) { | |
| 1235 | ✗ | t1t2[cur] = ushort(t); | |
| 1236 | } else { | ||
| 1237 | ✗ | result += tet_volume( | |
| 1238 | triangle_point_[t_origin], | ||
| 1239 | ✗ | triangle_point_[t1t2[0]], | |
| 1240 | ✗ | triangle_point_[t1t2[1]], | |
| 1241 | triangle_point_[t] | ||
| 1242 | ); | ||
| 1243 | ✗ | t1t2[1] = ushort(t); | |
| 1244 | } | ||
| 1245 | ✗ | ++cur; | |
| 1246 | index_t lv = triangle_find_vertex(t,v); | ||
| 1247 | ✗ | t = triangle_adjacent(t, (lv + 1)%3); | |
| 1248 | ++count; | ||
| 1249 | ✗ | geo_assert(count < 100000); | |
| 1250 | ✗ | } while(t != v2t_[v]); | |
| 1251 | } | ||
| 1252 | ✗ | return result; | |
| 1253 | } | ||
| 1254 | |||
| 1255 | 442 | vec3 ConvexCell::barycenter() const { | |
| 1256 | vec3 result; | ||
| 1257 | double m; | ||
| 1258 | 442 | compute_mg(m, result); | |
| 1259 |
1/2✓ Branch 0 taken 442 times.
✗ Branch 1 not taken.
|
442 | if(m != 0.0) { |
| 1260 | 442 | result.x /= m; | |
| 1261 | 442 | result.y /= m; | |
| 1262 | 442 | result.z /= m; | |
| 1263 | } | ||
| 1264 | 442 | return result; | |
| 1265 | } | ||
| 1266 | |||
| 1267 | 442 | void ConvexCell::compute_mg(double& m, vec3& result) const { | |
| 1268 | vbw_assert(!geometry_dirty_); | ||
| 1269 | 442 | result = make_vec3(0.0, 0.0, 0.0); | |
| 1270 | 442 | m = 0.0; | |
| 1271 | |||
| 1272 | ushort t_origin = END_OF_LIST; | ||
| 1273 |
2/2✓ Branch 0 taken 7354 times.
✓ Branch 1 taken 442 times.
|
7796 | for(index_t v=0; v<nb_v_; ++v) { |
| 1274 |
2/2✓ Branch 0 taken 1628 times.
✓ Branch 1 taken 5726 times.
|
7354 | if(v2t_[v] == END_OF_LIST) { |
| 1275 | 2070 | continue; | |
| 1276 | } | ||
| 1277 |
2/2✓ Branch 0 taken 442 times.
✓ Branch 1 taken 5284 times.
|
5726 | if(t_origin == END_OF_LIST) { |
| 1278 | t_origin = v2t_[v]; | ||
| 1279 | 442 | continue; | |
| 1280 | } | ||
| 1281 | ushort t1t2[2]; | ||
| 1282 | index_t cur=0; | ||
| 1283 | 5284 | index_t t = v2t_[v]; | |
| 1284 | index_t count = 0; | ||
| 1285 | do { | ||
| 1286 |
2/2✓ Branch 0 taken 10568 times.
✓ Branch 1 taken 16204 times.
|
26772 | if(cur < 2) { |
| 1287 | 10568 | t1t2[cur] = ushort(t); | |
| 1288 | } else { | ||
| 1289 | 16204 | vec3 p = triangle_point_[t_origin]; | |
| 1290 | 16204 | vec3 q = triangle_point_[t1t2[0]]; | |
| 1291 | 16204 | vec3 r = triangle_point_[t1t2[1]]; | |
| 1292 | 16204 | vec3 s = triangle_point_[t]; | |
| 1293 | 16204 | double cur_m = tet_volume(p,q,r,s); | |
| 1294 | 16204 | m += cur_m; | |
| 1295 | 16204 | result.x += cur_m*(p.x + q.x + r.x + s.x)/4.0; | |
| 1296 | 16204 | result.y += cur_m*(p.y + q.y + r.y + s.y)/4.0; | |
| 1297 | 16204 | result.z += cur_m*(p.z + q.z + r.z + s.z)/4.0; | |
| 1298 | 16204 | t1t2[1] = ushort(t); | |
| 1299 | } | ||
| 1300 |
2/2✓ Branch 0 taken 17898 times.
✓ Branch 1 taken 8874 times.
|
26772 | ++cur; |
| 1301 | index_t lv = triangle_find_vertex(t,v); | ||
| 1302 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 26772 times.
|
26772 | t = triangle_adjacent(t, (lv + 1)%3); |
| 1303 | ++count; | ||
| 1304 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 26772 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
26772 | geo_assert(count < 100000); |
| 1305 |
2/2✓ Branch 0 taken 21488 times.
✓ Branch 1 taken 5284 times.
|
26772 | } while(t != v2t_[v]); |
| 1306 | } | ||
| 1307 | 442 | } | |
| 1308 | |||
| 1309 | /***********************************************************************/ | ||
| 1310 | |||
| 1311 | |||
| 1312 | ✗ | double ConvexCell::squared_radius(vec3 center) const { | |
| 1313 | ✗ | double result = 0.0; | |
| 1314 | ✗ | index_t t = first_valid_; | |
| 1315 | ✗ | while(t != END_OF_LIST) { | |
| 1316 | TriangleWithFlags T = get_triangle_and_flags(t); | ||
| 1317 | ✗ | if(geometry_dirty_) { | |
| 1318 | ✗ | vec4 p4 = compute_triangle_point(t); | |
| 1319 | ✗ | vec3 p3 = make_vec3( | |
| 1320 | ✗ | p4.x/p4.w, p4.y/p4.w, p4.z/p4.w | |
| 1321 | ); | ||
| 1322 | ✗ | result = std::max(result, squared_distance(center,p3)); | |
| 1323 | } else { | ||
| 1324 | ✗ | vec3 p = triangle_point_[t]; | |
| 1325 | ✗ | result = std::max(result, squared_distance(center,p)); | |
| 1326 | } | ||
| 1327 | ✗ | t = index_t(T.flags); | |
| 1328 | } | ||
| 1329 | ✗ | return result; | |
| 1330 | } | ||
| 1331 | |||
| 1332 | ✗ | double ConvexCell::squared_inner_radius(vec3 center) const { | |
| 1333 | ✗ | double result = std::numeric_limits<double>::max(); | |
| 1334 | ✗ | for(index_t v=0; v<nb_v(); ++v) { | |
| 1335 | vec4 P = vertex_plane(v); | ||
| 1336 | // Ignore vertex at infinity. | ||
| 1337 | ✗ | if(P.x == 0.0 && P.y == 0.0 && P.z == 0.0) { | |
| 1338 | continue; | ||
| 1339 | } | ||
| 1340 | ✗ | result = std::min( | |
| 1341 | result, squared_point_plane_distance(center, P) | ||
| 1342 | ); | ||
| 1343 | } | ||
| 1344 | ✗ | return result; | |
| 1345 | } | ||
| 1346 | |||
| 1347 | |||
| 1348 |
2/2✓ Branch 0 taken 431 times.
✓ Branch 1 taken 11 times.
|
442 | void ConvexCell::connect_triangles() { |
| 1349 | |||
| 1350 | // create array that maps vertices pairs to triangles. | ||
| 1351 | // size of the array is nb_v squared. | ||
| 1352 | // If nb_v is small, allocate it on the stack, else | ||
| 1353 | // allocate it on the heap (allocating on the stack is | ||
| 1354 | // interesting for multithreading). | ||
| 1355 | |||
| 1356 | const index_t MAX_NV_ON_STACK = 50; | ||
| 1357 | index_t NV = nb_v(); | ||
| 1358 | ushort* vv2t = | ||
| 1359 |
2/2✓ Branch 0 taken 431 times.
✓ Branch 1 taken 11 times.
|
442 | (NV <= MAX_NV_ON_STACK) ? (ushort*)alloca(NV*NV*sizeof(ushort)) |
| 1360 | 11 | : new ushort[NV*NV] | |
| 1361 | ; | ||
| 1362 | |||
| 1363 | #ifdef GEO_DEBUG | ||
| 1364 | for(index_t i=0; i<NV*NV; ++i) { | ||
| 1365 | vv2t[i] = END_OF_LIST; | ||
| 1366 | } | ||
| 1367 | #endif | ||
| 1368 | |||
| 1369 | // For each triangle t, remember | ||
| 1370 | // that t is adjacent to the three | ||
| 1371 | // oriented edges (j,k), (k,i), (i,j) | ||
| 1372 | 442 | for( | |
| 1373 | ushort t = first_triangle(); | ||
| 1374 |
2/2✓ Branch 0 taken 9762 times.
✓ Branch 1 taken 442 times.
|
10204 | t != END_OF_LIST; |
| 1375 | t = next_triangle(t) | ||
| 1376 | ) { | ||
| 1377 | 9762 | Triangle T = t_[t]; | |
| 1378 | vbw_assert(T.i < nb_v()); | ||
| 1379 | vbw_assert(T.j < nb_v()); | ||
| 1380 | vbw_assert(T.k < nb_v()); | ||
| 1381 | 9762 | vv2t[NV*T.j + T.k] = t; | |
| 1382 | 9762 | vv2t[NV*T.k + T.i] = t; | |
| 1383 | 9762 | vv2t[NV*T.i + T.j] = t; | |
| 1384 | } | ||
| 1385 | |||
| 1386 | // For each triangle t, find the | ||
| 1387 | // three triangles adjacent to its | ||
| 1388 | // three edges (k,j), (i,k), (j,i) | ||
| 1389 | // (in reverse order because we | ||
| 1390 | // want to find the three triangles | ||
| 1391 | // on the other side of the edges) | ||
| 1392 | for( | ||
| 1393 | ushort t = first_triangle(); | ||
| 1394 |
2/2✓ Branch 0 taken 9762 times.
✓ Branch 1 taken 442 times.
|
10204 | t != END_OF_LIST; |
| 1395 | t = next_triangle(t) | ||
| 1396 | ) { | ||
| 1397 | 9762 | Triangle T = t_[t]; | |
| 1398 | vbw_assert(vv2t[NV*T.j + T.i] != END_OF_LIST); | ||
| 1399 | vbw_assert(vv2t[NV*T.k + T.j] != END_OF_LIST); | ||
| 1400 | vbw_assert(vv2t[NV*T.i + T.k] != END_OF_LIST); | ||
| 1401 | 9762 | t_adj_[t] = make_triangle( | |
| 1402 | 9762 | vv2t[NV*T.k + T.j], | |
| 1403 | 9762 | vv2t[NV*T.i + T.k], | |
| 1404 | 9762 | vv2t[NV*T.j + T.i] | |
| 1405 | ); | ||
| 1406 | } | ||
| 1407 | |||
| 1408 |
2/2✓ Branch 0 taken 11 times.
✓ Branch 1 taken 431 times.
|
442 | if(NV > MAX_NV_ON_STACK) { |
| 1409 | 11 | delete[] vv2t; | |
| 1410 | } | ||
| 1411 | 442 | } | |
| 1412 | |||
| 1413 | /************************************************************************/ | ||
| 1414 | |||
| 1415 | |||
| 1416 | } | ||
| 1417 |