Function intersectionPosition

include/CDT/include/Triangulation.hpp:448–471

// Precondition: ab and cd intersect normally
template <typename T>
V2d<T> intersectionPosition(
    const V2d<T>& a,
    const V2d<T>& b,
    const V2d<T>& c,
    const V2d<T>& d)
{
    using namespace predicates::adaptive;

    // note: for better accuracy we interpolate x and y separately
    // on a segment with the shortest x/y-projection correspondingly
    const T a_cd = orient2d(c.x, c.y, d.x, d.y, a.x, a.y);
    const T b_cd = orient2d(c.x, c.y, d.x, d.y, b.x, b.y);
    const T t_ab = a_cd / (a_cd - b_cd);

    const T c_ab = orient2d(a.x, a.y, b.x, b.y, c.x, c.y);
    const T d_ab = orient2d(a.x, a.y, b.x, b.y, d.x, d.y);
    const T t_cd = c_ab / (c_ab - d_ab);

    return V2d<T>(
        std::fabs(a.x - b.x) < std::fabs(c.x - d.x) ? lerp(a.x, b.x, t_ab)
                                                    : lerp(c.x, d.x, t_cd),
        std::fabs(a.y - b.y) < std::fabs(c.y - d.y) ? lerp(a.y, b.y, t_ab)
                                                    : lerp(c.y, d.y, t_cd));
}

} // namespace detail