15742553254

Committed 18 Jun 2025 08:07PM UTC coverage: 91.02% (+0.009%) from 91.011%

Build # 15742553254

Build Type

push

github

Committed by

mbotsch

Commit Message

make fairing more robust

Coverage Stats

28 of 31 new or added lines in 2 files covered. (90.32%)

5230 of 5746 relevant lines covered (91.02%)

642835.95 hits per line

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

91.22

/src/pmp/algorithms/hole_filling.cpp

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

#include "pmp/algorithms/hole_filling.h"

#include <Eigen/Dense>
#include <Eigen/Sparse>

#include <limits>
#include <vector>

#include "pmp/algorithms/fairing.h"
#include "pmp/algorithms/normals.h"

namespace pmp {
namespace {

class HoleFilling
{
public:
    explicit HoleFilling(SurfaceMesh& mesh);
    void fill_hole(Halfedge h);

private:
    struct Weight
    {
        Weight(Scalar angle = std::numeric_limits<Scalar>::max(),
               Scalar area = std::numeric_limits<Scalar>::max())
            : angle_(angle), area_(area)
        {
        }

        Weight operator+(const Weight& rhs) const
        {
            return Weight(std::max(angle_, rhs.angle_), area_ + rhs.area_);
        }

        bool operator<(const Weight& rhs) const
        {
            return (angle_ < rhs.angle_ ||
                    (angle_ == rhs.angle_ && area_ < rhs.area_));
        }

        Scalar angle_;
        Scalar area_;
    };

    // compute optimal triangulation of hole
    // throws InvalidInputException in case of a non-manifold hole.
    void triangulate_hole(Halfedge h);

    // compute the weight of the triangle (i,j,k).
    Weight compute_weight(int i, int j, int k) const;

    // refine triangulation (isotropic remeshing)
    void refine();
    void split_long_edges(const Scalar lmax);
    void collapse_short_edges(const Scalar lmin);
    void flip_edges();
    void relaxation();
    void remove_caps();
    void fairing();

    // return i'th vertex of hole
    Vertex hole_vertex(unsigned int i) const
    {
        assert(i < hole_.size());
        return mesh_.to_vertex(hole_[i]);
    }

    // return normal of face opposite to edge (i-1,i)
    Normal opposite_normal(unsigned int i) const
    {
        assert(i < hole_.size());
        return face_normal(mesh_,
                           mesh_.face(mesh_.opposite_halfedge(hole_[i])));
    }

    // does interior edge (_a,_b) exist already?
    bool is_interior_edge(Vertex a, Vertex b) const;

    // triangle area
    Scalar compute_area(Vertex a, Vertex b, Vertex c) const;

    // triangle normal
    Point compute_normal(Vertex a, Vertex b, Vertex c) const;

    // dihedral angle
    Scalar compute_angle(const Point& n1, const Point& n2) const;

    // mesh and properties
    SurfaceMesh& mesh_;
    VertexProperty<Point> points_;
    VertexProperty<bool> vlocked_;
    EdgeProperty<bool> elocked_;

    std::vector<Halfedge> hole_;

