GCC Code Coverage Report

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 <exploragram/hexdom/quad_dominant.h>
41			#include <exploragram/hexdom/meshcomesh.h>
42			#include <exploragram/hexdom/PGP.h>
43			#include <exploragram/hexdom/basic.h>
44			#include <exploragram/hexdom/extra_connectivity.h>
45			#include <exploragram/hexdom/mesh_inspector.h>
46			#include <exploragram/hexdom/intersect_tools.h>
47			#include <exploragram/hexdom/polygon.h>
48
49			#include <geogram/numerics/matrix_util.h>
50			#include <geogram/basic/permutation.h>
51			#include <geogram/mesh/triangle_intersection.h>
52			#include <geogram/mesh/mesh_tetrahedralize.h>
53			#include <geogram/delaunay/delaunay.h>
54			#include <geogram/points/nn_search.h>
55			#include <geogram/points/colocate.h>
56			#include <geogram/mesh/mesh_repair.h>
57			#include <geogram/mesh/mesh_fill_holes.h>
58			#include <geogram/mesh/mesh_geometry.h>
59
60
61			#include <algorithm>
62			#include <queue>
63			#include <stack>
64			#include <map>
65
66			namespace GEO {
67
68		✗	void export_boundary_with_uv(Mesh* m, Mesh* hex, const char* uv_name, const char* singular_name) {
69		✗	Attribute<GEO::vec2> uv(hex->facet_corners.attributes(), uv_name);
70		✗	Attribute<index_t> singtri(hex->facets.attributes(), singular_name);
71
72		✗	Attribute<bool> has_param(m->cell_facets.attributes(), "has_param");
73		✗	Attribute<vec3> UC(m->cell_corners.attributes(), "U");
74
75			// STEP 1: compute number of boundary vertices AND create a lookup table tet mesh vertex -> boundary surface vertex
76		✗	vector<index_t> tetV_to_facetV(m->vertices.nb(), NOT_AN_ID);
77		✗	index_t nb_boundary_V = 0;
78		✗	FOR(c, m->cells.nb()) FOR(cf, m->cells.nb_facets(c)) {
79		✗	if (NO_CELL != m->cells.adjacent(c, cf)) continue;
80		✗	FOR(cfv, m->cells.facet_nb_vertices(c, cf)) {
81		✗	index_t v = m->cells.facet_vertex(c, cf, cfv);
82		✗	if (NOT_AN_ID == tetV_to_facetV[v]) {
83		✗	tetV_to_facetV[v] = nb_boundary_V++;
84			}
85			}
86			}
87
88			// STEP 2: copy boundary vertices to the new mesh
89		✗	hex->vertices.create_vertices(nb_boundary_V);
90		✗	FOR(v, m->vertices.nb())
91		✗	if (NOT_AN_ID != tetV_to_facetV[v])
92		✗	hex->vertices.point(tetV_to_facetV[v]) = m->vertices.point(v);
93
94			// STEP 3: extract triangles and their parameterization (when possible)
95		✗	FOR(c, m->cells.nb()) FOR(cf, m->cells.nb_facets(c)) {
96		✗	if (m->cells.adjacent(c, cf) != NO_CELL) continue;
97			index_t tet_verts[3];
98		✗	FOR(cfv, 3) {
99		✗	tet_verts[cfv] = m->cells.facet_vertex(c, cf, cfv);
100			}
101
102			index_t tet_corners[3];
103		✗	FOR(cfv, 3) FOR(test, 4) if (m->cell_corners.vertex(m->cells.corner(c, test)) == tet_verts[cfv])
104		✗	tet_corners[cfv] = m->cells.corner(c, test);
105
106
107		✗	vec3 lX[3], lU[3];
108		✗	FOR(cfv, 3) {
109		✗	lX[cfv] = m->vertices.point(tet_verts[cfv]);
110		✗	lU[cfv] = UC[tet_corners[cfv]];
111			}
112
113			// STEP 2.2: non singular boundary triangles are easy to extract
114		✗	index_t fid = hex->facets.create_triangle(
115		✗	tetV_to_facetV[tet_verts[0]],
116		✗	tetV_to_facetV[tet_verts[1]],
117		✗	tetV_to_facetV[tet_verts[2]]
118			);
119
120		✗	bool has_valid_2d_param = false;
121		✗	if (has_param[m->cells.facet(c, cf)]) {
122		✗	FOR(dim, 3) { // we are looking for the param dimension that goes inside the volume
123		✗	if (lU[0][dim] != lU[1][dim] || lU[0][dim] != lU[2][dim]) continue;
124		✗	has_valid_2d_param = true;
125		✗	FOR(lv, 3) {
126		✗	uv[hex->facets.corner(fid, lv)] = vec2(lU[lv][(dim + 1) % 3], lU[lv][(dim + 2) % 3]);
127			}
128			}
129			#if 0
130			// mirror the copy when necessary
131			if (det(uv[hex->facets.corner(fid, 1)] - uv[hex->facets.corner(fid, 0)],
132			uv[hex->facets.corner(fid, 2)] - uv[hex->facets.corner(fid, 0)]) < 0)
133			FOR(lv, 3) uv[hex->facets.corner(fid, lv)][0] *= -1.;
134			#endif
135			}
136
137		✗	if (has_valid_2d_param) { // check param quality and invalidate it, if necessary
138		✗	TrglGradient grd(lX[0], lX[1], lX[2]);
139		✗	vec3 grduv[2];
140		✗	FOR(d, 2) grduv[d] = grd.gradient_3d(uv[hex->facets.corner(fid, 0)][d], uv[hex->facets.corner(fid, 1)][d], uv[hex->facets.corner(fid, 2)][d]);
141		✗	FOR(d, 2) FOR(dd, 3) if (Numeric::is_nan(grduv[d][dd])) has_valid_2d_param = false;
142			#if 0
143			if (grduv[0].length() > 10. * grduv[1].length()) has_valid_2d_param = false;
144			if (grduv[1].length() > 10. * grduv[0].length()) has_valid_2d_param = false;
145			if (std::abs(dot(normalize(grduv[0]), normalize(grduv[1]))) > cos(M_PI / 4.)) has_valid_2d_param = false;
146			#endif
147			}
148
149		✗	singtri[fid] = !has_valid_2d_param;
150			}
151		✗	}
152
153		✗	static bool triangulate_surface_preserve_attributes(Mesh* m) {
154		✗	bool result = true;
155		✗	plop("triangulating surface after embedding isoUV");
156		✗	vector<index_t> to_kill(m->facets.nb(), false);
157		✗	FOR(f, m->facets.nb()) {
158		✗	geo_assert(m->facets.nb_corners(f) >= 3);
159		✗	if (m->facets.nb_corners(f) == 3) continue;
160
161		✗	vector<vec3> pts;
162		✗	FOR(lc, m->facets.nb_corners(f)) {
163		✗	pts.push_back(X(m)[m->facets.vertex(f, lc)]);
164			}
165
166		✗	vector<index_t> triangles;
167		✗	bool success = Poly3d(pts).try_triangulate_minweight(triangles);
168		✗	result = result && success;
169
170		✗	geo_assert(0 == triangles.size() % 3);
171
172		✗	FOR(new_face, triangles.size() / 3) {
173		✗	index_t new_f = m->facets.create_triangle(
174		✗	m->facets.vertex(f, triangles[3 * new_face + 0]),
175		✗	m->facets.vertex(f, triangles[3 * new_face + 1]),
176		✗	m->facets.vertex(f, triangles[3 * new_face + 2]));
177
178		✗	to_kill.push_back(false);
179		✗	m->facets.attributes().copy_item(new_f, f);
180		✗	FOR(nfc,3) {
181		✗	m->facet_corners.attributes().copy_item(m->facets.corner(new_f, nfc), m->facets.corner(f, triangles[3 * new_face + nfc]));
182			}
183			}
184		✗	to_kill[f] = true;
185		✗	}
186		✗	m->facets.delete_elements(to_kill);
187		✗	return result;
188		✗	}
189
190
191			/**
192			* INPUT: facets with xyz geometry and uv coordinates s.t. no edge is of size 0
193			singular bool per triangle
194			* integer values of uv --- interpolated along edge 'e' --- matches on both side of 'e'
195			on presuppose que les U ont ete pre-snappes sur la grille entiere
196			* OUTPUT: facets with uv coordinates s.t.
197			* a new vertex is inserted (with interpolated xyz and uv) at every integer values of uv along edges
198			* nothing prevents some vertices around a facet to share the same uv's
199			ATTENTION: it modifies the parameterization when snaps grid corners on edges!
200			*/
201
202		✗	void split_edges_by_iso_uvs(Mesh* m, const char *uv_name, const char *singular_name) {
203		✗	Attribute<GEO::vec2> uv(m->facet_corners.attributes(), uv_name);
204		✗	Attribute<index_t> singular(m->facets.attributes(), singular_name);
205
206		✗	index_t nb_init_facets = m->facets.nb();
207		✗	Attribute<bool> added_vertices(m->vertices.attributes(), "added_vertices");
208		✗	Attribute<index_t> orig_tri_fid(m->facets.attributes(), "orig_tri_fid");
209		✗	Attribute<index_t> resp_facet(m->vertices.attributes(), "resp_facet");
210
211			typedef std::pair<index_t /*v_id*/, vec2 /*uv*/> NewCorner;
212		✗	vector<vector<NewCorner> > new_corners(m->facet_corners.nb());
213
214		✗	reach("create vertices and init arrays of NewCorner to insert on edges");
215			{
216		✗	FacetsExtraConnectivity fec(m);
217		✗	FOR(h, m->facet_corners.nb()) {
218		✗	index_t opp = fec.opposite(h);
219
220		✗	if (singular[fec.facet(h)]) continue; // the face is responsible for the edge if opp<h or if its opposite is singular
221		✗	if (NOT_AN_ID != opp && !singular[fec.facet(opp)] && opp > h) continue;
222
223		✗	vec2 lU[2] = { uv[h], uv[fec.next(h)] };
224
225			#if 0
226			bool trans_func_found = false;
227			index_t R;
228			vec2 T;
229			if (NOT_AN_ID != opp && !singular[fec.facet(opp)]) {
230			vec2 A = lU[1]-lU[0], B = uv[opp]-uv[fec.next(opp)], C = lU[0];
231			if (std::abs(A.length()-B.length())<1e-3) {
232			double mindist = std::numeric_limits<double>::max();
233			FOR(r,4) {
234			if ((A-B).length()<mindist) {
235			R = r;
236			mindist = (A-B).length();
237			}
238			A = vec2(-A[1], A[0]);
239			}
240			FOR(r, R) {
241			A = vec2(-A[1], A[0]);
242			C = vec2(-C[1], C[0]);
243			}
244			plop(mindist);
245			vec2 Ttmp = uv[fec.next(opp)] - C;
246			T = vec2(round(Ttmp[0]), round(Ttmp[1]));
247			trans_func_found = (mindist<1e-3 && (T-Ttmp).length()<1e-3);
248			} else {
249			plop("GNA?!");
250			}
251			if (!trans_func_found && NOT_AN_ID!=opp && !singular[fec.facet(opp)]) {
252			plop("ATTENTION, grids do not match");
253			singular[fec.facet(h)] = true;
254			singular[fec.facet(opp)] = true;
255			continue;
256			}
257			}
258			#endif
259		✗	vector<double> coeff;
260		✗	FOR(coord, 2) { // find barycentric coordinates of both u and v integer isos
261		✗	double v[2] = { lU[0][coord], lU[1][coord] };
262		✗	double from = floor(std::min(v[0], v[1])) + 1;
263		✗	double to = std::max(v[0], v[1]);
264		✗	if (to-from > 1000) continue;
265		✗	for (double iso = from; iso < to; iso += 1.) {
266		✗	double c = (iso - v[0]) / (v[1] - v[0]); // v[0] is far from v[1] (U was pre-snapped to integers with .05 tolerance)
267		✗	if (!Numeric::is_nan(c) && c > 0 && c < 1) {
268			// vec2 u = lU[0] + c*(lU[1] - lU[0]);
269			// u[coord] = iso;
270		✗	coeff.push_back(c);
271			}
272			}
273			}
274		✗	std::sort(coeff.begin(), coeff.end(), std::less<double>()); // we need to sort in order to merge close values
275
276		✗	vector<vec3> pts;
277		✗	vector<vec2> lu;
278		✗	vector<vec2> luopp;
279
280		✗	FOR(i, coeff.size()) {
281			// a + c*(b-a) returns exactly a when a==b and no NaNs involved
282			// no need to worry about the cases when c==0 and c==1, because of the presnapping of the parameterization
283		✗	vec3 pt = X(m)[fec.org(h)] + coeff[i] * (X(m)[fec.dest(h)] - X(m)[fec.org(h)]);
284		✗	vec2 u = lU[0] + coeff[i] * (lU[1] - lU[0]);
285		✗	vec2 uopp;
286		✗	if (NOT_AN_ID != opp) uopp = uv[fec.next(opp)] + coeff[i] * (uv[opp] - uv[fec.next(opp)]);
287
288		✗	FOR(d, 2) { // we must guarantee that new vertices have (at least) one integer component
289		✗	if (std::abs(u[d] - round(u[d])) < 1e-10) {
290		✗	u[d] = round(u[d]);
291			}
292
293		✗	if (NOT_AN_ID != opp && std::abs(uopp[d] - round(uopp[d])) < 1e-10) {
294		✗	uopp[d] = round(uopp[d]);
295			}
296			}
297
298		✗	pts.push_back(pt);
299		✗	lu.push_back(u);
300		✗	if (NOT_AN_ID != opp) luopp.push_back(uopp);
301			}
302
303			// create vertices
304		✗	index_t off = m->vertices.create_vertices(pts.size());
305		✗	FOR(i, pts.size()) {
306		✗	added_vertices[i+off] = true;
307		✗	resp_facet[i+off] = orig_tri_fid[fec.facet(h)];
308			}
309		✗	FOR(i, pts.size()) {
310		✗	m->vertices.point(off + i) = pts[i];
311		✗	new_corners[h].push_back(NewCorner(off + i, lu[i]));
312		✗	if (NOT_AN_ID != opp)
313		✗	new_corners[opp].push_back(NewCorner(off + i, luopp[i]));
314			}
315		✗	if (NOT_AN_ID != opp)
316		✗	std::reverse(new_corners[opp].begin(), new_corners[opp].end());
317		✗	}
318		✗	}
319
320			// now we have all the corners to insert, we create new facets for all the surface, old facets are to delete
321
322		✗	reach("split edges");
323		✗	FOR(f, nb_init_facets) {
324		✗	vector<index_t> polyV;
325		✗	vector<vec2> poly_uv;
326		✗	FOR(fc, m->facets.nb_corners(f)) {
327		✗	polyV.push_back(m->facets.vertex(f, fc));
328		✗	index_t c = m->facets.corner(f, fc);
329		✗	poly_uv.push_back(uv[c]);
330		✗	FOR(i, new_corners[c].size()) {
331		✗	polyV.push_back(new_corners[c][i].first);
332		✗	poly_uv.push_back(new_corners[c][i].second);
333			}
334			}
335		✗	index_t nf = m->facets.create_polygon(polyV);
336		✗	m->facets.attributes().copy_item(nf ,f);
337		✗	FOR(fc, m->facets.nb_corners(nf))
338		✗	uv[m->facets.corner(nf, fc)] = poly_uv[fc];
339		✗	}
340
341		✗	reach("kill facets");
342		✗	vector<index_t> to_kill(nb_init_facets, true); // kill old (pre-split) facets
343		✗	to_kill.resize(m->facets.nb(), false);
344		✗	m->facets.delete_elements(to_kill);
345		✗	}
346
347
348			/**
349			* INPUT: facets with uv coordinates, where many vertices have integer values in uv
350			* OUTPUT: facets with uv coordinates, inclusing edges that have one coordinate of uv that is constant and integer valued
351			*
352			* remark: some polylines of these new edges are likely to be pre-images of edges of the regular grid...
353			*/
354
355
356		✗	void find_degenerate_facets(Mesh* m, vector<index_t> &degenerate) {
357		✗	degenerate.clear();
358		✗	FOR(f, m->facets.nb()) {
359		✗	index_t nbv = m->facets.nb_corners(f);
360			// geo_assert(3 == m->facets.nb_corners(f));
361		✗	FOR (c1, nbv) {
362		✗	FOR (c2, nbv) {
363		✗	if (c1 == c2) continue;
364		✗	index_t v1 = m->facets.vertex(f, c1);
365		✗	index_t v2 = m->facets.vertex(f, c2);
366
367		✗	geo_assert( v1!= v2 );
368		✗	if (X(m)[v1][0] == X(m)[v2][0] && X(m)[v1][1] == X(m)[v2][1] && X(m)[v1][2] == X(m)[v2][2]) {
369		✗	degenerate.push_back(f);
370			}
371			}
372			}
373			}
374		✗	}
375
376		✗	void imprint(Mesh* m, const char *uv_name, const char *singular_name) {
377			{
378		✗	Attribute<index_t> singular(m->facets.attributes(), singular_name);
379		✗	Attribute<index_t> resp_facet(m->vertices.attributes(), "resp_facet");
380		✗	Attribute<index_t> orig_tri_fid(m->facets.attributes(), "orig_tri_fid");
381		✗	FOR(i, m->facets.nb() ) orig_tri_fid[i] = i;
382		✗	FOR(v, m->vertices.nb()) resp_facet[v] = NOT_AN_ID;
383
384		✗	check_no_intersecting_faces(m);
385		✗	}
386
387
388
389		✗	Mesh m_bak;
390		✗	m_bak.copy(*m);
391
392			for(;;) {
393		✗	split_edges_by_iso_uvs(m, uv_name, singular_name);
394		✗	facets_split(m, uv_name, singular_name);
395
396		✗	triangulate_surface_preserve_attributes(m);
397
398		✗	vector<index_t> fails;
399		✗	find_degenerate_facets(m, fails);
400
401		✗	vector<index_t> intersections;
402		✗	find_self_intersections(m, intersections);
403		✗	fails.insert(fails.end(), intersections.begin(), intersections.end()); // degenerate triangles or intersecting ones, all the same, I do not want them
404		✗	if (!fails.size()) break;
405
406		✗	plop("imprint fail, need to undo");
407			{
408		✗	Attribute<index_t> singular(m_bak.facets.attributes(), singular_name);
409		✗	Attribute<index_t> resp_facet(m->vertices.attributes(), "resp_facet");
410
411		✗	FOR(i, fails.size()) {
412		✗	FOR(j, 3) {
413		✗	index_t f = resp_facet[m->facets.vertex(fails[i], j)];
414		✗	if (NOT_AN_ID!=f) singular[f] = true;
415			}
416			}
417		✗	}
418		✗	m->copy(m_bak);
419		✗	}
420		✗	}
421
422		✗	void facets_split(Mesh* m, const char *uv_name, const char *singular_name) {
423		✗	Attribute<bool> added_vertices(m->vertices.attributes(), "added_vertices");
424		✗	Attribute<GEO::vec2> uv(m->facet_corners.attributes(), uv_name);
425		✗	Attribute<index_t> singular(m->facets.attributes(), singular_name);
426
427		✗	Attribute<index_t> orig_tri_fid(m->facets.attributes(), "orig_tri_fid");
428		✗	Attribute<index_t> resp_facet(m->vertices.attributes(), "resp_facet");
429
430		✗	vector<index_t> to_kill(m->facets.nb(), 0);
431		✗	FOR(f, m->facets.nb()) {
432		✗	if (singular[f]) continue;
433		✗	index_t nbc = m->facets.nb_corners(f);
434
435
436			// STEP 1: find a couple (coordinate, iso value) to split the facet
437		✗	index_t coord = index_t(-1);
438		✗	double iso = 0;
439		✗	FOR(test_coord, 2) {
440		✗	double min_v = 1e20;
441		✗	double max_v = -1e20;
442		✗	FOR(fc, nbc) {
443		✗	min_v = std::min(min_v, uv[m->facets.corner(f, fc)][test_coord]);
444		✗	max_v = std::max(max_v, uv[m->facets.corner(f, fc)][test_coord]);
445			}
446		✗	if (floor(min_v)+1 < max_v) { // floor(min_v)+1 is the first integer value strictly superior to min_v
447		✗	coord = test_coord;
448		✗	iso = floor(min_v)+1;
449		✗	break;
450			}
451			}
452		✗	if (coord == index_t(-1)) continue;
453
454			// STEP 2: if STEP 1 succedeed, compute the extremities of the new edge
455		✗	index_t cut[2] = { NOT_AN_ID, NOT_AN_ID };
456
457		✗	FOR(fc, nbc) {
458		✗	if (uv[m->facets.corner(f, fc)][coord] != iso) continue;
459		✗	if (cut[0] == NOT_AN_ID) {
460		✗	cut[0] = fc;
461		✗	} else if (fc != next_mod(cut[0], nbc) && next_mod(fc, nbc) != cut[0]) { // no biangles
462		✗	cut[1] = fc;
463			}
464			}
465		✗	if (cut[1] == NOT_AN_ID) continue;
466
467			{ // let us check that the cut separates the facet
468		✗	bool inf1=true, sup1=true, inf2=true, sup2=true;
469		✗	for (index_t fc = cut[0]+1; fc<cut[1]; fc++) {
470		✗	double tex = uv[m->facets.corner(f, fc)][coord];
471		✗	inf1 = inf1 && (tex<iso);
472		✗	sup1 = sup1 && (tex>iso);
473			}
474		✗	for (index_t fc = cut[1]+1; fc<cut[0]+nbc; fc++) {
475		✗	double tex = uv[m->facets.corner(f, fc%nbc)][coord];
476		✗	inf2 = inf2 && (tex<iso);
477		✗	sup2 = sup2 && (tex>iso);
478			}
479		✗	if (!( (inf1 && sup2) || (sup1 && inf2) )) {
480		✗	plop("WARNING: facet cannot be properly cut by the iso");
481		✗	continue;
482			}
483			}
484
485
486		✗	to_kill[f] = true;
487
488			// STEP 3: compute new vertices (iso-integer value of uv) on the new edge
489		✗	vector<vec3> nv_pts;
490		✗	vector<vec2> nv_uv;
491			{
492		✗	vec3 lX[2];
493		✗	FOR(i, 2) lX[i] = m->vertices.point(m->facets.vertex(f, cut[i]));
494		✗	vec2 lU[2];
495		✗	FOR(i, 2) lU[i] = uv[m->facets.corner(f, cut[i])];
496		✗	vector<double> coeff;
497		✗	double v[2] = { lU[0][(coord+1)%2], lU[1][(coord+1)%2] }; // recall that coord is the cutting dimension
498
499		✗	for (double cur_iso = ceil(std::min(v[0], v[1])); cur_iso < std::max(v[0], v[1]); cur_iso += 1.0) {
500		✗	double c = (cur_iso - v[0]) / (v[1] - v[0]);
501		✗	if (!Numeric::is_nan(c) && c > 0 && c < 1)
502		✗	coeff.push_back(c);
503			}
504
505		✗	std::sort(coeff.begin(), coeff.end(), std::less<double>());
506		✗	FOR(i, coeff.size()) {
507		✗	vec3 x = lX[0] + coeff[i] * (lX[1] - lX[0]); // it guarantees x==lX[0] when lX[0]==lX[1]
508		✗	vec2 u = lU[0] + coeff[i] * (lU[1] - lU[0]); // no need to worry about coeff[i]==0 and coeff[i]==1 because of the parameterization pre-snapping
509		✗	nv_pts.push_back(x);
510		✗	nv_uv.push_back(u);
511			}
512			// new vertices must have only integer values of uv --- remove possible numerical imprecision
513		✗	FOR(i, nv_pts.size()) {
514		✗	FOR(d, 2) {
515		✗	nv_uv[i][d] = round(nv_uv[i][d]);
516			}
517			}
518		✗	}
519
520			// STEP 4: create new vertices and new faces
521		✗	index_t off = m->vertices.create_vertices(nv_pts.size());
522		✗	FOR(i, nv_pts.size()) {
523		✗	resp_facet[off+i] = orig_tri_fid[f];
524		✗	added_vertices[off+i] = true;
525		✗	X(m)[off + i] = nv_pts[i];
526			}
527		✗	FOR(half, 2) {
528		✗	vector <index_t> lv;
529		✗	vector <vec2> luv;
530
531			// add original vertices
532		✗	index_t cir = cut[half];
533			do {
534		✗	lv.push_back(m->facets.vertex(f, cir));
535		✗	luv.push_back(uv[m->facets.corner(f, cir)]);
536		✗	cir = next_mod(cir, nbc);
537		✗	} while (cir != cut[(half + 1) % 2]);
538		✗	lv.push_back(m->facets.vertex(f, cir));
539		✗	luv.push_back(uv[m->facets.corner(f, cir)]);
540
541			// add new vertices
542		✗	FOR(i, nv_pts.size()) {
543		✗	index_t ind = i;
544		✗	if (half == 0) ind = nv_pts.size() - 1 - i;
545		✗	lv.push_back(off + ind);
546		✗	luv.push_back(nv_uv[ind]);
547			}
548
549		✗	index_t fid = m->facets.create_polygon(lv);
550		✗	m->facets.attributes().copy_item(fid, f);
551		✗	FOR(fc, m->facets.nb_corners(fid)) {
552		✗	uv[m->facets.corner(fid, fc)] = luv[fc];
553			}
554		✗	to_kill.push_back(false);
555		✗	}
556		✗	}
557		✗	m->facets.delete_elements(to_kill);
558		✗	}
559
560			/**
561			* INPUT: facets with uv coordinates
562			* OUTPUT: chart facet attribute s.t. the chart frontier is included in edges that are iso-integer value of 'uv'
563			*/
564
565		✗	inline index_t indir_root(index_t i, vector<index_t>& indir) {
566		✗	while (i != indir[i]) i = indir[i];
567		✗	return i;
568			}
569
570			// merges charts for adjacent facets under two conditions:
571			// 1) both facets are not singular
572			// 2) the shared edge is not integer iso in both facets
573		✗	void mark_charts(Mesh* m, const char *uv_name, const char *chart_name, const char *singular_name) { // 2-manifold surface is supposed
574		✗	Attribute<bool> isovalue(m->facet_corners.attributes(), "isovalue");
575		✗	Attribute<bool> quadelement(m->facets.attributes(), "quadelement");
576		✗	Attribute<bool> quadcorners(m->vertices.attributes(), "quadcorners");
577		✗	Attribute<bool> added_vertices(m->vertices.attributes(), "added_vertices");
578
579		✗	Attribute<GEO::vec2> uv(m->facet_corners.attributes(), uv_name);
580		✗	Attribute<index_t> singular(m->facets.attributes(), singular_name);
581		✗	Attribute<index_t> chart(m->facets.attributes(), chart_name);
582
583		✗	FacetsExtraConnectivity fec(m);
584
585		✗	vector<index_t> indir(m->facets.nb());
586			{ // fill isovalue edge attribute and merge quad candidate charts
587		✗	FOR(f, m->facets.nb()) {
588		✗	indir[f] = f;
589			}
590
591		✗	FOR(h, m->facet_corners.nb()) {
592		✗	index_t hopp = fec.opposite(h);
593		✗	index_t f = fec.facet(h);
594		✗	index_t fopp = NOT_AN_ID==hopp ? NOT_AN_ID : fec.facet(hopp);
595
596		✗	if (singular[f] || (NOT_AN_ID!=hopp && hopp>h)) continue;
597
598		✗	bool cut = false;
599		✗	index_t test[2] = { h, hopp };
600		✗	FOR(lh, 2) {
601		✗	if (lh && (NOT_AN_ID==hopp || singular[fopp])) break;
602		✗	FOR(coord, 2) {
603		✗	cut = cut || (uv[test[lh]][coord] == uv[fec.next(test[lh])][coord] && uv[test[lh]][coord] == round(uv[test[lh]][coord]));
604			}
605			}
606
607		✗	if (cut) {
608		✗	isovalue[h] = true;
609		✗	if (NOT_AN_ID!=hopp) isovalue[hopp] = true;
610			}
611
612		✗	if (!cut && NOT_AN_ID!=fopp && !singular[fopp]) {
613		✗	indir[indir_root(fopp, indir)] = indir_root(f, indir);
614			}
615			}
616
617		✗	FOR(f, m->facets.nb()) {
618		✗	chart[f] = indir_root(f, indir);
619			}
620			}
621
622
623			{ // determine quad corner vertices: fill quadcorners attribute
624		✗	FOR (h, m->facet_corners.nb()) {
625		✗	int cnt = 0;
626		✗	index_t cir = h;
627			do {
628		✗	if (NOT_AN_ID==cir) break;
629		✗	cnt += int(isovalue[cir]);
630		✗	cir = fec.next_around_vertex(cir);
631		✗	} while (cir != h);
632
633		✗	if ((NOT_AN_ID==cir && cnt>=2) || cnt>2) {
634		✗	quadcorners[fec.org(h)] = true;
635			}
636			}
637			}
638
639
640			{ // fill quadelement attribute
641
642			// this code verifies only if chart boundaries are marked as isovalues + if it has 4 quadcorners
643			// normally it is also necessary to check if the "quad" is a 2d disk (one boundary + Euler characterisic (imagine a torus with a hole))
644
645		✗	vector<bool> seen(m->facets.nb(), false);
646		✗	FOR(fseed, m->facets.nb()) {
647		✗	if (fseed != chart[fseed]) continue;
648		✗	int nb_quad_corners = 0;
649		✗	std::deque<index_t> Q;
650		✗	Q.push_back(fseed);
651		✗	seen[fseed] = true;
652		✗	bool iso_only_at_boundaries = true;
653		✗	std::vector<index_t> C;
654		✗	while (Q.size()) {
655		✗	index_t f = Q.front();
656		✗	if (singular[f]) {
657		✗	iso_only_at_boundaries = false;
658		✗	break;
659			}
660		✗	C.push_back(f);
661		✗	Q.pop_front();
662		✗	FOR(fc, m->facets.nb_corners(f)) {
663		✗	index_t h = m->facets.corner(f, fc);
664		✗	index_t hopp = fec.opposite(h);
665		✗	index_t fopp = NOT_AN_ID==hopp ? NOT_AN_ID : fec.facet(hopp);
666		✗	if (fopp==NOT_AN_ID) {
667		✗	if (isovalue[h] && quadcorners[fec.org(h)]) nb_quad_corners++;
668			} else {
669		✗	if (chart[fopp] != chart[fseed]) {
670		✗	if (isovalue[h]) {
671		✗	if (quadcorners[fec.org(h)])
672		✗	nb_quad_corners++;
673		✗	continue;
674			}
675		✗	Q.clear();
676		✗	iso_only_at_boundaries = false;
677		✗	break;
678			}
679			}
680		✗	if (NOT_AN_ID!=fopp && !seen[fopp]) {
681		✗	Q.push_back(fopp);
682		✗	seen[fopp] = true;
683			}
684			}
685			}
686		✗	if (iso_only_at_boundaries && nb_quad_corners == 4) {
687		✗	FOR(i, C.size()) {
688		✗	quadelement[C[i]] = true;
689			}
690			}
691		✗	}
692		✗	}
693
694			{ // charts = orig triangles in non-quad regions
695		✗	FOR(f, m->facets.nb()) {
696		✗	if (!quadelement[f]) indir[f] = f;
697			}
698
699		✗	Attribute<index_t> orig_tri_fid(m->facets.attributes(), "orig_tri_fid");
700		✗	FOR(h, m->facet_corners.nb()) {
701		✗	index_t hopp = fec.opposite(h);
702		✗	index_t f = fec.facet(h);
703		✗	index_t fopp = NOT_AN_ID==hopp ? NOT_AN_ID : fec.facet(hopp);
704		✗	if (NOT_AN_ID == fopp || quadelement[f] || quadelement[fopp] || orig_tri_fid[f]!=orig_tri_fid[fopp]) continue;
705		✗	indir[indir_root(fopp, indir)] = indir_root(f, indir);
706			}
707
708		✗	FOR(f, m->facets.nb()) {
709		✗	chart[f] = indir_root(f, indir);
710			}
711		✗	}
712
713		✗	vector<vector<index_t> > v2f = generate_v2f(m);
714			{ // mark the boundary between triangles and quads
715		✗	FOR(v, m->vertices.nb()) {
716		✗	TrFan fan = TrFan(v, m, v2f, chart);
717		✗	bool touches_a_quad = false;
718		✗	FOR(i, fan.nb_fan_triangles()) {
719		✗	touches_a_quad = touches_a_quad || quadelement[fan[i].f];
720			}
721		✗	if (quadcorners[v] || !touches_a_quad || fan.ncharts()<=2) continue;
722
723		✗	index_t first_triangle = NOT_AN_ID;
724		✗	FOR(i, fan.nb_fan_triangles()) {
725		✗	if (quadelement[fan[i].f]) continue;
726		✗	if (NOT_AN_ID==first_triangle) {
727		✗	first_triangle = fan[i].f;
728			}
729		✗	indir[indir_root(first_triangle, indir)] = indir_root(fan[i].f, indir);
730			}
731		✗	}
732		✗	FOR(f, m->facets.nb()) {
733		✗	chart[f] = indir_root(f, indir);
734			}
735			}
736
737			{ // fill verts to remove attribute
738		✗	Attribute<bool> verts_to_remove(m->vertices.attributes(), "verts_to_remove");
739		✗	FOR(v, m->vertices.nb()) {
740		✗	TrFan fan = TrFan(v, m, v2f, chart);
741		✗	bool touches_a_quad = false;
742		✗	FOR(i, fan.nb_fan_triangles()) {
743		✗	touches_a_quad = touches_a_quad || quadelement[fan[i].f];
744			}
745		✗	if (quadcorners[v]) continue;
746		✗	if ( fan.incomplete_ && fan.ncharts()!=1) continue;
747		✗	if (!fan.incomplete_ && fan.ncharts() >2) continue;
748		✗	if (!touches_a_quad && !added_vertices[v]) continue;
749		✗	verts_to_remove[v] = true;
750		✗	}
751		✗	}
752		✗	}
753
754			/****************************************************************************************************/
755		✗	static bool try_export_quadtri_from_charts(Mesh* m, vector<BBox>& locked_regions) {
756		✗	bool modified = false;
757		✗	Attribute<index_t> chart(m->facets.attributes(), "chart"); // TODO document this
758		✗	Attribute<bool> quadelement(m->facets.attributes(), "quadelement");
759
760		✗	vector<index_t> to_kill(m->facets.nb(), false);
761
762		✗	index_t max_chart_no = 0;
763		✗	FOR(f, m->facets.nb()) {
764		✗	max_chart_no = std::max(max_chart_no, chart[f] + 1);
765			}
766		✗	vector<int> nbtri_in_chart(max_chart_no, 0);
767		✗	FOR(f, m->facets.nb()) {
768		✗	nbtri_in_chart[chart[f]]++;
769			}
770
771		✗	FacetsExtraConnectivity fec(m);
772		✗	index_t nbf = m->facets.nb();
773		✗	FOR(f, nbf) {
774		✗	geo_assert(3 == m->facets.nb_corners(f));
775
776		✗	if (nbtri_in_chart[chart[f]] != 2 || !quadelement[f]) continue;
777		✗	FOR(ih, 3) {
778		✗	index_t h = m->facets.corner(f, ih);
779		✗	index_t hopp = fec.opposite(h);
780		✗	if (NOT_AN_ID==hopp) continue;
781		✗	index_t fopp = fec.facet(hopp);
782		✗	geo_assert(NOT_AN_ID != fopp);
783		✗	if (chart[fopp] != chart[f]) continue;
784		✗	if (f < fopp) break;
785
786		✗	vector<index_t> pts;
787		✗	pts.push_back(fec.dest(h));
788		✗	pts.push_back(fec.dest(fec.next(h)));
789		✗	pts.push_back(fec.org(h));
790		✗	pts.push_back(fec.dest(fec.next(hopp)));
791
792		✗	bool intersect_locked_region = false;
793		✗	BBox b;
794		✗	FOR(v, 4) {
795		✗	b.add(X(m)[pts[v]]);
796			}
797		✗	FOR(i, locked_regions.size()) {
798		✗	intersect_locked_region = intersect_locked_region || locked_regions[i].intersect(b);
799			}
800
801		✗	if (!intersect_locked_region) {
802		✗	index_t nf = m->facets.create_polygon(pts);
803		✗	modified = true;
804		✗	m->facets.attributes().copy_item(nf ,f);
805		✗	to_kill.push_back(false);
806		✗	to_kill[f] = true;
807		✗	to_kill[fopp] = true;
808			}
809		✗	}
810			}
811		✗	m->facets.delete_elements(to_kill, false);
812		✗	return modified;
813		✗	}
814
815
816
817
818		✗	static void sample_triangle(vec3 *ABC, double eps, vector<vec3>& samples) {
819		✗	double max_edge_length = 0;
820		✗	FOR(p, 3) max_equal(max_edge_length, (ABC[(p + 1) % 3] - ABC[p]).length());
821		✗	index_t nb_steps = 10;
822		✗	if (eps>0) min_equal(nb_steps, index_t(max_edge_length / eps + 2));
823
824		✗	FOR(i, nb_steps)FOR(j, nb_steps - i) {
825		✗	double u = double(i) / double(nb_steps - 1);
826		✗	double v = double(j) / double(nb_steps - 1);
827		✗	samples.push_back(ABC[0] + u*(ABC[1] - ABC[0]) + v*(ABC[2] - ABC[0]));
828			}
829		✗	}
830			//double upper_bound_min_dist2_to_triangles(vec3 P, vector<vec3>& triangles) {
831			// vec3 closest_point;
832			// double l0, l1, l2;
833			// double min_dist2 = 1e20;
834			// FOR(t, triangles.size() / 3)
835			// min_equal(min_dist2,
836			// Geom::point_triangle_squared_distance<vec3>(P,
837			// triangles[t * 3], triangles[t * 3 + 1], triangles[t * 3 + 2], closest_point, l0, l1, l2)
838			// );
839			// return min_dist2;
840			//}
841
842		✗	static vector<index_t> facets_having_a_point_further_than_eps(Mesh* m,Mesh* ref,double epsilon) {
843
844		✗	vector<index_t> res;
845
846		✗	vector<BBox> inboxes = facets_bbox(ref);
847		✗	DynamicHBoxes hb; hb.init(inboxes);
848
849		✗	FOR(f, m->facets.nb()) {
850		✗	bool fail = false;
851		✗	vector<vec3> samples;
852		✗	vec3 ABC[3];
853		✗	FOR(lv,3) ABC[lv] = X(m)[m->facets.vertex(f,lv)];
854		✗	sample_triangle(ABC, epsilon / 2., samples);
855		✗	if (m->facets.nb_vertices(f) == 4) {
856		✗	FOR(lv, 3) ABC[lv] = X(m)[m->facets.vertex(f, (lv+2)%3)];
857		✗	sample_triangle(ABC, epsilon / 2., samples);
858			}
859
860
861		✗	FOR(p, samples.size()) {
862		✗	vec3 P = samples[p];
863
864		✗	double min_dist2 = 1e20;
865
866		✗	BBox bbox; bbox.add(P); bbox.dilate(epsilon/2.);
867		✗	vector<index_t> prim;
868		✗	hb.intersect(bbox, prim);
869		✗	FOR(fid, prim.size()) {
870		✗	index_t other_f = prim[fid];
871		✗	vec3 closest_point;
872			double l0, l1, l2;
873		✗	min_equal(min_dist2, Geom::point_triangle_squared_distance<vec3>(P,
874		✗	X(ref)[ref->facets.vertex(other_f, 0)],
875		✗	X(ref)[ref->facets.vertex(other_f, 1)],
876		✗	X(ref)[ref->facets.vertex(other_f, 2)],
877			closest_point, l0, l1, l2));
878
879		✗	if (ref->facets.nb_vertices(other_f) == 4) {
880		✗	min_equal(min_dist2, Geom::point_triangle_squared_distance<vec3>(P,
881		✗	X(ref)[ref->facets.vertex(other_f, 0)],
882		✗	X(ref)[ref->facets.vertex(other_f, 2)],
883		✗	X(ref)[ref->facets.vertex(other_f, 3)],
884			closest_point, l0, l1, l2));
885			}
886			}
887
888		✗	if (::sqrt(min_dist2) > epsilon / 2.) {
889		✗	fail = true;
890		✗	break;
891			}
892		✗	}
893		✗	if (fail) res.push_back(f);
894		✗	}
895		✗	return res;
896		✗	}
897
898			// Attention, sub-functions of this function need to access to attributes "chart" and "singular"
899		✗	void simplify_quad_charts(Mesh* m) {
900		✗	std::string msg;
901		✗	if (!surface_is_manifold(m, msg)) plop(msg);
902
903		✗	Mesh m_bak;
904		✗	m_bak.copy(*m);
905		✗	double epsilon = 0;// .4*get_facet_average_edge_size(&m_bak);
906			//epsilon = 0;
907		✗	vector<BBox> locked_regions;
908			// int cnt = 0;
909			for(;;) {
910		✗	vector<index_t> invalid_m ;
911		✗	vector<index_t> invalid_bak ;
912		✗	vector<index_t> intersections;
913			{
914		✗	Attribute<index_t> chart(m->facets.attributes(), "chart");
915		✗	Attribute<index_t> undo(m->facets.attributes(), "undo");
916		✗	FOR(fid, m->facets.nb()) {
917		✗	undo[fid] = NOT_AN_ID;
918			}
919		✗	Attribute<bool> verts_to_remove(m->vertices.attributes(), "verts_to_remove");
920		✗	plop("try_simplify(m, chart, verts_to_remove, undo)");
921		✗	try_simplify(m, chart, verts_to_remove, undo);
922		✗	plop("try_export_quadtri_from_charts(m, locked_regions)");
923		✗	try_export_quadtri_from_charts(m, locked_regions);
924
925		✗	plop("check for intersections");
926		✗	plop(m->facets.nb());
927		✗	find_self_intersections(m, intersections);
928		✗	FOR(f, intersections.size()) {
929		✗	if (3 == m->facets.nb_vertices(f)) {
930		✗	if (undo[intersections[f]] != NOT_AN_ID) verts_to_remove[undo[intersections[f]]] = false;
931			// GEO::Logger::out("HexDom") << intersections[f] << " " << undo[intersections[f]] << std::endl;
932			} else {
933		✗	geo_assert(4 == m->facets.nb_vertices(f));
934		✗	BBox inbox;
935		✗	FOR(fv, 4) {
936		✗	inbox.add(X(m)[m->facets.vertex(intersections[f], fv)]);
937			}
938		✗	locked_regions.push_back(inbox);
939			}
940			// Attribute<bool> gna(m->facets.attributes(), "gna");
941			// gna[intersections[f]] = true;
942			}
943			//intersections.clear();// HACK (simulation de quadhex)
944
945		✗	if (epsilon > 0) {
946		✗	plop("check for Hausdorff distance --- dist to init mesh");
947		✗	invalid_m = facets_having_a_point_further_than_eps(m, &m_bak, epsilon);
948		✗	plop(invalid_m.size());
949		✗	FOR(i, invalid_m.size()) {
950		✗	if (3 == m->facets.nb_vertices(invalid_m[i])) {
951		✗	if (undo[invalid_m[i]] != NOT_AN_ID) verts_to_remove[undo[invalid_m[i]]] = false;
952			}
953			else {
954		✗	BBox inbox;
955		✗	FOR(fv, 4) {
956		✗	inbox.add(X(m)[m->facets.vertex(invalid_m[i], fv)]);
957			}
958		✗	locked_regions.push_back(inbox);
959			}
960			}
961
962		✗	plop("check for Hausdorff distance --- dist to new mesh");
963		✗	invalid_bak = facets_having_a_point_further_than_eps(&m_bak, m, epsilon);
964		✗	plop(invalid_bak.size());
965
966		✗	FOR(i, invalid_bak.size()) FOR(lv, 3)
967		✗	verts_to_remove[m_bak.facets.vertex(invalid_bak[i], lv)] = false;
968			}
969		✗	}
970
971
972		✗	if (intersections.empty() && invalid_m.empty() && invalid_bak.empty()) break;
973
974		✗	plop("conflict detected");
975			{
976		✗	Attribute<bool> m_verts_to_remove(m->vertices.attributes(), "verts_to_remove");
977		✗	Attribute<bool> m_bak_verts_to_remove(m_bak.vertices.attributes(), "verts_to_remove");
978		✗	geo_assert(m->vertices.nb() == m_bak.vertices.nb());
979		✗	FOR(v, m->vertices.nb()) m_bak_verts_to_remove[v] = m_verts_to_remove[v];
980		✗	}
981
982			// char filename[1024], filename2[1024];
983			// sprintf(filename, "/home/ssloy/tmp/hexdom_nightly/geogram/zdebug%i_bak.geogram", cnt);
984			// sprintf(filename2, "/home/ssloy/tmp/hexdom_nightly/geogram/zdebug%i_simp.geogram", cnt);
985			// cnt++;
986			// mesh_save(m_bak, filename);
987			// mesh_save(*m, filename2);
988
989		✗	m->copy(m_bak);
990		✗	}
991		✗	kill_isolated_vertices(m);
992		✗	}
993
994
995
996
997			}
998