triangle.h

/*****************************************************************************
 * VCLib                                                                     *
 * Visual Computing Library                                                  *
 *                                                                           *
 * Copyright(C) 2021-2025                                                    *
 * Visual Computing Lab                                                      *
 * ISTI - Italian National Research Council                                  *
 *                                                                           *
 * All rights reserved.                                                      *
 *                                                                           *
 * This program is free software; you can redistribute it and/or modify      *
 * it under the terms of the Mozilla Public License Version 2.0 as published *
 * by the Mozilla Foundation; either version 2 of the License, or            *
 * (at your option) any later version.                                       *
 *                                                                           *
 * This program is distributed in the hope that it will be useful,           *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the              *
 * Mozilla Public License Version 2.0                                        *
 * (https://www.mozilla.org/en-US/MPL/2.0/) for more details.                *
 ****************************************************************************/
 
#ifndef VCL_SPACE_CORE_TRIANGLE_H
#define VCL_SPACE_CORE_TRIANGLE_H
 
#include "point.h"
 
namespace vcl {
 
template<PointConcept PointT>
class Triangle
{
    std::array<PointT, 3> mPoints;
 
public:
    using ScalarType = PointT::ScalarType;
 
    using PointType = PointT;
 
    static const uint DIM = PointT::DIM;
 
    Triangle() = default;
 
    Triangle(const PointT& p0, const PointT& p1, const PointT& p2) :
            mPoints(p0, p1, p2)
    {
    }
 
    constexpr uint size() const { return 3; }
 
    const PointT& point(uint i) const { return mPoints.at(i); }
 
    PointT& point(uint i) { return mPoints.at(i); }
 
    PointT& point0() { return mPoints[0]; }
 
    const PointT& point0() const { return mPoints[0]; }
 
    PointT& point1() { return mPoints[1]; }
 
    const PointT& point1() const { return mPoints[1]; }
 
    PointT& point2() { return mPoints[2]; }
 
    const PointT& point2() const { return mPoints[2]; }
 
    ScalarType sideLength(uint i) const
    {
        return sideLength(mPoints[0], mPoints[1], mPoints[2], i);
    }
 
    ScalarType sideLength0() const { return sideLength(0); }
 
    ScalarType sideLength1() const { return sideLength(1); }
 
    ScalarType sideLength2() const { return sideLength(2); }
 
    PointT normal() const requires (PointT::DIM == 3)
    {
        return normal(mPoints[0], mPoints[1], mPoints[2]);
    }
 
    PointT barycenter() const
    {
        return barycenter(mPoints[0], mPoints[1], mPoints[2]);
    }
 
    PointT weightedBarycenter(ScalarType w0, ScalarType w1, ScalarType w2) const
    {
        return weightedBarycenter(
            mPoints[0], mPoints[1], mPoints[2], w0, w1, w2);
    }
 
    PointT weightedBarycenter(const Point3<ScalarType>& w) const
    {
        return weightedBarycenter(w(0), w(1), w(2));
    }
 
    PointT barycentricCoordinatePoint(
        ScalarType b0,
        ScalarType b1,
        ScalarType b2) const
    {
        return barycentricCoordinatePoint(
            mPoints[0], mPoints[1], mPoints[2], b0, b1, b2);
    }
 
    PointT barycentricCoordinatePoint(const Point3<ScalarType>& b) const
    {
        return barycentricCoordinatePoint(
            mPoints[0], mPoints[1], mPoints[2], b(0), b(1), b(2));
    }
 
    PointT circumcenter() const
    {
        return circumcenter(mPoints[0], mPoints[1], mPoints[2]);
    }
 
    ScalarType perimeter() const
    {
        return perimeter(mPoints[0], mPoints[1], mPoints[2]);
    }
 
    ScalarType area() const { return area(mPoints[0], mPoints[1], mPoints[2]); }
 
    ScalarType quality() const
    {
        return quality(mPoints[0], mPoints[1], mPoints[2]);
    }
 
    ScalarType qualityRadii() const
    {
        return qualityRadii(mPoints[0], mPoints[1], mPoints[2]);
    }
 
    ScalarType qualityMeanRatio() const
    {
        return qualityMeanRatio(mPoints[0], mPoints[1], mPoints[2]);
    }
 
    /* Static member functions */
 
    ScalarType sideLength(
        const PointT& p0,
        const PointT& p1,
        const PointT& p2,
        uint          i) const
    {
        switch (i % 3) {
        case 0: return p0.dist(p1);
        case 1: return p1.dist(p2);
        case 2: return p2.dist(p0);
        default: throw std::out_of_range("Invalid side index");
        }
    }
 
