OrcaSlicer/src/libslic3r/PrincipalComponents2D.cpp

#include "PrincipalComponents2D.hpp"
#include "Point.hpp"
#include <tuple>
#include <cstdlib>
#include "Polygon.hpp"
#include <cstddef>
#include "libslic3r.h"
#include <cmath>

namespace Slic3r {



// returns triangle area, first_moment_of_area_xy, second_moment_of_area_xy, second_moment_of_area_covariance
// none of the values is divided/normalized by area.
// The function computes intgeral over the area of the triangle, with function f(x,y) = x for first moments of area (y is analogous)
// f(x,y) = x^2 for second moment of area
// and f(x,y) = x*y for second moment of area covariance
std::tuple<float, Vec2f, Vec2f, float> compute_moments_of_area_of_triangle(const Vec2f &a, const Vec2f &b, const Vec2f &c)
{
    // based on the following guide:
    // Denote the vertices of S by a, b, c. Then the map
    //  g:(u,v)↦a+u(b−a)+v(c−a) ,
    //  which in coordinates appears as
    //  g:(u,v)↦{x(u,v)y(u,v)=a1+u(b1−a1)+v(c1−a1)=a2+u(b2−a2)+v(c2−a2) ,(1)
    //  obviously maps S′ bijectively onto S. Therefore the transformation formula for multiple integrals steps into action, and we obtain
    //  ∫Sf(x,y)d(x,y)=∫S′f(x(u,v),y(u,v))∣∣Jg(u,v)∣∣ d(u,v) .
    //  In the case at hand the Jacobian determinant is a constant: From (1) we obtain
    //  Jg(u,v)=det[xuyuxvyv]=(b1−a1)(c2−a2)−(c1−a1)(b2−a2) .
    //  Therefore we can write
    //  ∫Sf(x,y)d(x,y)=∣∣Jg∣∣∫10∫1−u0f~(u,v) dv du ,
    //  where f~ denotes the pullback of f to S′:
    //  f~(u,v):=f(x(u,v),y(u,v)) .
    //  Don't forget taking the absolute value of Jg!

    float jacobian_determinant_abs = std::abs((b.x() - a.x()) * (c.y() - a.y()) - (c.x() - a.x()) * (b.y() - a.y()));

    // coordinate transform: gx(u,v) = a.x + u * (b.x - a.x) + v * (c.x - a.x)
    // coordinate transform: gy(u,v) = a.y + u * (b.y - a.y) + v * (c.y - a.y)
    // second moment of area for x: f(x, y) = x^2;
    //              f(gx(u,v), gy(u,v)) = gx(u,v)^2 = ... (long expanded form)

    // result is Int_T func = jacobian_determinant_abs * Int_0^1 Int_0^1-u func(gx(u,v), gy(u,v)) dv du
    // integral_0^1 integral_0^(1 - u) (a + u (b - a) + v (c - a))^2 dv du = 1/12 (a^2 + a (b + c) + b^2 + b c + c^2)

    Vec2f second_moment_of_area_xy = jacobian_determinant_abs *
                                     (a.cwiseProduct(a) + b.cwiseProduct(b) + b.cwiseProduct(c) + c.cwiseProduct(c) +
                                      a.cwiseProduct(b + c)) /
                                     12.0f;
    // second moment of area covariance : f(x, y) = x*y;
    //              f(gx(u,v), gy(u,v)) = gx(u,v)*gy(u,v) = ... (long expanded form)
    //(a_1 + u * (b_1 - a_1) + v * (c_1 - a_1)) * (a_2 + u * (b_2 - a_2) + v * (c_2 - a_2))
    // ==    (a_1 + u (b_1 - a_1) + v (c_1 - a_1)) (a_2 + u (b_2 - a_2) + v (c_2 - a_2))

    // intermediate result: integral_0^(1 - u) (a_1 + u (b_1 - a_1) + v (c_1 - a_1)) (a_2 + u (b_2 - a_2) + v (c_2 - a_2)) dv =
    //  1/6 (u - 1) (-c_1 (u - 1) (a_2 (u - 1) - 3 b_2 u) - c_2 (u - 1) (a_1 (u - 1) - 3 b_1 u + 2 c_1 (u - 1)) + 3 b_1 u (a_2 (u - 1) - 2
    //  b_2 u) + a_1 (u - 1) (3 b_2 u - 2 a_2 (u - 1))) result = integral_0^1 1/6 (u - 1) (-c_1 (u - 1) (a_2 (u - 1) - 3 b_2 u) - c_2 (u -
    //  1) (a_1 (u - 1) - 3 b_1 u + 2 c_1 (u - 1)) + 3 b_1 u (a_2 (u - 1) - 2 b_2 u) + a_1 (u - 1) (3 b_2 u - 2 a_2 (u - 1))) du =
    //   1/24 (a_2 (b_1 + c_1) + a_1 (2 a_2 + b_2 + c_2) + b_2 c_1 + b_1 c_2 + 2 b_1 b_2 + 2 c_1 c_2)
    //  result is Int_T func = jacobian_determinant_abs * Int_0^1 Int_0^1-u func(gx(u,v), gy(u,v)) dv du
    float second_moment_of_area_covariance = jacobian_determinant_abs * (1.0f / 24.0f) *
                                             (a.y() * (b.x() + c.x()) + a.x() * (2.0f * a.y() + b.y() + c.y()) + b.y() * c.x() +
                                              b.x() * c.y() + 2.0f * b.x() * b.y() + 2.0f * c.x() * c.y());

    float area = jacobian_determinant_abs * 0.5f;

    Vec2f first_moment_of_area_xy = jacobian_determinant_abs * (a + b + c) / 6.0f;

    return {area, first_moment_of_area_xy, second_moment_of_area_xy, second_moment_of_area_covariance};
};

// returns two eigenvectors of the area covered by given polygons. The vectors are sorted by their corresponding eigenvalue, largest first
std::tuple<Vec2f, Vec2f> compute_principal_components(const Polygons &polys)
{
    Vec2f centroid_accumulator                         = Vec2f::Zero();
    Vec2f second_moment_of_area_accumulator            = Vec2f::Zero();
    float second_moment_of_area_covariance_accumulator = 0.0f;
    float area                                         = 0.0f;

    for (const Polygon &poly : polys) {
        Vec2f p0 = unscaled(poly.first_point()).cast<float>();
        for (size_t i = 2; i < poly.points.size(); i++) {
            Vec2f p1 = unscaled(poly.points[i - 1]).cast<float>();
            Vec2f p2 = unscaled(poly.points[i]).cast<float>();

            float sign = cross2(p1 - p0, p2 - p1) > 0 ? 1.0f : -1.0f;

            auto [triangle_area, first_moment_of_area, second_moment_area,
                  second_moment_of_area_covariance] = compute_moments_of_area_of_triangle(p0, p1, p2);
            area += sign * triangle_area;
            centroid_accumulator += sign * first_moment_of_area;
            second_moment_of_area_accumulator += sign * second_moment_area;
            second_moment_of_area_covariance_accumulator += sign * second_moment_of_area_covariance;
        }
    }

    if (area <= 0.0) {
        return {Vec2f::Zero(), Vec2f::Zero()};
    }

    Vec2f  centroid   = centroid_accumulator / area;
    Vec2f  variance   = second_moment_of_area_accumulator / area - centroid.cwiseProduct(centroid);
    double covariance = second_moment_of_area_covariance_accumulator / area - centroid.x() * centroid.y();
#if 0
        std::cout << "area : " << area << std::endl;
        std::cout << "variancex : " << variance.x() << std::endl;
        std::cout << "variancey : " << variance.y() << std::endl;
        std::cout << "covariance : " << covariance << std::endl;
#endif
    if (abs(covariance) < EPSILON) {
        std::tuple<Vec2f, Vec2f> result{Vec2f{variance.x(), 0.0}, Vec2f{0.0, variance.y()}};
        if (variance.y() > variance.x()) {
            return {std::get<1>(result), std::get<0>(result)};
        } else
            return result;
    }

    // now we find the first principal component of the covered area by computing max eigenvalue and the correspoding eigenvector of
    // covariance matrix
    //  covaraince matrix C is :  | VarX  Cov  |
    //                            | Cov   VarY |
    // Eigenvalues are solutions to det(C - lI) = 0, where l is the eigenvalue and I unit matrix
    // Eigenvector for eigenvalue l is any vector v such that Cv = lv

    float eigenvalue_a = 0.5f * (variance.x() + variance.y() +
                                 sqrt((variance.x() - variance.y()) * (variance.x() - variance.y()) + 4.0f * covariance * covariance));
    float eigenvalue_b = 0.5f * (variance.x() + variance.y() -
                                 sqrt((variance.x() - variance.y()) * (variance.x() - variance.y()) + 4.0f * covariance * covariance));
    Vec2f  eigenvector_a{(eigenvalue_a - variance.y()) / covariance, 1.0f};
    Vec2f  eigenvector_b{(eigenvalue_b - variance.y()) / covariance, 1.0f};

#if 0
        std::cout << "eigenvalue_a: " << eigenvalue_a << std::endl;
        std::cout << "eigenvalue_b: " << eigenvalue_b << std::endl;
        std::cout << "eigenvectorA: " << eigenvector_a.x() <<  "  " << eigenvector_a.y() << std::endl;
        std::cout << "eigenvectorB: " << eigenvector_b.x() <<  "  " << eigenvector_b.y() << std::endl;
#endif

    if (eigenvalue_a > eigenvalue_b) {
        return {eigenvector_a, eigenvector_b};
    } else {
        return {eigenvector_b, eigenvector_a};
    }
}

}