14956893126

Committed 11 May 2025 02:44PM UTC coverage: 91.033% (-0.1%) from 91.129%

Build # 14956893126

Build Type

push

github

Committed by

web-flow

Commit Message

Write normal index if vertex normals are written (#234)

Signed-off-by: Mohamed El Shorbagy <mohrizq895@gmail.com>

Run Details

5 of 9 new or added lines in 1 file covered. (55.56%)

9 existing lines in 4 files now uncovered.

5208 of 5721 relevant lines covered (91.03%)

645685.96 hits per line

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

92.65

/src/pmp/algorithms/curvature.cpp

// Copyright 2011-2020 the Polygon Mesh Processing Library developers.
// Distributed under a MIT-style license, see LICENSE.txt for details.

#include "pmp/algorithms/curvature.h"
#include "pmp/exceptions.h"
#include "pmp/algorithms/normals.h"
#include "pmp/algorithms/differential_geometry.h"
#include "pmp/algorithms/laplace.h"

#include <algorithm>
#include <numbers>

namespace pmp {

class CurvatureAnalyzer
{
public:
    //! construct with mesh to be analyzed
    CurvatureAnalyzer(SurfaceMesh& mesh);

    ~CurvatureAnalyzer()
    {
        mesh_.remove_vertex_property(min_curvature_);
        mesh_.remove_vertex_property(max_curvature_);
    }

    //! compute curvature information for each vertex, optionally followed
    //! by some smoothing iterations of the curvature values
    void analyze(unsigned int post_smoothing_steps = 0);

    //! compute curvature information for each vertex, optionally followed
    //! by some smoothing iterations of the curvature values
    void analyze_tensor(unsigned int post_smoothing_steps = 0,
                        bool two_ring_neighborhood = false);

    //! return mean curvature
    Scalar mean_curvature(Vertex v) const
    {
        return Scalar(0.5) * (min_curvature_[v] + max_curvature_[v]);
    }

    //! return Gaussian curvature
    Scalar gauss_curvature(Vertex v) const
    {
        return min_curvature_[v] * max_curvature_[v];
    }

    //! return minimum (signed) curvature
    Scalar min_curvature(Vertex v) const { return min_curvature_[v]; }

    //! return maximum (signed) curvature
    Scalar max_curvature(Vertex v) const { return max_curvature_[v]; }

    //! return maximum absolute curvature
    Scalar max_abs_curvature(Vertex v) const
    {
        return std::max(fabs(min_curvature_[v]), fabs(max_curvature_[v]));
    }

private:
    // determine curvature values on boundary from non-boundary neighbors
    void set_boundary_curvatures();

    // smooth curvature values
    void smooth_curvatures(unsigned int iterations);

    SurfaceMesh& mesh_;
    VertexProperty<Scalar> min_curvature_;
    VertexProperty<Scalar> max_curvature_;
};

CurvatureAnalyzer::CurvatureAnalyzer(SurfaceMesh& mesh) : mesh_(mesh)
{
    min_curvature_ = mesh_.add_vertex_property<Scalar>("curv:min");
    max_curvature_ = mesh_.add_vertex_property<Scalar>("curv:max");
}

void CurvatureAnalyzer::analyze(unsigned int post_smoothing_steps)
{
    Scalar kmin, kmax, mean, gauss;
    Scalar area, sum_angles;
    Point p0, p1, p2;

    // compute area-normalized Laplace
    SparseMatrix L;
    laplace_matrix(mesh_, L);
    DiagonalMatrix M;
    mass_matrix(mesh_, M);
    DenseMatrix X;
    coordinates_to_matrix(mesh_, X);
    DenseMatrix LX = L * X;

    // mean curvature as norm of Laplace
    // Gauss curvatures as angle deficit
    // min/max from mean/gauss
    for (auto v : mesh_.vertices())
    {
        kmin = kmax = 0.0;

        if (!mesh_.is_isolated(v) && !mesh_.is_boundary(v))
        {
            p0 = mesh_.position(v);

            // Voronoi area
            area = M.diagonal()[v.idx()];

            // angle sum
            sum_angles = 0.0;
            for (auto vh : mesh_.halfedges(v))
            {
                p1 = mesh_.position(mesh_.to_vertex(vh));
                p2 = mesh_.position(
                    mesh_.to_vertex(mesh_.ccw_rotated_halfedge(vh)));
                sum_angles += angle(p1 - p0, p2 - p0);
            }

            mean = 0.5 * LX.row(v.idx()).norm() / area;
            gauss = (2.0 * std::numbers::pi - sum_angles) / area;

            const Scalar s = sqrt(std::max(Scalar(0.0), mean * mean - gauss));
            kmin = mean - s;
            kmax = mean + s;
        }

        min_curvature_[v] = kmin;
        max_curvature_[v] = kmax;
    }

    // boundary vertices: interpolate from interior neighbors
    set_boundary_curvatures();

    // smooth curvatures values
    smooth_curvatures(post_smoothing_steps);
}

void CurvatureAnalyzer::analyze_tensor(unsigned int post_smoothing_steps,
                                       bool two_ring_neighborhood)
{
    auto area = mesh_.add_vertex_property<double>("curv:area", 0.0);
    auto normal = mesh_.add_face_property<dvec3>("curv:normal");
    auto evec = mesh_.add_edge_property<dvec3>("curv:evec", dvec3(0, 0, 0));
    auto angle = mesh_.add_edge_property<double>("curv:angle", 0.0);

    dvec3 n0, n1, ev;
    double l, A, beta, a1, a2, a3;
    dmat3 tensor;

    double eval1, eval2, eval3, kmin, kmax;
    dvec3 evec1, evec2, evec3;

    std::vector<Vertex> neighborhood;
    neighborhood.reserve(15);

    // precompute Voronoi area per vertex
    DiagonalMatrix M;
    mass_matrix(mesh_, M);
    for (auto v : mesh_.vertices())
    {
        area[v] = M.diagonal()[v.idx()];
    }

    // precompute face normals
    for (auto f : mesh_.faces())
    {
        normal[f] = (dvec3)face_normal(mesh_, f);
    }

    // precompute dihedralAngle*edge_length*edge per edge
    for (auto e : mesh_.edges())
    {
        auto h0 = mesh_.halfedge(e, 0);
        auto h1 = mesh_.halfedge(e, 1);
        auto f0 = mesh_.face(h0);
        auto f1 = mesh_.face(h1);
        if (f0.is_valid() && f1.is_valid())
        {
            n0 = normal[f0];
            n1 = normal[f1];
            ev = (dvec3)mesh_.position(mesh_.to_vertex(h0));
            ev -= (dvec3)mesh_.position(mesh_.to_vertex(h1));
            l = norm(ev);
            ev /= l;
            l *= 0.5; // only consider half of the edge (matching Voronoi area)
            angle[e] = atan2(dot(cross(n0, n1), ev), dot(n0, n1));
            evec[e] = sqrt(l) * ev;
        }
    }

    // compute curvature tensor for each vertex
    for (auto v : mesh_.vertices())
    {
        kmin = 0.0;
        kmax = 0.0;

        if (!mesh_.is_isolated(v) && !mesh_.is_boundary(v))
        {
            // one-ring or two-ring neighborhood?
            neighborhood.clear();
            neighborhood.push_back(v);
            if (two_ring_neighborhood)
            {
                for (auto vv : mesh_.vertices(v))
                    neighborhood.push_back(vv);
            }

            A = 0.0;
            tensor = dmat3(0.0);

            // compute tensor over vertex neighborhood stored in vertices
            for (auto nit : neighborhood)
            {
                if (mesh_.is_boundary(nit))
                    continue;

                // accumulate tensor from dihedral angles around vertices
                for (auto e : mesh_.edges(nit))
                {
                    ev = evec[e];
                    beta = angle[e];
                    for (int i = 0; i < 3; ++i)
                        for (int j = 0; j < 3; ++j)
                            tensor(i, j) += beta * ev[i] * ev[j];
                }

                // accumulate area
                A += area[nit];
            }

            // normalize tensor by accumulated
            tensor /= A;

            // Eigen-decomposition
            const bool ok = symmetric_eigendecomposition(
                tensor, eval1, eval2, eval3, evec1, evec2, evec3);
            if (ok)
            {
                // curvature values:
                //   normal vector -> eval with smallest absolute value
                //   evals are sorted in decreasing order
                a1 = fabs(eval1);
                a2 = fabs(eval2);
                a3 = fabs(eval3);
                if (a1 < a2)
                {
                    if (a1 < a3)
                    {
                        // e1 is normal
                        kmax = eval2;
                        kmin = eval3;
                    }
                    else
                    {
                        // e3 is normal
                        kmax = eval1;
                        kmin = eval2;
                    }
                }
                else
                {
                    if (a2 < a3)
                    {
                        // e2 is normal
                        kmax = eval1;
                        kmin = eval3;
                    }
                    else
                    {
                        // e3 is normal
                        kmax = eval1;
                        kmin = eval2;
                    }
                }
            }
        }

        assert(kmin <= kmax);

        min_curvature_[v] = kmin;
        max_curvature_[v] = kmax;
    }

    // clean-up properties
    mesh_.remove_vertex_property(area);
    mesh_.remove_edge_property(evec);
    mesh_.remove_edge_property(angle);
    mesh_.remove_face_property(normal);

    // boundary vertices: interpolate from interior neighbors
    set_boundary_curvatures();

    // smooth curvature values
    smooth_curvatures(post_smoothing_steps);
}

void CurvatureAnalyzer::set_boundary_curvatures()
{
    for (auto v : mesh_.vertices())
    {
        if (mesh_.is_boundary(v))
        {
            Scalar kmin(0.0), kmax(0.0), sum(0.0);
            for (auto vv : mesh_.vertices(v))
            {
                if (!mesh_.is_boundary(vv))
                {
                    sum += 1.0;
                    kmin += min_curvature_[vv];
                    kmax += max_curvature_[vv];
                }
            }

            if (sum)
            {
                kmin /= sum;
                kmax /= sum;
            }

            min_curvature_[v] = kmin;
            max_curvature_[v] = kmax;
        }
    }
}

void CurvatureAnalyzer::smooth_curvatures(unsigned int iterations)
{
    if (iterations == 0)
        return;

    // Laplace matrix (clamp negative cotan weights to zero)
    SparseMatrix L;
    laplace_matrix(mesh_, L, true);

    // normalize each row by sum of weights
    // scale by 0.5 to make it more robust
    // multiply by -1 to make it neg. definite again
    const DiagonalMatrix D = -0.5 * L.diagonal().asDiagonal().inverse();
    L = D * L;

    // copy vertex curvatures to matrix
    const int n = mesh_.n_vertices();
    DenseMatrix curv(n, 2);
    for (auto v : mesh_.vertices())
    {
        curv(v.idx(), 0) = min_curvature_[v];
        curv(v.idx(), 1) = max_curvature_[v];
    }

    // perform smoothing iterations
    for (unsigned int i = 0; i < iterations; ++i)
    {
        curv += L * curv;
    }

    // copy result to curvatures
    for (auto v : mesh_.vertices())
    {
        min_curvature_[v] = curv(v.idx(), 0);
        max_curvature_[v] = curv(v.idx(), 1);
    }
}

void curvature_to_texture_coordinates(SurfaceMesh& mesh)
{
    auto curvatures = mesh.get_vertex_property<Scalar>("v:curv");
    assert(curvatures);

    // sort curvature values
    std::vector<Scalar> values;
    values.reserve(mesh.n_vertices());
    for (auto v : mesh.vertices())
    {
        values.push_back(curvatures[v]);
    }
    std::ranges::sort(values);
    const unsigned int n = values.size() - 1;
    // std::cout << "curvatures in [" << values[0] << ", " << values[n] << "]\n";

    // clamp upper/lower 5%
    const unsigned int i = n / 20;
    const Scalar kmin = values[i];
    Scalar kmax = values[n - 1 - i];

    // generate 1D texture coordinates
    auto tex = mesh.vertex_property<TexCoord>("v:tex");
    if (kmin < 0.0) // signed
    {
        kmax = std::max(fabs(kmin), fabs(kmax));
        for (auto v : mesh.vertices())
        {
            tex[v] = TexCoord((0.5f * curvatures[v] / kmax) + 0.5f, 0.0);
        }
    }
    else // unsigned
    {
        for (auto v : mesh.vertices())
        {
            tex[v] = TexCoord((curvatures[v] - kmin) / (kmax - kmin), 0.0);
        }
    }
}

void curvature(SurfaceMesh& mesh, Curvature c, int smoothing_steps,
               bool use_tensor, bool use_two_ring)
{
    CurvatureAnalyzer analyzer(mesh);
    if (use_tensor)
        analyzer.analyze_tensor(smoothing_steps, use_two_ring);
    else
        analyzer.analyze(smoothing_steps);

    auto curvatures = mesh.vertex_property<Scalar>("v:curv");

    switch (c)
    {
        case Curvature::Min:
        {
            for (auto v : mesh.vertices())
                curvatures[v] = analyzer.min_curvature(v);
            break;
        }
        case Curvature::Max:
        {
            for (auto v : mesh.vertices())
                curvatures[v] = analyzer.max_curvature(v);
            break;
        }
        case Curvature::Mean:
        {
            for (auto v : mesh.vertices())
                curvatures[v] = fabs(analyzer.mean_curvature(v));
            break;
        }
        case Curvature::Gauss:
        {
            for (auto v : mesh.vertices())
                curvatures[v] = analyzer.gauss_curvature(v);
            break;
        }
        case Curvature::MaxAbs:
        {
            for (auto v : mesh.vertices())
                curvatures[v] = analyzer.max_abs_curvature(v);
            break;
        }
        default:
            throw InvalidInputException("Unknown Curvature type");
    }
}

} // namespace pmp

1	// Copyright 2011-2020 the Polygon Mesh Processing Library developers.
2	// Distributed under a MIT-style license, see LICENSE.txt for details.
3
4	#include "pmp/algorithms/curvature.h"
5	#include "pmp/exceptions.h"
6	#include "pmp/algorithms/normals.h"
7	#include "pmp/algorithms/differential_geometry.h"
8	#include "pmp/algorithms/laplace.h"
9
10	#include <algorithm>
11	#include <numbers>
12
13	namespace pmp {
14
15	class CurvatureAnalyzer
16	{
17	public:
18	//! construct with mesh to be analyzed
19	CurvatureAnalyzer(SurfaceMesh& mesh);
20
21	~CurvatureAnalyzer()	8✔
22	{
23	mesh_.remove_vertex_property(min_curvature_);	8✔
24	mesh_.remove_vertex_property(max_curvature_);	8✔
25	}	8✔
26
27	//! compute curvature information for each vertex, optionally followed
28	//! by some smoothing iterations of the curvature values
29	void analyze(unsigned int post_smoothing_steps = 0);
30
31	//! compute curvature information for each vertex, optionally followed
32	//! by some smoothing iterations of the curvature values
33	void analyze_tensor(unsigned int post_smoothing_steps = 0,
34	bool two_ring_neighborhood = false);
35
36	//! return mean curvature
37	Scalar mean_curvature(Vertex v) const	20,484✔
38	{
39	return Scalar(0.5) * (min_curvature_[v] + max_curvature_[v]);	20,484✔
40	}
41
42	//! return Gaussian curvature
43	Scalar gauss_curvature(Vertex v) const	10,242✔
44	{
45	return min_curvature_[v] * max_curvature_[v];	10,242✔
46	}
47
48	//! return minimum (signed) curvature
49	Scalar min_curvature(Vertex v) const { return min_curvature_[v]; }	10,242✔
50
51	//! return maximum (signed) curvature
52	Scalar max_curvature(Vertex v) const { return max_curvature_[v]; }	10,242✔
53
54	//! return maximum absolute curvature
55	Scalar max_abs_curvature(Vertex v) const	111✔
56	{
57	return std::max(fabs(min_curvature_[v]), fabs(max_curvature_[v]));	111✔
58	}
59
60	private:
61	// determine curvature values on boundary from non-boundary neighbors
62	void set_boundary_curvatures();
63
64	// smooth curvature values
65	void smooth_curvatures(unsigned int iterations);
66
67	SurfaceMesh& mesh_;
68	VertexProperty<Scalar> min_curvature_;
69	VertexProperty<Scalar> max_curvature_;
70	};
71
72	CurvatureAnalyzer::CurvatureAnalyzer(SurfaceMesh& mesh) : mesh_(mesh)	8✔
73	{
74	min_curvature_ = mesh_.add_vertex_property<Scalar>("curv:min");	16✔
75	max_curvature_ = mesh_.add_vertex_property<Scalar>("curv:max");	16✔
76	}	8✔
77
78	void CurvatureAnalyzer::analyze(unsigned int post_smoothing_steps)	5✔
79	{
80	Scalar kmin, kmax, mean, gauss;
81	Scalar area, sum_angles;
82	Point p0, p1, p2;
83
84	// compute area-normalized Laplace
85	SparseMatrix L;	5✔
86	laplace_matrix(mesh_, L);	5✔
87	DiagonalMatrix M;	5✔
88	mass_matrix(mesh_, M);	5✔
89	DenseMatrix X;	5✔
90	coordinates_to_matrix(mesh_, X);	5✔
91	DenseMatrix LX = L * X;	5✔
92
93	// mean curvature as norm of Laplace
94	// Gauss curvatures as angle deficit
95	// min/max from mean/gauss
96	for (auto v : mesh_.vertices())	51,215✔
97	{
98	kmin = kmax = 0.0;	51,210✔
99
100	if (!mesh_.is_isolated(v) && !mesh_.is_boundary(v))	51,210✔
101	{
102	p0 = mesh_.position(v);	51,210✔
103
104	// Voronoi area
105	area = M.diagonal()[v.idx()];	51,210✔
106
107	// angle sum
108	sum_angles = 0.0;	51,210✔
109	for (auto vh : mesh_.halfedges(v))	358,410✔
110	{
111	p1 = mesh_.position(mesh_.to_vertex(vh));	307,200✔
112	p2 = mesh_.position(	307,200✔
113	mesh_.to_vertex(mesh_.ccw_rotated_halfedge(vh)));	307,200✔
114	sum_angles += angle(p1 - p0, p2 - p0);	307,200✔
115	}
116
117	mean = 0.5 * LX.row(v.idx()).norm() / area;	51,210✔
118	gauss = (2.0 * std::numbers::pi - sum_angles) / area;	51,210✔
119
120	const Scalar s = sqrt(std::max(Scalar(0.0), mean * mean - gauss));	51,210✔
121	kmin = mean - s;	51,210✔
122	kmax = mean + s;	51,210✔
123	}
124
125	min_curvature_[v] = kmin;	51,210✔
126	max_curvature_[v] = kmax;	51,210✔
127	}
128
129	// boundary vertices: interpolate from interior neighbors
130	set_boundary_curvatures();	5✔
131
132	// smooth curvatures values
133	smooth_curvatures(post_smoothing_steps);	5✔
134	}	5✔
135
136	void CurvatureAnalyzer::analyze_tensor(unsigned int post_smoothing_steps,	3✔
137	bool two_ring_neighborhood)
138	{
139	auto area = mesh_.add_vertex_property<double>("curv:area", 0.0);	6✔
140	auto normal = mesh_.add_face_property<dvec3>("curv:normal");	12✔
141	auto evec = mesh_.add_edge_property<dvec3>("curv:evec", dvec3(0, 0, 0));	6✔
142	auto angle = mesh_.add_edge_property<double>("curv:angle", 0.0);	6✔
143
144	dvec3 n0, n1, ev;
145	double l, A, beta, a1, a2, a3;
146	dmat3 tensor;
147
148	double eval1, eval2, eval3, kmin, kmax;
149	dvec3 evec1, evec2, evec3;
150
151	std::vector<Vertex> neighborhood;	3✔
152	neighborhood.reserve(15);	3✔
153
154	// precompute Voronoi area per vertex
155	DiagonalMatrix M;	3✔
156	mass_matrix(mesh_, M);	3✔
157	for (auto v : mesh_.vertices())	114✔
158	{
159	area[v] = M.diagonal()[v.idx()];	111✔
160	}
161
162	// precompute face normals
163	for (auto f : mesh_.faces())	207✔
164	{
165	normal[f] = (dvec3)face_normal(mesh_, f);	204✔
166	}
167
168	// precompute dihedralAngleedge_lengthedge per edge
169	for (auto e : mesh_.edges())	313✔
170	{
171	auto h0 = mesh_.halfedge(e, 0);	310✔
172	auto h1 = mesh_.halfedge(e, 1);	310✔
173	auto f0 = mesh_.face(h0);	310✔
174	auto f1 = mesh_.face(h1);	310✔
175	if (f0.is_valid() && f1.is_valid())	310✔
176	{
177	n0 = normal[f0];	302✔
178	n1 = normal[f1];	302✔
179	ev = (dvec3)mesh_.position(mesh_.to_vertex(h0));	302✔
180	ev -= (dvec3)mesh_.position(mesh_.to_vertex(h1));	302✔
181	l = norm(ev);	302✔
182	ev /= l;	302✔
183	l *= 0.5; // only consider half of the edge (matching Voronoi area)	302✔
184	angle[e] = atan2(dot(cross(n0, n1), ev), dot(n0, n1));	302✔
185	evec[e] = sqrt(l) * ev;	302✔
186	}
187	}
188
189	// compute curvature tensor for each vertex
190	for (auto v : mesh_.vertices())	114✔
191	{
192	kmin = 0.0;	111✔
193	kmax = 0.0;	111✔
194
195	if (!mesh_.is_isolated(v) && !mesh_.is_boundary(v))	111✔
196	{
197	// one-ring or two-ring neighborhood?
198	neighborhood.clear();	103✔
199	neighborhood.push_back(v);	103✔
200	if (two_ring_neighborhood)	103✔
201	{
202	for (auto vv : mesh_.vertices(v))	×
203	neighborhood.push_back(vv);	×
204	}
205
206	A = 0.0;	103✔
207	tensor = dmat3(0.0);	103✔
208
209	// compute tensor over vertex neighborhood stored in vertices
210	for (auto nit : neighborhood)	206✔
211	{
212	if (mesh_.is_boundary(nit))	103✔
213	continue;	×
214
215	// accumulate tensor from dihedral angles around vertices
216	for (auto e : mesh_.edges(nit))	699✔
217	{
218	ev = evec[e];	596✔
219	beta = angle[e];	596✔
220	for (int i = 0; i < 3; ++i)	2,384✔
221	for (int j = 0; j < 3; ++j)	7,152✔
222	tensor(i, j) += beta * ev[i] * ev[j];	5,364✔
223	}
224
225	// accumulate area
226	A += area[nit];	103✔
227	}
228
229	// normalize tensor by accumulated
230	tensor /= A;	103✔
231
232	// Eigen-decomposition
233	const bool ok = symmetric_eigendecomposition(	103✔
234	tensor, eval1, eval2, eval3, evec1, evec2, evec3);
235	if (ok)	103✔
236	{
237	// curvature values:
238	// normal vector -> eval with smallest absolute value
239	// evals are sorted in decreasing order
240	a1 = fabs(eval1);	103✔
241	a2 = fabs(eval2);	103✔
242	a3 = fabs(eval3);	103✔
243	if (a1 < a2)	103✔
244	{
245	if (a1 < a3)	×
246	{
247	// e1 is normal
248	kmax = eval2;	×
249	kmin = eval3;	×
250	}
251	else
252	{
253	// e3 is normal
254	kmax = eval1;	×
255	kmin = eval2;	×
256	}
257	}
258	else
259	{
260	if (a2 < a3)	103✔
261	{
262	// e2 is normal
263	kmax = eval1;	×
264	kmin = eval3;	×
265	}
266	else
267	{
268	// e3 is normal
269	kmax = eval1;	103✔
270	kmin = eval2;	103✔
271	}
272	}
273	}
274	}
275
276	assert(kmin <= kmax);	111✔
277
278	min_curvature_[v] = kmin;	111✔
279	max_curvature_[v] = kmax;	111✔
280	}
281
282	// clean-up properties
283	mesh_.remove_vertex_property(area);	3✔
284	mesh_.remove_edge_property(evec);	3✔
285	mesh_.remove_edge_property(angle);	3✔
286	mesh_.remove_face_property(normal);	3✔
287
288	// boundary vertices: interpolate from interior neighbors
289	set_boundary_curvatures();	3✔
290
291	// smooth curvature values
292	smooth_curvatures(post_smoothing_steps);	3✔
293	}	3✔
294
295	void CurvatureAnalyzer::set_boundary_curvatures()	8✔
296	{
297	for (auto v : mesh_.vertices())	102,650✔
298	{
299	if (mesh_.is_boundary(v))	51,321✔
300	{
301	Scalar kmin(0.0), kmax(0.0), sum(0.0);	8✔
302	for (auto vv : mesh_.vertices(v))	56✔
303	{
304	if (!mesh_.is_boundary(vv))	24✔
305	{
306	sum += 1.0;	8✔
307	kmin += min_curvature_[vv];	8✔
308	kmax += max_curvature_[vv];	8✔
309	}
310	}
311
312	if (sum)	8✔
313	{
314	kmin /= sum;	8✔
315	kmax /= sum;	8✔
316	}
317
318	min_curvature_[v] = kmin;	8✔
319	max_curvature_[v] = kmax;	8✔
320	}
321	}
322	}	8✔
323
324	void CurvatureAnalyzer::smooth_curvatures(unsigned int iterations)	8✔
325	{
326	if (iterations == 0)	8✔
327	return;	3✔
328
329	// Laplace matrix (clamp negative cotan weights to zero)
330	SparseMatrix L;	5✔
331	laplace_matrix(mesh_, L, true);	5✔
332
333	// normalize each row by sum of weights
334	// scale by 0.5 to make it more robust
335	// multiply by -1 to make it neg. definite again
336	const DiagonalMatrix D = -0.5 * L.diagonal().asDiagonal().inverse();	5✔
337	L = D * L;	5✔
338
339	// copy vertex curvatures to matrix
340	const int n = mesh_.n_vertices();	5✔
341	DenseMatrix curv(n, 2);	5✔
342	for (auto v : mesh_.vertices())	51,215✔
343	{
344	curv(v.idx(), 0) = min_curvature_[v];	51,210✔
345	curv(v.idx(), 1) = max_curvature_[v];	51,210✔
346	}
347
348	// perform smoothing iterations
349	for (unsigned int i = 0; i < iterations; ++i)	10✔
350	{
351	curv += L * curv;	5✔
352	}
353
354	// copy result to curvatures
355	for (auto v : mesh_.vertices())	51,215✔
356	{
357	min_curvature_[v] = curv(v.idx(), 0);	51,210✔
358	max_curvature_[v] = curv(v.idx(), 1);	51,210✔
359	}
360	}	5✔
361
362	void curvature_to_texture_coordinates(SurfaceMesh& mesh)	1✔
363	{
364	auto curvatures = mesh.get_vertex_property<Scalar>("v:curv");	1✔
365	assert(curvatures);	1✔
366
367	// sort curvature values
368	std::vector<Scalar> values;	1✔
369	values.reserve(mesh.n_vertices());	1✔
370	for (auto v : mesh.vertices())	10,243✔
371	{
372	values.push_back(curvatures[v]);	10,242✔
373	}
374	std::ranges::sort(values);	1✔
375	const unsigned int n = values.size() - 1;	1✔
376	// std::cout << "curvatures in [" << values[0] << ", " << values[n] << "]\n";
377
378	// clamp upper/lower 5%
379	const unsigned int i = n / 20;	1✔
380	const Scalar kmin = values[i];	1✔
381	Scalar kmax = values[n - 1 - i];	1✔
382
383	// generate 1D texture coordinates
384	auto tex = mesh.vertex_property<TexCoord>("v:tex");	3✔
385	if (kmin < 0.0) // signed	1✔
386	{
387	kmax = std::max(fabs(kmin), fabs(kmax));	×
UNCOV 388	for (auto v : mesh.vertices())	×
389	{
UNCOV 390	tex[v] = TexCoord((0.5f * curvatures[v] / kmax) + 0.5f, 0.0);	×
391	}
392	}
393	else // unsigned
394	{
395	for (auto v : mesh.vertices())	10,243✔
396	{
397	tex[v] = TexCoord((curvatures[v] - kmin) / (kmax - kmin), 0.0);	10,242✔
398	}
399	}
400	}	1✔
401
402	void curvature(SurfaceMesh& mesh, Curvature c, int smoothing_steps,	8✔
403	bool use_tensor, bool use_two_ring)
404	{
405	CurvatureAnalyzer analyzer(mesh);	8✔
406	if (use_tensor)	8✔
407	analyzer.analyze_tensor(smoothing_steps, use_two_ring);	3✔
408	else
409	analyzer.analyze(smoothing_steps);	5✔
410
411	auto curvatures = mesh.vertex_property<Scalar>("v:curv");	8✔
412
413	switch (c)	8✔
414	{
415	case Curvature::Min:	1✔
416	{
417	for (auto v : mesh.vertices())	10,243✔
418	curvatures[v] = analyzer.min_curvature(v);	10,242✔
419	break;	1✔
420	}
421	case Curvature::Max:	1✔
422	{
423	for (auto v : mesh.vertices())	10,243✔
424	curvatures[v] = analyzer.max_curvature(v);	10,242✔
425	break;	1✔
426	}
427	case Curvature::Mean:	2✔
428	{
429	for (auto v : mesh.vertices())	20,486✔
430	curvatures[v] = fabs(analyzer.mean_curvature(v));	20,484✔
431	break;	2✔
432	}
433	case Curvature::Gauss:	1✔
434	{
435	for (auto v : mesh.vertices())	10,243✔
436	curvatures[v] = analyzer.gauss_curvature(v);	10,242✔
437	break;	1✔
438	}
439	case Curvature::MaxAbs:	3✔
440	{
441	for (auto v : mesh.vertices())	114✔
442	curvatures[v] = analyzer.max_abs_curvature(v);	111✔
443	break;	3✔
444	}
UNCOV 445	default:	×
UNCOV 446	throw InvalidInputException("Unknown Curvature type");	×
447	}
448	}	8✔
449
450	} // namespace pmp

pmp-library / pmp-library / 14956893126

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