predicates.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2000-2022 Inria
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *
 *  * Redistributions of source code must retain the above copyright notice,
 *  this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright notice,
 *  this list of conditions and the following disclaimer in the documentation
 *  and/or other materials provided with the distribution.
 *  * Neither the name of the ALICE Project-Team nor the names of its
 *  contributors may be used to endorse or promote products derived from this
 *  software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 *
 *  Contact: Bruno Levy
 *
 *     https://www.inria.fr/fr/bruno-levy
 *
 *     Inria,
 *     Domaine de Voluceau,
 *     78150 Le Chesnay - Rocquencourt
 *     FRANCE
 *
 */
 
#ifndef GEOGRAM_NUMERICS_PREDICATES
#define GEOGRAM_NUMERICS_PREDICATES
 
#include <geogram/basic/common.h>
#include <geogram/basic/numeric.h>
#include <geogram/basic/geometry.h>
#include <geogram/numerics/PCK.h>
 
namespace GEO {
 
    namespace PCK {
 
        enum SOSMode { SOS_ADDRESS, SOS_LEXICO };
 
        void GEOGRAM_API set_SOS_mode(SOSMode m);
 
        SOSMode GEOGRAM_API get_SOS_mode();
 
        Sign GEOGRAM_API side1_SOS(
            const double* p0, const double* p1,
            const double* q0,
            coord_index_t DIM
        );
 
        Sign GEOGRAM_API side2_SOS(
            const double* p0, const double* p1, const double* p2,
            const double* q0, const double* q1,
            coord_index_t DIM
        );
 
        Sign GEOGRAM_API side3_SOS(
            const double* p0, const double* p1,
            const double* p2, const double* p3,
            const double* q0, const double* q1, const double* q2,
            coord_index_t DIM
        );
 
        Sign GEOGRAM_API side3_3dlifted_SOS(
            const double* p0, const double* p1,
            const double* p2, const double* p3,
            double h0, double h1, double h2, double h3,
            const double* q0, const double* q1, const double* q2,
            bool SOS=true
        );
 
        Sign GEOGRAM_API side4_SOS(
            const double* p0,
            const double* p1, const double* p2,
            const double* p3, const double* p4,
            const double* q0, const double* q1,
            const double* q2, const double* q3,
            coord_index_t DIM
        );
 
 
        Sign GEOGRAM_API side4_3d(
            const double* p0,
            const double* p1, const double* p2,
            const double* p3, const double* p4
        );
 
        Sign GEOGRAM_API side4_3d_SOS(
            const double* p0, const double* p1,
            const double* p2, const double* p3, const double* p4
        );
 
        Sign GEOGRAM_API in_sphere_3d_SOS(
            const double* p0, const double* p1,
            const double* p2, const double* p3,
            const double* p4
        );
 
 
        Sign GEOGRAM_API in_circle_2d_SOS(
            const double* p0, const double* p1, const double* p2,
            const double* p3
        );
 
 
        Sign GEOGRAM_API in_circle_3d_SOS(
            const double* p0, const double* p1, const double* p2,
            const double* p3
        );
 
 
        Sign GEOGRAM_API in_circle_3dlifted_SOS(
            const double* p0, const double* p1, const double* p2,
            const double* p3,
            double h0, double h1, double h2, double h3,
            bool SOS=true
        );
 
        Sign GEOGRAM_API orient_2d(
            const double* p0, const double* p1, const double* p2
        );
 
 
#ifndef GEOGRAM_PSM
        inline Sign orient_2d(
            const vec2& p0, const vec2& p1, const vec2& p2
        ) {
            return orient_2d(p0.data(),p1.data(),p2.data());
        }
#endif
 
        Sign GEOGRAM_API orient_2dlifted_SOS(
            const double* p0, const double* p1,
            const double* p2, const double* p3,
            double h0, double h1, double h2, double h3
        );
 
 
        Sign GEOGRAM_API orient_3d(
            const double* p0, const double* p1,
            const double* p2, const double* p3
        );
 
        Sign GEOGRAM_API orient_3d_SOS(
            const double* p0, const double* p1,
            const double* p2, const double* p3
        );
 
#ifndef GEOGRAM_PSM
        inline Sign orient_3d(
            const vec3& p0, const vec3& p1,
            const vec3& p2, const vec3& p3
        ) {
            return orient_3d(p0.data(),p1.data(),p2.data(),p3.data());
        }
 
