| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * Copyright (c) 2000-2023 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/numerics/exact_geometry.h> | ||
| 41 | #include <geogram/numerics/interval_nt.h> | ||
| 42 | #include <geogram/numerics/predicates.h> | ||
| 43 | #include <geogram/numerics/PCK.h> | ||
| 44 | #include <geogram/basic/logger.h> | ||
| 45 | |||
| 46 | namespace GEO { | ||
| 47 | |||
| 48 | |||
| 49 | namespace PCK { | ||
| 50 | |||
| 51 | 640064 | Sign orient_2d( | |
| 52 | const vec2HE& p0, const vec2HE& p1, const vec2HE& p2 | ||
| 53 | ) { | ||
| 54 |
4/4✓ Branch 0 taken 35 times.
✓ Branch 1 taken 640029 times.
✓ Branch 3 taken 34 times.
✓ Branch 4 taken 1 times.
|
640064 | static PredicateStats stats("orient_2d(vec2HE)"); |
| 55 | stats.log_invoke(); | ||
| 56 | // Filter, using interval arithmetics | ||
| 57 | { | ||
| 58 | interval_nt::Rounding rounding; | ||
| 59 | 640064 | interval_nt Delta = det3x3( | |
| 60 | 640064 | interval_nt(p0.x),interval_nt(p0.y),interval_nt(p0.w), | |
| 61 | 640064 | interval_nt(p1.x),interval_nt(p1.y),interval_nt(p1.w), | |
| 62 | 640064 | interval_nt(p2.x),interval_nt(p2.y),interval_nt(p2.w) | |
| 63 | ); | ||
| 64 |
2/2✓ Branch 0 taken 7953 times.
✓ Branch 1 taken 632111 times.
|
640064 | interval_nt::Sign2 s = Delta.sign(); |
| 65 | if(interval_nt::sign_is_determined(s)) { | ||
| 66 | return Sign( | ||
| 67 |
1/2✓ Branch 1 taken 632111 times.
✗ Branch 2 not taken.
|
1264222 | interval_nt::convert_sign(s)* |
| 68 |
1/2✓ Branch 0 taken 632111 times.
✗ Branch 1 not taken.
|
632111 | p0.w.sign()*p1.w.sign()*p2.w.sign() |
| 69 | 632111 | ); | |
| 70 | } | ||
| 71 | } | ||
| 72 | stats.log_exact(); | ||
| 73 | #ifdef GEO_HAS_BIG_STACK | ||
| 74 | 7953 | const expansion& Delta = expansion_det3x3( | |
| 75 | p0.x.rep(), p0.y.rep(), p0.w.rep(), | ||
| 76 | p1.x.rep(), p1.y.rep(), p1.w.rep(), | ||
| 77 | p2.x.rep(), p2.y.rep(), p2.w.rep() | ||
| 78 | ); | ||
| 79 | #else | ||
| 80 | expansion_nt Delta = det3x3( | ||
| 81 | p0.x, p0.y, p0.w, | ||
| 82 | p1.x, p1.y, p1.w, | ||
| 83 | p2.x, p2.y, p2.w | ||
| 84 | ); | ||
| 85 | #endif | ||
| 86 | return Sign( | ||
| 87 |
1/2✓ Branch 0 taken 7953 times.
✗ Branch 1 not taken.
|
7953 | Delta.sign()* |
| 88 |
1/2✓ Branch 0 taken 7953 times.
✗ Branch 1 not taken.
|
7953 | p0.w.rep().sign()* |
| 89 | p1.w.rep().sign()* | ||
| 90 | p2.w.rep().sign() | ||
| 91 | 7953 | ); | |
| 92 | } | ||
| 93 | |||
| 94 | 23591 | Sign orient_3d( | |
| 95 | const vec3HE& p0, const vec3HE& p1, | ||
| 96 | const vec3HE& p2, const vec3HE& p3 | ||
| 97 | ) { | ||
| 98 |
3/4✓ Branch 0 taken 23 times.
✓ Branch 1 taken 23568 times.
✓ Branch 3 taken 23 times.
✗ Branch 4 not taken.
|
23591 | static PredicateStats stats("orient_3d(vec3HE)"); |
| 99 | stats.log_invoke(); | ||
| 100 | // Filter | ||
| 101 | { | ||
| 102 | interval_nt::Rounding rounding; | ||
| 103 | 23591 | vec3HI p0I(p0); | |
| 104 | 23591 | vec3HI U = vec3HI(p1)-p0I; | |
| 105 | 23591 | vec3HI V = vec3HI(p2)-p0I; | |
| 106 | 23591 | vec3HI W = vec3HI(p3)-p0I; | |
| 107 | 23591 | interval_nt::Sign2 s1 = U.w.sign(); | |
| 108 | 23591 | interval_nt::Sign2 s2 = V.w.sign(); | |
| 109 |
1/2✓ Branch 0 taken 23591 times.
✗ Branch 1 not taken.
|
23591 | interval_nt::Sign2 s3 = W.w.sign(); |
| 110 | if( | ||
| 111 |
1/2✓ Branch 0 taken 23591 times.
✗ Branch 1 not taken.
|
23591 | interval_nt::sign_is_non_zero(s1) && |
| 112 |
2/4✓ Branch 0 taken 23591 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 23591 times.
✗ Branch 3 not taken.
|
47182 | interval_nt::sign_is_non_zero(s2) && |
| 113 | interval_nt::sign_is_non_zero(s3) | ||
| 114 | ) { | ||
| 115 | 23591 | interval_nt Delta = det3x3( | |
| 116 | U.x, U.y, U.z, | ||
| 117 | V.x, V.y, V.z, | ||
| 118 | W.x, W.y, W.z | ||
| 119 | ); | ||
| 120 |
2/2✓ Branch 0 taken 16096 times.
✓ Branch 1 taken 7495 times.
|
23591 | interval_nt::Sign2 s = Delta.sign(); |
| 121 |
2/2✓ Branch 0 taken 16096 times.
✓ Branch 1 taken 7495 times.
|
23591 | if(interval_nt::sign_is_non_zero(s)) { |
| 122 | return Sign( | ||
| 123 |
1/2✓ Branch 1 taken 16096 times.
✗ Branch 2 not taken.
|
16096 | interval_nt::convert_sign(s)* |
| 124 |
1/2✓ Branch 1 taken 16096 times.
✗ Branch 2 not taken.
|
16096 | interval_nt::convert_sign(s1)* |
| 125 |
1/2✓ Branch 1 taken 16096 times.
✗ Branch 2 not taken.
|
16096 | interval_nt::convert_sign(s2)* |
| 126 |
1/2✓ Branch 1 taken 16096 times.
✗ Branch 2 not taken.
|
16096 | interval_nt::convert_sign(s3) |
| 127 | 16096 | ); | |
| 128 | } | ||
| 129 | } | ||
| 130 | } | ||
| 131 | |||
| 132 | stats.log_exact(); | ||
| 133 | |||
| 134 | 7495 | vec3HE U = p1-p0; | |
| 135 |
1/2✓ Branch 1 taken 7495 times.
✗ Branch 2 not taken.
|
7495 | vec3HE V = p2-p0; |
| 136 |
1/2✓ Branch 1 taken 7495 times.
✗ Branch 2 not taken.
|
7495 | vec3HE W = p3-p0; |
| 137 | |||
| 138 | // Here we do not use expansion_det3x3() (that | ||
| 139 | // allocates on the stack), because | ||
| 140 | // RadialSort uses generted points that | ||
| 141 | // can have very looonng expansions that | ||
| 142 | // can cause stack overflow. | ||
| 143 | expansion_nt Delta = det3x3( | ||
| 144 | U.x, U.y, U.z, | ||
| 145 | V.x, V.y, V.z, | ||
| 146 | W.x, W.y, W.z | ||
| 147 | ); | ||
| 148 | |||
| 149 | Sign result = Sign( | ||
| 150 |
1/2✓ Branch 0 taken 7495 times.
✗ Branch 1 not taken.
|
7495 | Delta.sign()* |
| 151 |
1/2✓ Branch 0 taken 7495 times.
✗ Branch 1 not taken.
|
7495 | U.w.rep().sign()* |
| 152 | V.w.rep().sign()* | ||
| 153 | W.w.rep().sign() | ||
| 154 | 7495 | ); | |
| 155 | |||
| 156 | return result; | ||
| 157 | 7495 | } | |
| 158 | |||
| 159 | 4180787 | Sign orient_2d_projected( | |
| 160 | const vec3HE& p0, const vec3HE& p1, const vec3HE& p2, | ||
| 161 | coord_index_t axis | ||
| 162 | ) { | ||
| 163 |
4/4✓ Branch 0 taken 47 times.
✓ Branch 1 taken 4180740 times.
✓ Branch 3 taken 34 times.
✓ Branch 4 taken 13 times.
|
4180787 | static PredicateStats stats("orient_2d_projected(vec3HE)"); |
| 164 | stats.log_invoke(); | ||
| 165 | |||
| 166 | 4180787 | coord_index_t u = coord_index_t((axis+1)%3); | |
| 167 | 4180787 | coord_index_t v = coord_index_t((axis+2)%3); | |
| 168 | |||
| 169 | // Filter, using interval arithmetics | ||
| 170 | { | ||
| 171 | interval_nt::Rounding rounding; | ||
| 172 | |||
| 173 | 4180787 | interval_nt Delta = det3x3( | |
| 174 | 4180787 | interval_nt(p0[u]),interval_nt(p0[v]),interval_nt(p0.w), | |
| 175 | 4180787 | interval_nt(p1[u]),interval_nt(p1[v]),interval_nt(p1.w), | |
| 176 | 4180787 | interval_nt(p2[u]),interval_nt(p2[v]),interval_nt(p2.w) | |
| 177 | ); | ||
| 178 |
2/2✓ Branch 0 taken 861190 times.
✓ Branch 1 taken 3319597 times.
|
4180787 | interval_nt::Sign2 s = Delta.sign(); |
| 179 | if(interval_nt::sign_is_determined(s)) { | ||
| 180 | return Sign( | ||
| 181 |
1/2✓ Branch 1 taken 3319597 times.
✗ Branch 2 not taken.
|
6639194 | interval_nt::convert_sign(s)* |
| 182 |
1/2✓ Branch 0 taken 3319597 times.
✗ Branch 1 not taken.
|
3319597 | p0.w.sign()*p1.w.sign()*p2.w.sign() |
| 183 | 3319597 | ); | |
| 184 | } | ||
| 185 | } | ||
| 186 | |||
| 187 | stats.log_exact(); | ||
| 188 | |||
| 189 | Sign result = ZERO; | ||
| 190 | { | ||
| 191 | #ifdef GEO_HAS_BIG_STACK | ||
| 192 | 861190 | const expansion& Delta = expansion_det3x3( | |
| 193 | p0[u].rep(), p0[v].rep(), p0.w.rep(), | ||
| 194 | p1[u].rep(), p1[v].rep(), p1.w.rep(), | ||
| 195 | p2[u].rep(), p2[v].rep(), p2.w.rep() | ||
| 196 | ); | ||
| 197 | #else | ||
| 198 | expansion_nt Delta = det3x3( | ||
| 199 | p0[u], p0[v], p0.w, | ||
| 200 | p1[u], p1[v], p1.w, | ||
| 201 | p2[u], p2[v], p2.w | ||
| 202 | ); | ||
| 203 | #endif | ||
| 204 | 861190 | result = Sign( | |
| 205 |
1/2✓ Branch 0 taken 861190 times.
✗ Branch 1 not taken.
|
861190 | Delta.sign()* |
| 206 |
1/2✓ Branch 0 taken 861190 times.
✗ Branch 1 not taken.
|
861190 | p0.w.rep().sign()* |
| 207 | p1.w.rep().sign()* | ||
| 208 | p2.w.rep().sign() | ||
| 209 | ); | ||
| 210 | } | ||
| 211 | 861190 | return result; | |
| 212 | } | ||
| 213 | |||
| 214 | ✗ | Sign dot_2d(const vec2HE& p0, const vec2HE& p1, const vec2HE& p2) { | |
| 215 | ✗ | static PredicateStats stats("dot_2d(vec2HE)"); | |
| 216 | stats.log_invoke(); | ||
| 217 | |||
| 218 | // TODO: filter | ||
| 219 | |||
| 220 | ✗ | vec2HE U = p1 - p0; | |
| 221 | ✗ | vec2HE V = p2 - p0; | |
| 222 | #ifdef GEO_HAS_BIG_STACK | ||
| 223 | ✗ | const expansion& x1x2 = expansion_product(U.x.rep(), V.x.rep()); | |
| 224 | ✗ | const expansion& y1y2 = expansion_product(U.y.rep(), V.y.rep()); | |
| 225 | ✗ | const expansion& S = expansion_sum(x1x2, y1y2); | |
| 226 | #else | ||
| 227 | expansion_nt S = U.x*V.x+U.y*V.y; | ||
| 228 | #endif | ||
| 229 | ✗ | return Sign(S.sign()*U.w.sign()*V.w.sign()); | |
| 230 | ✗ | } | |
| 231 | |||
| 232 | /******************************************************************************/ | ||
| 233 | |||
| 234 | /** | ||
| 235 | * \brief Computes the sign of | ||
| 236 | * det3x3(x1,y1,1,x2,y2,1,x3,y3,1) | ||
| 237 | * \param[in] p1 , p2 , p3 the three points in | ||
| 238 | * homogeneous exact coordiates | ||
| 239 | * \return the sign of the determinant | ||
| 240 | */ | ||
| 241 | 12735 | static inline Sign det3_111_sign( | |
| 242 | const vec2HE& p1, | ||
| 243 | const vec2HE& p2, | ||
| 244 | const vec2HE& p3 | ||
| 245 | ) { | ||
| 246 | 12735 | expansion_nt m1 = det2x2(p2.x, p2.y, p3.x, p3.y); | |
| 247 |
1/2✓ Branch 1 taken 12735 times.
✗ Branch 2 not taken.
|
12735 | expansion_nt m2 = det2x2(p1.x, p1.y, p3.x, p3.y); |
| 248 | expansion_nt m3 = det2x2(p1.x, p1.y, p2.x, p2.y); | ||
| 249 | m1.optimize(); m2.optimize(); m3.optimize(); | ||
| 250 |
8/22✓ Branch 1 taken 12735 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12735 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12735 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12735 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 12735 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 12735 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 12735 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 12735 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
|
50940 | expansion_nt D = p1.w*m1-p2.w*m2+p3.w*m3 ; |
| 251 |
3/6✓ Branch 0 taken 12735 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 12735 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 12735 times.
✗ Branch 5 not taken.
|
50940 | return Sign(p1.w.sign()*p2.w.sign()*p3.w.sign()*D.sign()); |
| 252 | } | ||
| 253 | |||
| 254 | 343754 | Sign incircle_2d_SOS_with_lengths( | |
| 255 | const vec2HE& p0, const vec2HE& p1, | ||
| 256 | const vec2HE& p2, const vec2HE& p3, | ||
| 257 | double l0, double l1, double l2, double l3 | ||
| 258 | ) { | ||
| 259 |
4/4✓ Branch 0 taken 39 times.
✓ Branch 1 taken 343715 times.
✓ Branch 3 taken 35 times.
✓ Branch 4 taken 4 times.
|
343754 | static PredicateStats stats("incircle_2d_SOS_with_lengths(vec2HE)"); |
| 260 | stats.log_invoke(); | ||
| 261 | |||
| 262 | Sign result = ZERO; | ||
| 263 | |||
| 264 | // "Documentation is a love letter that you write to your | ||
| 265 | // future self." - Damian Conway (or that you write to the | ||
| 266 | // poor guy who will have one day to dive again in this stuff, | ||
| 267 | // but it's probably me anyway) | ||
| 268 | // | ||
| 269 | // Determinant to compute: | ||
| 270 | // | x0 y0 l0 1 | | ||
| 271 | // | x1 y1 l1 1 | | ||
| 272 | // | x2 y2 l2 1 | | ||
| 273 | // | x3 y3 l3 1 | | ||
| 274 | // where li = xi^2 + yi^2 | ||
| 275 | // (positive if (p0,p1,p2) counterclockwise and p3 in circumcircle | ||
| 276 | // of (p0,p1,p2)). Sign changes if (p0,p1,p2) is clockwise). | ||
| 277 | // We suppose that the li's are *given* numbers (it is like | ||
| 278 | // perturbating a regular (weighted) triangulation instead of a | ||
| 279 | // Delaunay triangulation). | ||
| 280 | // It allows to use arithmetic expansions without | ||
| 281 | // overflowing/underflowing too soon. | ||
| 282 | // | ||
| 283 | // Subtract last row to first three rows | ||
| 284 | // (does not change determinant): | ||
| 285 | // | ||
| 286 | // | x0-x3 y0-y3 l0-l3 0 | | ||
| 287 | // | x1-x3 y1-y3 l1-l3 0 | | ||
| 288 | // | x2-x3 y2-y3 l2-l3 0 | | ||
| 289 | // | x3 y3 l3 1 | | ||
| 290 | // | ||
| 291 | // Develop along last column: | ||
| 292 | // | x0-x3 y0-y3 l0-l3 | | ||
| 293 | // | x1-x3 y1-y3 l1-l3 | | ||
| 294 | // | x2-x3 y2-y3 l2-l3 | | ||
| 295 | // | ||
| 296 | // let (Xi+1,Yi+1,Wi+1) = (xi,yi)-(x3,y3) in homogeneous coordinates | ||
| 297 | // let Li+1 = li-l3: | ||
| 298 | // | X1/W1 Y1/W1 L1 | | ||
| 299 | // | X2/W2 Y2/W2 L2 | | ||
| 300 | // | X3/W3 Y3/W3 L3 | | ||
| 301 | // | ||
| 302 | // Develop along last column, factor-out the Wi's | ||
| 303 | // | ||
| 304 | // | X2 Y2 | | X1 Y1 | | X1 Y1 | | ||
| 305 | // (L1/W2W3) | X3 Y3 | - (L2/W1W3) | X3 Y3 | + (L3/W1W2) | X2 Y2 | | ||
| 306 | // | ||
| 307 | // Multiply everything by W1W2W3: | ||
| 308 | // | ||
| 309 | // | X2 Y2 | | X1 Y1 | | X1 Y1 | | ||
| 310 | // sign( L1W1 | X3 Y3 | - L2W2 | X3 Y3 | + L3W3 | X2 Y2 | ) * | ||
| 311 | // sign(W1) * sign(W2) * sign(W3) | ||
| 312 | |||
| 313 | // The four approximated li's. It is OK since they will | ||
| 314 | // always have the same value for the same vertex. | ||
| 315 | // We do it like that because computing them exactly (and | ||
| 316 | // properly propagating the wi^2's) makes expansions | ||
| 317 | // overflow/underflow. | ||
| 318 | // It is like perturbating a regular (weighted) triangulation | ||
| 319 | // instead of a Delaunay triangulation. | ||
| 320 | // However, if incircle(p1,p2,p3,p4) is lower than 0, it does | ||
| 321 | // not imply that (p1,p2,p3,p4) forms a convex quadrilateral | ||
| 322 | // (needs to be tested in addition, it is what CDT2d does | ||
| 323 | // when exact_incircle_ is set to false). | ||
| 324 | |||
| 325 | // We could also compute them each time, as in incircle_2d_SOS() | ||
| 326 | // (but we are caching them in MeshSurfaceIntersection's temporary | ||
| 327 | // Vertex objects), this gains 20-25% performance so it is | ||
| 328 | // worth it. | ||
| 329 | |||
| 330 | // Filter | ||
| 331 | { | ||
| 332 | interval_nt::Rounding rounding; | ||
| 333 | interval_nt l3I(l3); | ||
| 334 | 343754 | interval_nt L1 = interval_nt(l0) - l3I; | |
| 335 | 343754 | interval_nt L2 = interval_nt(l1) - l3I; | |
| 336 | 343754 | interval_nt L3 = interval_nt(l2) - l3I; | |
| 337 | |||
| 338 | 343754 | vec2HI p3I(p3); | |
| 339 | 343754 | vec2HI P1 = vec2HI(p0) - p3I; | |
| 340 | 343754 | vec2HI P2 = vec2HI(p1) - p3I; | |
| 341 | 343754 | vec2HI P3 = vec2HI(p2) - p3I; | |
| 342 | |||
| 343 | 343754 | interval_nt::Sign2 s1 = P1.w.sign(); | |
| 344 | 343754 | interval_nt::Sign2 s2 = P2.w.sign(); | |
| 345 |
1/2✓ Branch 0 taken 343754 times.
✗ Branch 1 not taken.
|
343754 | interval_nt::Sign2 s3 = P3.w.sign(); |
| 346 | |||
| 347 | if( | ||
| 348 |
1/2✓ Branch 0 taken 343754 times.
✗ Branch 1 not taken.
|
343754 | interval_nt::sign_is_non_zero(s1) && |
| 349 |
2/4✓ Branch 0 taken 343754 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 343754 times.
✗ Branch 3 not taken.
|
687508 | interval_nt::sign_is_non_zero(s2) && |
| 350 | interval_nt::sign_is_non_zero(s3) | ||
| 351 | ) { | ||
| 352 | |||
| 353 | 343754 | interval_nt M1 = det2x2(P2.x, P2.y, P3.x, P3.y); | |
| 354 | 343754 | interval_nt M2 = det2x2(P1.x, P1.y, P3.x, P3.y); | |
| 355 | 343754 | interval_nt M3 = det2x2(P1.x, P1.y, P2.x, P2.y); | |
| 356 | |||
| 357 | 343754 | interval_nt D = L1*P1.w*M1 | |
| 358 | 343754 | - L2*P2.w*M2 | |
| 359 | 343754 | + L3*P3.w*M3 ; | |
| 360 | |||
| 361 |
2/2✓ Branch 0 taken 307691 times.
✓ Branch 1 taken 36063 times.
|
343754 | interval_nt::Sign2 s = D.sign(); |
| 362 |
2/2✓ Branch 0 taken 307691 times.
✓ Branch 1 taken 36063 times.
|
343754 | if(interval_nt::sign_is_non_zero(s)) { |
| 363 | return Sign( | ||
| 364 |
1/2✓ Branch 1 taken 307691 times.
✗ Branch 2 not taken.
|
307691 | interval_nt::convert_sign(s) * |
| 365 |
1/2✓ Branch 1 taken 307691 times.
✗ Branch 2 not taken.
|
307691 | interval_nt::convert_sign(s1) * |
| 366 |
1/2✓ Branch 1 taken 307691 times.
✗ Branch 2 not taken.
|
307691 | interval_nt::convert_sign(s2) * |
| 367 |
1/2✓ Branch 1 taken 307691 times.
✗ Branch 2 not taken.
|
307691 | interval_nt::convert_sign(s3) |
| 368 | 307691 | ); | |
| 369 | } | ||
| 370 | } | ||
| 371 | } | ||
| 372 | |||
| 373 | // Exact | ||
| 374 | stats.log_exact(); | ||
| 375 | { | ||
| 376 | // These ones can be computed on the stack even | ||
| 377 | // under MacOSX since they are at most of length 2 | ||
| 378 | 36063 | expansion_nt L1(expansion_nt::DIFF, l0, l3); | |
| 379 | 36063 | expansion_nt L2(expansion_nt::DIFF, l1, l3); | |
| 380 | 36063 | expansion_nt L3(expansion_nt::DIFF, l2, l3); | |
| 381 | L1.optimize(); L2.optimize(); L3.optimize(); | ||
| 382 | |||
| 383 |
1/2✓ Branch 1 taken 36063 times.
✗ Branch 2 not taken.
|
36063 | vec2HE P1 = p0 - p3; |
| 384 |
1/2✓ Branch 1 taken 36063 times.
✗ Branch 2 not taken.
|
36063 | vec2HE P2 = p1 - p3; |
| 385 |
1/2✓ Branch 1 taken 36063 times.
✗ Branch 2 not taken.
|
36063 | vec2HE P3 = p2 - p3; |
| 386 |
3/6✓ Branch 1 taken 36063 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36063 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36063 times.
✗ Branch 8 not taken.
|
36063 | P1.optimize(); P2.optimize(); P3.optimize(); |
| 387 | |||
| 388 | |||
| 389 | expansion_nt M1 = det2x2(P2.x, P2.y, P3.x, P3.y); | ||
| 390 | expansion_nt M2 = det2x2(P1.x, P1.y, P3.x, P3.y); | ||
| 391 | expansion_nt M3 = det2x2(P1.x, P1.y, P2.x, P2.y); | ||
| 392 | M1.optimize(); M2.optimize(); M3.optimize(); | ||
| 393 | |||
| 394 |
3/8✓ Branch 1 taken 36063 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36063 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 36063 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
72126 | expansion_nt D = L1*P1.w*M1 |
| 395 |
4/10✓ Branch 1 taken 36063 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36063 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36063 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 36063 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
108189 | - L2*P2.w*M2 |
| 396 |
4/10✓ Branch 1 taken 36063 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36063 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36063 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 36063 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
108189 | + L3*P3.w*M3 ; |
| 397 | |||
| 398 |
3/6✓ Branch 0 taken 36063 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 36063 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 36063 times.
✗ Branch 5 not taken.
|
108189 | result = Sign(D.sign()*P1.w.sign()*P2.w.sign()*P3.w.sign()); |
| 399 | 36063 | } | |
| 400 | |||
| 401 |
2/2✓ Branch 0 taken 12735 times.
✓ Branch 1 taken 23328 times.
|
36063 | if(result != ZERO) { |
| 402 | return result; | ||
| 403 | } | ||
| 404 | |||
| 405 | stats.log_SOS(); | ||
| 406 | |||
| 407 | // Symbolic perturbation. | ||
| 408 | // | ||
| 409 | // We use the simple form of the predicate: | ||
| 410 | // | x0 y0 (x0^2+y0^2+eps^i0) 1 | | ||
| 411 | // | x1 y1 (x1^2+y1^2+eps^i1) 1 | | ||
| 412 | // | x2 y2 (x2^2+y2^2+eps^i2) 1 | | ||
| 413 | // | x3 y3 (x3^2+y3^2+eps^i3) 1 | | ||
| 414 | // where i0,i1,i2,i3 denote the indices of the points (here, they | ||
| 415 | // are local indices, coming from geometric sorting, lexico order) | ||
| 416 | // Develop along the third row (keeping only the terms in epsilon): | ||
| 417 | // | x1 y1 1 | | x0 y0 1 | | ||
| 418 | // eps^i0 | x2 y2 1 | - eps^i1 | x2 y2 1 | | ||
| 419 | // | x3 y3 1 | | x3 y3 1 | | ||
| 420 | // | ||
| 421 | // | x0 y0 1 | | x0 y0 1 | | ||
| 422 | // + eps^i2 | x1 y1 1 | - eps^i3 | x1 y1 1 | | ||
| 423 | // | x3 y3 1 | | x2 y2 1 | | ||
| 424 | 12735 | return SOS( | |
| 425 | vec2HgLexicoCompare<expansion_nt>(), | ||
| 426 | 5656 | p0, SOS_result( det3_111_sign(p1,p2,p3)), | |
| 427 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 2242 times.
|
2242 | p1, SOS_result(-det3_111_sign(p0,p2,p3)), |
| 428 | 3407 | p2, SOS_result( det3_111_sign(p0,p1,p3)), | |
| 429 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1430 times.
|
1430 | p3, SOS_result(-det3_111_sign(p0,p1,p2)) |
| 430 | 12735 | ); | |
| 431 | } | ||
| 432 | |||
| 433 | 724351 | coord_index_t triangle_normal_axis( | |
| 434 | const vec3& p1, const vec3& p2, const vec3& p3 | ||
| 435 | ) { | ||
| 436 |
4/4✓ Branch 0 taken 37 times.
✓ Branch 1 taken 724314 times.
✓ Branch 3 taken 35 times.
✓ Branch 4 taken 2 times.
|
724351 | static PredicateStats stats("triangle_normal_axis"); |
| 437 | stats.log_invoke(); | ||
| 438 | |||
| 439 | // Filter using interval arithmetics | ||
| 440 | { | ||
| 441 | interval_nt::Rounding rounding; | ||
| 442 | vec3I p1I(p1); | ||
| 443 | 724351 | vec3I U = vec3I(p2) - p1I; | |
| 444 | 724351 | vec3I V = vec3I(p3) - p1I; | |
| 445 | 724351 | vec3I N = cross(U,V); | |
| 446 | 724351 | interval_nt::Sign2 sx = N.x.sign(); | |
| 447 | 724351 | interval_nt::Sign2 sy = N.y.sign(); | |
| 448 |
2/2✓ Branch 0 taken 197019 times.
✓ Branch 1 taken 527332 times.
|
724351 | interval_nt::Sign2 sz = N.z.sign(); |
| 449 | if( | ||
| 450 | !interval_nt::sign_is_determined(sx) || | ||
| 451 | !interval_nt::sign_is_determined(sy) || | ||
| 452 | !interval_nt::sign_is_determined(sz) | ||
| 453 | ) { | ||
| 454 | 577082 | goto exact; // Yes, goto, why not ? | |
| 455 | } | ||
| 456 | |||
| 457 |
2/2✓ Branch 1 taken 74848 times.
✓ Branch 2 taken 72421 times.
|
147269 | if(interval_nt::convert_sign(sx) != POSITIVE) { |
| 458 | N.x.negate(); | ||
| 459 | } | ||
| 460 |
2/2✓ Branch 1 taken 72309 times.
✓ Branch 2 taken 74960 times.
|
147269 | if(interval_nt::convert_sign(sy) != POSITIVE) { |
| 461 | N.y.negate(); | ||
| 462 | } | ||
| 463 |
2/2✓ Branch 1 taken 65653 times.
✓ Branch 2 taken 81616 times.
|
147269 | if(interval_nt::convert_sign(sz) != POSITIVE) { |
| 464 | N.z.negate(); | ||
| 465 | } | ||
| 466 | 147269 | interval_nt::Sign2 sxy = (N.x - N.y).sign(); | |
| 467 |
2/2✓ Branch 0 taken 147245 times.
✓ Branch 1 taken 24 times.
|
147269 | interval_nt::Sign2 sxz = (N.x - N.z).sign(); |
| 468 | if( | ||
| 469 | !interval_nt::sign_is_determined(sxy) || | ||
| 470 | !interval_nt::sign_is_determined(sxz) | ||
| 471 | ) { | ||
| 472 | 146 | goto exact; // Ahaha, another one !!! | |
| 473 | } | ||
| 474 | if( | ||
| 475 |
4/4✓ Branch 1 taken 73644 times.
✓ Branch 2 taken 73479 times.
✓ Branch 3 taken 23719 times.
✓ Branch 4 taken 49925 times.
|
220767 | interval_nt::convert_sign(sxy) >= 0 && |
| 476 | 73644 | interval_nt::convert_sign(sxz) >= 0 | |
| 477 | ) { | ||
| 478 | return 0; | ||
| 479 | } | ||
| 480 |
1/2✓ Branch 0 taken 97198 times.
✗ Branch 1 not taken.
|
97198 | interval_nt::Sign2 syz = (N.y - N.z).sign(); |
| 481 | if(!interval_nt::sign_is_determined(syz)) { | ||
| 482 | ✗ | goto exact; // The last one (for now !) | |
| 483 | } | ||
| 484 |
2/2✓ Branch 1 taken 50208 times.
✓ Branch 2 taken 46990 times.
|
97198 | if(interval_nt::convert_sign(syz) >=0 ) { |
| 485 | 50208 | return 1; | |
| 486 | } | ||
| 487 | return 2; | ||
| 488 | } | ||
| 489 | |||
| 490 | // Exact computation, using low-level expansion API | ||
| 491 | // (expansions allocated on the stack, better for | ||
| 492 | // multithreading) | ||
| 493 | exact: | ||
| 494 | stats.log_exact(); | ||
| 495 | |||
| 496 | // These ones can be computed on the stack even | ||
| 497 | // under MacOSX since they are at most of length 2 | ||
| 498 | |||
| 499 | 577228 | const expansion& Ux = expansion_diff(p2.x, p1.x); | |
| 500 | 577228 | const expansion& Uy = expansion_diff(p2.y, p1.y); | |
| 501 | 577228 | const expansion& Uz = expansion_diff(p2.z, p1.z); | |
| 502 | |||
| 503 | 577228 | const expansion& Vx = expansion_diff(p3.x, p1.x); | |
| 504 | 577228 | const expansion& Vy = expansion_diff(p3.y, p1.y); | |
| 505 | 577228 | const expansion& Vz = expansion_diff(p3.z, p1.z); | |
| 506 | |||
| 507 | 577228 | expansion& Nx = expansion_det2x2(Uy,Vy,Uz,Vz); | |
| 508 | 577228 | expansion& Ny = expansion_det2x2(Uz,Vz,Ux,Vx); | |
| 509 | 577228 | expansion& Nz = expansion_det2x2(Ux,Vx,Uy,Vy); | |
| 510 | |||
| 511 |
2/2✓ Branch 0 taken 552600 times.
✓ Branch 1 taken 24628 times.
|
577228 | if(Nx.sign() != POSITIVE) { |
| 512 | Nx.negate(); | ||
| 513 | } | ||
| 514 | |||
| 515 |
2/2✓ Branch 0 taken 550599 times.
✓ Branch 1 taken 26629 times.
|
577228 | if(Ny.sign() != POSITIVE) { |
| 516 | Ny.negate(); | ||
| 517 | } | ||
| 518 | |||
| 519 |
2/2✓ Branch 0 taken 422643 times.
✓ Branch 1 taken 154585 times.
|
577228 | if(Nz.sign() != POSITIVE) { |
| 520 | Nz.negate(); | ||
| 521 | } | ||
| 522 | |||
| 523 |
4/4✓ Branch 1 taken 535494 times.
✓ Branch 2 taken 41734 times.
✓ Branch 4 taken 494203 times.
✓ Branch 5 taken 41291 times.
|
577228 | if(Nx.compare(Ny) >= 0 && Nx.compare(Nz) >= 0) { |
| 524 | geo_debug_assert(Nx.sign() != ZERO); | ||
| 525 | return 0; | ||
| 526 | } | ||
| 527 |
2/2✓ Branch 1 taken 497968 times.
✓ Branch 2 taken 37969 times.
|
535937 | if(Ny.compare(Nz) >= 0) { |
| 528 | geo_debug_assert(Ny.sign() != ZERO); | ||
| 529 | return 1; | ||
| 530 | } | ||
| 531 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 497968 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
497968 | geo_assert(Nz.sign() != ZERO); |
| 532 | return 2; | ||
| 533 | } | ||
| 534 | |||
| 535 | ✗ | bool aligned_3d( | |
| 536 | const vec3HE& p0, const vec3HE& p1, const vec3HE& p2 | ||
| 537 | ) { | ||
| 538 | // TODO: filter if need be | ||
| 539 | ✗ | vec3HE U = p1-p0; | |
| 540 | ✗ | vec3HE V = p2-p0; | |
| 541 | return ( | ||
| 542 | ✗ | det2x2(U.x,V.x,U.y,V.y).sign() == ZERO && | |
| 543 | ✗ | det2x2(U.y,V.y,U.z,V.z).sign() == ZERO && | |
| 544 | ✗ | det2x2(U.z,V.z,U.x,V.x).sign() == ZERO | |
| 545 | ✗ | ); | |
| 546 | ✗ | } | |
| 547 | |||
| 548 | 246918 | bool on_segment_3d( | |
| 549 | const vec3HE& p, const vec3HE& q1, const vec3HE& q2 | ||
| 550 | ) { | ||
| 551 | // TODO: filter if need be | ||
| 552 | 246918 | vec3HE U = p-q1; | |
| 553 |
1/2✓ Branch 1 taken 246918 times.
✗ Branch 2 not taken.
|
246918 | vec3HE V = p-q2; |
| 554 | 246918 | if ( | |
| 555 |
4/4✓ Branch 0 taken 11583 times.
✓ Branch 1 taken 8379 times.
✓ Branch 2 taken 8439 times.
✓ Branch 3 taken 238479 times.
|
266880 | det2x2(U.x,V.x,U.y,V.y).sign() != ZERO || |
| 556 |
5/8✓ Branch 0 taken 19962 times.
✓ Branch 1 taken 226956 times.
✓ Branch 2 taken 8439 times.
✓ Branch 3 taken 3144 times.
✓ Branch 4 taken 246918 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
505419 | det2x2(U.y,V.y,U.z,V.z).sign() != ZERO || |
| 557 |
2/4✓ Branch 0 taken 19962 times.
✓ Branch 1 taken 226956 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
246918 | det2x2(U.z,V.z,U.x,V.x).sign() != ZERO |
| 558 | ) { | ||
| 559 | return false; | ||
| 560 | } | ||
| 561 | |||
| 562 | return ( | ||
| 563 | ( | ||
| 564 |
8/20✓ Branch 1 taken 8439 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8439 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8439 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 8439 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 8439 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 8439 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 8439 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 8439 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
42195 | (U.x*V.x + U.y*V.y + U.z*V.z).sign() * |
| 565 | 8439 | U.w.sign() * V.w.sign() | |
| 566 | ) <= ZERO | ||
| 567 |
1/2✓ Branch 0 taken 8439 times.
✗ Branch 1 not taken.
|
8439 | ); |
| 568 | 246918 | } | |
| 569 | |||
| 570 | |||
| 571 | 24896 | vec3 approximate(const vec3HE& p) { | |
| 572 | // TODO: find a way of computing the round to nearest approxomation. | ||
| 573 | // see division operation for expansions, | ||
| 574 | // here: | ||
| 575 | // https://www.jucs.org/jucs_5_6/division_of_floating_point/Daumas_M.pdf | ||
| 576 | double w = p.w.estimate(); | ||
| 577 | 49792 | return vec3(p.x.estimate()/w, p.y.estimate()/w, p.z.estimate()/w); | |
| 578 | } | ||
| 579 | |||
| 580 | 81 | vec2 approximate(const vec2HE& p) { | |
| 581 | // TODO: find a way of computing the round to nearest approxomation. | ||
| 582 | // see division operation for expansions, | ||
| 583 | // here: | ||
| 584 | // https://www.jucs.org/jucs_5_6/division_of_floating_point/Daumas_M.pdf | ||
| 585 | double w = p.w.estimate(); | ||
| 586 | 81 | return vec2(p.x.estimate()/w, p.y.estimate()/w); | |
| 587 | } | ||
| 588 | |||
| 589 | } | ||
| 590 | |||
| 591 | /*****************************************************************/ | ||
| 592 | |||
| 593 | // Under Linux we got 10 Mb of stack (!) Then some operations can be | ||
| 594 | // made faster by using the low-level expansion API (that allocates | ||
| 595 | // intermediary multiprecision values on stack rather than in the heap). | ||
| 596 | // These optimized functions are written as template specializations | ||
| 597 | // (used automatically). | ||
| 598 | |||
| 599 | #ifdef GEO_HAS_BIG_STACK | ||
| 600 | |||
| 601 | 310486 | template<> expansion_nt det(const vec2E& v1, const vec2E& v2) { | |
| 602 | expansion* result = expansion::new_expansion_on_heap( | ||
| 603 | expansion::det2x2_capacity( | ||
| 604 | v1.x.rep(), v1.y.rep(), | ||
| 605 | v2.x.rep(), v2.y.rep() | ||
| 606 | ) | ||
| 607 | ); | ||
| 608 | 310486 | result->assign_det2x2( | |
| 609 | v1.x.rep(), v1.y.rep(), | ||
| 610 | v2.x.rep(), v2.y.rep() | ||
| 611 | ); | ||
| 612 | 310486 | return expansion_nt(result); | |
| 613 | } | ||
| 614 | |||
| 615 | ✗ | template<> expansion_nt dot(const vec2E& v1, const vec2E& v2) { | |
| 616 | ✗ | const expansion& m1 = expansion_product(v1.x.rep(), v2.x.rep()); | |
| 617 | ✗ | const expansion& m2 = expansion_product(v1.y.rep(), v2.y.rep()); | |
| 618 | ✗ | return expansion_nt(expansion_nt::SUM, m1, m2); | |
| 619 | } | ||
| 620 | |||
| 621 | 180328 | template<> expansion_nt dot(const vec3E& v1, const vec3E& v2) { | |
| 622 | 180328 | const expansion& m1 = expansion_product(v1.x.rep(), v2.x.rep()); | |
| 623 | 180328 | const expansion& m2 = expansion_product(v1.y.rep(), v2.y.rep()); | |
| 624 | 180328 | const expansion& m3 = expansion_product(v1.z.rep(), v2.z.rep()); | |
| 625 | 180328 | return expansion_nt(expansion_nt::SUM,m1,m2,m3); | |
| 626 | } | ||
| 627 | |||
| 628 | /*********************************************/ | ||
| 629 | |||
| 630 | 222801 | template<> vec3Hg<expansion_nt> mix( | |
| 631 | const rationalg<expansion_nt>& t, const vec3& p1, const vec3& p2 | ||
| 632 | ) { | ||
| 633 | expansion& st_d = const_cast<expansion&>(t.denom().rep()); | ||
| 634 | 222801 | st_d.optimize(); | |
| 635 | expansion& t_n = const_cast<expansion&>(t.num().rep()); | ||
| 636 | 222801 | t_n.optimize(); | |
| 637 | 222801 | const expansion& s_n = expansion_diff(st_d, t_n); | |
| 638 | 222801 | const expansion& sx = expansion_product(s_n, p1.x); | |
| 639 | 222801 | const expansion& tx = expansion_product(t_n, p2.x); | |
| 640 | 222801 | const expansion& sy = expansion_product(s_n, p1.y); | |
| 641 | 222801 | const expansion& ty = expansion_product(t_n, p2.y); | |
| 642 | 222801 | const expansion& sz = expansion_product(s_n, p1.z); | |
| 643 | 222801 | const expansion& tz = expansion_product(t_n, p2.z); | |
| 644 | return vec3HE( | ||
| 645 |
2/6✓ Branch 1 taken 222801 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 222801 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
445602 | expansion_nt(expansion_nt::SUM, sx,tx), |
| 646 |
2/6✓ Branch 1 taken 222801 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 222801 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
445602 | expansion_nt(expansion_nt::SUM, sy,ty), |
| 647 |
2/6✓ Branch 1 taken 222801 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 222801 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
445602 | expansion_nt(expansion_nt::SUM, sz,tz), |
| 648 | 222801 | expansion_nt(st_d) | |
| 649 | 445602 | ); | |
| 650 | } | ||
| 651 | |||
| 652 | ✗ | template<> vec2Hg<expansion_nt> mix( | |
| 653 | const rationalg<expansion_nt>& t, const vec2& p1, const vec2& p2 | ||
| 654 | ) { | ||
| 655 | expansion& st_d = const_cast<expansion&>(t.denom().rep()); | ||
| 656 | ✗ | st_d.optimize(); | |
| 657 | expansion& t_n = const_cast<expansion&>(t.num().rep()); | ||
| 658 | ✗ | t_n.optimize(); | |
| 659 | ✗ | const expansion& s_n = expansion_diff(st_d, t_n); | |
| 660 | ✗ | const expansion& sx = expansion_product(s_n, p1.x); | |
| 661 | ✗ | const expansion& tx = expansion_product(t_n, p2.x); | |
| 662 | ✗ | const expansion& sy = expansion_product(s_n, p1.y); | |
| 663 | ✗ | const expansion& ty = expansion_product(t_n, p2.y); | |
| 664 | return vec2HE( | ||
| 665 | ✗ | expansion_nt(expansion_nt::SUM, sx,tx), | |
| 666 | ✗ | expansion_nt(expansion_nt::SUM, sy,ty), | |
| 667 | ✗ | expansion_nt(st_d) | |
| 668 | ✗ | ); | |
| 669 | } | ||
| 670 | |||
| 671 | /*********************************************/ | ||
| 672 | |||
| 673 | 17052 | template <> vec3E triangle_normal<vec3E>( | |
| 674 | const vec3& p1, const vec3& p2, const vec3& p3 | ||
| 675 | ) { | ||
| 676 | 17052 | const expansion& Ux = expansion_diff(p2.x,p1.x); | |
| 677 | 17052 | const expansion& Uy = expansion_diff(p2.y,p1.y); | |
| 678 | 17052 | const expansion& Uz = expansion_diff(p2.z,p1.z); | |
| 679 | 17052 | const expansion& Vx = expansion_diff(p3.x,p1.x); | |
| 680 | 17052 | const expansion& Vy = expansion_diff(p3.y,p1.y); | |
| 681 | 17052 | const expansion& Vz = expansion_diff(p3.z,p1.z); | |
| 682 | expansion* Nx = expansion::new_expansion_on_heap( | ||
| 683 | expansion::det2x2_capacity(Uy,Uz,Vy,Vz) | ||
| 684 | ); | ||
| 685 | 17052 | Nx->assign_det2x2(Uy,Uz,Vy,Vz); | |
| 686 | expansion* Ny = expansion::new_expansion_on_heap( | ||
| 687 | expansion::det2x2_capacity(Uz,Ux,Vz,Vx) | ||
| 688 | ); | ||
| 689 | 17052 | Ny->assign_det2x2(Uz,Ux,Vz,Vx); | |
| 690 | expansion* Nz = expansion::new_expansion_on_heap( | ||
| 691 | expansion::det2x2_capacity(Ux,Uy,Vx,Vy) | ||
| 692 | ); | ||
| 693 | 17052 | Nz->assign_det2x2(Ux,Uy,Vx,Vy); | |
| 694 | 17052 | return vec3E(expansion_nt(Nx), expansion_nt(Ny), expansion_nt(Nz)); | |
| 695 | } | ||
| 696 | #endif | ||
| 697 | } | ||
| 698 | |||
| 699 | /**************************************************************************/ | ||
| 700 | |||
| 701 | #ifndef GEOGRAM_PSM | ||
| 702 | namespace GEO { | ||
| 703 | namespace PCK { | ||
| 704 | ✗ | Sign orient_3d_SOS( | |
| 705 | const exact::vec3h& p0, const exact::vec3h& p1, | ||
| 706 | const exact::vec3h& p2, const exact::vec3h& p3 | ||
| 707 | ) { | ||
| 708 | |||
| 709 | struct SOS { | ||
| 710 | enum {W=3}; | ||
| 711 | |||
| 712 | ✗ | SOS( | |
| 713 | const exact::vec3h& p0, const exact::vec3h& p1, | ||
| 714 | const exact::vec3h& p2, const exact::vec3h& p3 | ||
| 715 | ✗ | ) : p_orig{&p0, &p1, &p2, &p3}, p_sort{&p0, &p1, &p2, &p3} { | |
| 716 | ✗ | std::sort( | |
| 717 | ✗ | p_sort, p_sort+4, | |
| 718 | []( | ||
| 719 | const exact::vec3h* pp1, const exact::vec3h* pp2 | ||
| 720 | )->bool { | ||
| 721 | vec3HgLexicoCompare<exact::scalar> cmp; | ||
| 722 | ✗ | return cmp(*pp1,*pp2); | |
| 723 | } | ||
| 724 | ); | ||
| 725 | ✗ | parity = Permutation::permutation_is_odd(p_orig, p_sort, 4) | |
| 726 | ✗ | ? NEGATIVE : POSITIVE; | |
| 727 | ✗ | } | |
| 728 | |||
| 729 | ✗ | Sign orient_1d(index_t i, index_t j, index_t axis) const { | |
| 730 | coord_index_t ax = coord_index_t(axis); | ||
| 731 | Sign s = geo_cmp( | ||
| 732 | ✗ | exact::rational((*p_sort[i])[ax], (*p_sort[i])[W]), | |
| 733 | ✗ | exact::rational((*p_sort[j])[ax], (*p_sort[j])[W]) | |
| 734 | ); | ||
| 735 | ✗ | return Sign(s*parity); | |
| 736 | } | ||
| 737 | |||
| 738 | ✗ | Sign orient_2d( | |
| 739 | index_t i, index_t j, index_t k, index_t axis1, index_t axis2 | ||
| 740 | ) const { | ||
| 741 | coord_index_t ax1 = coord_index_t(axis1); | ||
| 742 | coord_index_t ax2 = coord_index_t(axis2); | ||
| 743 | exact::vec2h pi{ | ||
| 744 | ✗ | (*p_sort[i])[ax1], (*p_sort[i])[ax2], (*p_sort[i])[W] | |
| 745 | ✗ | }; | |
| 746 | exact::vec2h pj{ | ||
| 747 | ✗ | (*p_sort[j])[ax1], (*p_sort[j])[ax2], (*p_sort[j])[W] | |
| 748 | ✗ | }; | |
| 749 | exact::vec2h pk{ | ||
| 750 | ✗ | (*p_sort[k])[ax1], (*p_sort[k])[ax2], (*p_sort[k])[W] | |
| 751 | ✗ | }; | |
| 752 | ✗ | Sign s = PCK::orient_2d(pi,pj,pk); | |
| 753 | ✗ | return Sign(s*parity); | |
| 754 | ✗ | } | |
| 755 | |||
| 756 | const exact::vec3h* p_orig[4]; | ||
| 757 | const exact::vec3h* p_sort[4]; | ||
| 758 | Sign parity; | ||
| 759 | }; | ||
| 760 | |||
| 761 | ✗ | return orient_3d_SOS_impl<exact::vec3h, SOS>(p0,p1,p2,p3); | |
| 762 | } | ||
| 763 | } | ||
| 764 | } | ||
| 765 | #endif | ||
| 766 |