    // data for computing optimal triangulation
    std::vector<std::vector<Weight>> weight_;
    std::vector<std::vector<int>> index_;
};

HoleFilling::HoleFilling(SurfaceMesh& mesh) : mesh_(mesh)
{
    points_ = mesh_.vertex_property<Point>("v:point");
}

bool HoleFilling::is_interior_edge(Vertex a, Vertex b) const
{
    const Halfedge h = mesh_.find_halfedge(a, b);
    if (!h.is_valid())
        return false; // edge does not exist
    return (!mesh_.is_boundary(h) &&
            !mesh_.is_boundary(mesh_.opposite_halfedge(h)));
}

Scalar HoleFilling::compute_area(Vertex a, Vertex b, Vertex c) const
{
    return sqrnorm(cross(points_[b] - points_[a], points_[c] - points_[a]));
}

Point HoleFilling::compute_normal(Vertex a, Vertex b, Vertex c) const
{
    return normalize(cross(points_[b] - points_[a], points_[c] - points_[a]));
}

Scalar HoleFilling::compute_angle(const Point& n1, const Point& n2) const
{
    return (1.0 - dot(n1, n2));
}

void HoleFilling::fill_hole(Halfedge h)
{
    if (!h.is_valid())
    {
        throw InvalidInputException("HoleFilling: Invalid halfedge.");
    }

    // is it really a hole?
    if (!mesh_.is_boundary(h))
    {
        auto what = std::string{__func__} + ": Not a boundary halfedge.";
        throw InvalidInputException(what);
    }

    // lock vertices/edge that already exist, to be later able to
    // identify the filled-in vertices/edges
    vlocked_ = mesh_.add_vertex_property<bool>("HoleFilling:vlocked", false);
    elocked_ = mesh_.add_edge_property<bool>("HoleFilling:elocked", false);
    for (auto v : mesh_.vertices())
        vlocked_[v] = true;
    for (auto e : mesh_.edges())
        elocked_[e] = true;

    try
    {
        triangulate_hole(h); // do minimal triangulation
        refine();            // refine filled-in edges
    }
    catch (InvalidInputException& e)
    {
        // clean up
        hole_.clear();
        mesh_.remove_vertex_property(vlocked_);
        mesh_.remove_edge_property(elocked_);

        throw e;
    }

    // clean up
    hole_.clear();
    mesh_.remove_vertex_property(vlocked_);
    mesh_.remove_edge_property(elocked_);
}

void HoleFilling::triangulate_hole(Halfedge h)
{
    // trace hole
    hole_.clear();
    Halfedge hit = h;
    do
    {
        // check for manifoldness
        if (!mesh_.is_manifold(mesh_.to_vertex(hit)))
        {
            auto what = std::string{__func__} + ": Non-manifold hole.";
            throw InvalidInputException(what);
        }

        hole_.emplace_back(hit);
    } while ((hit = mesh_.next_halfedge(hit)) != h);
    const int n = hole_.size();

    // compute minimal triangulation by dynamic programming
    weight_.clear();
    weight_.resize(n, std::vector<Weight>(n, Weight()));
    index_.clear();
    index_.resize(n, std::vector<int>(n, 0));

    int i, j, m, k, imin;
    Weight w, wmin;

    // initialize 2-gons
    for (i = 0; i < n - 1; ++i)
    {
        weight_[i][i + 1] = Weight(0, 0);
        index_[i][i + 1] = -1;
    }

    // n-gons with n>2
    for (j = 2; j < n; ++j)
    {
        // for all n-gons [i,i+j]
        for (i = 0; i < n - j; ++i)
        {
            k = i + j;
            wmin = Weight();
            imin = -1;

            // find best split i < m < i+j
            for (m = i + 1; m < k; ++m)
            {
                w = weight_[i][m] + compute_weight(i, m, k) + weight_[m][k];
                if (w < wmin)
                {
                    wmin = w;
                    imin = m;
                }
            }

            weight_[i][k] = wmin;
            index_[i][k] = imin;
        }
    }

    // now add triangles to mesh
    std::vector<ivec2> todo;
    todo.reserve(n);
    todo.emplace_back(0, n - 1);
    while (!todo.empty())
    {
        ivec2 tri = todo.back();
        todo.pop_back();
        const int start = tri[0];
        const int end = tri[1];
        if (end - start < 2)
            continue;
        const int split = index_[start][end];

        mesh_.add_triangle(hole_vertex(start), hole_vertex(split),
                           hole_vertex(end));

        todo.emplace_back(start, split);
        todo.emplace_back(split, end);
    }

    // clean up
    weight_.clear();
    index_.clear();
}

HoleFilling::Weight HoleFilling::compute_weight(int i, int j, int k) const
{
    const Vertex a = hole_vertex(i);
    const Vertex b = hole_vertex(j);
    const Vertex c = hole_vertex(k);

    // if one of the potential edges already exists, this would result
    // in an invalid triangulation -> prevent by giving infinite weight
    if (is_interior_edge(a, b) || is_interior_edge(b, c) ||
        is_interior_edge(c, a))
    {
        return {};
    }

    // compute area
    const Scalar area = compute_area(a, b, c);

    // compute dihedral angles with...
    Scalar angle(0);
    const Point n = compute_normal(a, b, c);
    Point n2;

    // ...neighbor to (i,j)
    n2 = (i + 1 == j) ? opposite_normal(j)
                      : compute_normal(a, hole_vertex(index_[i][j]), b);
    angle = std::max(angle, compute_angle(n, n2));

    // ...neighbor to (j,k)
    n2 = (j + 1 == k) ? opposite_normal(k)
                      : compute_normal(b, hole_vertex(index_[j][k]), c);
    angle = std::max(angle, compute_angle(n, n2));

    // ...neighbor to (k,i) if (k,i)==(n-1, 0)
    if (i == 0 && k + 1 == (int)hole_.size())
    {
        n2 = opposite_normal(0);
        angle = std::max(angle, compute_angle(n, n2));
    }

    return {angle, area};
}

void HoleFilling::refine()
{
    // compute target edge length
    const int n = hole_.size();
    Scalar mean_length(0);
    for (int i = 0; i < n; ++i)
    {
        mean_length += distance(points_[hole_vertex(i)],
                                points_[hole_vertex((i + 1) % n)]);
    }
    mean_length /= (Scalar)n;

    // do some iterations
    const auto min_length = Scalar{0.7} * mean_length;
    const auto max_length = Scalar{1.5} * mean_length;
    for (int iter = 0; iter < 10; ++iter)
    {
        split_long_edges(max_length);
        collapse_short_edges(min_length);
        flip_edges();
        relaxation();
    }
    remove_caps();
    fairing();
}

void HoleFilling::split_long_edges(const Scalar lmax)
{
    bool ok;
    int i;

    for (ok = false, i = 0; !ok && i < 10; ++i)
    {
        ok = true;

        for (auto e : mesh_.edges())
        {
            if (!elocked_[e])
            {
                const Halfedge h10 = mesh_.halfedge(e, 0);
                const Halfedge h01 = mesh_.halfedge(e, 1);
                const Point& p0 = points_[mesh_.to_vertex(h10)];
                const Point& p1 = points_[mesh_.to_vertex(h01)];

                if (distance(p0, p1) > lmax)
                {
                    mesh_.split(e, 0.5 * (p0 + p1));
                    ok = false;
                }
            }
        }
    }
}

void HoleFilling::collapse_short_edges(const Scalar _lmin)
{
    bool ok;
    int i;

    for (ok = false, i = 0; !ok && i < 10; ++i)
    {
        ok = true;

        for (auto e : mesh_.edges())
        {
            if (!mesh_.is_deleted(e) && !elocked_[e])
            {
                const Halfedge h10 = mesh_.halfedge(e, 0);
                const Halfedge h01 = mesh_.halfedge(e, 1);
                const Vertex v0 = mesh_.to_vertex(h10);
                const Vertex v1 = mesh_.to_vertex(h01);
                const Point& p0 = points_[v0];
                const Point& p1 = points_[v1];

                // edge too short?
                if (distance(p0, p1) < _lmin)
                {
                    Halfedge h;
                    if (!vlocked_[v0])
                        h = h01;
                    else if (!vlocked_[v1])
                        h = h10;

                    if (h.is_valid() && mesh_.is_collapse_ok(h))
                    {
                        mesh_.collapse(h);
                        ok = false;
                    }
                }
            }
        }
    }

    mesh_.garbage_collection();
}

void HoleFilling::flip_edges()
{
    Vertex v0, v1, v2, v3;
    Halfedge h;
    int val0, val1, val2, val3;
    int val_opt0, val_opt1, val_opt2, val_opt3;
    int ve0, ve1, ve2, ve3, ve_before, ve_after;
    bool ok;
    int i;

    for (ok = false, i = 0; !ok && i < 10; ++i)
    {
        ok = true;

        for (auto e : mesh_.edges())
        {
            if (!elocked_[e])
            {
                h = mesh_.halfedge(e, 0);
                v0 = mesh_.to_vertex(h);
                v2 = mesh_.to_vertex(mesh_.next_halfedge(h));
                h = mesh_.halfedge(e, 1);
                v1 = mesh_.to_vertex(h);
                v3 = mesh_.to_vertex(mesh_.next_halfedge(h));

                val0 = mesh_.valence(v0);
                val1 = mesh_.valence(v1);
                val2 = mesh_.valence(v2);
                val3 = mesh_.valence(v3);

                val_opt0 = (mesh_.is_boundary(v0) ? 4 : 6);
                val_opt1 = (mesh_.is_boundary(v1) ? 4 : 6);
                val_opt2 = (mesh_.is_boundary(v2) ? 4 : 6);
                val_opt3 = (mesh_.is_boundary(v3) ? 4 : 6);

                ve0 = (val0 - val_opt0);
                ve1 = (val1 - val_opt1);
                ve2 = (val2 - val_opt2);
                ve3 = (val3 - val_opt3);

                ve_before = ve0 * ve0 + ve1 * ve1 + ve2 * ve2 + ve3 * ve3;

                --val0;
                --val1;
                ++val2;
                ++val3;

                ve0 = (val0 - val_opt0);
                ve1 = (val1 - val_opt1);
                ve2 = (val2 - val_opt2);
                ve3 = (val3 - val_opt3);

                ve_after = ve0 * ve0 + ve1 * ve1 + ve2 * ve2 + ve3 * ve3;

                if (ve_before > ve_after && mesh_.is_flip_ok(e))
                {
                    mesh_.flip(e);
                    ok = false;
                }
            }
        }
    }
}

void HoleFilling::relaxation()
{
    // properties
    VertexProperty<int> idx =
        mesh_.add_vertex_property<int>("HoleFilling:idx", -1);

    // collect free vertices
    std::vector<Vertex> vertices;
    vertices.reserve(mesh_.n_vertices());
    for (auto v : mesh_.vertices())
    {
        if (!vlocked_[v])
        {
            idx[v] = vertices.size();
            vertices.emplace_back(v);
        }
    }
    const int n = vertices.size();

    // setup matrix & rhs
    Eigen::MatrixXd B(n, 3);
    using Triplet = Eigen::Triplet<double>;
    std::vector<Triplet> triplets;
    for (int i = 0; i < n; ++i)
    {
        const Vertex v = vertices[i];
        Point b(0, 0, 0);
        Scalar c(0);

        for (auto vv : mesh_.vertices(v))
        {
            if (vlocked_[vv])
                b += points_[vv];
            else
                triplets.emplace_back(i, idx[vv], -1.0);
            ++c;
        }

        if (vlocked_[v])
            b -= c * points_[v];
        else
            triplets.emplace_back(i, idx[v], c);

        B.row(i) = (Eigen::Vector3d)b;
    }

    // solve least squares system
    using SparseMatrix = Eigen::SparseMatrix<double>;
    SparseMatrix A(n, n);
    A.setFromTriplets(triplets.begin(), triplets.end());
    const Eigen::SimplicialLDLT<SparseMatrix> solver(A);
    Eigen::MatrixXd X = solver.solve(B);

    if (solver.info() != Eigen::Success)
    {
        // clean up
        mesh_.remove_vertex_property(idx);
        auto what = std::string{__func__} + ": Failed to solve linear system.";
        throw SolverException(what);
    }

    // copy solution to mesh vertices
    for (int i = 0; i < n; ++i)
    {
        points_[vertices[i]] = X.row(i);
    }

    // clean up
    mesh_.remove_vertex_property(idx);
}

void HoleFilling::remove_caps()
{
    Halfedge h;
    Vertex v, vb, vd;
    Scalar a0, a1, amin;
    const Scalar aa(::cos(170.0 * std::numbers::pi / 180.0));
    Point a, b, c, d;

    for (auto e : mesh_.edges())
    {
        if (mesh_.is_flip_ok(e))
        {
            h = mesh_.halfedge(e, 0);
            a = points_[mesh_.to_vertex(h)];
            h = mesh_.next_halfedge(h);
            b = points_[vb = mesh_.to_vertex(h)];
            h = mesh_.halfedge(e, 1);
            c = points_[mesh_.to_vertex(h)];
            h = mesh_.next_halfedge(h);
            d = points_[vd = mesh_.to_vertex(h)];

            a0 = dot(normalize(a - b), normalize(c - b));
            a1 = dot(normalize(a - d), normalize(c - d));

            if (a0 < a1)
            {
                amin = a0;
                v = vb;
            }
            else
            {
                amin = a1;
                v = vd;
            }

            // is it a cap?
            if (amin < aa)
            {
                // flip
                mesh_.flip(e);
            }
        }
    }
}

void HoleFilling::fairing()
{
    // did the refinement insert new vertices?
    // if yes, then trigger fairing; otherwise don't.
    bool new_vertices = false;
    for (auto v : mesh_.vertices())
        if (!vlocked_[v])
            new_vertices = true;
    if (!new_vertices)
        return;

    // convert non-locked into selection
    bool add_vsel = !mesh_.has_vertex_property("v:selected");
    auto vsel = mesh_.vertex_property<bool>("v:selected");
    for (auto v : mesh_.vertices())
        vsel[v] = !vlocked_[v];

    try
    {
        // fair new vertices
        minimize_curvature(mesh_);
    }
    catch (SolverException& e)
    {
        // clean up
        if (add_vsel)
            mesh_.remove_vertex_property(vsel);
        throw e;
    }

    // clean up
    if (add_vsel)
        mesh_.remove_vertex_property(vsel);
}

} // namespace

void fill_hole(SurfaceMesh& mesh, Halfedge h)
{
    HoleFilling(mesh).fill_hole(h);
}
} // 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/hole_filling.h"
5
6	#include <Eigen/Dense>
7	#include <Eigen/Sparse>
8
9	#include <limits>
10	#include <vector>
11
12	#include "pmp/algorithms/fairing.h"
13	#include "pmp/algorithms/normals.h"
14
15	namespace pmp {
16	namespace {
17
18	class HoleFilling
19	{
20	public:
21	explicit HoleFilling(SurfaceMesh& mesh);
22	void fill_hole(Halfedge h);
23
24	private:
25	struct Weight
26	{
27	Weight(Scalar angle = std::numeric_limits<Scalar>::max(),	199✔
28	Scalar area = std::numeric_limits<Scalar>::max())
29	: angle_(angle), area_(area)	199✔
30	{
31	}	199✔
32
33	Weight operator+(const Weight& rhs) const	112✔
34	{
35	return Weight(std::max(angle_, rhs.angle_), area_ + rhs.area_);	112✔
36	}
37
38	bool operator<(const Weight& rhs) const	56✔
39	{
40	return (angle_ < rhs.angle_ \|\|	91✔
41	(angle_ == rhs.angle_ && area_ < rhs.area_));	91✔
42	}
43
44	Scalar angle_;
45	Scalar area_;
46	};
47
48	// compute optimal triangulation of hole
49	// throws InvalidInputException in case of a non-manifold hole.
50	void triangulate_hole(Halfedge h);
51
52	// compute the weight of the triangle (i,j,k).
53	Weight compute_weight(int i, int j, int k) const;
54
55	// refine triangulation (isotropic remeshing)
56	void refine();
57	void split_long_edges(const Scalar lmax);
58	void collapse_short_edges(const Scalar lmin);
59	void flip_edges();
60	void relaxation();
61	void remove_caps();
62	void fairing();
63
64	// return i'th vertex of hole
65	Vertex hole_vertex(unsigned int i) const	272✔
66	{
67	assert(i < hole_.size());	272✔
68	return mesh_.to_vertex(hole_[i]);	272✔
69	}
70
71	// return normal of face opposite to edge (i-1,i)
72	Normal opposite_normal(unsigned int i) const	48✔
73	{
74	assert(i < hole_.size());	48✔
75	return face_normal(mesh_,	48✔
76	mesh_.face(mesh_.opposite_halfedge(hole_[i])));	96✔
77	}
78
79	// does interior edge (_a,_b) exist already?
80	bool is_interior_edge(Vertex a, Vertex b) const;
81
82	// triangle area
83	Scalar compute_area(Vertex a, Vertex b, Vertex c) const;
84
85	// triangle normal
86	Point compute_normal(Vertex a, Vertex b, Vertex c) const;
87
88	// dihedral angle
89	Scalar compute_angle(const Point& n1, const Point& n2) const;
90
91	// mesh and properties
92	SurfaceMesh& mesh_;
93	VertexProperty<Point> points_;
94	VertexProperty<bool> vlocked_;
95	EdgeProperty<bool> elocked_;
96
97	std::vector<Halfedge> hole_;
98
99	// data for computing optimal triangulation
100	std::vector<std::vector<Weight>> weight_;
101	std::vector<std::vector<int>> index_;
102	};
103
104	HoleFilling::HoleFilling(SurfaceMesh& mesh) : mesh_(mesh)	1✔
105	{
106	points_ = mesh_.vertex_property<Point>("v:point");	4✔
107	}	1✔
108
109	bool HoleFilling::is_interior_edge(Vertex a, Vertex b) const	168✔
110	{
111	const Halfedge h = mesh_.find_halfedge(a, b);	168✔
112	if (!h.is_valid())	168✔
113	return false; // edge does not exist	120✔
114	return (!mesh_.is_boundary(h) &&	48✔
115	!mesh_.is_boundary(mesh_.opposite_halfedge(h)));	48✔
116	}
117
118	Scalar HoleFilling::compute_area(Vertex a, Vertex b, Vertex c) const	56✔
119	{
120	return sqrnorm(cross(points_[b] - points_[a], points_[c] - points_[a]));	56✔
121	}
122
123	Point HoleFilling::compute_normal(Vertex a, Vertex b, Vertex c) const	126✔
124	{
125	return normalize(cross(points_[b] - points_[a], points_[c] - points_[a]));	126✔
126	}
127
128	Scalar HoleFilling::compute_angle(const Point& n1, const Point& n2) const	118✔
129	{
130	return (1.0 - dot(n1, n2));	118✔
131	}
132
133	void HoleFilling::fill_hole(Halfedge h)	1✔
134	{
135	if (!h.is_valid())	1✔
136	{
137	throw InvalidInputException("HoleFilling: Invalid halfedge.");	×
138	}
139
140	// is it really a hole?
141	if (!mesh_.is_boundary(h))	1✔
142	{
143	auto what = std::string{__func__} + ": Not a boundary halfedge.";	×
144	throw InvalidInputException(what);	×
145	}	×
146
147	// lock vertices/edge that already exist, to be later able to
148	// identify the filled-in vertices/edges
149	vlocked_ = mesh_.add_vertex_property<bool>("HoleFilling:vlocked", false);	2✔
150	elocked_ = mesh_.add_edge_property<bool>("HoleFilling:elocked", false);	2✔
151	for (auto v : mesh_.vertices())	10✔
152	vlocked_[v] = true;	9✔
153	for (auto e : mesh_.edges())	17✔
154	elocked_[e] = true;	16✔
155
156	try
157	{
158	triangulate_hole(h); // do minimal triangulation	1✔
159	refine(); // refine filled-in edges	1✔
160	}
161	catch (InvalidInputException& e)	×
162	{
163	// clean up
164	hole_.clear();	×
165	mesh_.remove_vertex_property(vlocked_);	×
166	mesh_.remove_edge_property(elocked_);	×
167
168	throw e;	×
169	}	×
170
171	// clean up
172	hole_.clear();	1✔
173	mesh_.remove_vertex_property(vlocked_);	1✔
174	mesh_.remove_edge_property(elocked_);	1✔
175	}	1✔
176
177	void HoleFilling::triangulate_hole(Halfedge h)	1✔
178	{
179	// trace hole
180	hole_.clear();	1✔
181	Halfedge hit = h;	1✔
182	do
183	{
184	// check for manifoldness
185	if (!mesh_.is_manifold(mesh_.to_vertex(hit)))	8✔
186	{
187	auto what = std::string{__func__} + ": Non-manifold hole.";	×
188	throw InvalidInputException(what);	×
189	}	×
190
191	hole_.emplace_back(hit);	8✔
192	} while ((hit = mesh_.next_halfedge(hit)) != h);	8✔
193	const int n = hole_.size();	1✔
194
195	// compute minimal triangulation by dynamic programming
196	weight_.clear();	1✔
197	weight_.resize(n, std::vector<Weight>(n, Weight()));	2✔
198	index_.clear();	1✔
199	index_.resize(n, std::vector<int>(n, 0));	2✔
200
201	int i, j, m, k, imin;
202	Weight w, wmin;	1✔
203
204	// initialize 2-gons
205	for (i = 0; i < n - 1; ++i)	8✔
206	{
207	weight_[i][i + 1] = Weight(0, 0);	7✔
208	index_[i][i + 1] = -1;	7✔
209	}
210
211	// n-gons with n>2
212	for (j = 2; j < n; ++j)	7✔
213	{
214	// for all n-gons [i,i+j]
215	for (i = 0; i < n - j; ++i)	27✔
216	{
217	k = i + j;	21✔
218	wmin = Weight();	21✔
219	imin = -1;	21✔
220
221	// find best split i < m < i+j
222	for (m = i + 1; m < k; ++m)	77✔
223	{
224	w = weight_[i][m] + compute_weight(i, m, k) + weight_[m][k];	56✔
225	if (w < wmin)	56✔
226	{
227	wmin = w;	29✔
228	imin = m;	29✔
229	}
230	}
231
232	weight_[i][k] = wmin;	21✔
233	index_[i][k] = imin;	21✔
234	}
235	}
236
237	// now add triangles to mesh
238	std::vector<ivec2> todo;	1✔
239	todo.reserve(n);	1✔
240	todo.emplace_back(0, n - 1);	1✔
241	while (!todo.empty())	14✔
242	{
243	ivec2 tri = todo.back();	13✔
244	todo.pop_back();	13✔
245	const int start = tri[0];	13✔
246	const int end = tri[1];	13✔
247	if (end - start < 2)	13✔
248	continue;	7✔
249	const int split = index_[start][end];	6✔
250
251	mesh_.add_triangle(hole_vertex(start), hole_vertex(split),	6✔
252	hole_vertex(end));
253
254	todo.emplace_back(start, split);	6✔
255	todo.emplace_back(split, end);	6✔
256	}
257
258	// clean up
259	weight_.clear();	1✔
260	index_.clear();	1✔
261	}	1✔
262
263	HoleFilling::Weight HoleFilling::compute_weight(int i, int j, int k) const	56✔
264	{
265	const Vertex a = hole_vertex(i);	56✔
266	const Vertex b = hole_vertex(j);	56✔
267	const Vertex c = hole_vertex(k);	56✔
268
269	// if one of the potential edges already exists, this would result
270	// in an invalid triangulation -> prevent by giving infinite weight
271	if (is_interior_edge(a, b) \|\| is_interior_edge(b, c) \|\|	112✔
272	is_interior_edge(c, a))	56✔
273	{
274	return {};	×
275	}
276
277	// compute area
278	const Scalar area = compute_area(a, b, c);	56✔
279
280	// compute dihedral angles with...
281	Scalar angle(0);	56✔
282	const Point n = compute_normal(a, b, c);	56✔
283	Point n2;
284
285	// ...neighbor to (i,j)
286	n2 = (i + 1 == j) ? opposite_normal(j)	56✔
287	: compute_normal(a, hole_vertex(index_[i][j]), b);	35✔
288	angle = std::max(angle, compute_angle(n, n2));	56✔
289
290	// ...neighbor to (j,k)
291	n2 = (j + 1 == k) ? opposite_normal(k)	56✔
292	: compute_normal(b, hole_vertex(index_[j][k]), c);	35✔
293	angle = std::max(angle, compute_angle(n, n2));	56✔
294
295	// ...neighbor to (k,i) if (k,i)==(n-1, 0)
296	if (i == 0 && k + 1 == (int)hole_.size())	56✔
297	{
298	n2 = opposite_normal(0);	6✔
299	angle = std::max(angle, compute_angle(n, n2));	6✔
300	}
301
302	return {angle, area};	56✔
303	}
304
305	void HoleFilling::refine()	1✔
306	{
307	// compute target edge length
308	const int n = hole_.size();	1✔
309	Scalar mean_length(0);	1✔
310	for (int i = 0; i < n; ++i)	9✔
311	{
312	mean_length += distance(points_[hole_vertex(i)],	8✔
313	points_[hole_vertex((i + 1) % n)]);	8✔
314	}
315	mean_length /= (Scalar)n;	1✔
316
317	// do some iterations
318	const auto min_length = Scalar{0.7} * mean_length;	1✔
319	const auto max_length = Scalar{1.5} * mean_length;	1✔
320	for (int iter = 0; iter < 10; ++iter)	11✔
321	{
322	split_long_edges(max_length);	10✔
323	collapse_short_edges(min_length);	10✔
324	flip_edges();	10✔
325	relaxation();	10✔
326	}
327	remove_caps();	1✔
328	fairing();	1✔
329	}	1✔
330
331	void HoleFilling::split_long_edges(const Scalar lmax)	10✔
332	{
333	bool ok;
334	int i;
335
336	for (ok = false, i = 0; !ok && i < 10; ++i)	26✔
337	{
338	ok = true;	16✔
339
340	for (auto e : mesh_.edges())	952✔
341	{
342	if (!elocked_[e])	468✔
343	{
344	const Halfedge h10 = mesh_.halfedge(e, 0);	212✔
345	const Halfedge h01 = mesh_.halfedge(e, 1);	212✔
346	const Point& p0 = points_[mesh_.to_vertex(h10)];	212✔
347	const Point& p1 = points_[mesh_.to_vertex(h01)];	212✔
348
349	if (distance(p0, p1) > lmax)	212✔
350	{
351	mesh_.split(e, 0.5 * (p0 + p1));	15✔
352	ok = false;	15✔
353	}
354	}
355	}
356	}
357	}	10✔
358
359	void HoleFilling::collapse_short_edges(const Scalar _lmin)	10✔
360	{
361	bool ok;
362	int i;
363
364	for (ok = false, i = 0; !ok && i < 10; ++i)	24✔
365	{
366	ok = true;	14✔
367
368	for (auto e : mesh_.edges())	822✔
369	{
370	if (!mesh_.is_deleted(e) && !elocked_[e])	404✔
371	{
372	const Halfedge h10 = mesh_.halfedge(e, 0);	180✔
373	const Halfedge h01 = mesh_.halfedge(e, 1);	180✔
374	const Vertex v0 = mesh_.to_vertex(h10);	180✔
375	const Vertex v1 = mesh_.to_vertex(h01);	180✔
376	const Point& p0 = points_[v0];	180✔
377	const Point& p1 = points_[v1];	180✔
378
379	// edge too short?
380	if (distance(p0, p1) < _lmin)	180✔
381	{
382	Halfedge h;	13✔
383	if (!vlocked_[v0])	13✔
384	h = h01;	4✔
385	else if (!vlocked_[v1])	9✔
386	h = h10;	9✔
387
388	if (h.is_valid() && mesh_.is_collapse_ok(h))	13✔
389	{
390	mesh_.collapse(h);	13✔
391	ok = false;	13✔
392	}
393	}
394	}
395	}
396	}
397
398	mesh_.garbage_collection();	10✔
399	}	10✔
400
401	void HoleFilling::flip_edges()	10✔
402	{
403	Vertex v0, v1, v2, v3;	10✔
404	Halfedge h;	10✔
405	int val0, val1, val2, val3;
406	int val_opt0, val_opt1, val_opt2, val_opt3;
407	int ve0, ve1, ve2, ve3, ve_before, ve_after;
408	bool ok;
409	int i;
410
411	for (ok = false, i = 0; !ok && i < 10; ++i)	24✔
412	{
413	ok = true;	14✔
414
415	for (auto e : mesh_.edges())	746✔
416	{
417	if (!elocked_[e])	366✔
418	{
419	h = mesh_.halfedge(e, 0);	142✔
420	v0 = mesh_.to_vertex(h);	142✔
421	v2 = mesh_.to_vertex(mesh_.next_halfedge(h));	142✔
422	h = mesh_.halfedge(e, 1);	142✔
423	v1 = mesh_.to_vertex(h);	142✔
424	v3 = mesh_.to_vertex(mesh_.next_halfedge(h));	142✔
425
426	val0 = mesh_.valence(v0);	142✔
427	val1 = mesh_.valence(v1);	142✔
428	val2 = mesh_.valence(v2);	142✔
429	val3 = mesh_.valence(v3);	142✔
430
431	val_opt0 = (mesh_.is_boundary(v0) ? 4 : 6);	142✔
432	val_opt1 = (mesh_.is_boundary(v1) ? 4 : 6);	142✔
433	val_opt2 = (mesh_.is_boundary(v2) ? 4 : 6);	142✔
434	val_opt3 = (mesh_.is_boundary(v3) ? 4 : 6);	142✔
435
436	ve0 = (val0 - val_opt0);	142✔
437	ve1 = (val1 - val_opt1);	142✔
438	ve2 = (val2 - val_opt2);	142✔
439	ve3 = (val3 - val_opt3);	142✔
440
441	ve_before = ve0 * ve0 + ve1 * ve1 + ve2 * ve2 + ve3 * ve3;	142✔
442
443	--val0;	142✔
444	--val1;	142✔
445	++val2;	142✔
446	++val3;	142✔
447
448	ve0 = (val0 - val_opt0);	142✔
449	ve1 = (val1 - val_opt1);	142✔
450	ve2 = (val2 - val_opt2);	142✔
451	ve3 = (val3 - val_opt3);	142✔
452
453	ve_after = ve0 * ve0 + ve1 * ve1 + ve2 * ve2 + ve3 * ve3;	142✔
454
455	if (ve_before > ve_after && mesh_.is_flip_ok(e))	142✔
456	{
457	mesh_.flip(e);	8✔
458	ok = false;	8✔
459	}
460	}
461	}
462	}
463	}	10✔
464
465	void HoleFilling::relaxation()	10✔
466	{
467	// properties
468	VertexProperty<int> idx =
469	mesh_.add_vertex_property<int>("HoleFilling:idx", -1);	20✔
470
471	// collect free vertices
472	std::vector<Vertex> vertices;	10✔
473	vertices.reserve(mesh_.n_vertices());	10✔
474	for (auto v : mesh_.vertices())	118✔
475	{
476	if (!vlocked_[v])	108✔
477	{
478	idx[v] = vertices.size();	18✔
479	vertices.emplace_back(v);	18✔
480	}
481	}
482	const int n = vertices.size();	10✔
483
484	// setup matrix & rhs
485	Eigen::MatrixXd B(n, 3);	10✔
486	using Triplet = Eigen::Triplet<double>;
487	std::vector<Triplet> triplets;	10✔
488	for (int i = 0; i < n; ++i)	28✔
489	{
490	const Vertex v = vertices[i];	18✔
491	Point b(0, 0, 0);	18✔
492	Scalar c(0);	18✔
493
494	for (auto vv : mesh_.vertices(v))	124✔
495	{
496	if (vlocked_[vv])	106✔
497	b += points_[vv];	90✔
498	else
499	triplets.emplace_back(i, idx[vv], -1.0);	16✔
500	++c;	106✔
501	}
502
503	if (vlocked_[v])	18✔
504	b -= c * points_[v];	×
505	else
506	triplets.emplace_back(i, idx[v], c);	18✔
507
508	B.row(i) = (Eigen::Vector3d)b;	18✔
509	}
510
511	// solve least squares system
512	using SparseMatrix = Eigen::SparseMatrix<double>;
513	SparseMatrix A(n, n);	10✔
514	A.setFromTriplets(triplets.begin(), triplets.end());	10✔
515	const Eigen::SimplicialLDLT<SparseMatrix> solver(A);	10✔
516	Eigen::MatrixXd X = solver.solve(B);	10✔
517
518	if (solver.info() != Eigen::Success)	10✔
519	{
520	// clean up
521	mesh_.remove_vertex_property(idx);	×
522	auto what = std::string{__func__} + ": Failed to solve linear system.";	×
523	throw SolverException(what);	×
524	}	×
525
526	// copy solution to mesh vertices
527	for (int i = 0; i < n; ++i)	28✔
528	{
529	points_[vertices[i]] = X.row(i);	18✔
530	}
531
532	// clean up
533	mesh_.remove_vertex_property(idx);	10✔
534	}	10✔
535
536	void HoleFilling::remove_caps()	1✔
537	{
538	Halfedge h;	1✔
539	Vertex v, vb, vd;	1✔
540	Scalar a0, a1, amin;
541	const Scalar aa(::cos(170.0 * std::numbers::pi / 180.0));	1✔
542	Point a, b, c, d;
543
544	for (auto e : mesh_.edges())	55✔
545	{
546	if (mesh_.is_flip_ok(e))	27✔
547	{
548	h = mesh_.halfedge(e, 0);	27✔
549	a = points_[mesh_.to_vertex(h)];	27✔
550	h = mesh_.next_halfedge(h);	27✔
551	b = points_[vb = mesh_.to_vertex(h)];	27✔
552	h = mesh_.halfedge(e, 1);	27✔
553	c = points_[mesh_.to_vertex(h)];	27✔
554	h = mesh_.next_halfedge(h);	27✔
555	d = points_[vd = mesh_.to_vertex(h)];	27✔
556
557	a0 = dot(normalize(a - b), normalize(c - b));	27✔
558	a1 = dot(normalize(a - d), normalize(c - d));	27✔
559
560	if (a0 < a1)	27✔
561	{
562	amin = a0;	6✔
563	v = vb;	6✔
564	}
565	else
566	{
567	amin = a1;	21✔
568	v = vd;	21✔
569	}
570
571	// is it a cap?
572	if (amin < aa)	27✔
573	{
574	// flip
NEW 575	mesh_.flip(e);	×
576	}
577	}
578	}
579	}	1✔
580
581	void HoleFilling::fairing()	1✔
582	{
583	// did the refinement insert new vertices?
584	// if yes, then trigger fairing; otherwise don't.
585	bool new_vertices = false;	1✔
586	for (auto v : mesh_.vertices())	12✔
587	if (!vlocked_[v])	11✔
588	new_vertices = true;	2✔
589	if (!new_vertices)	1✔
590	return;	×
591
592	// convert non-locked into selection
593	bool add_vsel = !mesh_.has_vertex_property("v:selected");	2✔
594	auto vsel = mesh_.vertex_property<bool>("v:selected");	2✔
595	for (auto v : mesh_.vertices())	12✔
596	vsel[v] = !vlocked_[v];	11✔
597
598	try
599	{
600	// fair new vertices
601	minimize_curvature(mesh_);	1✔
602	}
603	catch (SolverException& e)	×
604	{
605	// clean up
NEW 606	if (add_vsel)	×
NEW 607	mesh_.remove_vertex_property(vsel);	×
608	throw e;	×
609	}	×
610
611	// clean up
612	if (add_vsel)	1✔
613	mesh_.remove_vertex_property(vsel);	1✔
614	}
615
616	} // namespace
617
618	void fill_hole(SurfaceMesh& mesh, Halfedge h)	1✔
619	{
620	HoleFilling(mesh).fill_hole(h);	1✔
621	}	1✔
622	} // namespace pmp

pmp-library / pmp-library / 15742553254

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