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 <geogram/delaunay/LFS.h>
41			#include <geogram/basic/geometry.h>
42			#include <geogram/basic/logger.h>
43			#include <geogram/basic/command_line.h>
44
45			namespace {
46
47			using namespace GEO;
48
49			const index_t pid_size = 6;
50			int pid[pid_size][4] = {
51			{0, 1, 2, 3},
52			{0, 2, 1, 3},
53			{0, 3, 1, 2},
54			{1, 2, 0, 3},
55			{1, 3, 2, 0},
56			{2, 3, 0, 1}
57			};
58
59			/**
60			* \brief Computes the circumcenter and squared radius of a tetrahedron.
61			* \param[in] p first vertex of the tetrahedron
62			* \param[in] q second vertex of the tetrahedron
63			* \param[in] r third vertex of the tetrahedron
64			* \param[in] s fourth vertex of the tetrahedron
65			* \param[out] t_circumcenter computed circumcenter of the tetrahedron
66			* \param[out] squared_radius computed squared radius of the tetrahedron
67			* \param[out] dihedral_angles if non-nullptr,
68			* computed dihedral angles of the tetrahedron
69			* \return true if successful, false otherwise
70			* (for instance, if the tetrahedron is flat).
71			*/
72		✗	bool tetra_circumcenter_squaredradius(
73			const vec3& p, const vec3& q,
74			const vec3& r, const vec3& s,
75			vec3& t_circumcenter, double& squared_radius,
76			double* dihedral_angles = nullptr
77			) {
78
79		✗	vec3 qp = q - p;
80		✗	vec3 rp = r - p;
81		✗	vec3 sp = s - p;
82
83		✗	double tet_vol_x_12 = 2.0 * dot(qp, cross(rp, sp));
84
85			// it is not a safe check, one should avoid
86			// to input a degenerated tetra.
87		✗	if(tet_vol_x_12 == 0.0) {
88		✗	return false;
89			}
90
91			vec3 d_t_circumcenter =
92		✗	length2(qp) * cross(rp, sp) +
93		✗	length2(rp) * cross(sp, qp) +
94		✗	length2(sp) * cross(qp, rp);
95
96		✗	squared_radius = length2(
97		✗	d_t_circumcenter) / (tet_vol_x_12 * tet_vol_x_12
98			);
99		✗	d_t_circumcenter = (1.0 / tet_vol_x_12) * d_t_circumcenter;
100		✗	t_circumcenter = p + d_t_circumcenter;
101
102		✗	if(dihedral_angles) {
103			// compute Dihedral angle directly
104			// the following computation is not optimized.
105		✗	vec3 point[4];
106		✗	point[0] = p;
107		✗	point[1] = q;
108		✗	point[2] = r;
109		✗	point[3] = s;
110		✗	for(index_t i = 0; i < pid_size; i++) {
111		✗	vec3 b1 = point[pid[i][0]] - point[pid[i][2]];
112		✗	vec3 b2 = point[pid[i][1]] - point[pid[i][0]];
113		✗	vec3 b3 = point[pid[i][3]] - point[pid[i][1]];
114		✗	vec3 b2b3 = cross(b2, b3);
115		✗	dihedral_angles[i] = ::fabs(
116		✗	::atan2(
117		✗	length(b2) * dot(b1, b2b3),
118		✗	dot(cross(b1, b2), b2b3)
119			)
120			);
121			}
122			}
123		✗	return true;
124			}
125
126			/**
127			* \brief Gets a Delaunay vertex by global index.
128			* \param[in] delaunay the Delaunay triangulation
129			* \param[in] v the index of the vertex
130			* \return a const reference to the vertex, as a vec3
131			*/
132		✗	inline const vec3& delaunay_vertex(Delaunay* delaunay, index_t v) {
133		✗	return *(const vec3*) delaunay->vertex_ptr(v);
134			}
135
136			/**
137			* \brief Gets a Delaunay vertex by tetrahedron index
138			* and local vertex index.
139			* \param[in] delaunay the Delaunay triangulation
140			* \param[in] c the index of the tetrahedron
141			* \param[in] lv the local index of the vertex (0,1,2 or 3)
142			* in tetrahedron \p c.
143			* \return a const reference to the vertex, as a vec3.
144			*/
145		✗	inline const vec3& delaunay_tet_vertex(
146			Delaunay* delaunay, index_t c, index_t lv
147			) {
148		✗	return delaunay_vertex(delaunay, index_t(delaunay->cell_vertex(c, lv)));
149			}
150
151			index_t tet_facet_vertex[4][3] = {
152			{1, 2, 3},
153			{0, 3, 2},
154			{3, 0, 1},
155			{2, 1, 0}
156			};
157
158			/**
159			* \brief Computes the normal to a tetrahedron facet.
160			* \param[in] delaunay the Delaunay triangulation
161			* \param[in] t the index of the tetrahedron
162			* \param[in] f the local index (0,1,2 or 3) of
163			* the facet in the tetrahedron \p t
164			* \return the normal to the facet \p f or tetrahedron \p t
165			*/
166		✗	vec3 delaunay_facet_normal(
167			Delaunay* delaunay, index_t t, index_t f
168			) {
169		✗	geo_debug_assert(f < 4);
170		✗	index_t v1 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][0]));
171		✗	index_t v2 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][1]));
172		✗	index_t v3 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][2]));
173		✗	const vec3& p1 = delaunay_vertex(delaunay, v1);
174		✗	const vec3& p2 = delaunay_vertex(delaunay, v2);
175		✗	const vec3& p3 = delaunay_vertex(delaunay, v3);
176		✗	return cross(p2 - p1, p3 - p1);
177			}
178			}
179
180			/****************************************************************************/
181
182			namespace GEO {
183
184		✗	void LocalFeatureSize::init(index_t nb_pts, const double* pts) {
185			// Note: I need PDEL here instead of ANN/BNN since I need
186			// to access the cells (and ANN/BNN only give me the
187			// neighbors)
188
189		✗	Logger::out("LFS") << "Delaunay" << std::endl;
190		✗	Delaunay_var delaunay = Delaunay::create(3, "PDEL");
191		✗	delaunay->set_vertices(nb_pts, pts);
192		✗	Logger::out("LFS") << "Done Delaunay" << std::endl;
193
194		✗	vector<vec3> circumcenter(delaunay->nb_cells());
195			vector<bool> voronoi_cell_is_infinite(
196			delaunay->nb_vertices(), false
197		✗	);
198		✗	vector<bool> is_sliver(delaunay->nb_cells());
199
200			static constexpr index_t NO_POLE = NO_INDEX;
201			// Index of incident tet whose circumcenter is the pole
202		✗	vector<index_t> positive_pole(delaunay->nb_vertices(),NO_POLE);
203		✗	vector<index_t> negative_pole(delaunay->nb_vertices(),NO_POLE);
204
205			vector<vec3> avg_infinite_dir(
206		✗	delaunay->nb_vertices(), vec3(0, 0, 0)
207		✗	);
208		✗	vector<double> dist(delaunay->nb_vertices(), 0.0);
209
210		✗	Logger::out("LFS") << "(1) Circumcenters and slivers" << std::endl;
211			// Step 1: compute circumcenters and check for slivers
212		✗	const double sliver_quality = sliver_angle_threshold_ / 180.0 * M_PI;
213		✗	for(index_t t = 0; t < delaunay->nb_cells(); t++) {
214		✗	const vec3& p = delaunay_tet_vertex(delaunay, t, 0);
215		✗	const vec3& q = delaunay_tet_vertex(delaunay, t, 1);
216		✗	const vec3& r = delaunay_tet_vertex(delaunay, t, 2);
217		✗	const vec3& s = delaunay_tet_vertex(delaunay, t, 3);
218
219			double dihedral_angle[6];
220			double tet_radius;
221		✗	bool ok = tetra_circumcenter_squaredradius(
222			p, q, r, s,
223		✗	circumcenter[t], tet_radius, dihedral_angle
224			);
225
226		✗	for(index_t a = 0; a < 6; a++) {
227		✗	if(!ok) {
228		✗	break;
229			}
230		✗	ok = ok &&
231		✗	dihedral_angle[a] >= sliver_quality &&
232		✗	dihedral_angle[a] <= M_PI - sliver_quality;
233			}
234		✗	is_sliver[t] = !ok;
235			}
236
237		✗	Logger::out("LFS") << "(2) Positive poles" << std::endl;
238			// Step 2: compute positive poles
239		✗	for(index_t t = 0; t < delaunay->nb_cells(); t++) {
240		✗	if(is_sliver[t]) {
241		✗	continue;
242			}
243		✗	index_t f_inf = NO_INDEX;
244		✗	for(index_t f = 0; f < 4; f++) {
245		✗	if(delaunay->cell_adjacent(t, f) == NO_INDEX) {
246		✗	f_inf = index_t(f);
247		✗	break;
248			}
249			}
250		✗	if(f_inf == NO_INDEX) {
251			// tet t does not have facet on border
252		✗	const vec3& c = circumcenter[t];
253		✗	for(index_t lv = 0; lv < 4; lv++) {
254		✗	index_t iv = index_t(delaunay->cell_vertex(t, lv));
255		✗	const vec3& v = delaunay_vertex(delaunay, iv);
256		✗	double d = length2(c - v);
257		✗	if(d > dist[iv]) {
258		✗	dist[iv] = d;
259		✗	positive_pole[iv] = t;
260			}
261			}
262			} else {
263			// tet t has facet f_inf on border
264			// positive pole = average direction of infinite Voronoi edges
265		✗	index_t v = index_t(delaunay->cell_vertex(t, index_t(f_inf)));
266		✗	voronoi_cell_is_infinite[v] = true;
267		✗	avg_infinite_dir[v] +=
268		✗	delaunay_facet_normal(delaunay, t, index_t(f_inf));
269			}
270			}
271
272		✗	Logger::out("LFS") << "(3) Negative poles" << std::endl;
273			// Step 3: compute negative poles
274		✗	std::fill(dist.begin(), dist.end(), 0.0);
275		✗	for(index_t t = 0; t < delaunay->nb_cells(); t++) {
276		✗	if(is_sliver[t]) {
277		✗	continue;
278			}
279		✗	bool t_is_infinite = false;
280		✗	for(index_t f = 0; f < 4; f++) {
281		✗	if(delaunay->cell_adjacent(t, f) == NO_INDEX) {
282		✗	t_is_infinite = true;
283		✗	break;
284			}
285			}
286		✗	if(t_is_infinite) {
287		✗	continue;
288			}
289		✗	const vec3& c = circumcenter[t];
290		✗	for(index_t lv = 0; lv < 4; lv++) {
291		✗	index_t iv = index_t(delaunay->cell_vertex(t, lv));
292		✗	const vec3& v = delaunay_vertex(delaunay, iv);
293		✗	vec3 N;
294		✗	if(voronoi_cell_is_infinite[iv]) {
295		✗	N = -avg_infinite_dir[iv];
296			} else {
297		✗	N = v - circumcenter[positive_pole[iv]];
298			}
299		✗	double d = dot(c - v, N);
300		✗	if(d > dist[iv]) {
301		✗	dist[iv] = d;
302		✗	negative_pole[iv] = t;
303			}
304			}
305			}
306
307		✗	Logger::out("LFS") << "(4) Kd-tree" << std::endl;
308			// Step 4: create search structure
309		✗	vector<bool> is_pole(delaunay->nb_cells(), false);
310		✗	for(index_t iv = 0; iv < delaunay->nb_vertices(); iv++) {
311		✗	if(negative_pole[iv] != NO_POLE) {
312		✗	is_pole[index_t(negative_pole[iv])] = true;
313			}
314		✗	if(!voronoi_cell_is_infinite[iv]) {
315		✗	if(positive_pole[iv] != NO_POLE) {
316		✗	is_pole[index_t(positive_pole[iv])] = true;
317			}
318			}
319			}
320
321		✗	index_t nb_poles = 0;
322		✗	for(index_t t = 0; t < delaunay->nb_cells(); t++) {
323		✗	if(is_pole[t]) {
324		✗	nb_poles++;
325			}
326			}
327
328		✗	poles_.reserve(nb_poles * 3);
329		✗	for(index_t t = 0; t < delaunay->nb_cells(); t++) {
330		✗	if(is_pole[t]) {
331		✗	poles_.push_back(circumcenter[t].x);
332		✗	poles_.push_back(circumcenter[t].y);
333		✗	poles_.push_back(circumcenter[t].z);
334			}
335			}
336
337		✗	spatial_search_ = Delaunay::create(3, "NN");
338		✗	spatial_search_->set_vertices(poles_.size() / 3, poles_.data());
339
340		✗	Logger::out("LFS") << "Done init." << std::endl;
341		✗	}
342			}
343