Method compute_triangle_voronoi_area

src/Attributes/FaceVoronoiAreaAttribute.cpp:26–66

}

VectorF FaceVoronoiAreaAttribute::compute_triangle_voronoi_area(
        Mesh& mesh, size_t face_idx) {
    size_t dim = mesh.get_dim();
    VectorI face = mesh.get_face(face_idx);

    Vector3F v0 = Vector3F::Zero();
    Vector3F v1 = Vector3F::Zero();
    Vector3F v2 = Vector3F::Zero();

    v0.segment(0, dim) = mesh.get_vertex(face[0]);
    v1.segment(0, dim) = mesh.get_vertex(face[1]);
    v2.segment(0, dim) = mesh.get_vertex(face[2]);

    Vector3F e0 = v2 - v1;
    Vector3F e1 = v0 - v2;
    Vector3F e2 = v1 - v0;

    Float double_area = e2.cross(-e1).norm();

    Float cot_0 = e2.dot(-e1) / e2.cross(-e1).norm();
    Float cot_1 = e0.dot(-e2) / e0.cross(-e2).norm();
    Float cot_2 = e1.dot(-e0) / e1.cross(-e0).norm();

    // Handle obtuse triangles.
    if (cot_0 < 0.0) {
        return Vector3F(0.5, 0.25, 0.25) * double_area * 0.5;
    } else if (cot_1 < 0.0) {
        return Vector3F(0.25, 0.5, 0.25) * double_area * 0.5;
    } else if (cot_2 < 0.0) {
        return Vector3F(0.25, 0.25, 0.5) * double_area * 0.5;
    }

    Float sq_l0 = e0.squaredNorm();
    Float sq_l1 = e1.squaredNorm();
    Float sq_l2 = e2.squaredNorm();

    return Vector3F(
            sq_l1 * cot_1 + sq_l2 * cot_2,
            sq_l0 * cot_0 + sq_l2 * cot_2,
            sq_l1 * cot_1 + sq_l0 * cot_0) / 8.0;
}