14956887918

Committed 11 May 2025 02:43PM UTC coverage: 91.097% (-0.03%) from 91.129%

Build # 14956887918

Build Type

push

github

Committed by

dsieger

Commit Message

Skip curvature smoothing iterations == 0

Run Details

2 of 2 new or added lines in 1 file covered. (100.0%)

2 existing lines in 1 file now uncovered.

5208 of 5717 relevant lines covered (91.1%)

646137.74 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

96.22

/src/pmp/algorithms/laplace.cpp

// Copyright 2011-2023 the Polygon Mesh Processing Library developers.
// Copyright 2020 Astrid Bunge, Philipp Herholz, Misha Kazhdan, Mario Botsch.
// Distributed under a MIT-style license, see LICENSE.txt for details.

#include "pmp/algorithms/laplace.h"

namespace pmp {

namespace {

// compute triangle area (in double)
double triarea(const Eigen::Vector3d& p0, const Eigen::Vector3d& p1,
               const Eigen::Vector3d& p2)
{
    return 0.5 * (p1 - p0).cross(p2 - p0).norm();
}

// compute virtual vertex per polygon, represented by affine weights,
// such that the resulting triangle fan minimizes the sum of squared triangle areas
void compute_virtual_vertex(const DenseMatrix& poly, Eigen::VectorXd& weights)
{
    const int n = poly.rows();

    // setup array of positions and edges
    std::vector<dvec3> x(n), d(n);
    for (int i = 0; i < n; ++i)
        x[i] = poly.row(i);
    for (int i = 0; i < n; ++i)
        d[i] = x[(i + 1) % n] - x[i];

    // setup matrix A and rhs b
    // see Equation (38) of "Polygon Laplacian made simple", Eurographics 2020
    DenseMatrix A(n + 1, n);
    Eigen::VectorXd b(n + 1);
    for (int j = 0; j < n; ++j)
    {
        for (int i = j; i < n; ++i)
        {
            double Aij(0.0), bi(0.0);
            for (int k = 0; k < n; ++k)
            {
                Aij += dot(cross(x[j], d[k]), cross(x[i], d[k]));
                bi += dot(cross(x[i], d[k]), cross(x[k], d[k]));
            }
            A(i, j) = A(j, i) = Aij;
            b(i) = bi;
        }
    }
    for (int j = 0; j < n; ++j)
    {
        A(n, j) = 1.0;
    }
    b(n) = 1.0;

    weights = A.completeOrthogonalDecomposition().solve(b).topRows(n);
}

void triangle_mass_matrix(const Eigen::Vector3d& p0, const Eigen::Vector3d& p1,
                          const Eigen::Vector3d& p2, DiagonalMatrix& Mtri)
{
    // three vertex positions
    const std::array<dvec3, 3> p = {p0, p1, p2};

    // edge vectors
    std::array<dvec3, 3> e;
    for (int i = 0; i < 3; ++i)
        e[i] = p[(i + 1) % 3] - p[i];

    // compute and check (twice the) triangle area
    const auto tri_area = norm(cross(e[0], e[1]));
    if (tri_area <= std::numeric_limits<double>::min())
    {
        Mtri.setZero(3);
        return;
    }

    // dot products for each corner (of its two emanating edge vectors)
    std::array<double, 3> d;
    for (int i = 0; i < 3; ++i)
        d[i] = -dot(e[i], e[(i + 2) % 3]);

    // cotangents for each corner: cot = cos/sin = dot(A,B)/norm(cross(A,B))
    std::array<double, 3> cot;
    for (int i = 0; i < 3; ++i)
        cot[i] = d[i] / tri_area;

    // compute area for each corner
    Eigen::Vector3d area;
    for (int i = 0; i < 3; ++i)
    {
        // angle at corner is obtuse
        if (d[i] < 0.0)
        {
            area[i] = 0.25 * tri_area;
        }
        // angle at some other corner is obtuse
        else if (d[(i + 1) % 3] < 0.0 || d[(i + 2) % 3] < 0.0)
        {
            area[i] = 0.125 * tri_area;
        }
        // no obtuse angles
        else
        {
            area[i] = 0.125 * (sqrnorm(e[i]) * cot[(i + 2) % 3] +
                               sqrnorm(e[(i + 2) % 3]) * cot[(i + 1) % 3]);
        }
    }

    Mtri = area.asDiagonal();
}

void polygon_mass_matrix(const DenseMatrix& polygon, DiagonalMatrix& Mpoly)
{
    const int n = (int)polygon.rows();

    // shortcut for triangles
    if (n == 3)
    {
        triangle_mass_matrix(polygon.row(0), polygon.row(1), polygon.row(2),
                             Mpoly);
        return;
    }

    // compute position of virtual vertex
    Eigen::VectorXd vweights;
    compute_virtual_vertex(polygon, vweights);
    const Eigen::Vector3d vvertex = polygon.transpose() * vweights;

    // laplace matrix of refined triangle fan
    DenseMatrix Mfan = DenseMatrix::Zero(n + 1, n + 1);
    DiagonalMatrix Mtri;
    for (int i = 0; i < n; ++i)
    {
        const int j = (i + 1) % n;

        // build laplace matrix of one triangle
        triangle_mass_matrix(polygon.row(i), polygon.row(j), vvertex, Mtri);

        // assemble to laplace matrix for refined triangle fan
        // (we are dealing with diagonal matrices)
        Mfan.diagonal()[i] += Mtri.diagonal()[0];
        Mfan.diagonal()[j] += Mtri.diagonal()[1];
        Mfan.diagonal()[n] += Mtri.diagonal()[2];
    }

    // build prolongation matrix
    DenseMatrix P(n + 1, n);
    P.setIdentity();
    P.row(n) = vweights;

    // build polygon laplace matrix by sandwiching
    DenseMatrix PMP = P.transpose() * Mfan * P;
    Mpoly = PMP.rowwise().sum().asDiagonal();
}

void triangle_laplace_matrix(const Eigen::Vector3d& p0,
                             const Eigen::Vector3d& p1,
                             const Eigen::Vector3d& p2, DenseMatrix& Ltri)
{
    std::array<double, 3> l, l2, cot;

    // squared edge lengths
    l2[0] = (p1 - p2).squaredNorm();
    l2[1] = (p0 - p2).squaredNorm();
    l2[2] = (p0 - p1).squaredNorm();

    // edge lengths
    l[0] = sqrt(l2[0]);
    l[1] = sqrt(l2[1]);
    l[2] = sqrt(l2[2]);

    // triangle area
    const double arg = (l[0] + (l[1] + l[2])) * (l[2] - (l[0] - l[1])) *
                       (l[2] + (l[0] - l[1])) * (l[0] + (l[1] - l[2]));
    const double area = 0.5 * sqrt(arg);

    if (area > std::numeric_limits<double>::min())
    {
        // cotangents of angles
        cot[0] = 0.25 * (l2[1] + l2[2] - l2[0]) / area;
        cot[1] = 0.25 * (l2[2] + l2[0] - l2[1]) / area;
        cot[2] = 0.25 * (l2[0] + l2[1] - l2[2]) / area;

        Ltri(0, 0) = cot[1] + cot[2];
        Ltri(1, 1) = cot[0] + cot[2];
        Ltri(2, 2) = cot[0] + cot[1];
        Ltri(1, 0) = Ltri(0, 1) = -cot[2];
        Ltri(2, 0) = Ltri(0, 2) = -cot[1];
        Ltri(2, 1) = Ltri(1, 2) = -cot[0];
    }
}

void polygon_laplace_matrix(const DenseMatrix& polygon, DenseMatrix& Lpoly)
{
    const int n = (int)polygon.rows();
    Lpoly = DenseMatrix::Zero(n, n);

    // shortcut for triangles
    if (n == 3)
    {
        triangle_laplace_matrix(polygon.row(0), polygon.row(1), polygon.row(2),
                                Lpoly);
        return;
    }

    // compute position of virtual vertex
    Eigen::VectorXd vweights;
    compute_virtual_vertex(polygon, vweights);
    const Eigen::Vector3d vvertex = polygon.transpose() * vweights;

    // laplace matrix of refined triangle fan
    DenseMatrix Lfan = DenseMatrix::Zero(n + 1, n + 1);
    DenseMatrix Ltri(3, 3);
    for (int i = 0; i < n; ++i)
    {
        const int j = (i + 1) % n;

        // build laplace matrix of one triangle
        triangle_laplace_matrix(polygon.row(i), polygon.row(j), vvertex, Ltri);

        // assemble to laplace matrix for refined triangle fan
        Lfan(i, i) += Ltri(0, 0);
        Lfan(i, j) += Ltri(0, 1);
        Lfan(i, n) += Ltri(0, 2);
        Lfan(j, i) += Ltri(1, 0);
        Lfan(j, j) += Ltri(1, 1);
        Lfan(j, n) += Ltri(1, 2);
        Lfan(n, i) += Ltri(2, 0);
        Lfan(n, j) += Ltri(2, 1);
        Lfan(n, n) += Ltri(2, 2);
    }

    // build prolongation matrix
    DenseMatrix P(n + 1, n);
    P.setIdentity();
    P.row(n) = vweights;

    // build polygon laplace matrix by sandwiching
    Lpoly = P.transpose() * Lfan * P;
}

void triangle_gradient_matrix(const Eigen::Vector3d& p0,
                              const Eigen::Vector3d& p1,
                              const Eigen::Vector3d& p2, DenseMatrix& G)
{
    G.resize(3, 3);
    Eigen::Vector3d n = (p1 - p0).cross(p2 - p0);
    n /= n.squaredNorm();
    G.col(0) = n.cross(p2 - p1);
    G.col(1) = n.cross(p0 - p2);
    G.col(2) = n.cross(p1 - p0);
}

void polygon_gradient_matrix(const DenseMatrix& polygon, DenseMatrix& Gpoly)
{
    const int n = (int)polygon.rows();

    // compute position of virtual vertex
    Eigen::VectorXd vweights;
    compute_virtual_vertex(polygon, vweights);
    const Eigen::Vector3d vvertex = polygon.transpose() * vweights;

    DenseMatrix Gfan = DenseMatrix::Zero(3 * n, n + 1);
    DenseMatrix Gtri(3, 3);
    int row = 0;
    for (int i = 0; i < n; ++i)
    {
        const int j = (i + 1) % n;

        // build laplace matrix of one triangle
        triangle_gradient_matrix(polygon.row(i), polygon.row(j), vvertex, Gtri);

        // assemble to matrix for triangle fan
        for (int k = 0; k < 3; ++k)
        {
            Gfan(row + k, i) += Gtri(k, 0);
            Gfan(row + k, j) += Gtri(k, 1);
            Gfan(row + k, n) += Gtri(k, 2);
        }

        row += 3;
    }

    // build prolongation matrix
    DenseMatrix P(n + 1, n);
    P.setIdentity();
    P.row(n) = vweights;

    // build polygon gradient matrix by sandwiching (from left only)
    Gpoly = Gfan * P;
}

void divmass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)
{
    // how many virtual triangles will we have after refinement?
    unsigned int nt = 0;
    for (auto f : mesh.faces())
    {
        nt += mesh.valence(f);
    }

    // initialize global matrix
    M.resize(3 * nt);
    auto& diag = M.diagonal();

    std::vector<Vertex> vertices; // polygon vertices
    DenseMatrix polygon;          // positions of polygon vertices

    unsigned int idx = 0;

    for (const auto f : mesh.faces())
    {
        // collect polygon vertices
        vertices.clear();
        for (const auto v : mesh.vertices(f))
        {
            vertices.push_back(v);
        }
        const int n = vertices.size();

        // collect their positions
        polygon.resize(n, 3);
        for (int i = 0; i < n; ++i)
        {
            polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);
        }

        // compute position of virtual vertex
        Eigen::VectorXd vweights;
        compute_virtual_vertex(polygon, vweights);
        const Eigen::Vector3d vvertex = polygon.transpose() * vweights;

        for (int i = 0; i < n; ++i)
        {
            const double area =
                triarea(polygon.row(i), polygon.row((i + 1) % n), vvertex);

            diag[idx++] = area;
            diag[idx++] = area;
            diag[idx++] = area;
        }
    }

