25581337924

Committed 08 May 2026 09:51PM UTC coverage: 89.809% (-1.2%) from 90.99%

Build # 25581337924

Build Type

push

github

Committed by

dsieger

Commit Message

Simplify decimation test, disable unsupported test

Closes #248

Coverage Stats

5173 of 5760 relevant lines covered (89.81%)

603378.12 hits per line

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

83.7

/src/pmp/algorithms/subdivision.cpp

// Copyright 2011-2022 the Polygon Mesh Processing Library developers.
// SPDX-License-Identifier: MIT

#include "pmp/algorithms/subdivision.h"
#include "pmp/algorithms/differential_geometry.h"

#include <numbers>

namespace pmp {

void catmull_clark_subdivision(SurfaceMesh& mesh,
                               BoundaryHandling boundary_handling)
{
    auto points = mesh.vertex_property<Point>("v:point");
    auto vfeature = mesh.get_vertex_property<bool>("v:feature");
    auto efeature = mesh.get_edge_property<bool>("e:feature");

    // reserve memory
    const size_t nv = mesh.n_vertices();
    const size_t ne = mesh.n_edges();
    const size_t nf = mesh.n_faces();
    mesh.reserve(nv + ne + nf, 2 * ne + 4 * nf, 4 * nf);

    // get properties
    auto vpoint = mesh.add_vertex_property<Point>("catmull:vpoint");
    auto epoint = mesh.add_edge_property<Point>("catmull:epoint");
    auto fpoint = mesh.add_face_property<Point>("catmull:fpoint");

    // compute face vertices
    for (auto f : mesh.faces())
    {
        fpoint[f] = centroid(mesh, f);
    }

    // compute edge vertices
    for (auto e : mesh.edges())
    {
        // boundary or feature edge?
        if (mesh.is_boundary(e) || (efeature && efeature[e]))
        {
            epoint[e] =
                0.5f * (points[mesh.vertex(e, 0)] + points[mesh.vertex(e, 1)]);
        }

        // interior edge
        else
        {
            Point p(0, 0, 0);
            p += points[mesh.vertex(e, 0)];
            p += points[mesh.vertex(e, 1)];
            p += fpoint[mesh.face(e, 0)];
            p += fpoint[mesh.face(e, 1)];
            p *= 0.25f;
            epoint[e] = p;
        }
    }

    // compute new positions for old vertices
    for (auto v : mesh.vertices())
    {
        // isolated vertex?
        if (mesh.is_isolated(v))
        {
            vpoint[v] = points[v];
        }

        // boundary vertex?
        else if (mesh.is_boundary(v))
        {
            if (boundary_handling == BoundaryHandling::Preserve)
            {
                vpoint[v] = points[v];
            }
            else
            {
                auto h1 = mesh.halfedge(v);
                auto h0 = mesh.prev_halfedge(h1);

                Point p = points[v];
                p *= 6.0;
                p += points[mesh.to_vertex(h1)];
                p += points[mesh.from_vertex(h0)];
                p *= 0.125;

                vpoint[v] = p;
            }
        }

        // interior feature vertex?
        else if (vfeature && vfeature[v])
        {
            Point p = points[v];
            p *= 6.0;
            int count(0);

            for (auto h : mesh.halfedges(v))
            {
                if (efeature[mesh.edge(h)])
                {
                    p += points[mesh.to_vertex(h)];
                    ++count;
                }
            }

            if (count == 2) // vertex is on feature edge
            {
                p *= 0.125;
                vpoint[v] = p;
            }
            else // keep fixed
            {
                vpoint[v] = points[v];
            }
        }

        // interior vertex
        else
        {
            // weights from SIGGRAPH paper "Subdivision Surfaces in Character Animation"

            const Scalar k = mesh.valence(v);
            Point p(0, 0, 0);

            for (auto vv : mesh.vertices(v))
                p += points[vv];

            for (auto f : mesh.faces(v))
                p += fpoint[f];

            p /= (k * k);

            p += ((k - 2.0f) / k) * points[v];

            vpoint[v] = p;
        }
    }

    // assign new positions to old vertices
    for (auto v : mesh.vertices())
    {
        points[v] = vpoint[v];
    }

    // split edges
    for (auto e : mesh.edges())
    {
        // feature edge?
        if (efeature && efeature[e])
        {
            auto h = mesh.insert_vertex(e, epoint[e]);
            auto v = mesh.to_vertex(h);
            auto e0 = mesh.edge(h);
            auto e1 = mesh.edge(mesh.next_halfedge(h));

            vfeature[v] = true;
            efeature[e0] = true;
            efeature[e1] = true;
        }

        // normal edge
        else
        {
            mesh.insert_vertex(e, epoint[e]);
        }
    }

    // split faces
    for (auto f : mesh.faces())
    {
        auto h0 = mesh.halfedge(f);
        mesh.insert_edge(h0, mesh.next_halfedge(mesh.next_halfedge(h0)));

        auto h1 = mesh.next_halfedge(h0);
        mesh.insert_vertex(mesh.edge(h1), fpoint[f]);

        auto h = mesh.next_halfedge(mesh.next_halfedge(mesh.next_halfedge(h1)));
        while (h != h0)
        {
            mesh.insert_edge(h1, h);
            h = mesh.next_halfedge(mesh.next_halfedge(mesh.next_halfedge(h1)));
        }
    }

    // clean-up properties
    mesh.remove_vertex_property(vpoint);
    mesh.remove_edge_property(epoint);
    mesh.remove_face_property(fpoint);
}

void loop_subdivision(SurfaceMesh& mesh, BoundaryHandling boundary_handling)
{
    auto points = mesh.vertex_property<Point>("v:point");
    auto vfeature = mesh.get_vertex_property<bool>("v:feature");
    auto efeature = mesh.get_edge_property<bool>("e:feature");

    if (!mesh.is_triangle_mesh())
    {
        auto what = std::string{__func__} + ": Not a triangle mesh.";
        throw InvalidInputException(what);
    }

    // reserve memory
    const size_t nv = mesh.n_vertices();
    const size_t ne = mesh.n_edges();
    const size_t nf = mesh.n_faces();
    mesh.reserve(nv + ne, 2 * ne + 3 * nf, 4 * nf);

    // add properties
    auto vpoint = mesh.add_vertex_property<Point>("loop:vpoint");
    auto epoint = mesh.add_edge_property<Point>("loop:epoint");

    // compute vertex positions
    for (auto v : mesh.vertices())
    {
        // isolated vertex?
        if (mesh.is_isolated(v))
        {
            vpoint[v] = points[v];
        }

        // boundary vertex?
        else if (mesh.is_boundary(v))
        {
            if (boundary_handling == BoundaryHandling::Preserve)
            {
                vpoint[v] = points[v];
            }
            else
            {
                auto h1 = mesh.halfedge(v);
                auto h0 = mesh.prev_halfedge(h1);

                Point p = points[v];
                p *= 6.0;
                p += points[mesh.to_vertex(h1)];
                p += points[mesh.from_vertex(h0)];
                p *= 0.125;
                vpoint[v] = p;
            }
        }

        // interior feature vertex?
        else if (vfeature && vfeature[v])
        {
            Point p = points[v];
            p *= 6.0;
            int count(0);

            for (auto h : mesh.halfedges(v))
            {
                if (efeature[mesh.edge(h)])
                {
                    p += points[mesh.to_vertex(h)];
                    ++count;
                }
            }

            if (count == 2) // vertex is on feature edge
            {
                p *= 0.125;
                vpoint[v] = p;
            }
            else // keep fixed
            {
                vpoint[v] = points[v];
            }
        }

        // interior vertex
        else
        {
            Point p(0, 0, 0);
            Scalar k(0);

            for (auto vv : mesh.vertices(v))
            {
                p += points[vv];
                ++k;
            }
            p /= k;

            const Scalar beta =
                (0.625 -
                 pow(0.375 + 0.25 * std::cos(2.0 * std::numbers::pi / k), 2.0));

            vpoint[v] = points[v] * (Scalar)(1.0 - beta) + beta * p;
        }
    }

    // compute edge positions
    for (auto e : mesh.edges())
    {
        // boundary or feature edge?
        if (mesh.is_boundary(e) || (efeature && efeature[e]))
        {
            epoint[e] =
                (points[mesh.vertex(e, 0)] + points[mesh.vertex(e, 1)]) *
                Scalar(0.5);
        }

        // interior edge
        else
        {
            auto h0 = mesh.halfedge(e, 0);
            auto h1 = mesh.halfedge(e, 1);
            Point p = points[mesh.to_vertex(h0)];
            p += points[mesh.to_vertex(h1)];
            p *= 3.0;
            p += points[mesh.to_vertex(mesh.next_halfedge(h0))];
            p += points[mesh.to_vertex(mesh.next_halfedge(h1))];
            p *= 0.125;
            epoint[e] = p;
        }
    }

    // set new vertex positions
    for (auto v : mesh.vertices())
    {
        points[v] = vpoint[v];
    }

    // insert new vertices on edges
    for (auto e : mesh.edges())
    {
        // feature edge?
        if (efeature && efeature[e])
        {
            auto h = mesh.insert_vertex(e, epoint[e]);
            auto v = mesh.to_vertex(h);
            auto e0 = mesh.edge(h);
            auto e1 = mesh.edge(mesh.next_halfedge(h));

            vfeature[v] = true;
            efeature[e0] = true;
            efeature[e1] = true;
        }

        // normal edge
        else
        {
            mesh.insert_vertex(e, epoint[e]);
        }
    }

    // split faces
    Halfedge h;
    for (auto f : mesh.faces())
    {
        h = mesh.halfedge(f);
        mesh.insert_edge(h, mesh.next_halfedge(mesh.next_halfedge(h)));
        h = mesh.next_halfedge(h);
        mesh.insert_edge(h, mesh.next_halfedge(mesh.next_halfedge(h)));
        h = mesh.next_halfedge(h);
        mesh.insert_edge(h, mesh.next_halfedge(mesh.next_halfedge(h)));
    }

    // clean-up properties
    mesh.remove_vertex_property(vpoint);
    mesh.remove_edge_property(epoint);
}

void quad_tri_subdivision(SurfaceMesh& mesh, BoundaryHandling boundary_handling)
{
    auto points = mesh.vertex_property<Point>("v:point");

    // split each edge evenly into two parts
    for (auto e : mesh.edges())
    {
        mesh.insert_vertex(
            e, 0.5f * (points[mesh.vertex(e, 0)] + points[mesh.vertex(e, 1)]));
    }

    // subdivide faces without repositioning
    for (auto f : mesh.faces())
    {
        const size_t f_val = mesh.valence(f) / 2;
        if (f_val == 3)
        {
            // face was a triangle
            Halfedge h0 = mesh.halfedge(f);
            Halfedge h1 = mesh.next_halfedge(mesh.next_halfedge(h0));
            mesh.insert_edge(h0, h1);

            h0 = mesh.next_halfedge(h0);
            h1 = mesh.next_halfedge(mesh.next_halfedge(h0));
            mesh.insert_edge(h0, h1);

            h0 = mesh.next_halfedge(h0);
            h1 = mesh.next_halfedge(mesh.next_halfedge(h0));
            mesh.insert_edge(h0, h1);
        }
        else
        {
            // quadrangulate the rest
            const Halfedge h0 = mesh.halfedge(f);
            Halfedge h1 = mesh.next_halfedge(mesh.next_halfedge(h0));
            //NOTE: It's important to calculate the centroid before inserting the new edge
            auto cen = centroid(mesh, f);
            h1 = mesh.insert_edge(h0, h1);
            mesh.insert_vertex(mesh.edge(h1), cen);

            auto h =
                mesh.next_halfedge(mesh.next_halfedge(mesh.next_halfedge(h1)));
            while (h != h0)
            {
                mesh.insert_edge(h1, h);
                h = mesh.next_halfedge(
                    mesh.next_halfedge(mesh.next_halfedge(h1)));
            }
        }
    }

    auto new_pos =
        mesh.add_vertex_property<Point>("quad_tri:new_position", Point(0));

    for (auto v : mesh.vertices())
    {
        if (mesh.is_boundary(v))
        {
            if (boundary_handling == BoundaryHandling::Preserve)
            {
                new_pos[v] = points[v];
            }
            else
            {
                new_pos[v] = 0.5 * points[v];

                // add neighboring vertices on boundary
                for (auto vv : mesh.vertices(v))
                {
                    if (mesh.is_boundary(vv))
                    {
                        new_pos[v] += 0.25 * points[vv];
                    }
                }
            }
        }
        else
        {
            // count the number of faces and quads surrounding the vertex
            int n_faces = 0;
            int n_quads = 0;
            for (auto f : mesh.faces(v))
            {
                n_faces++;
                if (mesh.valence(f) == 4)
                    n_quads++;
            }

            if (n_quads == 0)
            {
                // vertex is surrounded only by triangles
                const double a =
                    2.0 *
                    pow(3.0 / 8.0 +
                            (std::cos(2.0 * std::numbers::pi / n_faces) - 1.0) /
                                4.0,
                        2.0);
                const double b = (1.0 - a) / n_faces;

                new_pos[v] = a * points[v];
                for (auto vv : mesh.vertices(v))
                {
                    new_pos[v] += b * points[vv];
                }
            }
            else if (n_quads == n_faces)
            {
                // vertex is surrounded only by quads
                const double c = (n_faces - 3.0) / n_faces;
                const double d = 2.0 / pow(n_faces, 2.0);
                const double e = 1.0 / pow(n_faces, 2.0);

                new_pos[v] = c * points[v];
                for (auto h : mesh.halfedges(v))
                {
                    new_pos[v] += d * points[mesh.to_vertex(h)];
                    new_pos[v] +=
                        e * points[mesh.to_vertex(mesh.next_halfedge(h))];
                }
            }
            else
            {
                // vertex is surrounded by triangles and quads
                const double alpha =
                    1.0 / (1.0 + 0.5 * n_faces + 0.25 * n_quads);
                const double beta = 0.5 * alpha;
                const double gamma = 0.25 * alpha;

                new_pos[v] = alpha * points[v];
                for (auto h : mesh.halfedges(v))
                {
                    new_pos[v] += beta * points[mesh.to_vertex(h)];
                    if (mesh.valence(mesh.face(h)) == 4)
                    {
                        new_pos[v] +=
                            gamma *
                            points[mesh.to_vertex(mesh.next_halfedge(h))];
                    }
                }
            }
        }
    }

    // apply new positions to the mesh
    for (auto v : mesh.vertices())
    {
        points[v] = new_pos[v];
    }

    mesh.remove_vertex_property(new_pos);
}

void linear_subdivision(SurfaceMesh& mesh)
{
    auto points = mesh.vertex_property<Point>("v:point");

    // linear subdivision of edges
    for (auto e : mesh.edges())
    {
        mesh.insert_vertex(
            e, 0.5f * (points[mesh.vertex(e, 0)] + points[mesh.vertex(e, 1)]));
    }

    // subdivide faces
    for (auto f : mesh.faces())
    {
        const size_t f_val = mesh.valence(f) / 2;

        if (f_val == 3) // triangle
        {
            Halfedge h0 = mesh.halfedge(f);
            Halfedge h1 = mesh.next_halfedge(mesh.next_halfedge(h0));
            mesh.insert_edge(h0, h1);

            h0 = mesh.next_halfedge(h0);
            h1 = mesh.next_halfedge(mesh.next_halfedge(h0));
            mesh.insert_edge(h0, h1);

            h0 = mesh.next_halfedge(h0);
            h1 = mesh.next_halfedge(mesh.next_halfedge(h0));
            mesh.insert_edge(h0, h1);
        }
        else // quadrangulate other faces
        {
            const Halfedge h0 = mesh.halfedge(f);
            Halfedge h1 = mesh.next_halfedge(mesh.next_halfedge(h0));

            // NOTE: It's important to calculate the centroid before inserting the new edge
            auto cen = centroid(mesh, f);
            h1 = mesh.insert_edge(h0, h1);
            mesh.insert_vertex(mesh.edge(h1), cen);

            auto h =
                mesh.next_halfedge(mesh.next_halfedge(mesh.next_halfedge(h1)));
            while (h != h0)
            {
                mesh.insert_edge(h1, h);
                h = mesh.next_halfedge(
                    mesh.next_halfedge(mesh.next_halfedge(h1)));
            }
        }
    }
}

} // namespace pmp

