GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh_subdivision.cpp
Date: 2026-09-07 02:37:58
Exec Total Coverage
Lines: 140 263 53.2%
Functions: 7 10 70.0%
Branches: 135 498 27.1%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/mesh/mesh_subdivision.h>
41 #include <geogram/mesh/mesh.h>
42 #include <geogram/bibliography/bibliography.h>
43
44 namespace GEO {
45
46 22 MeshSplitCallbacks::MeshSplitCallbacks(Mesh* mesh) : mesh_(mesh) {
47 22 }
48
49 44 MeshSplitCallbacks::~MeshSplitCallbacks() {
50 44 }
51
52 ✗ index_t MeshSplitCallbacks::create_vertex() {
53 ✗ index_t result = mesh_->vertices.create_vertex();
54 ✗ mesh_->vertices.attributes().zero_item(result);
55 ✗ return result;
56 }
57
58 5120 void MeshSplitCallbacks::scale_vertex(index_t v, double s) {
59 5120 mesh_->vertices.attributes().scale_item(v,s);
60 5120 }
61
62 209840 void MeshSplitCallbacks::zero_vertex(index_t v) {
63 209840 mesh_->vertices.attributes().zero_item(v);
64 209840 }
65
66 429900 void MeshSplitCallbacks::madd_vertex(index_t v1, double s, index_t v2) {
67 429900 mesh_->vertices.attributes().madd_item(v1,s,v2);
68 429900 }
69
70 /*************************************************************************/
71
72 17 void mesh_split_triangles(
73 Mesh& M, index_t facets_begin, index_t facets_end,
74 MeshSplitCallbacks* cb
75 ) {
76
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 17 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 17 geo_assert(M.facets.are_simplices());
77
78
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 17 MeshSplitCallbacks default_cb(&M);
79
1/2
✓ Branch 0 taken 17 times.
✗ Branch 1 not taken.
 17 if(cb == nullptr) {
80 17 cb = &default_cb;
81 }
82
83
1/2
✓ Branch 0 taken 17 times.
✗ Branch 1 not taken.
 17 if(facets_end == NO_INDEX) {
84 17 facets_end = M.facets.nb();
85 }
86
87 17 index_t nv0 = M.vertices.nb();
88 17 index_t nf0 = M.facets.nb();
89
90 // Compute corner to new vertex mapping
91
1/2
✓ Branch 2 taken 17 times.
✗ Branch 3 not taken.
 17 vector<index_t> ctov(M.facet_corners.nb(), NO_VERTEX);
92 17 index_t nbnewv=0;
93
2/2
✓ Branch 0 taken 61420 times.
✓ Branch 1 taken 17 times.
 61437 for(index_t f=facets_begin; f<facets_end; ++f) {
94
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 245680 for(index_t c=M.facets.corners_begin(f);
95
3/4
✓ Branch 1 taken 245680 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 184260 times.
✓ Branch 4 taken 61420 times.
 245680 c<M.facets.corners_end(f); ++c
96 ) {
97
3/4
✓ Branch 1 taken 184260 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 92130 times.
✓ Branch 4 taken 92130 times.
 184260 if(ctov[c] == NO_INDEX) {
98
1/2
✓ Branch 1 taken 92130 times.
✗ Branch 2 not taken.
 92130 ctov[c] = nbnewv;
99
1/2
✓ Branch 1 taken 92130 times.
✗ Branch 2 not taken.
 92130 index_t f2 = M.facet_corners.adjacent_facet(c);
100
1/2
✓ Branch 0 taken 92130 times.
✗ Branch 1 not taken.
 92130 if(f2 != NO_FACET) {
101
1/2
✓ Branch 1 taken 92130 times.
✗ Branch 2 not taken.
 182847 for(index_t c2=M.facets.corners_begin(f2);
102
2/4
✓ Branch 1 taken 182847 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 182847 times.
✗ Branch 4 not taken.
 182847 c2!=M.facets.corners_end(f2); ++c2) {
103
3/4
✓ Branch 1 taken 182847 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 92130 times.
✓ Branch 4 taken 90717 times.
 182847 if(M.facet_corners.adjacent_facet(c2) == f) {
104
1/2
✓ Branch 1 taken 92130 times.
✗ Branch 2 not taken.
 92130 ctov[c2] = nbnewv;
105 92130 break;
106 }
107 }
108 }
109 92130 ++nbnewv;
110 }
111 }
112 }
113
114 // Create vertices
115
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 17 M.vertices.create_vertices(nbnewv);
116
2/2
✓ Branch 0 taken 61420 times.
✓ Branch 1 taken 17 times.
 61437 for(index_t f=facets_begin; f<facets_end; ++f) {
117
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 245680 for(index_t c1=M.facets.corners_begin(f);
118
3/4
✓ Branch 1 taken 245680 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 184260 times.
✓ Branch 4 taken 61420 times.
 245680 c1<M.facets.corners_end(f); ++c1
119 ) {
120
1/2
✓ Branch 1 taken 184260 times.
✗ Branch 2 not taken.
 184260 index_t c2 = M.facets.next_corner_around_facet(f,c1);
121
1/2
✓ Branch 1 taken 184260 times.
✗ Branch 2 not taken.
 184260 index_t v1 = M.facet_corners.vertex(c1);
122
1/2
✓ Branch 1 taken 184260 times.
✗ Branch 2 not taken.
 184260 index_t v2 = M.facet_corners.vertex(c2);
123
1/2
✓ Branch 1 taken 184260 times.
✗ Branch 2 not taken.
 184260 index_t v12 = ctov[c1] + nv0;
124
1/2
✓ Branch 1 taken 184260 times.
✗ Branch 2 not taken.
 184260 cb->zero_vertex(v12);
125
1/2
✓ Branch 1 taken 184260 times.
✗ Branch 2 not taken.
 184260 cb->madd_vertex(v12, 0.5, v1);
126
1/2
✓ Branch 1 taken 184260 times.
✗ Branch 2 not taken.
 184260 cb->madd_vertex(v12, 0.5, v2);
127 }
128 }
129
130 // Create facets
131
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 17 M.facets.create_triangles(3*(facets_end - facets_begin));
132
2/2
✓ Branch 0 taken 61420 times.
✓ Branch 1 taken 17 times.
 61437 for(index_t f=facets_begin; f<facets_end; ++f) {
133
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 index_t v1 = M.facets.vertex(f,0);
134
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 index_t v2 = M.facets.vertex(f,1);
135
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 index_t v3 = M.facets.vertex(f,2);
136
2/4
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61420 times.
✗ Branch 5 not taken.
 61420 index_t v12 = ctov[M.facets.corners_begin(f) ] + nv0;
137
2/4
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61420 times.
✗ Branch 5 not taken.
 61420 index_t v23 = ctov[M.facets.corners_begin(f) + 1] + nv0;
138
2/4
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61420 times.
✗ Branch 5 not taken.
 61420 index_t v31 = ctov[M.facets.corners_begin(f) + 2] + nv0;
139
140
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(f,0,v31);
141
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(f,1,v12);
142
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(f,2,v23);
143
144
1/2
✓ Branch 2 taken 61420 times.
✗ Branch 3 not taken.
 61420 M.facets.attributes().copy_item(nf0+3*(f-facets_begin),f);
145
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(nf0+3*(f-facets_begin),0,v1);
146
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(nf0+3*(f-facets_begin),1,v12);
147
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(nf0+3*(f-facets_begin),2,v31);
148
149
1/2
✓ Branch 2 taken 61420 times.
✗ Branch 3 not taken.
 61420 M.facets.attributes().copy_item(nf0+3*(f-facets_begin)+1,f);
150
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(nf0+3*(f-facets_begin)+1,0,v12);
151
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(nf0+3*(f-facets_begin)+1,1,v2);
152
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(nf0+3*(f-facets_begin)+1,2,v23);
153
154
1/2
✓ Branch 2 taken 61420 times.
✗ Branch 3 not taken.
 61420 M.facets.attributes().copy_item(nf0+3*(f-facets_begin)+2,f);
155
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(nf0+3*(f-facets_begin)+2,0,v31);
156
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(nf0+3*(f-facets_begin)+2,1,v23);
157
1/2
✓ Branch 1 taken 61420 times.
✗ Branch 2 not taken.
 61420 M.facets.set_vertex(nf0+3*(f-facets_begin)+2,2,v3);
158 }
159
1/2
✓ Branch 2 taken 17 times.
✗ Branch 3 not taken.
 17 M.facets.connect(facets_begin, M.facets.nb());
160 17 }
161
162
163 5 void mesh_split_quads(
164 Mesh& M, index_t facets_begin, index_t facets_end,
165 MeshSplitCallbacks* cb
166 ) {
167
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 MeshSplitCallbacks default_cb(&M);
168
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if(cb == nullptr) {
169 5 cb = &default_cb;
170 }
171
172
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if(facets_end == NO_INDEX) {
173 5 facets_end = M.facets.nb();
174 }
175
176 5 index_t nv0 = M.vertices.nb();
177 5 index_t nf0 = M.facets.nb();
178
179 // Compute corner to new vertex and facet to new vertex
180 // mappings.
181
182
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 vector<index_t> ctov(M.facet_corners.nb(), NO_VERTEX);
183
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 vector<index_t> ftov(M.facets.nb(), NO_VERTEX);
184
185 5 index_t nbnewv=0;
186 5 index_t nbnewf=0;
187
188
2/2
✓ Branch 0 taken 5120 times.
✓ Branch 1 taken 5 times.
 5125 for(index_t f=facets_begin; f<facets_end; ++f) {
189
1/2
✓ Branch 1 taken 5120 times.
✗ Branch 2 not taken.
 5120 ftov[f] = nbnewv;
190 5120 ++nbnewv;
191
1/2
✓ Branch 1 taken 5120 times.
✗ Branch 2 not taken.
 25580 for(index_t c=M.facets.corners_begin(f);
192
3/4
✓ Branch 1 taken 25580 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 20460 times.
✓ Branch 4 taken 5120 times.
 25580 c<M.facets.corners_end(f); ++c
193 ) {
194 20460 ++nbnewf;
195
3/4
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10230 times.
✓ Branch 4 taken 10230 times.
 20460 if(ctov[c] == NO_INDEX) {
196
1/2
✓ Branch 1 taken 10230 times.
✗ Branch 2 not taken.
 10230 ctov[c] = nbnewv;
197
1/2
✓ Branch 1 taken 10230 times.
✗ Branch 2 not taken.
 10230 index_t f2 = M.facet_corners.adjacent_facet(c);
198
1/2
✓ Branch 0 taken 10230 times.
✗ Branch 1 not taken.
 10230 if(f2 != NO_FACET) {
199
1/2
✓ Branch 1 taken 10230 times.
✗ Branch 2 not taken.
 26821 for(index_t c2=M.facets.corners_begin(f2);
200
2/4
✓ Branch 1 taken 26821 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26821 times.
✗ Branch 4 not taken.
 26821 c2!=M.facets.corners_end(f2); ++c2) {
201
3/4
✓ Branch 1 taken 26821 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10230 times.
✓ Branch 4 taken 16591 times.
 26821 if(M.facet_corners.adjacent_facet(c2) == f) {
202
1/2
✓ Branch 1 taken 10230 times.
✗ Branch 2 not taken.
 10230 ctov[c2] = nbnewv;
203 10230 break;
204 }
205 }
206 }
207 10230 ++nbnewv;
208 }
209 }
210 }
211
212 // Create vertices
213
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 M.vertices.create_vertices(nbnewv);
214
2/2
✓ Branch 0 taken 5120 times.
✓ Branch 1 taken 5 times.
 5125 for(index_t f=facets_begin; f<facets_end; ++f) {
215
2/4
✓ Branch 1 taken 5120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5120 times.
✗ Branch 5 not taken.
 5120 cb->zero_vertex(ftov[f] + nv0);
216
1/2
✓ Branch 1 taken 5120 times.
✗ Branch 2 not taken.
 25580 for(index_t c1=M.facets.corners_begin(f);
217
3/4
✓ Branch 1 taken 25580 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 20460 times.
✓ Branch 4 taken 5120 times.
 25580 c1<M.facets.corners_end(f); ++c1
218 ) {
219
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 index_t c2 = M.facets.next_corner_around_facet(f,c1);
220
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 index_t v1 = M.facet_corners.vertex(c1);
221
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 index_t v2 = M.facet_corners.vertex(c2);
222
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 index_t v12 = ctov[c1] + nv0;
223
224
2/4
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20460 times.
✗ Branch 5 not taken.
 20460 cb->madd_vertex(ftov[f] + nv0, 1.0, v1);
225
226
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 cb->zero_vertex(v12);
227
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 cb->madd_vertex(v12,0.5,v1);
228
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 cb->madd_vertex(v12,0.5,v2);
229 }
230
1/2
✓ Branch 1 taken 5120 times.
✗ Branch 2 not taken.
 5120 double s = 1.0 / double(M.facets.nb_vertices(f));
231
2/4
✓ Branch 1 taken 5120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5120 times.
✗ Branch 5 not taken.
 5120 cb->scale_vertex(ftov[f]+nv0, s);
232 }
233
234 // Create facets
235
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 M.facets.create_quads(nbnewf);
236 5 index_t cur_f = 0;
237
2/2
✓ Branch 0 taken 5120 times.
✓ Branch 1 taken 5 times.
 5125 for(index_t f=facets_begin; f<facets_end; ++f) {
238 25580 for(
239
1/2
✓ Branch 1 taken 5120 times.
✗ Branch 2 not taken.
 5120 index_t c1=M.facets.corners_begin(f);
240
3/4
✓ Branch 1 taken 25580 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 20460 times.
✓ Branch 4 taken 5120 times.
 25580 c1<M.facets.corners_end(f); ++c1) {
241
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 index_t c0 = M.facets.prev_corner_around_facet(f,c1);
242
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 index_t v = M.facet_corners.vertex(c1);
243
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 index_t v1 = ctov[c0] + nv0;
244 20460 index_t v2 = v;
245
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 index_t v3 = ctov[c1] + nv0;
246
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 index_t v4 = ftov[f] + nv0;
247
1/2
✓ Branch 2 taken 20460 times.
✗ Branch 3 not taken.
 20460 M.facets.attributes().copy_item(cur_f + nf0, f);
248
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 M.facets.set_vertex(cur_f + nf0, 0, v1);
249
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 M.facets.set_vertex(cur_f + nf0, 1, v2);
250
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 M.facets.set_vertex(cur_f + nf0, 2, v3);
251
1/2
✓ Branch 1 taken 20460 times.
✗ Branch 2 not taken.
 20460 M.facets.set_vertex(cur_f + nf0, 3, v4);
252 20460 ++cur_f;
253 }
254 }
255
256
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 vector<index_t> to_delete(M.facets.nb(),0);
257
2/2
✓ Branch 0 taken 5120 times.
✓ Branch 1 taken 5 times.
 5125 for(index_t f=facets_begin; f<facets_end; ++f) {
258
1/2
✓ Branch 1 taken 5120 times.
✗ Branch 2 not taken.
 5120 to_delete[f] = 1;
259 }
260
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 M.facets.delete_elements(to_delete);
261
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 M.facets.connect();
262 5 }
263
264
265 ✗ void mesh_triangulate_center_vertex(
266 Mesh& M, index_t facets_begin, index_t facets_end,
267 MeshSplitCallbacks* cb
268 ) {
269 ✗ MeshSplitCallbacks default_cb(&M);
270 ✗ if(cb == nullptr) {
271 ✗ cb = &default_cb;
272 }
273
274 ✗ if(facets_end == NO_INDEX) {
275 ✗ facets_end = M.facets.nb();
276 }
277
278 ✗ index_t nv0 = M.vertices.nb();
279 ✗ index_t nf0 = M.facets.nb();
280
281 // Compute corner to new vertex and facet to new vertex
282 // mappings.
283
284 ✗ vector<index_t> ftov(M.facets.nb(), NO_VERTEX);
285
286 ✗ index_t nbnewv=0;
287 ✗ index_t nbnewf=0;
288 ✗ for(index_t f=facets_begin; f<facets_end; ++f) {
289 ✗ ftov[f] = nbnewv;
290 ✗ ++nbnewv;
291 ✗ nbnewf += M.facets.nb_vertices(f);
292 }
293
294 // Create vertices
295 ✗ M.vertices.create_vertices(nbnewv);
296 ✗ for(index_t f=facets_begin; f<facets_end; ++f) {
297 ✗ cb->zero_vertex(ftov[f] + nv0);
298 ✗ for(index_t c1=M.facets.corners_begin(f);
299 ✗ c1<M.facets.corners_end(f); ++c1
300 ) {
301 ✗ index_t v1 = M.facet_corners.vertex(c1);
302 ✗ cb->madd_vertex(ftov[f]+nv0, 1.0, v1);
303 }
304 ✗ double s = 1.0 / double(M.facets.nb_vertices(f));
305 ✗ cb->scale_vertex(ftov[f]+nv0, s);
306 }
307
308 // Create facets
309 ✗ M.facets.create_triangles(nbnewf);
310 ✗ index_t cur_f = 0;
311 ✗ for(index_t f=facets_begin; f<facets_end; ++f) {
312 ✗ for(
313 ✗ index_t c1=M.facets.corners_begin(f);
314 ✗ c1<M.facets.corners_end(f); ++c1
315 ) {
316 ✗ index_t c2 = M.facets.next_corner_around_facet(f,c1);
317 ✗ index_t v1 = M.facet_corners.vertex(c1);
318 ✗ index_t v2 = M.facet_corners.vertex(c2);
319 ✗ M.facets.attributes().copy_item(cur_f + nf0, f);
320 ✗ M.facets.set_vertex(cur_f + nf0, 0, v1);
321 ✗ M.facets.set_vertex(cur_f + nf0, 1, v2);
322 ✗ M.facets.set_vertex(cur_f + nf0, 2, nv0+ftov[f]);
323 ✗ ++cur_f;
324 }
325 }
326
327 ✗ vector<index_t> to_delete(M.facets.nb(),0);
328 ✗ for(index_t f=facets_begin; f<facets_end; ++f) {
329 ✗ to_delete[f] = 1;
330 }
331 ✗ M.facets.delete_elements(to_delete);
332 ✗ M.facets.connect();
333 ✗ }
334
335 ✗ void mesh_split_catmull_clark(Mesh& M, MeshSplitCallbacks* cb) {
336
337 ✗ geo_cite("journals/CAD/CatmullRGB");
338
339 ✗ MeshSplitCallbacks default_cb(&M);
340 ✗ if(cb == nullptr) {
341 ✗ cb = &default_cb;
342 }
343
344 ✗ vector<index_t> vertex_degree(M.vertices.nb(),0);
345 ✗ vector<index_t> corner_vertex(M.facet_corners.nb(),NO_VERTEX);
346 ✗ vector<index_t> facet_vertex(M.facets.nb(), NO_VERTEX);
347 ✗ std::vector<bool> v_on_border(M.vertices.nb(),false);
348
349 ✗ index_t nb_v_orig = M.vertices.nb();
350 ✗ index_t nb_f_orig = M.facets.nb();
351
352 // Create edge and facet vertices
353 ✗ for(index_t f1: M.facets) {
354 ✗ facet_vertex[f1] = cb->create_vertex();
355 ✗ for(index_t c1: M.facets.corners(f1)) {
356 ✗ index_t v = M.facet_corners.vertex(c1);
357 ✗ ++vertex_degree[v];
358 ✗ index_t f2 = M.facet_corners.adjacent_facet(c1);
359
360 ✗ if(f1 < f2 || f2 == NO_FACET) {
361 ✗ corner_vertex[c1] = cb->create_vertex();
362 ✗ if(f2 != NO_FACET) {
363 ✗ index_t cn = M.facets.next_corner_around_facet(f1,c1);
364 ✗ index_t v2 = M.facet_corners.vertex(cn);
365 ✗ index_t c2 = NO_CORNER;
366 ✗ for(c2=M.facets.corners_begin(f2);
367 ✗ c2 != M.facets.corners_end(f2); ++c2) {
368 ✗ if(M.facet_corners.vertex(c2) == v2) {
369 ✗ break;
370 }
371 }
372 ✗ geo_assert(M.facet_corners.vertex(c2) == v2);
373 ✗ corner_vertex[c2] = corner_vertex[c1];
374 }
375 }
376 }
377 }
378
379 // Compute facet vertices
380 ✗ for(index_t f: M.facets) {
381 ✗ double f_degree = double(M.facets.nb_vertices(f));
382 ✗ for(index_t c: M.facets.corners(f)) {
383 ✗ index_t v = M.facet_corners.vertex(c);
384 ✗ cb->madd_vertex(facet_vertex[f], 1.0 / f_degree, v);
385 ✗ if(M.facet_corners.adjacent_facet(c) == NO_FACET) {
386 ✗ v_on_border[v] = true;
387 }
388 }
389 }
390
391 // Compute edge vertices
392 ✗ for(index_t f: M.facets) {
393 ✗ for(index_t c: M.facets.corners(f)) {
394 ✗ index_t v = M.facet_corners.vertex(c);
395 ✗ if(M.facet_corners.adjacent_facet(c) == NO_FACET) {
396 ✗ cb->madd_vertex(corner_vertex[c], 1.0/2.0, v);
397 ✗ index_t c2 = M.facets.next_corner_around_facet(f,c);
398 ✗ index_t v2 = M.facet_corners.vertex(c2);
399 ✗ cb->madd_vertex(corner_vertex[c], 1.0/2.0, v2);
400 } else {
401 ✗ cb->madd_vertex(corner_vertex[c], 0.25, v);
402 ✗ cb->madd_vertex(corner_vertex[c], 0.25, facet_vertex[f]);
403 }
404 }
405 }
406
407 // Compute new position of original vertices
408
409 ✗ FOR(v, nb_v_orig) {
410 ✗ if(v_on_border[v]) {
411 ✗ continue;
412 }
413 ✗ double n = double(vertex_degree[v]);
414 ✗ if(n != 0.0) {
415 ✗ if(!v_on_border[v]) {
416 ✗ cb->scale_vertex(v, (n - 3.0) / n);
417 }
418 }
419 }
420
421 ✗ for(index_t f: M.facets) {
422 ✗ for(index_t c: M.facets.corners(f)) {
423 ✗ index_t v = M.facet_corners.vertex(c);
424 ✗ double n = double(vertex_degree[v]);
425
426 // As compared to original Catmull-Clark documentation:
427 // add 4.0 times edge vertex
428 // then remove contribution of facet vertices
429 // (this retrieves the original edges barycenters without
430 // needing intermediary storage).
431
432 ✗ if(!v_on_border[v]) {
433 ✗ index_t f2 = M.facet_corners.adjacent_facet(c);
434 ✗ cb->madd_vertex(v, 4.0 / (n*n), corner_vertex[c]);
435 ✗ cb->madd_vertex(v, -1.0 / (n*n), facet_vertex[f2]);
436 }
437 }
438 }
439
440 // Create new facets
441 ✗ FOR(f, nb_f_orig) {
442 ✗ for(index_t c: M.facets.corners(f)) {
443 ✗ index_t v = M.facet_corners.vertex(c);
444 ✗ index_t c2 = M.facets.prev_corner_around_facet(f,c);
445 ✗ index_t new_f = M.facets.create_quad(
446 ✗ corner_vertex[c2],
447 v,
448 ✗ corner_vertex[c],
449 ✗ facet_vertex[f]
450 );
451 ✗ M.facets.attributes().copy_item(new_f, f);
452 }
453 }
454
455 // Delete old facets
456 ✗ vector<index_t> delete_f(nb_f_orig, 1);
457 ✗ delete_f.resize(M.facets.nb(),0);
458 ✗ M.facets.delete_elements(delete_f);
459 ✗ M.facets.connect();
460 ✗ }
461 }
462