    assert(idx == 3 * nt);
}

} // anonymous namespace

void uniform_mass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)
{
    const unsigned int n = mesh.n_vertices();
    Eigen::VectorXd diag(n);
    for (auto v : mesh.vertices())
        diag[v.idx()] = mesh.valence(v);
    M = diag.asDiagonal();
}

void mass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)
{
    const int nv = mesh.n_vertices();
    std::vector<Vertex> vertices; // polygon vertices
    DenseMatrix polygon;          // positions of polygon vertices
    DiagonalMatrix Mpoly;         // local mass matrix

    M.setZero(nv);

    for (const auto f : mesh.faces())
    {
        // collect polygon vertices
        vertices.clear();
        for (const auto v : mesh.vertices(f))
        {
            vertices.push_back(v);
        }
        const int n = vertices.size();

        // collect their positions
        polygon.resize(n, 3);
        for (int i = 0; i < n; ++i)
        {
            polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);
        }

        // setup local mass matrix
        polygon_mass_matrix(polygon, Mpoly);

        // assemble to global mass matrix
        for (int k = 0; k < n; ++k)
        {
            M.diagonal()[vertices[k].idx()] += Mpoly.diagonal()[k];
        }
    }
}

void uniform_laplace_matrix(const SurfaceMesh& mesh, SparseMatrix& L)
{
    const unsigned int n = mesh.n_vertices();

    std::vector<Triplet> triplets;
    triplets.reserve(8 * n); // conservative estimate for triangle meshes

    for (auto vi : mesh.vertices())
    {
        Scalar sum_weights = 0.0;
        for (auto vj : mesh.vertices(vi))
        {
            sum_weights += 1.0;
            triplets.emplace_back(vi.idx(), vj.idx(), 1.0);
        }
        triplets.emplace_back(vi.idx(), vi.idx(), -sum_weights);
    }

    L.resize(n, n);
    L.setFromTriplets(triplets.begin(), triplets.end());
}