        inline Sign GEOGRAM_API orient_3d_SOS(
            const vec3& p0, const vec3& p1,
            const vec3& p2, const vec3& p3
        ) {
            return orient_3d_SOS(p0.data(),p1.data(),p2.data(),p3.data());
        }
 
#endif
 
        Sign GEOGRAM_API orient_3dlifted(
            const double* p0, const double* p1,
            const double* p2, const double* p3, const double* p4,
            double h0, double h1, double h2, double h3, double h4
        );
 
 
        Sign GEOGRAM_API orient_3dlifted_SOS(
            const double* p0, const double* p1,
            const double* p2, const double* p3, const double* p4,
            double h0, double h1, double h2, double h3, double h4
        );
 
 
        Sign GEOGRAM_API det_3d(
            const double* p0, const double* p1, const double* p2
        );
 
#ifndef GEOGRAM_PSM
        inline Sign det_3d(
            const vec3& p0, const vec3& p1, const vec3& p2
        ) {
            return det_3d(p0.data(), p1.data(), p2.data());
        }
#endif
 
        Sign GEOGRAM_API det_4d(
            const double* p0, const double* p1,
            const double* p2, const double* p3
        );
 
#ifndef GEOGRAM_PSM
        inline Sign det_4d(
            const vec4& p0, const vec4& p1,
            const vec4& p2, const vec4& p3
        ) {
            return det_4d(p0.data(), p1.data(), p2.data(), p3.data());
        }
#endif
 
        Sign GEOGRAM_API det_compare_4d(
            const double* p0, const double* p1,
            const double* p2, const double* p3,
            const double* p4
        );
 
        bool GEOGRAM_API aligned_3d(
            const double* p0, const double* p1, const double* p2
        );
 
        Sign GEOGRAM_API dot_3d(
            const double* p0, const double* p1, const double* p2
        );
 
#ifndef GEOGRAM_PSM
 
        inline bool aligned_3d(
            const vec3& p0, const vec3& p1, const vec3& p2
        ) {
            return aligned_3d(p0.data(), p1.data(), p2.data());
        }
 
        inline Sign dot_3d(
            const vec3& p0, const vec3& p1, const vec3& p2
        ) {
            return dot_3d(p0.data(), p1.data(), p2.data());
        }
#endif
 
        Sign GEOGRAM_API dot_compare_3d(
            const double* v0, const double* v1, const double* v2
        );
 
        bool points_are_identical_2d(
            const double* p1,
            const double* p2
        );
 
        bool GEOGRAM_API points_are_identical_3d(
            const double* p1,
            const double* p2
        );
 
        bool GEOGRAM_API points_are_colinear_3d(
            const double* p1,
            const double* p2,
            const double* p3
        );
 
        inline Sign orient_3d_inexact(
            const double* p0, const double* p1,
            const double* p2, const double* p3
        ) {
            double a11 = p1[0] - p0[0] ;
            double a12 = p1[1] - p0[1] ;
            double a13 = p1[2] - p0[2] ;
 
            double a21 = p2[0] - p0[0] ;
            double a22 = p2[1] - p0[1] ;
            double a23 = p2[2] - p0[2] ;
 
            double a31 = p3[0] - p0[0] ;
            double a32 = p3[1] - p0[1] ;
            double a33 = p3[2] - p0[2] ;
 
            double Delta = det3x3(
                a11,a12,a13,
                a21,a22,a23,
                a31,a32,a33
            );
 
            return geo_sgn(Delta);
        }
 
        void GEOGRAM_API show_stats();
 
        void GEOGRAM_API initialize();
 
        void GEOGRAM_API terminate();
    }
}
 
/**************************************************************************/
 
