exact_geometry.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2000-2023 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_EXACT_GEOMETRY
#define GEOGRAM_NUMERICS_EXACT_GEOMETRY
 
#include <geogram/basic/common.h>
#include <geogram/basic/geometry.h>
#include <geogram/basic/vechg.h>
#include <geogram/numerics/expansion_nt.h>
#include <geogram/numerics/interval_nt.h>
 
#include <geogram/numerics/predicates.h>
#include <geogram/numerics/exact_geometry.h>
 
#ifdef GEOGRAM_WITH_GEOGRAMPLUS
#include <geogram/geogramplus/numerics/exact_geometry.h>
#endif
 
// If Tessael's geogramplus is available, use exact_nt coordinates,
// else use expansion_nt coordinates.
// exact_nt coordinates makes the algorithm  10x to 20x faster
// and have no risk of underflow / overflow.
#ifdef GEOGRAM_WITH_GEOGRAMPLUS
#define GEOGRAM_USE_EXACT_NT  
#endif
 
namespace GEO {
 
    typedef vecng<2,expansion_nt> vec2E;
 
    typedef vecng<3,expansion_nt> vec3E;    
 
    typedef vecng<2,interval_nt> vec2I;    
    
    typedef vecng<3,interval_nt> vec3I;    
 
    typedef vec2Hg<expansion_nt> vec2HE;
 
    typedef vec3Hg<expansion_nt> vec3HE;
 
    typedef vec2Hg<interval_nt> vec2HI;
 
    typedef vec3Hg<interval_nt> vec3HI;
 
    /***********************************************************************/
 
    template <class VEC3 = vec3>
    inline VEC3 make_vec3(const vec3& p1, const vec3& p2) {
        typedef typename VEC3::value_type value_type;
        return VEC3(
            value_type(p2.x) - value_type(p1.x),
            value_type(p2.y) - value_type(p1.y),
            value_type(p2.z) - value_type(p1.z)
        );
    }
 
    template <class VEC2>
    inline VEC2 make_vec2(
        const vec2& p1, const vec2& p2
    ) {
        typedef typename VEC2::value_type value_type;        
        return VEC2(
            value_type(p2.x) - value_type(p1.x),
            value_type(p2.y) - value_type(p1.y)
        );
    }
    
    template <class VEC3>
    inline VEC3 triangle_normal(
        const vec3& p1, const vec3& p2, const vec3& p3
    ) {
        return cross(
            make_vec3<VEC3>(p1,p2),
            make_vec3<VEC3>(p1,p3)
        );
    }
 
    /***********************************************************************/
 
    namespace PCK {
 
        Sign GEOGRAM_API orient_2d(
            const vec2HE& p0, const vec2HE& p1, const vec2HE& p2
        );
 
        Sign GEOGRAM_API orient_2d_projected(
            const vec3HE& p0, const vec3HE& p1, const vec3HE& p2,
            coord_index_t axis
        );
 
        Sign GEOGRAM_API orient_3d(
            const vec3HE& p0, const vec3HE& p1,
            const vec3HE& p2, const vec3HE& p3
        );
 
        Sign GEOGRAM_API dot_2d(
            const vec2HE& p0, const vec2HE& p1, const vec2HE& p2
        );
 
        Sign GEOGRAM_API incircle_2d_SOS_with_lengths(
            const vec2HE& p0, const vec2HE& p1,
            const vec2HE& p2, const vec2HE& p3,
            double l0, double l1, double l2, double l3
        );
 
        Sign GEOGRAM_API incircle_2d_SOS_with_lengths(
            const vec2HE& p0, const vec2HE& p1,
            const vec2HE& p2, const vec2HE& p3,
            double l0, double l1, double l2, double l3
        );
        
        inline Sign incircle_2d_SOS(
            const vec2HE& p0, const vec2HE& p1,
            const vec2HE& p2, const vec2HE& p3
        ) {
            double l0 = (geo_sqr(p0.x) + geo_sqr(p0.y)).estimate() /
                         geo_sqr(p0.w).estimate();
            double l1 = (geo_sqr(p1.x) + geo_sqr(p1.y)).estimate() /
                         geo_sqr(p1.w).estimate();
            double l2 = (geo_sqr(p2.x) + geo_sqr(p2.y)).estimate() /
                         geo_sqr(p2.w).estimate();
            double l3 = (geo_sqr(p3.x) + geo_sqr(p3.y)).estimate() /
                         geo_sqr(p3.w).estimate();
            return incircle_2d_SOS_with_lengths(p0,p1,p2,p3,l0,l1,l2,l3);
        }
 
        coord_index_t GEOGRAM_API triangle_normal_axis(
            const vec3& p1, const vec3& p2, const vec3& p3
        );
 
        bool GEOGRAM_API aligned_3d(
            const vec3HE& p0, const vec3HE& p1, const vec3HE& p2
        );
 
        bool GEOGRAM_API on_segment_3d(
            const vec3HE& p, const vec3HE& q1, const vec3HE& q2
        );
        
        vec3 GEOGRAM_API approximate(const vec3HE& p);
 
        vec2 GEOGRAM_API approximate(const vec2HE& p);
        
    }
 
    /************************************************************************/
 
    template <>
    inline vec2E make_vec2<vec2E>(const vec2& p1, const vec2& p2) {
        return vec2E(
            expansion_nt(expansion_nt::DIFF, p2.x, p1.x),
            expansion_nt(expansion_nt::DIFF, p2.y, p1.y)
        );
    }
    
    template <>
    inline vec3E make_vec3<vec3E>(const vec3& p1, const vec3& p2) {
        return vec3E(
            expansion_nt(expansion_nt::DIFF, p2.x, p1.x),
            expansion_nt(expansion_nt::DIFF, p2.y, p1.y),
            expansion_nt(expansion_nt::DIFF, p2.z, p1.z)            
        );
    }
    
// Under Linux we got 10 Mb of stack (!) Then some operations can be
// made faster by using the low-level expansion API (that allocates
// intermediary multiprecision values on stack rather than in the heap).
// These optimized functions are written as template specializations
// (used automatically).    
    
#ifdef GEO_HAS_BIG_STACK
 
    template<> expansion_nt GEOGRAM_API det(const vec2E& v1, const vec2E& v2);
 
    template<> expansion_nt GEOGRAM_API dot(const vec2E& v1, const vec2E& v2);
 
    template<> expansion_nt GEOGRAM_API dot(const vec3E& v1, const vec3E& v2);
    
    template<> vec2Hg<expansion_nt> GEOGRAM_API mix(
        const rationalg<expansion_nt>& t,
        const vecng<2,double>& p1, const vecng<2,double>& p2
    );
 
    template<> vec3Hg<expansion_nt> GEOGRAM_API mix(
        const rationalg<expansion_nt>& t,
        const vecng<3,double>& p1, const vecng<3,double>& p2
    );
    
    template <> GEOGRAM_API vec3E triangle_normal<vec3E>(
        const vec3& p1, const vec3& p2, const vec3& p3
    );
    
#endif
    
    /************************************************************************/
 
    namespace exact {
#ifdef GEOGRAM_USE_EXACT_NT
        typedef exact_nt scalar;     
#else
        typedef expansion_nt scalar; 
#endif
        typedef vecng<2,scalar> vec2; 
        typedef vecng<3,scalar> vec3; 
        typedef vec2Hg<scalar> vec2h;
        
        typedef vec3Hg<scalar> vec3h;
 
        typedef rationalg<scalar> rational;
    }
}
 
#endif