void laplace_matrix(const SurfaceMesh& mesh, SparseMatrix& L, bool clamp)
{
    const int nv = mesh.n_vertices();
    const int nf = mesh.n_faces();
    std::vector<Vertex> vertices; // polygon vertices
    DenseMatrix polygon;          // positions of polygon vertices
    DenseMatrix Lpoly;            // local laplace matrix

    std::vector<Triplet> triplets;
    triplets.reserve(9 * nf); // estimate for triangle meshes

    for (const auto f : mesh.faces())
    {
        // collect polygon vertices
        vertices.clear();
        for (const auto v : mesh.vertices(f))
        {
            vertices.push_back(v);
        }
        const int n = vertices.size();

        // collect their positions
        polygon.resize(n, 3);
        for (int i = 0; i < n; ++i)
        {
            polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);
        }

        // setup local laplace matrix
        polygon_laplace_matrix(polygon, Lpoly);

        // assemble to global laplace matrix
        for (int j = 0; j < n; ++j)
        {
            for (int k = 0; k < n; ++k)
            {
                triplets.emplace_back(vertices[k].idx(), vertices[j].idx(),
                                      -Lpoly(k, j));
            }
        }
    }

    // build sparse matrix from triplets
    L.resize(nv, nv);
    L.setFromTriplets(triplets.begin(), triplets.end());

    // clamp negative off-diagonal entries to zero
    if (clamp)
    {
        for (unsigned int k = 0; k < L.outerSize(); k++)
        {
            double diag_offset(0.0);

            for (SparseMatrix::InnerIterator iter(L, k); iter; ++iter)
            {
                if (iter.row() != iter.col() && iter.value() < 0.0)
                {
                    diag_offset += -iter.value();
                    iter.valueRef() = 0.0;
                }
            }
            for (SparseMatrix::InnerIterator iter(L, k); iter; ++iter)
            {
                if (iter.row() == iter.col() && iter.value() < 0.0)
                    iter.valueRef() -= diag_offset;
            }
        }
    }
}