1	// Copyright 2011-2022 the Polygon Mesh Processing Library developers.
2	// SPDX-License-Identifier: MIT
3
4	#include "pmp/algorithms/subdivision.h"
5	#include "pmp/algorithms/differential_geometry.h"
6
7	#include <numbers>
8
9	namespace pmp {
10
11	void catmull_clark_subdivision(SurfaceMesh& mesh,	15✔
12	BoundaryHandling boundary_handling)
13	{
14	auto points = mesh.vertex_property<Point>("v:point");	75✔
15	auto vfeature = mesh.get_vertex_property<bool>("v:feature");	30✔
16	auto efeature = mesh.get_edge_property<bool>("e:feature");	15✔
17
18	// reserve memory
19	const size_t nv = mesh.n_vertices();	15✔
20	const size_t ne = mesh.n_edges();	15✔
21	const size_t nf = mesh.n_faces();	15✔
22	mesh.reserve(nv + ne + nf, 2 * ne + 4 * nf, 4 * nf);	15✔
23
24	// get properties
25	auto vpoint = mesh.add_vertex_property<Point>("catmull:vpoint");	60✔
26	auto epoint = mesh.add_edge_property<Point>("catmull:epoint");	60✔
27	auto fpoint = mesh.add_face_property<Point>("catmull:fpoint");	60✔
28
29	// compute face vertices
30	for (auto f : mesh.faces())	532✔
31	{
32	fpoint[f] = centroid(mesh, f);	517✔
33	}
34
35	// compute edge vertices
36	for (auto e : mesh.edges())	1,051✔
37	{
38	// boundary or feature edge?
39	if (mesh.is_boundary(e) \|\| (efeature && efeature[e]))	1,036✔
40	{
41	epoint[e] =	16✔
42	0.5f * (points[mesh.vertex(e, 0)] + points[mesh.vertex(e, 1)]);	32✔
43	}
44
45	// interior edge
46	else
47	{
48	Point p(0, 0, 0);	1,020✔
49	p += points[mesh.vertex(e, 0)];	1,020✔
50	p += points[mesh.vertex(e, 1)];	1,020✔
51	p += fpoint[mesh.face(e, 0)];	1,020✔
52	p += fpoint[mesh.face(e, 1)];	1,020✔
53	p *= 0.25f;	1,020✔
54	epoint[e] = p;	1,020✔
55	}
56	}
57
58	// compute new positions for old vertices
59	for (auto v : mesh.vertices())	563✔
60	{
61	// isolated vertex?
62	if (mesh.is_isolated(v))	548✔
63	{
64	vpoint[v] = points[v];	×
65	}
66
67	// boundary vertex?
68	else if (mesh.is_boundary(v))	548✔
69	{
70	if (boundary_handling == BoundaryHandling::Preserve)	4✔
71	{
72	vpoint[v] = points[v];	×
73	}
74	else
75	{
76	auto h1 = mesh.halfedge(v);	4✔
77	auto h0 = mesh.prev_halfedge(h1);	4✔
78
79	Point p = points[v];	4✔
80	p *= 6.0;	4✔
81	p += points[mesh.to_vertex(h1)];	4✔
82	p += points[mesh.from_vertex(h0)];	4✔
83	p *= 0.125;	4✔
84
85	vpoint[v] = p;	4✔
86	}
87	}
88
89	// interior feature vertex?
90	else if (vfeature && vfeature[v])	544✔
91	{
92	Point p = points[v];	8✔
93	p *= 6.0;	8✔
94	int count(0);	8✔
95
96	for (auto h : mesh.halfedges(v))	32✔
97	{
98	if (efeature[mesh.edge(h)])	24✔
99	{
100	p += points[mesh.to_vertex(h)];	24✔
101	++count;	24✔
102	}
103	}
104
105	if (count == 2) // vertex is on feature edge	8✔
106	{
107	p *= 0.125;	×
108	vpoint[v] = p;	×
109	}
110	else // keep fixed
111	{
112	vpoint[v] = points[v];	8✔
113	}
114	}
115
116	// interior vertex
117	else
118	{
119	// weights from SIGGRAPH paper "Subdivision Surfaces in Character Animation"
120
121	const Scalar k = mesh.valence(v);	536✔
122	Point p(0, 0, 0);	536✔
123
124	for (auto vv : mesh.vertices(v))	2,576✔
125	p += points[vv];	2,040✔
126
127	for (auto f : mesh.faces(v))	2,576✔
128	p += fpoint[f];	2,040✔
129
130	p /= (k * k);	536✔
131
132	p += ((k - 2.0f) / k) * points[v];	536✔
133
134	vpoint[v] = p;	536✔
135	}
136	}
137
138	// assign new positions to old vertices
139	for (auto v : mesh.vertices())	563✔
140	{
141	points[v] = vpoint[v];	548✔
142	}
143
144	// split edges
145	for (auto e : mesh.edges())	1,051✔
146	{
147	// feature edge?
148	if (efeature && efeature[e])	1,036✔
149	{
150	auto h = mesh.insert_vertex(e, epoint[e]);	12✔
151	auto v = mesh.to_vertex(h);	12✔
152	auto e0 = mesh.edge(h);	12✔
153	auto e1 = mesh.edge(mesh.next_halfedge(h));	12✔
154
155	vfeature[v] = true;	12✔
156	efeature[e0] = true;	12✔
157	efeature[e1] = true;	12✔
158	}
159
160	// normal edge
161	else
162	{
163	mesh.insert_vertex(e, epoint[e]);	1,024✔
164	}
165	}
166
167	// split faces
168	for (auto f : mesh.faces())	532✔
169	{
170	auto h0 = mesh.halfedge(f);	517✔
171	mesh.insert_edge(h0, mesh.next_halfedge(mesh.next_halfedge(h0)));	517✔
172
173	auto h1 = mesh.next_halfedge(h0);	517✔
174	mesh.insert_vertex(mesh.edge(h1), fpoint[f]);	517✔
175
176	auto h = mesh.next_halfedge(mesh.next_halfedge(mesh.next_halfedge(h1)));	517✔
177	while (h != h0)	1,551✔
178	{
179	mesh.insert_edge(h1, h);	1,034✔
180	h = mesh.next_halfedge(mesh.next_halfedge(mesh.next_halfedge(h1)));	1,034✔
181	}
182	}
183
184	// clean-up properties
185	mesh.remove_vertex_property(vpoint);	15✔
186	mesh.remove_edge_property(epoint);	15✔
187	mesh.remove_face_property(fpoint);	15✔
188	}	15✔
189
190	void loop_subdivision(SurfaceMesh& mesh, BoundaryHandling boundary_handling)	675✔
191	{
192	auto points = mesh.vertex_property<Point>("v:point");	3,375✔
193	auto vfeature = mesh.get_vertex_property<bool>("v:feature");	1,350✔
194	auto efeature = mesh.get_edge_property<bool>("e:feature");	675✔
195
196	if (!mesh.is_triangle_mesh())	675✔
197	{
198	auto what = std::string{__func__} + ": Not a triangle mesh.";	×
199	throw InvalidInputException(what);	×
200	}	×
201
202	// reserve memory
203	const size_t nv = mesh.n_vertices();	675✔
204	const size_t ne = mesh.n_edges();	675✔
205	const size_t nf = mesh.n_faces();	675✔
206	mesh.reserve(nv + ne, 2 * ne + 3 * nf, 4 * nf);	675✔
207
208	// add properties
209	auto vpoint = mesh.add_vertex_property<Point>("loop:vpoint");	2,700✔
210	auto epoint = mesh.add_edge_property<Point>("loop:epoint");	2,700✔
211
212	// compute vertex positions
213	for (auto v : mesh.vertices())	446,640✔
214	{
215	// isolated vertex?
216	if (mesh.is_isolated(v))	445,965✔
217	{
218	vpoint[v] = points[v];	×
219	}
220
221	// boundary vertex?
222	else if (mesh.is_boundary(v))	445,965✔
223	{
224	if (boundary_handling == BoundaryHandling::Preserve)	22✔
225	{
226	vpoint[v] = points[v];	×
227	}
228	else
229	{
230	auto h1 = mesh.halfedge(v);	22✔
231	auto h0 = mesh.prev_halfedge(h1);	22✔
232
233	Point p = points[v];	22✔
234	p *= 6.0;	22✔
235	p += points[mesh.to_vertex(h1)];	22✔
236	p += points[mesh.from_vertex(h0)];	22✔
237	p *= 0.125;	22✔
238	vpoint[v] = p;	22✔
239	}
240	}
241
242	// interior feature vertex?
243	else if (vfeature && vfeature[v])	445,943✔
244	{
245	Point p = points[v];	20✔
246	p *= 6.0;	20✔
247	int count(0);	20✔
248
249	for (auto h : mesh.halfedges(v))	116✔
250	{
251	if (efeature[mesh.edge(h)])	96✔
252	{
253	p += points[mesh.to_vertex(h)];	84✔
254	++count;	84✔
255	}
256	}
257
258	if (count == 2) // vertex is on feature edge	20✔
259	{
260	p *= 0.125;	×
261	vpoint[v] = p;	×
262	}
263	else // keep fixed
264	{
265	vpoint[v] = points[v];	20✔
266	}
267	}
268
269	// interior vertex
270	else
271	{
272	Point p(0, 0, 0);	445,923✔
273	Scalar k(0);	445,923✔
274
275	for (auto vv : mesh.vertices(v))	3,113,425✔
276	{
277	p += points[vv];	2,667,502✔
278	++k;	2,667,502✔
279	}
280	p /= k;	445,923✔
281
282	const Scalar beta =
283	(0.625 -	445,923✔
284	pow(0.375 + 0.25 * std::cos(2.0 * std::numbers::pi / k), 2.0));	445,923✔
285
286	vpoint[v] = points[v] * (Scalar)(1.0 - beta) + beta * p;	445,923✔
287	}
288	}
289
290	// compute edge positions
291	for (auto e : mesh.edges())	1,334,507✔
292	{
293	// boundary or feature edge?
294	if (mesh.is_boundary(e) \|\| (efeature && efeature[e]))	1,333,832✔
295	{
296	epoint[e] =	64✔
297	(points[mesh.vertex(e, 0)] + points[mesh.vertex(e, 1)]) *	128✔
298	Scalar(0.5);	64✔
299	}
300
301	// interior edge
302	else
303	{
304	auto h0 = mesh.halfedge(e, 0);	1,333,768✔
305	auto h1 = mesh.halfedge(e, 1);	1,333,768✔
306	Point p = points[mesh.to_vertex(h0)];	1,333,768✔
307	p += points[mesh.to_vertex(h1)];	1,333,768✔
308	p *= 3.0;	1,333,768✔
309	p += points[mesh.to_vertex(mesh.next_halfedge(h0))];	1,333,768✔
310	p += points[mesh.to_vertex(mesh.next_halfedge(h1))];	1,333,768✔
311	p *= 0.125;	1,333,768✔
312	epoint[e] = p;	1,333,768✔
313	}
314	}
315
316	// set new vertex positions
317	for (auto v : mesh.vertices())	446,640✔
318	{
319	points[v] = vpoint[v];	445,965✔
320	}
321
322	// insert new vertices on edges
323	for (auto e : mesh.edges())	1,334,507✔
324	{
325	// feature edge?
326	if (efeature && efeature[e])	1,333,832✔
327	{
328	auto h = mesh.insert_vertex(e, epoint[e]);	42✔
329	auto v = mesh.to_vertex(h);	42✔
330	auto e0 = mesh.edge(h);	42✔
331	auto e1 = mesh.edge(mesh.next_halfedge(h));	42✔
332
333	vfeature[v] = true;	42✔
334	efeature[e0] = true;	42✔
335	efeature[e1] = true;	42✔
336	}
337
338	// normal edge
339	else
340	{
341	mesh.insert_vertex(e, epoint[e]);	1,333,790✔
342	}
343	}
344
345	// split faces
346	Halfedge h;	675✔
347	for (auto f : mesh.faces())	889,889✔
348	{
349	h = mesh.halfedge(f);	889,214✔
350	mesh.insert_edge(h, mesh.next_halfedge(mesh.next_halfedge(h)));	889,214✔
351	h = mesh.next_halfedge(h);	889,214✔
352	mesh.insert_edge(h, mesh.next_halfedge(mesh.next_halfedge(h)));	889,214✔
353	h = mesh.next_halfedge(h);	889,214✔
354	mesh.insert_edge(h, mesh.next_halfedge(mesh.next_halfedge(h)));	889,214✔
355	}
356
357	// clean-up properties
358	mesh.remove_vertex_property(vpoint);	675✔
359	mesh.remove_edge_property(epoint);	675✔
360	}	675✔
361
362	void quad_tri_subdivision(SurfaceMesh& mesh, BoundaryHandling boundary_handling)	3✔
363	{
364	auto points = mesh.vertex_property<Point>("v:point");	12✔
365
366	// split each edge evenly into two parts
367	for (auto e : mesh.edges())	29✔
368	{
369	mesh.insert_vertex(	26✔
370	e, 0.5f * (points[mesh.vertex(e, 0)] + points[mesh.vertex(e, 1)]));	52✔
371	}
372
373	// subdivide faces without repositioning
374	for (auto f : mesh.faces())	18✔
375	{
376	const size_t f_val = mesh.valence(f) / 2;	15✔
377	if (f_val == 3)	15✔
378	{
379	// face was a triangle
380	Halfedge h0 = mesh.halfedge(f);	8✔
381	Halfedge h1 = mesh.next_halfedge(mesh.next_halfedge(h0));	8✔
382	mesh.insert_edge(h0, h1);	8✔
383
384	h0 = mesh.next_halfedge(h0);	8✔
385	h1 = mesh.next_halfedge(mesh.next_halfedge(h0));	8✔
386	mesh.insert_edge(h0, h1);	8✔
387
388	h0 = mesh.next_halfedge(h0);	8✔
389	h1 = mesh.next_halfedge(mesh.next_halfedge(h0));	8✔
390	mesh.insert_edge(h0, h1);	8✔
391	}
392	else
393	{
394	// quadrangulate the rest
395	const Halfedge h0 = mesh.halfedge(f);	7✔
396	Halfedge h1 = mesh.next_halfedge(mesh.next_halfedge(h0));	7✔
397	//NOTE: It's important to calculate the centroid before inserting the new edge
398	auto cen = centroid(mesh, f);	7✔
399	h1 = mesh.insert_edge(h0, h1);	7✔
400	mesh.insert_vertex(mesh.edge(h1), cen);	7✔
401
402	auto h =
403	mesh.next_halfedge(mesh.next_halfedge(mesh.next_halfedge(h1)));	7✔
404	while (h != h0)	21✔
405	{
406	mesh.insert_edge(h1, h);	14✔
407	h = mesh.next_halfedge(	14✔
408	mesh.next_halfedge(mesh.next_halfedge(h1)));
409	}
410	}
411	}
412
413	auto new_pos =
414	mesh.add_vertex_property<Point>("quad_tri:new_position", Point(0));	6✔
415
416	for (auto v : mesh.vertices())	53✔
417	{
418	if (mesh.is_boundary(v))	50✔
419	{
420	if (boundary_handling == BoundaryHandling::Preserve)	×
421	{
422	new_pos[v] = points[v];	×
423	}
424	else
425	{
426	new_pos[v] = 0.5 * points[v];	×
427
428	// add neighboring vertices on boundary
429	for (auto vv : mesh.vertices(v))	×
430	{
431	if (mesh.is_boundary(vv))	×
432	{
433	new_pos[v] += 0.25 * points[vv];	×
434	}
435	}
436	}
437	}
438	else
439	{
440	// count the number of faces and quads surrounding the vertex
441	int n_faces = 0;	50✔
442	int n_quads = 0;	50✔
443	for (auto f : mesh.faces(v))	258✔
444	{
445	n_faces++;	208✔
446	if (mesh.valence(f) == 4)	208✔
447	n_quads++;	112✔
448	}
449
450	if (n_quads == 0)	50✔
451	{
452	// vertex is surrounded only by triangles
453	const double a =
454	2.0 *
455	pow(3.0 / 8.0 +	15✔
456	(std::cos(2.0 * std::numbers::pi / n_faces) - 1.0) /	15✔
457	4.0,
458	2.0);	15✔
459	const double b = (1.0 - a) / n_faces;	15✔
460
461	new_pos[v] = a * points[v];	15✔
462	for (auto vv : mesh.vertices(v))	91✔
463	{
464	new_pos[v] += b * points[vv];	76✔
465	}
466	}
467	else if (n_quads == n_faces)	35✔
468	{
469	// vertex is surrounded only by quads
470	const double c = (n_faces - 3.0) / n_faces;	27✔
471	const double d = 2.0 / pow(n_faces, 2.0);	27✔
472	const double e = 1.0 / pow(n_faces, 2.0);	27✔
473
474	new_pos[v] = c * points[v];	27✔
475	for (auto h : mesh.halfedges(v))	127✔
476	{
477	new_pos[v] += d * points[mesh.to_vertex(h)];	100✔
478	new_pos[v] +=	100✔
479	e * points[mesh.to_vertex(mesh.next_halfedge(h))];	200✔
480	}
481	}
482	else
483	{
484	// vertex is surrounded by triangles and quads
485	const double alpha =	8✔
486	1.0 / (1.0 + 0.5 * n_faces + 0.25 * n_quads);	8✔
487	const double beta = 0.5 * alpha;	8✔
488	const double gamma = 0.25 * alpha;	8✔
489
490	new_pos[v] = alpha * points[v];	8✔
491	for (auto h : mesh.halfedges(v))	40✔
492	{
493	new_pos[v] += beta * points[mesh.to_vertex(h)];	32✔
494	if (mesh.valence(mesh.face(h)) == 4)	32✔
495	{
496	new_pos[v] +=	12✔
497	gamma *	12✔
498	points[mesh.to_vertex(mesh.next_halfedge(h))];	24✔
499	}
500	}
501	}
502	}
503	}
504
505	// apply new positions to the mesh
506	for (auto v : mesh.vertices())	53✔
507	{
508	points[v] = new_pos[v];	50✔
509	}
510
511	mesh.remove_vertex_property(new_pos);	3✔
512	}	3✔
513
514	void linear_subdivision(SurfaceMesh& mesh)	×
515	{
516	auto points = mesh.vertex_property<Point>("v:point");	×
517
518	// linear subdivision of edges
519	for (auto e : mesh.edges())	×
520	{
521	mesh.insert_vertex(	×
522	e, 0.5f * (points[mesh.vertex(e, 0)] + points[mesh.vertex(e, 1)]));	×
523	}
524
525	// subdivide faces
526	for (auto f : mesh.faces())	×
527	{
528	const size_t f_val = mesh.valence(f) / 2;	×
529
530	if (f_val == 3) // triangle	×
531	{
532	Halfedge h0 = mesh.halfedge(f);	×
533	Halfedge h1 = mesh.next_halfedge(mesh.next_halfedge(h0));	×
534	mesh.insert_edge(h0, h1);	×
535
536	h0 = mesh.next_halfedge(h0);	×
537	h1 = mesh.next_halfedge(mesh.next_halfedge(h0));	×
538	mesh.insert_edge(h0, h1);	×
539
540	h0 = mesh.next_halfedge(h0);	×
541	h1 = mesh.next_halfedge(mesh.next_halfedge(h0));	×
542	mesh.insert_edge(h0, h1);	×
543	}
544	else // quadrangulate other faces
545	{
546	const Halfedge h0 = mesh.halfedge(f);	×
547	Halfedge h1 = mesh.next_halfedge(mesh.next_halfedge(h0));	×
548
549	// NOTE: It's important to calculate the centroid before inserting the new edge
550	auto cen = centroid(mesh, f);	×
551	h1 = mesh.insert_edge(h0, h1);	×
552	mesh.insert_vertex(mesh.edge(h1), cen);	×
553
554	auto h =
555	mesh.next_halfedge(mesh.next_halfedge(mesh.next_halfedge(h1)));	×
556	while (h != h0)	×
557	{
558	mesh.insert_edge(h1, h);	×
559	h = mesh.next_halfedge(	×
560	mesh.next_halfedge(mesh.next_halfedge(h1)));
561	}
562	}
563	}
564	}	×
565
566	} // namespace pmp

pmp-library / pmp-library / 25581337924

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