| 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 | #ifndef GEOGRAM_NUMERICS_PREDICATES | ||
| 41 | #define GEOGRAM_NUMERICS_PREDICATES | ||
| 42 | |||
| 43 | #include <geogram/basic/common.h> | ||
| 44 | #include <geogram/basic/numeric.h> | ||
| 45 | #include <geogram/basic/geometry.h> | ||
| 46 | #include <geogram/numerics/PCK.h> | ||
| 47 | |||
| 48 | /** | ||
| 49 | * \file geogram/numerics/predicates.h | ||
| 50 | * \brief Filtered exact predicates for restricted Voronoi diagrams. | ||
| 51 | */ | ||
| 52 | |||
| 53 | |||
| 54 | namespace GEO { | ||
| 55 | |||
| 56 | /** | ||
| 57 | * \brief PCK (Predicate Construction Kit) implements a set of | ||
| 58 | * geometric predicates. PCK uses arithmetic filters (Meyer and Pion), | ||
| 59 | * expansion arithmetics (Shewchuk) and simulation of simplicity | ||
| 60 | * (Edelsbrunner). | ||
| 61 | */ | ||
| 62 | namespace PCK { | ||
| 63 | |||
| 64 | /** | ||
| 65 | * \brief Mode for symbolic perturbations. | ||
| 66 | */ | ||
| 67 | enum SOSMode { SOS_ADDRESS, SOS_LEXICO }; | ||
| 68 | |||
| 69 | /** | ||
| 70 | * \brief Sets the current mode for handling symbolic perturbations | ||
| 71 | * (SOS for Simulation Of Simplicity). | ||
| 72 | * \param[in] m one of SOS_ADDRESS, SOS_LEXICO | ||
| 73 | * \details If SOS_ADDRESS mode is used, then points are supposed | ||
| 74 | * to be allocated in a fixed array, and the same point always | ||
| 75 | * designated by the same address. If SOS_LEXICO is used then points | ||
| 76 | * are sorted in lexicographic order for computing the symbolic | ||
| 77 | * perturbation. SOS_LEXICO works for points that are generated | ||
| 78 | * dynamically (with no fixed address). | ||
| 79 | */ | ||
| 80 | void GEOGRAM_API set_SOS_mode(SOSMode m); | ||
| 81 | |||
| 82 | /** | ||
| 83 | * \brief Gets the current mode for handling symbolic perturbations. | ||
| 84 | * \return one of SOS_ADDRESS, SOS_LEXICO | ||
| 85 | * \see set_SOS_mode() | ||
| 86 | */ | ||
| 87 | SOSMode GEOGRAM_API get_SOS_mode(); | ||
| 88 | |||
| 89 | /** | ||
| 90 | * \brief Computes the side of a point (given directly) | ||
| 91 | * relative to a bisector. | ||
| 92 | * \details Computes the side of \f$ q0 \f$ relative to | ||
| 93 | * \f$ \Pi(p0,p1) \f$. | ||
| 94 | * Symbolic perturbation is applied whenever equality holds. | ||
| 95 | * \param[in] p0 , p1 extremities of the bisector | ||
| 96 | * \param[in] q0 point to be tested | ||
| 97 | * \param[in] DIM number of coordinates of the point | ||
| 98 | * \retval POSITIVE if d(p0,q0) < d(p1,q0) | ||
| 99 | * \retval NEGATIVE if d(p0,q0) > d(p1,q1) | ||
| 100 | * \retval perturb() if f(p0,q0) = d(p1,q1), | ||
| 101 | * where \c perturb() denotes a globally | ||
| 102 | * consistent perturbation, that returns either POSITIVE or NEGATIVE | ||
| 103 | * \note Only some specific dimensions are implemented (3,4,6 and 7) | ||
| 104 | */ | ||
| 105 | Sign GEOGRAM_API side1_SOS( | ||
| 106 | const double* p0, const double* p1, | ||
| 107 | const double* q0, | ||
| 108 | coord_index_t DIM | ||
| 109 | ); | ||
| 110 | |||
| 111 | /** | ||
| 112 | * \brief Computes the side of a point (given as the intersection | ||
| 113 | * between a segment and a bisector) relative to another bisector. | ||
| 114 | * \details Computes the side of \f$ q = \Pi(p0,p1) \cap [q0,q1] \f$ | ||
| 115 | * relative to \f$ \Pi(p0,p2) \f$. | ||
| 116 | * Symbolic perturbation is applied whenever equality holds. | ||
| 117 | * \param[in] p0 first extremity of the bisectors | ||
| 118 | * \param[in] p1 second extremity of the first bisector | ||
| 119 | * (that defines the intersection q) | ||
| 120 | * \param[in] p2 second extremity of the second bisector | ||
| 121 | * (against which orientation is tested) | ||
| 122 | * \param[in] q0 , q1 extremities of the segment | ||
| 123 | * (that defines the intersection q) | ||
| 124 | * \retval POSITIVE if d(p0,q) < d(p2,q) | ||
| 125 | * \retval NEGATIVE if d(p0,q) > d(p2,q) | ||
| 126 | * \retval perturb() if d(p0,q) = d(p2,q), | ||
| 127 | * where \c perturb() denotes a globally | ||
| 128 | * consistent perturbation, that returns either POSITIVE or NEGATIVE | ||
| 129 | * \note Only some specific dimensions are implemented (3,4,6 and 7) | ||
| 130 | */ | ||
| 131 | Sign GEOGRAM_API side2_SOS( | ||
| 132 | const double* p0, const double* p1, const double* p2, | ||
| 133 | const double* q0, const double* q1, | ||
| 134 | coord_index_t DIM | ||
| 135 | ); | ||
| 136 | |||
| 137 | /** | ||
| 138 | * \brief Computes the side of a point (given as the intersection | ||
| 139 | * between a facet and two bisectors) relative to another bisector. | ||
| 140 | * \details Computes the side of | ||
| 141 | * \f$ q = \Pi(p0,p1) \cap Pi(p0,p2) \cap \Delta[q0,q1,q2] \f$ | ||
| 142 | * relative to \f$ \Pi(p0,p3) \f$. | ||
| 143 | * Symbolic perturbation is applied whenever equality holds. | ||
| 144 | * \param[in] p0 first extremity of the bisectors | ||
| 145 | * \param[in] p1 second extremity of the first bisector | ||
| 146 | * (that defines the intersection q) | ||
| 147 | * \param[in] p2 second extremity of the second bisector | ||
| 148 | * (that defines the intersection q) | ||
| 149 | * \param[in] p3 second extremity of the third bisector | ||
| 150 | * (against which orientation is tested) | ||
| 151 | * \param[in] q0 , q1 , q2 vertices of the triangle | ||
| 152 | * (that defines the intersection q) | ||
| 153 | * \retval POSITIVE if d(p0,q) < d(p3,q) | ||
| 154 | * \retval NEGATIVE if d(p0,q) > d(p3,q) | ||
| 155 | * \retval perturb() if d(p0,q) = d(p3,q), | ||
| 156 | * where \c perturb() denotes a globally | ||
| 157 | * consistent perturbation, that returns either POSITIVE or NEGATIVE | ||
| 158 | * \note Only some specific dimensions are implemented (3,4,6 and 7) | ||
| 159 | */ | ||
| 160 | Sign GEOGRAM_API side3_SOS( | ||
| 161 | const double* p0, const double* p1, | ||
| 162 | const double* p2, const double* p3, | ||
| 163 | const double* q0, const double* q1, const double* q2, | ||
| 164 | coord_index_t DIM | ||
| 165 | ); | ||
| 166 | |||
| 167 | /** | ||
| 168 | * \brief Computes the side of a point (given as the intersection | ||
| 169 | * between a facet and two bisectors) relative to another bisector. | ||
| 170 | * \details Computes the side of | ||
| 171 | * \f$ q = \Pi(p0 h0,p1 h1) \cap Pi(p0 h0,p2 h2) \cap \Delta[q0, q1, q2] \f$ | ||
| 172 | * relative to \f$ \Pi(p0 hp0,p3 hp3) \f$. | ||
| 173 | * Symbolic perturbation is applied whenever equality holds. | ||
| 174 | * \param[in] p0 first extremity of the bisectors | ||
| 175 | * \param[in] p1 second extremity of the first bisector | ||
| 176 | * (that defines the intersection q) | ||
| 177 | * \param[in] p2 second extremity of the second bisector | ||
| 178 | * (that defines the intersection q) | ||
| 179 | * \param[in] p3 second extremity of the third bisector | ||
| 180 | * (against which orientation is tested) | ||
| 181 | * \param h0 , h1 , h2 , h3 lifted coordinates of \p p0, \p p1, \p p2 | ||
| 182 | * and \p p3 | ||
| 183 | * \param[in] q0 , q1 , q2 vertices of the triangle | ||
| 184 | * (that defines the intersection q) | ||
| 185 | * \param[in] SOS if true, do the symbolic perturbation in the | ||
| 186 | * degenerate case | ||
| 187 | * \retval POSITIVE if d(p0 hp0,q) < d(p3 hp3, q) | ||
| 188 | * \retval NEGATIVE if d(p0 hp0,q) > d(p3 hp3, q) | ||
| 189 | * \retval perturb() if d(p0 hp0,q) = d(p3 hp3, q), | ||
| 190 | * where \c perturb() denotes a globally | ||
| 191 | * consistent perturbation, that returns either POSITIVE or NEGATIVE | ||
| 192 | */ | ||
| 193 | Sign GEOGRAM_API side3_3dlifted_SOS( | ||
| 194 | const double* p0, const double* p1, | ||
| 195 | const double* p2, const double* p3, | ||
| 196 | double h0, double h1, double h2, double h3, | ||
| 197 | const double* q0, const double* q1, const double* q2, | ||
| 198 | bool SOS=true | ||
| 199 | ); | ||
| 200 | |||
| 201 | /** | ||
| 202 | * \brief Computes the side of a point (given as the intersection | ||
| 203 | * between a tetrahedron and three bisectors) relative to | ||
| 204 | * another bisector. | ||
| 205 | * \details Computes the side of | ||
| 206 | * \f$ q = \Pi(p0,p1) \cap Pi(p0,p2) \cap Pi(p0,p3) | ||
| 207 | * \cap \Delta[q0,q1,q2,q3] \f$ relative to \f$ \Pi(p0,p4) \f$. | ||
| 208 | * Symbolic perturbation is applied whenever equality holds. | ||
| 209 | * \param[in] p0 first extremity of the bisectors | ||
| 210 | * \param[in] p1 second extremity of the first bisector | ||
| 211 | * (that defines the intersection q) | ||
| 212 | * \param[in] p2 second extremity of the second bisector | ||
| 213 | * (that defines the intersection q) | ||
| 214 | * \param[in] p3 second extremity of the third bisector | ||
| 215 | * (that defines the intersection q) | ||
| 216 | * \param[in] p4 second extremity of the fourth bisector | ||
| 217 | * (against which orientation is tested) | ||
| 218 | * \param[in] q0 , q1 , q2 , q3 vertices of the tetrahedron | ||
| 219 | * (that defines the intersection q) | ||
| 220 | * (that defines the intersection q) | ||
| 221 | * \retval POSITIVE if d(p0,q) < d(p4,q) | ||
| 222 | * \retval NEGATIVE if d(p0,q) > d(p4,q) | ||
| 223 | * \retval perturb() if d(p0,q) = d(p4,q), | ||
| 224 | * where \c perturb() denotes a globally | ||
| 225 | * consistent perturbation, that returns either POSITIVE or NEGATIVE | ||
| 226 | * \note Only some specific dimensions are implemented (3,4,6 and 7) | ||
| 227 | */ | ||
| 228 | Sign GEOGRAM_API side4_SOS( | ||
| 229 | const double* p0, | ||
| 230 | const double* p1, const double* p2, | ||
| 231 | const double* p3, const double* p4, | ||
| 232 | const double* q0, const double* q1, | ||
| 233 | const double* q2, const double* q3, | ||
| 234 | coord_index_t DIM | ||
| 235 | ); | ||
| 236 | |||
| 237 | |||
| 238 | /** | ||
| 239 | * \brief Computes the side of a point (given as the intersection | ||
| 240 | * between three bisectors) relative to another bisector. | ||
| 241 | * \details Computes the side of | ||
| 242 | * \f$ q = \Pi(p0,p1) \cap \Pi(p0,p2) \cap \Pi(p0,p3) \f$ | ||
| 243 | * relative to \f$ Pi(p0,p4) \f$. | ||
| 244 | * This version does not apply symbolic perturbation when equality | ||
| 245 | * holds. | ||
| 246 | * side4_3d() is a special case of side4(), where the ambient and | ||
| 247 | * intrinsic dimensions coincide (therefore no embedding tetrahedron | ||
| 248 | * is needed). | ||
| 249 | * \param[in] p0 first extremity of the bisectors | ||
| 250 | * \param[in] p1 second extremity of the first bisector | ||
| 251 | * (that defines the intersection q) | ||
| 252 | * \param[in] p2 second extremity of the second bisector | ||
| 253 | * (that defines the intersection q) | ||
| 254 | * \param[in] p3 second extremity of the third bisector | ||
| 255 | * (that defines the intersection q) | ||
| 256 | * \param[in] p4 second extremity of the fourth bisector | ||
| 257 | * (against which orientation is tested) | ||
| 258 | * \retval POSITIVE if d(p0,q) < d(p4,q) | ||
| 259 | * \retval NEGATIVE if d(p0,q) > d(p4,q) | ||
| 260 | * \retval ZERO if d(p0,q) = d(p4,q), | ||
| 261 | */ | ||
| 262 | Sign GEOGRAM_API side4_3d( | ||
| 263 | const double* p0, | ||
| 264 | const double* p1, const double* p2, | ||
| 265 | const double* p3, const double* p4 | ||
| 266 | ); | ||
| 267 | |||
| 268 | /** | ||
| 269 | * \brief Computes the side of a point (given as the intersection | ||
| 270 | * between three bisectors) relative to another bisector. | ||
| 271 | * \details Computes the side of | ||
| 272 | * \f$ q = \Pi(p0,p1) \cap \Pi(p0,p2) \cap \Pi(p0,p3) \f$ | ||
| 273 | * relative to \f$ Pi(p0,p4) \f$. | ||
| 274 | * Symbolic perturbation is applied whenever equality holds. | ||
| 275 | * side4_3d() is a special case of side4(), where the ambient and | ||
| 276 | * intrinsic dimensions coincide (therefore no embedding tetrahedron | ||
| 277 | * is needed). | ||
| 278 | * \param[in] p0 first extremity of the bisectors | ||
| 279 | * \param[in] p1 second extremity of the first bisector | ||
| 280 | * (that defines the intersection q) | ||
| 281 | * \param[in] p2 second extremity of the second bisector | ||
| 282 | * (that defines the intersection q) | ||
| 283 | * \param[in] p3 second extremity of the third bisector | ||
| 284 | * (that defines the intersection q) | ||
| 285 | * \param[in] p4 second extremity of the fourth bisector | ||
| 286 | * (against which orientation is tested) | ||
| 287 | * \retval POSITIVE if d(p0,q) < d(p4,q) | ||
| 288 | * \retval NEGATIVE if d(p0,q) > d(p4,q) | ||
| 289 | * \retval perturb() if d(p0,q) = d(p4,q), | ||
| 290 | * where \c perturb() denotes a globally | ||
| 291 | * consistent perturbation, that returns either POSITIVE or NEGATIVE | ||
| 292 | */ | ||
| 293 | Sign GEOGRAM_API side4_3d_SOS( | ||
| 294 | const double* p0, const double* p1, | ||
| 295 | const double* p2, const double* p3, const double* p4 | ||
| 296 | ); | ||
| 297 | |||
| 298 | /** | ||
| 299 | * \brief Tests whether a 3d point is inside the circumscribed | ||
| 300 | * sphere of a 3d tetrahedron. | ||
| 301 | * \param[in] p0 , p1 , p2 , p3 the four vertices of the tetrahedron | ||
| 302 | * \param[in] p4 the point | ||
| 303 | * \retval POSITIVE whenever \p p4 is inside the circumscribed sphere | ||
| 304 | * of the tetrahedron \p p0, \p p1, \p p2, \p p3 | ||
| 305 | * \retval NEGATIVE whenever \p p4 is outside the circumscribed sphere | ||
| 306 | * of the tetrahedron \p p0, \p p1, \p p2, \p p3 | ||
| 307 | * \retval perturb() if \p p4 is exactly on the circumscribed sphere | ||
| 308 | * of the tetrahedron \p p0, \p p1, \p p2, \p p3, where \c perturb() | ||
| 309 | * denotes a globally consistent perturbation, that returns | ||
| 310 | * either POSITIVE or NEGATIVE | ||
| 311 | * \pre orient_3d(p0,p1,p2,p3) > 0 | ||
| 312 | */ | ||
| 313 | Sign GEOGRAM_API in_sphere_3d_SOS( | ||
| 314 | const double* p0, const double* p1, | ||
| 315 | const double* p2, const double* p3, | ||
| 316 | const double* p4 | ||
| 317 | ); | ||
| 318 | |||
| 319 | |||
| 320 | /** | ||
| 321 | * \brief Tests whether a 2d point is inside the | ||
| 322 | * circumscribed circle of a 3d triangle. | ||
| 323 | * \param[in] p0 , p1 , p2 vertices of the triangle | ||
| 324 | * \param[in] p3 the point to be tested | ||
| 325 | * \retval POSITIVE whenever \p p3 is inside the circumscribed circle | ||
| 326 | * of the triangle \p p0, \p p1, \p p2 | ||
| 327 | * \retval NEGATIVE whenever \p p2 is outside the circumscribed circle | ||
| 328 | * of the triangle \p p0, \p p1, \p p2 | ||
| 329 | * \retval perturb() if \p p3 is exactly on the circumscribed circle | ||
| 330 | * of the triangle \p p0, \p p1, \p p2, where \c perturb() | ||
| 331 | * denotes a globally consistent perturbation, that returns | ||
| 332 | * either POSITIVE or NEGATIVE | ||
| 333 | * \pre \p p3 belongs to the plane yielded by \p p0, \p p1 and \p p2 | ||
| 334 | */ | ||
| 335 | Sign GEOGRAM_API in_circle_2d_SOS( | ||
| 336 | const double* p0, const double* p1, const double* p2, | ||
| 337 | const double* p3 | ||
| 338 | ); | ||
| 339 | |||
| 340 | |||
| 341 | /** | ||
| 342 | * \brief Tests whether a 3d point is inside the | ||
| 343 | * circumscribed circle of a 3d triangle. | ||
| 344 | * \param[in] p0 , p1 , p2 vertices of the triangle | ||
| 345 | * \param[in] p3 the point to be tested | ||
| 346 | * \retval POSITIVE whenever \p p3 is inside the circumscribed circle | ||
| 347 | * of the triangle \p p0, \p p1, \p p2 | ||
| 348 | * \retval NEGATIVE whenever \p p2 is outside the circumscribed circle | ||
| 349 | * of the triangle \p p0, \p p1, \p p2 | ||
| 350 | * \retval perturb() if \p p3 is exactly on the circumscribed circle | ||
| 351 | * of the triangle \p p0, \p p1, \p p2, where \c perturb() | ||
| 352 | * denotes a globally consistent perturbation, that returns | ||
| 353 | * either POSITIVE or NEGATIVE | ||
| 354 | * \pre \p p3 belongs to the plane yielded by \p p0, \p p1 and \p p2 | ||
| 355 | */ | ||
| 356 | Sign GEOGRAM_API in_circle_3d_SOS( | ||
| 357 | const double* p0, const double* p1, const double* p2, | ||
| 358 | const double* p3 | ||
| 359 | ); | ||
| 360 | |||
| 361 | |||
| 362 | /** | ||
| 363 | * \brief Tests whether a lifted 3d point is inside the | ||
| 364 | * circumscribed circle of a lifted 3d triangle. | ||
| 365 | * \param[in] p0 , p1 , p2 vertices of the triangle | ||
| 366 | * \param[in] p3 the point to be tested | ||
| 367 | * \param[in] h0 , h1 , h2 lifted coordinate of the triangle vertices | ||
| 368 | * \param[in] h3 lifted coordinate of the point to be tested | ||
| 369 | * \param[in] SOS if true, do the symbolic perturbation in the degenerate | ||
| 370 | * cases | ||
| 371 | * \retval POSITIVE whenever (\p p3, \p h3) is inside the | ||
| 372 | * circumscribed circle of the triangle (\p p0,\p h0) (\p p1,\p h1), | ||
| 373 | * (\p p2, \p h2) | ||
| 374 | * \retval NEGATIVE whenever (\p p3, \p h3) is outside the | ||
| 375 | * circumscribed circle | ||
| 376 | * of the triangle (\p p0,\p h0) (\p p1,\p h1), (\p p2, \p h2) | ||
| 377 | * \retval perturb() if (\p p3, \p h3) is exactly | ||
| 378 | * on the circumscribed circle | ||
| 379 | * of the triangle (\p p0,\p h0) (\p p1,\p h1), (\p p2, \p h2) | ||
| 380 | * where \c perturb() denotes a globally consistent perturbation, | ||
| 381 | * that returns either POSITIVE or NEGATIVE | ||
| 382 | * \pre (\p p3, \p h3) belongs to the hyperplane yielded by | ||
| 383 | * (\p p0, \p h0), (\p p1, \p h1) and (\p p2, \p h2) | ||
| 384 | */ | ||
| 385 | Sign GEOGRAM_API in_circle_3dlifted_SOS( | ||
| 386 | const double* p0, const double* p1, const double* p2, | ||
| 387 | const double* p3, | ||
| 388 | double h0, double h1, double h2, double h3, | ||
| 389 | bool SOS=true | ||
| 390 | ); | ||
| 391 | |||
| 392 | /** | ||
| 393 | * \brief Computes the orientation predicate in 2d. | ||
| 394 | * \details Computes the sign of the signed area of | ||
| 395 | * the triangle p0, p1, p2. | ||
| 396 | * \param[in] p0 , p1 , p2 vertices of the triangle | ||
| 397 | * \retval POSITIVE if the triangle is oriented counter-clockwise | ||
| 398 | * \retval ZERO if the triangle is flat | ||
| 399 | * \retval NEGATIVE if the triangle is oriented clockwise | ||
| 400 | */ | ||
| 401 | Sign GEOGRAM_API orient_2d( | ||
| 402 | const double* p0, const double* p1, const double* p2 | ||
| 403 | ); | ||
| 404 | |||
| 405 | |||
| 406 | #ifndef GEOGRAM_PSM | ||
| 407 | /** | ||
| 408 | * \brief Computes the orientation predicate in 2d. | ||
| 409 | * \details Computes the sign of the signed area of | ||
| 410 | * the triangle p0, p1, p2. | ||
| 411 | * \param[in] p0 , p1 , p2 vertices of the triangle | ||
| 412 | * \retval POSITIVE if the triangle is oriented counter-clockwise | ||
| 413 | * \retval ZERO if the triangle is flat | ||
| 414 | * \retval NEGATIVE if the triangle is oriented clockwise | ||
| 415 | */ | ||
| 416 | inline Sign orient_2d( | ||
| 417 | const vec2& p0, const vec2& p1, const vec2& p2 | ||
| 418 | ) { | ||
| 419 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
146168 | return orient_2d(p0.data(),p1.data(),p2.data()); |
| 420 | } | ||
| 421 | #endif | ||
| 422 | |||
| 423 | /** | ||
| 424 | * \brief Computes the 3d orientation test with lifted points. | ||
| 425 | * \details Given three lifted points p0', p1', p2' in | ||
| 426 | * R^2, tests if the lifted point p3' in R^3 lies below or above | ||
| 427 | * the plane passing through the three points | ||
| 428 | * p0', p1', p2'. | ||
| 429 | * The first two coordinates and the | ||
| 430 | * third one are specified in separate arguments for each vertex. | ||
| 431 | * \param[in] p0 , p1 , p2 , p3 first 2 coordinates | ||
| 432 | * of the vertices of the 3-simplex | ||
| 433 | * \param[in] h0 , h1 , h2 , h3 heights of the vertices of | ||
| 434 | * the 3-simplex | ||
| 435 | * \retval POSITIVE if p3' lies below the plane | ||
| 436 | * \retval NEGATIVE if p3' lies above the plane | ||
| 437 | * \retval perturb() if p3' lies exactly on the hyperplane | ||
| 438 | * where \c perturb() denotes a globally | ||
| 439 | * consistent perturbation, that returns either POSITIVE or NEGATIVE | ||
| 440 | */ | ||
| 441 | Sign GEOGRAM_API orient_2dlifted_SOS( | ||
| 442 | const double* p0, const double* p1, | ||
| 443 | const double* p2, const double* p3, | ||
| 444 | double h0, double h1, double h2, double h3 | ||
| 445 | ); | ||
| 446 | |||
| 447 | |||
| 448 | /** | ||
| 449 | * \brief Computes the orientation predicate in 3d. | ||
| 450 | * \details Computes the sign of the signed volume of | ||
| 451 | * the tetrahedron p0, p1, p2, p3. | ||
| 452 | * \param[in] p0 , p1 , p2 , p3 vertices of the tetrahedron | ||
| 453 | * \retval POSITIVE if the tetrahedron is oriented positively | ||
| 454 | * \retval ZERO if the tetrahedron is flat | ||
| 455 | * \retval NEGATIVE if the tetrahedron is oriented negatively | ||
| 456 | */ | ||
| 457 | Sign GEOGRAM_API orient_3d( | ||
| 458 | const double* p0, const double* p1, | ||
| 459 | const double* p2, const double* p3 | ||
| 460 | ); | ||
| 461 | |||
| 462 | /** | ||
| 463 | * \brief Computes the orientation predicate in 3d. | ||
| 464 | * \details Computes the sign of the signed volume of | ||
| 465 | * the tetrahedron p0, p1, p2, p3. | ||
| 466 | * \param[in] p0 , p1 , p2 , p3 vertices of the tetrahedron | ||
| 467 | * \retval POSITIVE if the tetrahedron is oriented positively | ||
| 468 | * \retval NEGATIVE if the tetrahedron is oriented negatively | ||
| 469 | * \retval perturb() if the tetrahedron is flat, | ||
| 470 | * where \c perturb() denotes a globally | ||
| 471 | * consistent perturbation, that returns either POSITIVE or NEGATIVE | ||
| 472 | */ | ||
| 473 | Sign GEOGRAM_API orient_3d_SOS( | ||
| 474 | const double* p0, const double* p1, | ||
| 475 | const double* p2, const double* p3 | ||
| 476 | ); | ||
| 477 | |||
| 478 | #ifndef GEOGRAM_PSM | ||
| 479 | /** | ||
| 480 | * \brief Computes the orientation predicate in 3d. | ||
| 481 | * \details Computes the sign of the signed volume of | ||
| 482 | * the tetrahedron p0, p1, p2, p3. | ||
| 483 | * \param[in] p0 , p1 , p2 , p3 vertices of the tetrahedron | ||
| 484 | * \retval POSITIVE if the tetrahedron is oriented positively | ||
| 485 | * \retval ZERO if the tetrahedron is flat | ||
| 486 | * \retval NEGATIVE if the tetrahedron is oriented negatively | ||
| 487 | */ | ||
| 488 | inline Sign orient_3d( | ||
| 489 | const vec3& p0, const vec3& p1, | ||
| 490 | const vec3& p2, const vec3& p3 | ||
| 491 | ) { | ||
| 492 |
3/8✓ Branch 1 taken 60502 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22272 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 128 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
4878109 | return orient_3d(p0.data(),p1.data(),p2.data(),p3.data()); |
| 493 | } | ||
| 494 | |||
| 495 | /** | ||
| 496 | * \brief Computes the orientation predicate in 3d. | ||
| 497 | * \details Computes the sign of the signed volume of | ||
| 498 | * the tetrahedron p0, p1, p2, p3. | ||
| 499 | * \param[in] p0 , p1 , p2 , p3 vertices of the tetrahedron | ||
| 500 | * \retval POSITIVE if the tetrahedron is oriented positively | ||
| 501 | * \retval NEGATIVE if the tetrahedron is oriented negatively | ||
| 502 | * \retval perturb() if the tetrahedron is flat, | ||
| 503 | * where \c perturb() denotes a globally | ||
| 504 | * consistent perturbation, that returns either POSITIVE or NEGATIVE | ||
| 505 | */ | ||
| 506 | inline Sign GEOGRAM_API orient_3d_SOS( | ||
| 507 | const vec3& p0, const vec3& p1, | ||
| 508 | const vec3& p2, const vec3& p3 | ||
| 509 | ) { | ||
| 510 | 3163037 | return orient_3d_SOS(p0.data(),p1.data(),p2.data(),p3.data()); | |
| 511 | } | ||
| 512 | |||
| 513 | #endif | ||
| 514 | |||
| 515 | /** | ||
| 516 | * \brief Computes the 4d orientation test. | ||
| 517 | * \details Given four lifted points p0', p1', p2', and p3' in | ||
| 518 | * R^4, tests if the lifted point p4' in R^4 lies below or above | ||
| 519 | * the hyperplance passing through the four points | ||
| 520 | * p0', p1', p2', and p3'. | ||
| 521 | * This version does not apply symbolic perturbation. | ||
| 522 | * The first three coordinates and the | ||
| 523 | * fourth one are specified in separate arguments for each vertex. | ||
| 524 | * \param[in] p0 , p1 , p2 , p3 , p4 first 3 coordinates | ||
| 525 | * of the vertices of the 4-simplex | ||
| 526 | * \param[in] h0 , h1 , h2 , h3 , h4 heights of the vertices of | ||
| 527 | * the 4-simplex | ||
| 528 | * \retval POSITIVE if p4' lies below the hyperplane | ||
| 529 | * \retval NEGATIVE if p4' lies above the hyperplane | ||
| 530 | * \retval ZERO if p4' lies exactly on the hyperplane | ||
| 531 | */ | ||
| 532 | Sign GEOGRAM_API orient_3dlifted( | ||
| 533 | const double* p0, const double* p1, | ||
| 534 | const double* p2, const double* p3, const double* p4, | ||
| 535 | double h0, double h1, double h2, double h3, double h4 | ||
| 536 | ); | ||
| 537 | |||
| 538 | |||
| 539 | /** | ||
| 540 | * \brief Computes the 4d orientation test with symbolic perturbation. | ||
| 541 | * \details Given four lifted points p0', p1', p2', and p3' in | ||
| 542 | * R^4, tests if the lifted point p4' in R^4 lies below or above | ||
| 543 | * the hyperplance passing through the four | ||
| 544 | * points p0', p1', p2', and p3'. | ||
| 545 | * Symbolic perturbation is applied whenever the 5 vertices are | ||
| 546 | * not linearly independent. The first three coordinates and the | ||
| 547 | * fourth one are specified in separate arguments for each vertex. | ||
| 548 | * \param[in] p0 , p1 , p2 , p3 , p4 first 3 coordinates | ||
| 549 | * of the vertices of the 4-simplex | ||
| 550 | * \param[in] h0 , h1 , h2 , h3 , h4 heights of the vertices of | ||
| 551 | * the 4-simplex | ||
| 552 | * \retval POSITIVE if p4' lies below the hyperplane | ||
| 553 | * \retval NEGATIVE if p4' lies above the hyperplane | ||
| 554 | * \retval perturb() if p4' lies exactly on the hyperplane | ||
| 555 | * where \c perturb() denotes a globally | ||
| 556 | * consistent perturbation, that returns either POSITIVE or NEGATIVE | ||
| 557 | */ | ||
| 558 | Sign GEOGRAM_API orient_3dlifted_SOS( | ||
| 559 | const double* p0, const double* p1, | ||
| 560 | const double* p2, const double* p3, const double* p4, | ||
| 561 | double h0, double h1, double h2, double h3, double h4 | ||
| 562 | ); | ||
| 563 | |||
| 564 | |||
| 565 | /** | ||
| 566 | * \brief Computes the sign of the determinant of a 3x3 | ||
| 567 | * matrix formed by three 3d points. | ||
| 568 | * \param[in] p0 , p1 , p2 the three points | ||
| 569 | * \return the sign of the determinant of the matrix. | ||
| 570 | */ | ||
| 571 | Sign GEOGRAM_API det_3d( | ||
| 572 | const double* p0, const double* p1, const double* p2 | ||
| 573 | ); | ||
| 574 | |||
| 575 | #ifndef GEOGRAM_PSM | ||
| 576 | /** | ||
| 577 | * \brief Computes the sign of the determinant of a 3x3 | ||
| 578 | * matrix formed by three 3d points. | ||
| 579 | * \param[in] p0 , p1 , p2 the three points | ||
| 580 | * \return the sign of the determinant of the matrix. | ||
| 581 | */ | ||
| 582 | inline Sign det_3d( | ||
| 583 | const vec3& p0, const vec3& p1, const vec3& p2 | ||
| 584 | ) { | ||
| 585 | ✗ | return det_3d(p0.data(), p1.data(), p2.data()); | |
| 586 | } | ||
| 587 | #endif | ||
| 588 | |||
| 589 | /** | ||
| 590 | * \brief Computes the sign of the determinant of a 4x4 | ||
| 591 | * matrix formed by four 4d points. | ||
| 592 | * \param[in] p0 , p1 , p2 , p3 the four points | ||
| 593 | * \return the sign of the determinant of the matrix. | ||
| 594 | */ | ||
| 595 | Sign GEOGRAM_API det_4d( | ||
| 596 | const double* p0, const double* p1, | ||
| 597 | const double* p2, const double* p3 | ||
| 598 | ); | ||
| 599 | |||
| 600 | #ifndef GEOGRAM_PSM | ||
| 601 | /** | ||
| 602 | * \brief Computes the sign of the determinant of a 4x4 | ||
| 603 | * matrix formed by four 4d points. | ||
| 604 | * \param[in] p0 , p1 , p2 , p3 the four points | ||
| 605 | * \return the sign of the determinant of the matrix. | ||
| 606 | */ | ||
| 607 | inline Sign det_4d( | ||
| 608 | const vec4& p0, const vec4& p1, | ||
| 609 | const vec4& p2, const vec4& p3 | ||
| 610 | ) { | ||
| 611 | ✗ | return det_4d(p0.data(), p1.data(), p2.data(), p3.data()); | |
| 612 | } | ||
| 613 | #endif | ||
| 614 | |||
| 615 | /** | ||
| 616 | * \brief Computes the sign of the determinant of a | ||
| 617 | * 4x4 matrix formed by three 4d points and the | ||
| 618 | * difference of two 4d points. | ||
| 619 | * \param[in] p0 , p1 , p2 , p3 , p4 the four points | ||
| 620 | * \return the sign of the determinant of the matrix | ||
| 621 | * p0 p1 p2 p4-p3 | ||
| 622 | */ | ||
| 623 | Sign GEOGRAM_API det_compare_4d( | ||
| 624 | const double* p0, const double* p1, | ||
| 625 | const double* p2, const double* p3, | ||
| 626 | const double* p4 | ||
| 627 | ); | ||
| 628 | |||
| 629 | /** | ||
| 630 | * \brief Tests whether three points are aligned. | ||
| 631 | * \param[in] p0 , p1 , p2 the three points | ||
| 632 | * \retval true if the three points are aligned. | ||
| 633 | * \retval false otherwise. | ||
| 634 | * \details Function to be tested, use points_are_colinear_3d() | ||
| 635 | * instead. | ||
| 636 | */ | ||
| 637 | bool GEOGRAM_API aligned_3d( | ||
| 638 | const double* p0, const double* p1, const double* p2 | ||
| 639 | ); | ||
| 640 | |||
| 641 | /** | ||
| 642 | * \brief Computes the sign of the dot product between two | ||
| 643 | * vectors. | ||
| 644 | * \param[in] p0 , p1 , p2 three 3d points. | ||
| 645 | * \return the sign of the dot product between the vectors | ||
| 646 | * p0p1 and p0p2. | ||
| 647 | */ | ||
| 648 | Sign GEOGRAM_API dot_3d( | ||
| 649 | const double* p0, const double* p1, const double* p2 | ||
| 650 | ); | ||
| 651 | |||
| 652 | #ifndef GEOGRAM_PSM | ||
| 653 | |||
| 654 | /** | ||
| 655 | * \brief Tests whether three points are aligned. | ||
| 656 | * \param[in] p0 , p1 , p2 the three points | ||
| 657 | * \retval true if the three points are aligned. | ||
| 658 | * \retval false otherwise. | ||
| 659 | * \details Function to be tested, use points_are_colinear_3d() | ||
| 660 | * instead. | ||
| 661 | */ | ||
| 662 | inline bool aligned_3d( | ||
| 663 | const vec3& p0, const vec3& p1, const vec3& p2 | ||
| 664 | ) { | ||
| 665 |
0/4✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
114536 | return aligned_3d(p0.data(), p1.data(), p2.data()); |
| 666 | } | ||
| 667 | |||
| 668 | /** | ||
| 669 | * \brief Computes the sign of the dot product between two | ||
| 670 | * vectors. | ||
| 671 | * \param[in] p0 , p1 , p2 three 3d points. | ||
| 672 | * \return the sign of the dot product between the vectors | ||
| 673 | * p0p1 and p0p2. | ||
| 674 | */ | ||
| 675 | inline Sign dot_3d( | ||
| 676 | const vec3& p0, const vec3& p1, const vec3& p2 | ||
| 677 | ) { | ||
| 678 | return dot_3d(p0.data(), p1.data(), p2.data()); | ||
| 679 | } | ||
| 680 | #endif | ||
| 681 | |||
| 682 | /** | ||
| 683 | * \brief Compares two dot products. | ||
| 684 | * \param[in] v0 , v1 , v2 three vectors. | ||
| 685 | * \return the sign of v0.v1 - v0.v2 | ||
| 686 | */ | ||
| 687 | Sign GEOGRAM_API dot_compare_3d( | ||
| 688 | const double* v0, const double* v1, const double* v2 | ||
| 689 | ); | ||
| 690 | |||
| 691 | /** | ||
| 692 | * \brief Tests whether two 2d points are identical. | ||
| 693 | * \param[in] p1 first point | ||
| 694 | * \param[in] p2 second point | ||
| 695 | * \retval true if \p p1 and \p p2 have exactly the same | ||
| 696 | * coordinates | ||
| 697 | * \retval false otherwise | ||
| 698 | */ | ||
| 699 | bool points_are_identical_2d( | ||
| 700 | const double* p1, | ||
| 701 | const double* p2 | ||
| 702 | ); | ||
| 703 | |||
| 704 | /** | ||
| 705 | * \brief Tests whether two 3d points are identical. | ||
| 706 | * \param[in] p1 first point | ||
| 707 | * \param[in] p2 second point | ||
| 708 | * \retval true if \p p1 and \p p2 have exactly the same | ||
| 709 | * coordinates | ||
| 710 | * \retval false otherwise | ||
| 711 | */ | ||
| 712 | bool GEOGRAM_API points_are_identical_3d( | ||
| 713 | const double* p1, | ||
| 714 | const double* p2 | ||
| 715 | ); | ||
| 716 | |||
| 717 | /** | ||
| 718 | * \brief Tests whether three 3d points are colinear. | ||
| 719 | * \param[in] p1 first point | ||
| 720 | * \param[in] p2 second point | ||
| 721 | * \param[in] p3 third point | ||
| 722 | * \retval true if \p p1, \p p2 and \p p3 are colinear | ||
| 723 | * \retbal false otherwise | ||
| 724 | */ | ||
| 725 | bool GEOGRAM_API points_are_colinear_3d( | ||
| 726 | const double* p1, | ||
| 727 | const double* p2, | ||
| 728 | const double* p3 | ||
| 729 | ); | ||
| 730 | |||
| 731 | /** | ||
| 732 | * \brief Computes the (approximate) orientation predicate in 3d. | ||
| 733 | * \details Computes the sign of the (approximate) signed volume of | ||
| 734 | * the tetrahedron p0, p1, p2, p3. | ||
| 735 | * \param[in] p0 first vertex of the tetrahedron | ||
| 736 | * \param[in] p1 second vertex of the tetrahedron | ||
| 737 | * \param[in] p2 third vertex of the tetrahedron | ||
| 738 | * \param[in] p3 fourth vertex of the tetrahedron | ||
| 739 | * \retval POSITIVE if the tetrahedron is oriented positively | ||
| 740 | * \retval ZERO if the tetrahedron is flat | ||
| 741 | * \retval NEGATIVE if the tetrahedron is oriented negatively | ||
| 742 | */ | ||
| 743 | 63522 | inline Sign orient_3d_inexact( | |
| 744 | const double* p0, const double* p1, | ||
| 745 | const double* p2, const double* p3 | ||
| 746 | ) { | ||
| 747 | 63522 | double a11 = p1[0] - p0[0] ; | |
| 748 | 63522 | double a12 = p1[1] - p0[1] ; | |
| 749 | 63522 | double a13 = p1[2] - p0[2] ; | |
| 750 | |||
| 751 | 63522 | double a21 = p2[0] - p0[0] ; | |
| 752 | 63522 | double a22 = p2[1] - p0[1] ; | |
| 753 | 63522 | double a23 = p2[2] - p0[2] ; | |
| 754 | |||
| 755 | 63522 | double a31 = p3[0] - p0[0] ; | |
| 756 | 63522 | double a32 = p3[1] - p0[1] ; | |
| 757 | 63522 | double a33 = p3[2] - p0[2] ; | |
| 758 | |||
| 759 | double Delta = det3x3( | ||
| 760 | a11,a12,a13, | ||
| 761 | a21,a22,a23, | ||
| 762 | a31,a32,a33 | ||
| 763 | ); | ||
| 764 | |||
| 765 | 63522 | return geo_sgn(Delta); | |
| 766 | } | ||
| 767 | |||
| 768 | /** | ||
| 769 | * \brief Displays some statistics about predicates, | ||
| 770 | * including the number of calls, the number of exact arithmetics | ||
| 771 | * calls, and the number of Simulation of Simplicity calls. | ||
| 772 | */ | ||
| 773 | void GEOGRAM_API show_stats(); | ||
| 774 | |||
| 775 | /** | ||
| 776 | * \brief Needs to be called before using any predicate. | ||
| 777 | */ | ||
| 778 | void GEOGRAM_API initialize(); | ||
| 779 | |||
| 780 | /** | ||
| 781 | * \brief Needs to be called at the end of the program. | ||
| 782 | */ | ||
| 783 | void GEOGRAM_API terminate(); | ||
| 784 | } | ||
| 785 | } | ||
| 786 | |||
| 787 | /**************************************************************************/ | ||
| 788 | |||
| 789 | namespace GEO { | ||
| 790 | |||
| 791 | |||
| 792 | /** | ||
| 793 | * \brief Implementation of the perturbed orient_3d predicate | ||
| 794 | * \param[in] p0 , p1 , p2 , p3 the four points, as const references to T | ||
| 795 | * \tparam T point class | ||
| 796 | * \tparam SOS a class with: | ||
| 797 | * - constructor that takes the four points as const references to T | ||
| 798 | * - Sign orient_1d(v1, v2, axis) | ||
| 799 | * - Sign orient_2d(v1, v2, v3, axis1, axis2) | ||
| 800 | * where v1, v2, v3 are in {0,1,2,3} and axis, axis1, axis2 in {0,1,2}. | ||
| 801 | * The returned sign is multiplied by the parity of the order of the | ||
| 802 | * four points. | ||
| 803 | * \details This is the result of a discussion with Marc Alexa (01/2026), | ||
| 804 | * see also their article: A practical algorithm for weighted k-hulls, | ||
| 805 | * Look, Meyer, Alexa, SGP 2026 | ||
| 806 | */ | ||
| 807 | 6044551 | template<class T, class SOS> inline Sign orient_3d_SOS_impl( | |
| 808 | const T& p0, const T& p1, const T& p2, const T& p3 | ||
| 809 | ) { | ||
| 810 | constexpr coord_index_t X = 0, Y = 1, Z = 2; | ||
| 811 | 6044551 | Sign s = ::GEO::PCK::orient_3d(p0, p1, p2, p3); | |
| 812 |
2/2✓ Branch 0 taken 555358 times.
✓ Branch 1 taken 5489193 times.
|
6044551 | if(s != ZERO) { |
| 813 | return s; | ||
| 814 | } | ||
| 815 | |||
| 816 | // The perturbed determinant is as follows: | ||
| 817 | // | x1+eps y1+eps^2 z1+eps^4 1 | | ||
| 818 | // | x2+eps^8 y2+eps^16 z2+eps^32 1 | | ||
| 819 | // | x3+eps^64 y3+eps^128 z3+eps^256 1 | | ||
| 820 | // | x4+eps^512 y4+eps^1024 z4+eps^2048 1 | | ||
| 821 | // | ||
| 822 | // By developping and sorting by exponents of eps | ||
| 823 | // one gets the perturbations. Did it with TinyCAS: | ||
| 824 | // https://github.com/BrunoLevy/Experiment/blob/main/algo/tiny_cas.h | ||
| 825 | // | ||
| 826 | // | a b 1 | | ||
| 827 | // - The minors | c d 1 | correspond to orient_2d((a,b), (c,d), (e,f)) | ||
| 828 | // | e f 1 | | ||
| 829 | // | ||
| 830 | // - The other terms are just difference of coordinates (orient_1d) | ||
| 831 | |||
| 832 | // Static array that encodes all the terms of the expansion. | ||
| 833 | static const struct SOSInfo { | ||
| 834 | index_t dim; // 0: constant, 1: orient_1d, 2: orient_2d | ||
| 835 | index_t v1, v2, v3; // local indices of the two or three vertices | ||
| 836 | index_t ax1, ax2; // one or two projection axes | ||
| 837 | Sign sign; // sign of the term | ||
| 838 | } sosInfo[] = { | ||
| 839 | {2, 1, 2, 3, Y, Z, POSITIVE}, // eps | ||
| 840 | {2, 1, 2, 3, X, Z, NEGATIVE}, // eps^2 | ||
| 841 | {2, 1, 2, 3, X, Y, POSITIVE}, // eps^4 | ||
| 842 | {2, 0, 2, 3, Y, Z, NEGATIVE}, // eps^8 | ||
| 843 | {1, 3, 2, NO_INDEX, Z, NO_INDEX, POSITIVE}, // eps^10 | ||
| 844 | {1, 2, 3, NO_INDEX, Y, NO_INDEX, POSITIVE}, // eps^12 | ||
| 845 | {2, 0, 2, 3, X, Z, POSITIVE}, // eps^16 | ||
| 846 | // z2-z3 = -term in eps^10, already seen // eps^17 | ||
| 847 | {1, 3, 2, NO_INDEX, X, NO_INDEX, POSITIVE}, // eps^20 | ||
| 848 | {2, 0, 2, 3, X, Y, NEGATIVE}, // eps^32 | ||
| 849 | // y3-y2 = -term in eps^12, already seen // eps^33 | ||
| 850 | // x2-x3 = -term in eps^20, already seen // eps^34 | ||
| 851 | {2, 0, 1, 3, Y, Z, POSITIVE}, // eps^64 | ||
| 852 | {1, 1, 3, NO_INDEX, Z, NO_INDEX, POSITIVE}, // eps^66 | ||
| 853 | {1, 3, 1, NO_INDEX, Y, NO_INDEX, POSITIVE}, // eps^68 | ||
| 854 | {1, 3, 0, NO_INDEX, Z, NO_INDEX, POSITIVE}, // eps^80 | ||
| 855 | {0,NO_INDEX,NO_INDEX,NO_INDEX,NO_INDEX,NO_INDEX, NEGATIVE} // eps^84 | ||
| 856 | // There are more terms (up to eps^2184) but we do not need them, | ||
| 857 | // since we got a (constant) non-zero coefficient for eps^84 | ||
| 858 | }; | ||
| 859 | |||
| 860 | 555358 | SOS sos(p0, p1, p2, p3); | |
| 861 | |||
| 862 | 699835 | for(index_t k=0; ;++k) { | |
| 863 | const SOSInfo& I = sosInfo[k]; | ||
| 864 |
3/4✓ Branch 0 taken 56 times.
✓ Branch 1 taken 77545 times.
✓ Branch 2 taken 1177592 times.
✗ Branch 3 not taken.
|
1255193 | switch(I.dim) { |
| 865 | 56 | case 0: { | |
| 866 | 56 | return I.sign; | |
| 867 | } break; | ||
| 868 | 77545 | case 1: { | |
| 869 |
2/2✓ Branch 0 taken 31675 times.
✓ Branch 1 taken 45870 times.
|
77545 | s = sos.orient_1d(I.v1, I.v2, I.ax1); |
| 870 |
2/2✓ Branch 0 taken 31675 times.
✓ Branch 1 taken 45870 times.
|
77545 | if(s != ZERO) { |
| 871 | 31675 | return Sign(I.sign*s); | |
| 872 | } | ||
| 873 | } break; | ||
| 874 | 1177592 | case 2: { | |
| 875 | 1177592 | s = sos.orient_2d(I.v1, I.v2, I.v3, I.ax1, I.ax2); | |
| 876 |
2/2✓ Branch 0 taken 523627 times.
✓ Branch 1 taken 653965 times.
|
1177592 | if(s != ZERO) { |
| 877 | 523627 | return Sign(I.sign*s); | |
| 878 | } | ||
| 879 | } break; | ||
| 880 | default: | ||
| 881 | ✗ | geo_assert_not_reached; | |
| 882 | } | ||
| 883 | } | ||
| 884 | geo_assert_not_reached; | ||
| 885 | } | ||
| 886 | |||
| 887 | namespace Permutation { | ||
| 888 | /** | ||
| 889 | * \brief Computes the parity of a permutation | ||
| 890 | * \param[in] orig an array of pointers, typically nD points | ||
| 891 | * \param[in] perm a permutation of \p orig | ||
| 892 | * \param[in] n size of \p orig and \p perm, 64 max | ||
| 893 | * \retval true if \p perm is an odd permutation of \p orig | ||
| 894 | * \retval false if \p perm is an even permutation of \p orig | ||
| 895 | * \pre \p n <= 64 and \p perm is a permutation of \p orig | ||
| 896 | * \details operates in O(n^2) (only use for small arrays) | ||
| 897 | */ | ||
| 898 | 555358 | template <class T> inline bool permutation_is_odd( | |
| 899 | const T** orig, const T** perm, index_t n | ||
| 900 | ) { | ||
| 901 | geo_debug_assert(n <= 64); | ||
| 902 | Numeric::uint64 visited = 0; | ||
| 903 | bool odd = false; | ||
| 904 |
2/2✓ Branch 0 taken 2221432 times.
✓ Branch 1 taken 555358 times.
|
2776790 | for (index_t i = 0; i < n; ++i) { |
| 905 |
2/2✓ Branch 0 taken 1125617 times.
✓ Branch 1 taken 1095815 times.
|
2221432 | if ((visited >> i) & 1) { |
| 906 | 1095815 | continue; | |
| 907 | } | ||
| 908 | // Compute the length of the cycle starting from perm[i] | ||
| 909 | index_t len = 0; | ||
| 910 |
2/2✓ Branch 0 taken 2221432 times.
✓ Branch 1 taken 1125617 times.
|
3347049 | for (index_t j = i; !((visited >> j) & 1); ) { |
| 911 | 2221432 | visited |= (Numeric::uint64(1) << j); | |
| 912 | 2221432 | ++len; | |
| 913 | 4442864 | j = index_t(std::find(orig, orig + n, perm[j]) - orig); | |
| 914 | } | ||
| 915 | // even-length cycle contributes odd parity | ||
| 916 |
2/2✓ Branch 0 taken 363961 times.
✓ Branch 1 taken 761656 times.
|
1125617 | if (len % 2 == 0) odd = !odd; |
| 917 | } | ||
| 918 | 555358 | return odd; | |
| 919 | } | ||
| 920 | } | ||
| 921 | } | ||
| 922 | |||
| 923 | #endif | ||
| 924 |