namespace GEO {
 
 
    template<class T, class SOS> inline Sign orient_3d_SOS_impl(
        const T& p0, const T& p1, const T& p2, const T& p3
    ) {
        constexpr coord_index_t X = 0, Y = 1, Z = 2;
        Sign s = ::GEO::PCK::orient_3d(p0, p1, p2, p3);
        if(s != ZERO) {
            return s;
        }
 
        // The perturbed determinant is as follows:
        // | x1+eps     y1+eps^2    z1+eps^4    1 |
        // | x2+eps^8   y2+eps^16   z2+eps^32   1 |
        // | x3+eps^64  y3+eps^128  z3+eps^256  1 |
        // | x4+eps^512 y4+eps^1024 z4+eps^2048 1 |
        //
        // By developping and sorting by exponents of eps
        // one gets the perturbations. Did it with TinyCAS:
        // https://github.com/BrunoLevy/Experiment/blob/main/algo/tiny_cas.h
        //
        //              | a b 1 |
        // - The minors | c d 1 | correspond to orient_2d((a,b), (c,d), (e,f))
        //              | e f 1 |
        //
        // - The other terms are just difference of coordinates (orient_1d)
 
        // Static array that encodes all the terms of the expansion.
        static const struct SOSInfo {
            index_t dim;        // 0: constant, 1: orient_1d, 2: orient_2d
            index_t v1, v2, v3; // local indices of the two or three vertices
            index_t ax1, ax2;   // one or two projection axes
            Sign sign;          // sign of the term
        } sosInfo[] = {
            {2,  1, 2, 3,          Y, Z,         POSITIVE}, // eps
            {2,  1, 2, 3,          X, Z,         NEGATIVE}, // eps^2
            {2,  1, 2, 3,          X, Y,         POSITIVE}, // eps^4
            {2,  0, 2, 3,          Y, Z,         NEGATIVE}, // eps^8
            {1,  3, 2, NO_INDEX,   Z, NO_INDEX,  POSITIVE}, // eps^10
            {1,  2, 3, NO_INDEX,   Y, NO_INDEX,  POSITIVE}, // eps^12
            {2,  0, 2, 3,          X, Z,         POSITIVE}, // eps^16
            // z2-z3 = -term in eps^10, already seen        // eps^17
            {1,  3, 2, NO_INDEX,   X, NO_INDEX,  POSITIVE}, // eps^20
            {2,  0, 2, 3,          X, Y,         NEGATIVE}, // eps^32
            // y3-y2 = -term in eps^12, already seen        // eps^33
            // x2-x3 = -term in eps^20, already seen        // eps^34
            {2,  0, 1, 3,          Y, Z,         POSITIVE}, // eps^64
            {1,  1, 3, NO_INDEX,   Z, NO_INDEX,  POSITIVE}, // eps^66
            {1,  3, 1, NO_INDEX,   Y, NO_INDEX,  POSITIVE}, // eps^68
            {1,  3, 0, NO_INDEX,   Z, NO_INDEX,  POSITIVE}, // eps^80
            {0,NO_INDEX,NO_INDEX,NO_INDEX,NO_INDEX,NO_INDEX, NEGATIVE} // eps^84
            // There are more terms (up to eps^2184) but we do not need them,
            // since we got a (constant) non-zero coefficient for eps^84
        };
 
        SOS sos(p0, p1, p2, p3);
 
        for(index_t k=0; ;++k) {
            const SOSInfo& I = sosInfo[k];
            switch(I.dim) {
            case 0: {
                return I.sign;
            } break;
            case 1: {
                s = sos.orient_1d(I.v1, I.v2, I.ax1);
                if(s != ZERO) {
                    return Sign(I.sign*s);
                }
            } break;
            case 2: {
                s = sos.orient_2d(I.v1, I.v2, I.v3, I.ax1, I.ax2);
                if(s != ZERO) {
                    return Sign(I.sign*s);
                }
            } break;
            default:
                geo_assert_not_reached;
            }
        }
        geo_assert_not_reached;
    }
 
    namespace Permutation {
        template <class T> inline bool permutation_is_odd(
            const T** orig, const T** perm, index_t n
        ) {
            geo_debug_assert(n <= 64);
            Numeric::uint64 visited = 0;
            bool odd = false;
            for (index_t i = 0; i < n; ++i) {
                if ((visited >> i) & 1) {
                    continue;
                }
                // Compute the length of the cycle starting from perm[i]
                index_t len = 0;
                for (index_t j = i; !((visited >> j) & 1); ) {
                    visited |= (Numeric::uint64(1) << j);
                    ++len;
                    j = index_t(std::find(orig, orig + n, perm[j]) - orig);
                }
                // even-length cycle contributes odd parity
                if (len % 2 == 0) odd = !odd;
            }
            return odd;
        }
    }
}
 
#endif