void gradient_matrix(const SurfaceMesh& mesh, SparseMatrix& G)
{
    const int nv = mesh.n_vertices();

    // how many virtual triangles will we have after refinement?
    // triangles are not refined, other polygons are.
    unsigned int nt = 0;
    for (auto f : mesh.faces())
    {
        nt += mesh.valence(f);
    }

    std::vector<Vertex> vertices; // polygon vertices
    DenseMatrix polygon;          // positions of polygon vertices
    DenseMatrix Gpoly;            // local gradient matrix

    std::vector<Triplet> triplets;
    triplets.reserve(9 * nt);

    unsigned int n_rows = 0;

    for (const auto f : mesh.faces())
    {
        // collect polygon vertices
        vertices.clear();
        for (const auto v : mesh.vertices(f))
        {
            vertices.push_back(v);
        }
        const int n = vertices.size();

        // collect their positions
        polygon.resize(n, 3);
        for (int i = 0; i < n; ++i)
        {
            polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);
        }

        // setup local element matrix
        polygon_gradient_matrix(polygon, Gpoly);

        // assemble to global matrix
        for (int j = 0; j < Gpoly.cols(); ++j)
        {
            for (int i = 0; i < Gpoly.rows(); ++i)
            {
                triplets.emplace_back(n_rows + i, vertices[j].idx(),
                                      Gpoly(i, j));
            }
        }

        n_rows += Gpoly.rows();
    }
    assert(n_rows == 3 * nt);

    // build sparse matrix from triplets
    G.resize(n_rows, nv);
    G.setFromTriplets(triplets.begin(), triplets.end());
}

void divergence_matrix(const SurfaceMesh& mesh, SparseMatrix& D)
{
    SparseMatrix G;
    gradient_matrix(mesh, G);
    DiagonalMatrix M;
    divmass_matrix(mesh, M);
    D = -G.transpose() * M;
}

void coordinates_to_matrix(const SurfaceMesh& mesh, DenseMatrix& X)
{
    X.resize(mesh.n_vertices(), 3);
    for (auto v : mesh.vertices())
        X.row(v.idx()) = static_cast<Eigen::Vector3d>(mesh.position(v));
}

void matrix_to_coordinates(const DenseMatrix& X, SurfaceMesh& mesh)
{
    assert((size_t)X.rows() == mesh.n_vertices() && X.cols() == 3);
    for (auto v : mesh.vertices())
        mesh.position(v) = X.row(v.idx());
}

} // namespace pmp