    static PointT normal(const PointT& p0, const PointT& p1, const PointT& p2)
        requires (PointT::DIM == 3)
    {
        return (p1 - p0).cross(p2 - p0);
    }
 
    static PointT barycenter(
        const PointT& p0,
        const PointT& p1,
        const PointT& p2)
    {
        return (p0 + p1 + p2) / 3;
    }
 
    static PointT weightedBarycenter(
        const PointT& p0,
        const PointT& p1,
        const PointT& p2,
        ScalarType    w0,
        ScalarType    w1,
        ScalarType    w2)
    {
        return (p0 * w0 + p1 * w1 + p2 * w2) / (w0 + w1 + w2);
    }
 
    static PointT barycentricCoordinatePoint(
        const PointT& p0,
        const PointT& p1,
        const PointT& p2,
        ScalarType    b0,
        ScalarType    b1,
        ScalarType    b2)
    {
        return p0 * b0 + p1 * b1 + p2 * b2;
    }
 
    static PointT circumcenter(
        const PointT& p0,
        const PointT& p1,
        const PointT& p2)
    {
        ScalarType a2 = p1.squaredDist(p2);
        ScalarType b2 = p2.squaredDist(p0);
        ScalarType c2 = p0.squaredDist(p1);
 
        PointType c = p0 * a2 * (-a2 + b2 + c2) + p1 * b2 * (a2 - b2 + c2) +
                      p2 * c2 * (a2 + b2 - c2);
 
        c /= 2 * (a2 * b2 + a2 * c2 + b2 * c2) - a2 * a2 - b2 * b2 - c2 * c2;
 
        return c;
    }
 
    static ScalarType perimeter(
        const PointT& p0,
        const PointT& p1,
        const PointT& p2)
    {
        return p0.dist(p1) + p1.dist(p2) + p2.dist(p0);
    }
 
    static ScalarType area(const PointT& p0, const PointT& p1, const PointT& p2)
    {
        if constexpr (DIM == 2) {
            return ((p1[0] - p0[0]) * (p2[1] - p0[1]) -
                    (p2[0] - p0[0]) * (p1[1] - p0[1])) /
                   2;
        }
        if constexpr (DIM == 3) {
            return normal(p0, p1, p2).norm() / 2;
        }
        else {
            // heron's formula
            ScalarType s = perimeter() / 2;
            return std::sqrt(
                s * (s - p0.dist(p1)) * (s - p1.dist(p2)) * (s - p2.dist(p0)));
        }
    }
 
    static ScalarType quality(
        const PointT& p0,
        const PointT& p1,
        const PointT& p2)
    {
        ScalarType a = area();
        if (a == 0)
            return 0; // Area zero triangles have surely quality==0;
 
        PointType d10 = p1 - p0;
        PointType d20 = p2 - p0;
        PointType d12 = p1 - p2;
        double    b =
            std::min({d10.squaredNorm(), d20.squaredNorm(), d12.squaredNorm()});
        if (b == 0)
            return 0; // Again: area zero triangles have surely quality==0;
        return (2 * a) / b;
    }
 
    static ScalarType qualityRadii(
        const PointT& p0,
        const PointT& p1,
        const PointT& p2)
    {
        double a = p0.dist(p1);
        double b = p2.dist(p0);
        double c = p1.dist(p2);
 
        double sum   = (a + b + c) * 0.5;
        double area2 = sum * (a + b - sum) * (a + c - sum) * (b + c - sum);
 
        if (area2 <= 0)
            return 0;
 
        return (8 * area2) / (a * b * c * sum);
    }
 
    static ScalarType qualityMeanRatio(
        const PointT& p0,
        const PointT& p1,
        const PointT& p2)
    {
        double a = p0.dist(p1);
        double b = p2.dist(p0);
        double c = p1.dist(p2);
 
        double sum   = (a + b + c) * 0.5; // semiperimeter
        double area2 = sum * (a + b - sum) * (a + c - sum) * (b + c - sum);
        if (area2 <= 0)
            return 0;
        return (4.0 * std::sqrt(3.0) * std::sqrt(area2)) /
               (a * a + b * b + c * c);
    }
};
 
/* Specialization Aliases */
 
template<typename Scalar>
using Triangle2 = Triangle<Point2<Scalar>>;
 
using Triangle2f = Triangle<Point2f>;
using Triangle2d = Triangle<Point2d>;
 
template<typename Scalar>
using Triangle3 = Triangle<Point3<Scalar>>;
 
using Triangle3f = Triangle<Point3f>;
using Triangle3d = Triangle<Point3d>;
 
} // namespace vcl
 
#endif // VCL_SPACE_CORE_TRIANGLE_H