1	// Copyright 2011-2023 the Polygon Mesh Processing Library developers.
2	// Copyright 2020 Astrid Bunge, Philipp Herholz, Misha Kazhdan, Mario Botsch.
3	// Distributed under a MIT-style license, see LICENSE.txt for details.
4
5	#include "pmp/algorithms/laplace.h"
6
7	namespace pmp {
8
9	namespace {
10
11	// compute triangle area (in double)
12	double triarea(const Eigen::Vector3d& p0, const Eigen::Vector3d& p1,	10,134✔
13	const Eigen::Vector3d& p2)
14	{
15	return 0.5 * (p1 - p0).cross(p2 - p0).norm();	10,134✔
16	}
17
18	// compute virtual vertex per polygon, represented by affine weights,
19	// such that the resulting triangle fan minimizes the sum of squared triangle areas
20	void compute_virtual_vertex(const DenseMatrix& poly, Eigen::VectorXd& weights)	10,745✔
21	{
22	const int n = poly.rows();	10,745✔
23
24	// setup array of positions and edges
25	std::vector<dvec3> x(n), d(n);	32,235✔
26	for (int i = 0; i < n; ++i)	45,907✔
27	x[i] = poly.row(i);	35,162✔
28	for (int i = 0; i < n; ++i)	45,907✔
29	d[i] = x[(i + 1) % n] - x[i];	35,162✔
30
31	// setup matrix A and rhs b
32	// see Equation (38) of "Polygon Laplacian made simple", Eurographics 2020
33	DenseMatrix A(n + 1, n);	10,745✔
34	Eigen::VectorXd b(n + 1);	10,745✔
35	for (int j = 0; j < n; ++j)	45,907✔
36	{
37	for (int i = j; i < n; ++i)	111,340✔
38	{
39	double Aij(0.0), bi(0.0);	76,178✔
40	for (int k = 0; k < n; ++k)	333,982✔
41	{
42	Aij += dot(cross(x[j], d[k]), cross(x[i], d[k]));	257,804✔
43	bi += dot(cross(x[i], d[k]), cross(x[k], d[k]));	257,804✔
44	}
45	A(i, j) = A(j, i) = Aij;	76,178✔
46	b(i) = bi;	76,178✔
47	}
48	}
49	for (int j = 0; j < n; ++j)	45,907✔
50	{
51	A(n, j) = 1.0;	35,162✔
52	}
53	b(n) = 1.0;	10,745✔
54
55	weights = A.completeOrthogonalDecomposition().solve(b).topRows(n);	10,745✔
56	}	10,745✔
57
58	void triangle_mass_matrix(const Eigen::Vector3d& p0, const Eigen::Vector3d& p1,	107,624✔
59	const Eigen::Vector3d& p2, DiagonalMatrix& Mtri)
60	{
61	// three vertex positions
62	const std::array<dvec3, 3> p = {p0, p1, p2};	107,624✔
63
64	// edge vectors
65	std::array<dvec3, 3> e;
66	for (int i = 0; i < 3; ++i)	430,496✔
67	e[i] = p[(i + 1) % 3] - p[i];	322,872✔
68
69	// compute and check (twice the) triangle area
70	const auto tri_area = norm(cross(e[0], e[1]));	107,624✔
71	if (tri_area <= std::numeric_limits<double>::min())	107,624✔
72	{
73	Mtri.setZero(3);	×
74	return;	×
75	}
76
77	// dot products for each corner (of its two emanating edge vectors)
78	std::array<double, 3> d;
79	for (int i = 0; i < 3; ++i)	430,496✔
80	d[i] = -dot(e[i], e[(i + 2) % 3]);	322,872✔
81
82	// cotangents for each corner: cot = cos/sin = dot(A,B)/norm(cross(A,B))
83	std::array<double, 3> cot;
84	for (int i = 0; i < 3; ++i)	430,496✔
85	cot[i] = d[i] / tri_area;	322,872✔
86
87	// compute area for each corner
88	Eigen::Vector3d area;	107,624✔
89	for (int i = 0; i < 3; ++i)	430,496✔
90	{
91	// angle at corner is obtuse
92	if (d[i] < 0.0)	322,872✔
93	{
94	area[i] = 0.25 * tri_area;	1,214✔
95	}
96	// angle at some other corner is obtuse
97	else if (d[(i + 1) % 3] < 0.0 \|\| d[(i + 2) % 3] < 0.0)	321,658✔
98	{
99	area[i] = 0.125 * tri_area;	2,428✔
100	}
101	// no obtuse angles
102	else
103	{
104	area[i] = 0.125 * (sqrnorm(e[i]) * cot[(i + 2) % 3] +	638,460✔
105	sqrnorm(e[(i + 2) % 3]) * cot[(i + 1) % 3]);	319,230✔
106	}
107	}
108
109	Mtri = area.asDiagonal();	107,624✔
110	}
111
112	void polygon_mass_matrix(const DenseMatrix& polygon, DiagonalMatrix& Mpoly)	105,887✔
113	{
114	const int n = (int)polygon.rows();	105,887✔
115
116	// shortcut for triangles
117	if (n == 3)	105,887✔
118	{
119	triangle_mass_matrix(polygon.row(0), polygon.row(1), polygon.row(2),	105,308✔
120	Mpoly);
121	return;	105,308✔
122	}
123
124	// compute position of virtual vertex
125	Eigen::VectorXd vweights;	579✔
126	compute_virtual_vertex(polygon, vweights);	579✔
127	const Eigen::Vector3d vvertex = polygon.transpose() * vweights;	579✔
128
129	// laplace matrix of refined triangle fan
130	DenseMatrix Mfan = DenseMatrix::Zero(n + 1, n + 1);	579✔
131	DiagonalMatrix Mtri;	579✔
132	for (int i = 0; i < n; ++i)	2,895✔
133	{
134	const int j = (i + 1) % n;	2,316✔
135
136	// build laplace matrix of one triangle
137	triangle_mass_matrix(polygon.row(i), polygon.row(j), vvertex, Mtri);	2,316✔
138
139	// assemble to laplace matrix for refined triangle fan
140	// (we are dealing with diagonal matrices)
141	Mfan.diagonal()[i] += Mtri.diagonal()[0];	2,316✔
142	Mfan.diagonal()[j] += Mtri.diagonal()[1];	2,316✔
143	Mfan.diagonal()[n] += Mtri.diagonal()[2];	2,316✔
144	}
145
146	// build prolongation matrix
147	DenseMatrix P(n + 1, n);	579✔
148	P.setIdentity();	579✔
149	P.row(n) = vweights;	579✔
150
151	// build polygon laplace matrix by sandwiching
152	DenseMatrix PMP = P.transpose() * Mfan * P;	579✔
153	Mpoly = PMP.rowwise().sum().asDiagonal();	579✔
154	}	579✔
155
156	void triangle_laplace_matrix(const Eigen::Vector3d& p0,	209,932✔
157	const Eigen::Vector3d& p1,
158	const Eigen::Vector3d& p2, DenseMatrix& Ltri)
159	{
160	std::array<double, 3> l, l2, cot;
161
162	// squared edge lengths
163	l2[0] = (p1 - p2).squaredNorm();	209,932✔
164	l2[1] = (p0 - p2).squaredNorm();	209,932✔
165	l2[2] = (p0 - p1).squaredNorm();	209,932✔
166
167	// edge lengths
168	l[0] = sqrt(l2[0]);	209,932✔
169	l[1] = sqrt(l2[1]);	209,932✔
170	l[2] = sqrt(l2[2]);	209,932✔
171
172	// triangle area
173	const double arg = (l[0] + (l[1] + l[2])) * (l[2] - (l[0] - l[1])) *	209,932✔
174	(l[2] + (l[0] - l[1])) * (l[0] + (l[1] - l[2]));	209,932✔
175	const double area = 0.5 * sqrt(arg);	209,932✔
176
177	if (area > std::numeric_limits<double>::min())	209,932✔
178	{
179	// cotangents of angles
180	cot[0] = 0.25 * (l2[1] + l2[2] - l2[0]) / area;	209,932✔
181	cot[1] = 0.25 * (l2[2] + l2[0] - l2[1]) / area;	209,932✔
182	cot[2] = 0.25 * (l2[0] + l2[1] - l2[2]) / area;	209,932✔
183
184	Ltri(0, 0) = cot[1] + cot[2];	209,932✔
185	Ltri(1, 1) = cot[0] + cot[2];	209,932✔
186	Ltri(2, 2) = cot[0] + cot[1];	209,932✔
187	Ltri(1, 0) = Ltri(0, 1) = -cot[2];	209,932✔
188	Ltri(2, 0) = Ltri(0, 2) = -cot[1];	209,932✔
189	Ltri(2, 1) = Ltri(1, 2) = -cot[0];	209,932✔
190	}
191	}	209,932✔
192
193	void polygon_laplace_matrix(const DenseMatrix& polygon, DenseMatrix& Lpoly)	208,147✔
194	{
195	const int n = (int)polygon.rows();	208,147✔
196	Lpoly = DenseMatrix::Zero(n, n);	208,147✔
197
198	// shortcut for triangles
199	if (n == 3)	208,147✔
200	{
201	triangle_laplace_matrix(polygon.row(0), polygon.row(1), polygon.row(2),	207,552✔
202	Lpoly);
203	return;	207,552✔
204	}
205
206	// compute position of virtual vertex
207	Eigen::VectorXd vweights;	595✔
208	compute_virtual_vertex(polygon, vweights);	595✔
209	const Eigen::Vector3d vvertex = polygon.transpose() * vweights;	595✔
210
211	// laplace matrix of refined triangle fan
212	DenseMatrix Lfan = DenseMatrix::Zero(n + 1, n + 1);	595✔
213	DenseMatrix Ltri(3, 3);	595✔
214	for (int i = 0; i < n; ++i)	2,975✔
215	{
216	const int j = (i + 1) % n;	2,380✔
217
218	// build laplace matrix of one triangle
219	triangle_laplace_matrix(polygon.row(i), polygon.row(j), vvertex, Ltri);	2,380✔
220
221	// assemble to laplace matrix for refined triangle fan
222	Lfan(i, i) += Ltri(0, 0);	2,380✔
223	Lfan(i, j) += Ltri(0, 1);	2,380✔
224	Lfan(i, n) += Ltri(0, 2);	2,380✔
225	Lfan(j, i) += Ltri(1, 0);	2,380✔
226	Lfan(j, j) += Ltri(1, 1);	2,380✔
227	Lfan(j, n) += Ltri(1, 2);	2,380✔
228	Lfan(n, i) += Ltri(2, 0);	2,380✔
229	Lfan(n, j) += Ltri(2, 1);	2,380✔
230	Lfan(n, n) += Ltri(2, 2);	2,380✔
231	}
232
233	// build prolongation matrix
234	DenseMatrix P(n + 1, n);	595✔
235	P.setIdentity();	595✔
236	P.row(n) = vweights;	595✔
237
238	// build polygon laplace matrix by sandwiching
239	Lpoly = P.transpose() * Lfan * P;	595✔
240	}	595✔
241
242	void triangle_gradient_matrix(const Eigen::Vector3d& p0,	20,332✔
243	const Eigen::Vector3d& p1,
244	const Eigen::Vector3d& p2, DenseMatrix& G)
245	{
246	G.resize(3, 3);	20,332✔
247	Eigen::Vector3d n = (p1 - p0).cross(p2 - p0);	20,332✔
248	n /= n.squaredNorm();	20,332✔
249	G.col(0) = n.cross(p2 - p1);	20,332✔
250	G.col(1) = n.cross(p0 - p2);	20,332✔
251	G.col(2) = n.cross(p1 - p0);	20,332✔
252	}	20,332✔
253
254	void polygon_gradient_matrix(const DenseMatrix& polygon, DenseMatrix& Gpoly)	6,386✔
255	{
256	const int n = (int)polygon.rows();	6,386✔
257
258	// compute position of virtual vertex
259	Eigen::VectorXd vweights;	6,386✔
260	compute_virtual_vertex(polygon, vweights);	6,386✔
261	const Eigen::Vector3d vvertex = polygon.transpose() * vweights;	6,386✔
262
263	DenseMatrix Gfan = DenseMatrix::Zero(3 * n, n + 1);	6,386✔
264	DenseMatrix Gtri(3, 3);	6,386✔
265	int row = 0;	6,386✔
266	for (int i = 0; i < n; ++i)	26,718✔
267	{
268	const int j = (i + 1) % n;	20,332✔
269
270	// build laplace matrix of one triangle
271	triangle_gradient_matrix(polygon.row(i), polygon.row(j), vvertex, Gtri);	20,332✔
272
273	// assemble to matrix for triangle fan
274	for (int k = 0; k < 3; ++k)	81,328✔
275	{
276	Gfan(row + k, i) += Gtri(k, 0);	60,996✔
277	Gfan(row + k, j) += Gtri(k, 1);	60,996✔
278	Gfan(row + k, n) += Gtri(k, 2);	60,996✔
279	}
280
281	row += 3;	20,332✔
282	}
283
284	// build prolongation matrix
285	DenseMatrix P(n + 1, n);	6,386✔
286	P.setIdentity();	6,386✔
287	P.row(n) = vweights;	6,386✔
288
289	// build polygon gradient matrix by sandwiching (from left only)
290	Gpoly = Gfan * P;	6,386✔
291	}	6,386✔
292
293	void divmass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)	6✔
294	{
295	// how many virtual triangles will we have after refinement?
296	unsigned int nt = 0;	6✔
297	for (auto f : mesh.faces())	3,191✔
298	{
299	nt += mesh.valence(f);	3,185✔
300	}
301
302	// initialize global matrix
303	M.resize(3 * nt);	6✔
304	auto& diag = M.diagonal();	6✔
305
306	std::vector<Vertex> vertices; // polygon vertices	6✔
307	DenseMatrix polygon; // positions of polygon vertices	6✔
308
309	unsigned int idx = 0;	6✔
310
311	for (const auto f : mesh.faces())	6,376✔
312	{
313	// collect polygon vertices
314	vertices.clear();	3,185✔
315	for (const auto v : mesh.vertices(f))	23,453✔
316	{
317	vertices.push_back(v);	10,134✔
318	}
319	const int n = vertices.size();	3,185✔
320
321	// collect their positions
322	polygon.resize(n, 3);	3,185✔
323	for (int i = 0; i < n; ++i)	13,319✔
324	{
325	polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);	10,134✔
326	}
327
328	// compute position of virtual vertex
329	Eigen::VectorXd vweights;	3,185✔
330	compute_virtual_vertex(polygon, vweights);	3,185✔
331	const Eigen::Vector3d vvertex = polygon.transpose() * vweights;	3,185✔
332
333	for (int i = 0; i < n; ++i)	13,319✔
334	{
335	const double area =
336	triarea(polygon.row(i), polygon.row((i + 1) % n), vvertex);	10,134✔
337
338	diag[idx++] = area;	10,134✔
339	diag[idx++] = area;	10,134✔
340	diag[idx++] = area;	10,134✔
341	}
342	}	3,185✔
343
344	assert(idx == 3 * nt);	6✔
345	}	6✔
346
347	} // anonymous namespace
348
349	void uniform_mass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)	×
350	{
351	const unsigned int n = mesh.n_vertices();	×
352	Eigen::VectorXd diag(n);	×
353	for (auto v : mesh.vertices())	×
354	diag[v.idx()] = mesh.valence(v);	×
355	M = diag.asDiagonal();	×
356	}	×
357
358	void mass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)	19✔
359	{
360	const int nv = mesh.n_vertices();	19✔
361	std::vector<Vertex> vertices; // polygon vertices	19✔
362	DenseMatrix polygon; // positions of polygon vertices	19✔
363	DiagonalMatrix Mpoly; // local mass matrix	19✔
364
365	M.setZero(nv);	19✔
366
367	for (const auto f : mesh.faces())	211,793✔
368	{
369	// collect polygon vertices
370	vertices.clear();	105,887✔
371	for (const auto v : mesh.vertices(f))	742,367✔
372	{
373	vertices.push_back(v);	318,240✔
374	}
375	const int n = vertices.size();	105,887✔
376
377	// collect their positions
378	polygon.resize(n, 3);	105,887✔
379	for (int i = 0; i < n; ++i)	424,127✔
380	{
381	polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);	318,240✔
382	}
383
384	// setup local mass matrix
385	polygon_mass_matrix(polygon, Mpoly);	105,887✔
386
387	// assemble to global mass matrix
388	for (int k = 0; k < n; ++k)	424,127✔
389	{
390	M.diagonal()[vertices[k].idx()] += Mpoly.diagonal()[k];	318,240✔
391	}
392	}
393	}	19✔
394
395	void uniform_laplace_matrix(const SurfaceMesh& mesh, SparseMatrix& L)	2✔
396	{
397	const unsigned int n = mesh.n_vertices();	2✔
398
399	std::vector<Triplet> triplets;	2✔
400	triplets.reserve(8 * n); // conservative estimate for triangle meshes	2✔
401
402	for (auto vi : mesh.vertices())	20✔
403	{
404	Scalar sum_weights = 0.0;	18✔
405	for (auto vj : mesh.vertices(vi))	82✔
406	{
407	sum_weights += 1.0;	64✔
408	triplets.emplace_back(vi.idx(), vj.idx(), 1.0);	64✔
409	}
410	triplets.emplace_back(vi.idx(), vi.idx(), -sum_weights);	18✔
411	}
412
413	L.resize(n, n);	2✔
414	L.setFromTriplets(triplets.begin(), triplets.end());	2✔
415	}	2✔
416
417	void laplace_matrix(const SurfaceMesh& mesh, SparseMatrix& L, bool clamp)	25✔
418	{
419	const int nv = mesh.n_vertices();	25✔
420	const int nf = mesh.n_faces();	25✔
421	std::vector<Vertex> vertices; // polygon vertices	25✔
422	DenseMatrix polygon; // positions of polygon vertices	25✔
423	DenseMatrix Lpoly; // local laplace matrix	25✔
424
425	std::vector<Triplet> triplets;	25✔
426	triplets.reserve(9 * nf); // estimate for triangle meshes	25✔
427
428	for (const auto f : mesh.faces())	416,319✔
429	{
430	// collect polygon vertices
431	vertices.clear();	208,147✔
432	for (const auto v : mesh.vertices(f))	1,458,219✔
433	{
434	vertices.push_back(v);	625,036✔
435	}
436	const int n = vertices.size();	208,147✔
437
438	// collect their positions
439	polygon.resize(n, 3);	208,147✔
440	for (int i = 0; i < n; ++i)	833,183✔
441	{
442	polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);	625,036✔
443	}
444
445	// setup local laplace matrix
446	polygon_laplace_matrix(polygon, Lpoly);	208,147✔
447
448	// assemble to global laplace matrix
449	for (int j = 0; j < n; ++j)	833,183✔
450	{
451	for (int k = 0; k < n; ++k)	2,502,524✔
452	{
453	triplets.emplace_back(vertices[k].idx(), vertices[j].idx(),	1,877,488✔
454	-Lpoly(k, j));	3,754,976✔
455	}
456	}
457	}
458
459	// build sparse matrix from triplets
460	L.resize(nv, nv);	25✔
461	L.setFromTriplets(triplets.begin(), triplets.end());	25✔
462
463	// clamp negative off-diagonal entries to zero
464	if (clamp)	25✔
465	{
466	for (unsigned int k = 0; k < L.outerSize(); k++)	51,299✔
467	{
468	double diag_offset(0.0);	51,289✔
469
470	for (SparseMatrix::InnerIterator iter(L, k); iter; ++iter)	410,120✔
471	{
472	if (iter.row() != iter.col() && iter.value() < 0.0)	358,831✔
473	{
UNCOV 474	diag_offset += -iter.value();	×
UNCOV 475	iter.valueRef() = 0.0;	×
476	}
477	}
478	for (SparseMatrix::InnerIterator iter(L, k); iter; ++iter)	410,120✔
479	{
480	if (iter.row() == iter.col() && iter.value() < 0.0)	358,831✔
481	iter.valueRef() -= diag_offset;	51,289✔
482	}
483	}
484	}
485	}	25✔
486
487	void gradient_matrix(const SurfaceMesh& mesh, SparseMatrix& G)	13✔
488	{
489	const int nv = mesh.n_vertices();	13✔
490
491	// how many virtual triangles will we have after refinement?
492	// triangles are not refined, other polygons are.
493	unsigned int nt = 0;	13✔
494	for (auto f : mesh.faces())	6,399✔
495	{
496	nt += mesh.valence(f);	6,386✔
497	}
498
499	std::vector<Vertex> vertices; // polygon vertices	13✔
500	DenseMatrix polygon; // positions of polygon vertices	13✔
501	DenseMatrix Gpoly; // local gradient matrix	13✔
502
503	std::vector<Triplet> triplets;	13✔
504	triplets.reserve(9 * nt);	13✔
505
506	unsigned int n_rows = 0;	13✔
507
508	for (const auto f : mesh.faces())	6,399✔
509	{
510	// collect polygon vertices
511	vertices.clear();	6,386✔
512	for (const auto v : mesh.vertices(f))	47,050✔
513	{
514	vertices.push_back(v);	20,332✔
515	}
516	const int n = vertices.size();	6,386✔
517
518	// collect their positions
519	polygon.resize(n, 3);	6,386✔
520	for (int i = 0; i < n; ++i)	26,718✔
521	{
522	polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);	20,332✔
523	}
524
525	// setup local element matrix
526	polygon_gradient_matrix(polygon, Gpoly);	6,386✔
527
528	// assemble to global matrix
529	for (int j = 0; j < Gpoly.cols(); ++j)	26,718✔
530	{
531	for (int i = 0; i < Gpoly.rows(); ++i)	217,408✔
532	{
533	triplets.emplace_back(n_rows + i, vertices[j].idx(),	197,076✔
534	Gpoly(i, j));	197,076✔
535	}
536	}
537
538	n_rows += Gpoly.rows();	6,386✔
539	}
540	assert(n_rows == 3 * nt);	13✔
541
542	// build sparse matrix from triplets
543	G.resize(n_rows, nv);	13✔
544	G.setFromTriplets(triplets.begin(), triplets.end());	13✔
545	}	13✔
546
547	void divergence_matrix(const SurfaceMesh& mesh, SparseMatrix& D)	6✔
548	{
549	SparseMatrix G;	6✔
550	gradient_matrix(mesh, G);	6✔
551	DiagonalMatrix M;	6✔
552	divmass_matrix(mesh, M);	6✔
553	D = -G.transpose() * M;	6✔
554	}	6✔
555
556	void coordinates_to_matrix(const SurfaceMesh& mesh, DenseMatrix& X)	15✔
557	{
558	X.resize(mesh.n_vertices(), 3);	15✔
559	for (auto v : mesh.vertices())	51,397✔
560	X.row(v.idx()) = static_cast<Eigen::Vector3d>(mesh.position(v));	51,382✔
561	}	15✔
562
563	void matrix_to_coordinates(const DenseMatrix& X, SurfaceMesh& mesh)	7✔
564	{
565	assert((size_t)X.rows() == mesh.n_vertices() && X.cols() == 3);	7✔
566	for (auto v : mesh.vertices())	104✔
567	mesh.position(v) = X.row(v.idx());	97✔
568	}	7✔
569
570	} // namespace pmp

pmp-library / pmp-library / 14956887918

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous