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/voronoi/convex_cell.h>
41
42			#ifndef STANDALONE_CONVEX_CELL
43			#include <geogram/numerics/predicates.h>
44			#include <geogram/mesh/mesh.h>
45			#endif
46
47			#include <iostream>
48			#include <fstream>
49			#include <vector>
50			#include <cmath>
51			#include <limits>
52			#include <stack>
53
54
55			namespace {
56			using namespace VBW;
57
58			/**
59			* \brief a class for a stack of ushorts allocated on the stack.
60			* \details Used by clip_by_plane_fast() internally. I do not want
61			* clip_by_plane_fast() to do dynamic allocation.
62			*/
63			class SmallStack_ushort {
64			public:
65
66			/**
67			* \brief SmallStack constructor.
68			* \param[in] buffer a buffer. Callers keeps ownership
69			* (in most cases it will be created by alloca()).
70			* \param[in] capacity the number of uints that can
71			* be stored in the buffer.
72			*/
73			SmallStack_ushort(
74			ushort* buffer, int capacity
75			) : buffer_(buffer),
76			index_(-1),
77			capacity_(capacity) {
78			}
79
80			/**
81			* \brief Tests whether this stack is empty.
82			* \retval true if this stack is empty.
83			* \retval false otherwise.
84			*/
85			bool empty() const {
86			return (index_ == -1);
87			}
88
89			/**
90			* \brief Pushes an item on the stack.
91			* \param[in] val the item to be pushed.
92			*/
93			void push(ushort val) {
94		✗	++index_;
95			vbw_assert(index_ < capacity_);
96			(void)capacity_; // To silence a warning.
97		✗	buffer_[index_] = val;
98		✗	}
99
100			/**
101			* \brief Pops an item from the stack.
102			*/
103			void pop() {
104			vbw_assert(!empty());
105		✗	--index_;
106			}
107
108			/**
109			* \brief Gets the item on the top of the stack.
110			* \return the item.
111			*/
112			ushort top() const {
113			vbw_assert(!empty());
114		✗	return buffer_[index_];
115			}
116			private:
117			ushort* buffer_;
118			int index_;
119			int capacity_;
120			};
121			}
122
123			/*************************************************************/
124
125			namespace VBW {
126
127		✗	ConvexCell::ConvexCell(ConvexCellFlags flags) :
128		✗	max_t_(64),
129		✗	max_v_(32),
130			t_(max_t_),
131		✗	t_adj_(max_t_),
132		✗	plane_eqn_(max_v_),
133		✗	v2t_(max_v_),
134		✗	v2e_(max_v_)
135			{
136			#ifndef STANDALONE_CONVEX_CELL
137		✗	use_exact_predicates_ = true;
138			#endif
139		✗	nb_t_ = 0;
140		✗	nb_v_ = 0;
141		✗	first_free_ = END_OF_LIST;
142		✗	first_valid_ = END_OF_LIST;
143		✗	geometry_dirty_ = true;
144		✗	has_vglobal_ = ((flags & WithVGlobal) != 0);
145		✗	if(has_vglobal_) {
146		✗	vglobal_.assign(max_v_,index_t(-1));
147			}
148		✗	has_tflags_ = ((flags & WithTFlags) != 0);
149		✗	if(has_tflags_) {
150		✗	tflags_.assign(max_t_,0);
151			}
152		✗	v2t_.assign(max_v_,ushort(-1));
153		✗	v2e_.assign(max_v_,uchar(-1));
154		✗	}
155
156			/***********************************************************************/
157
158		✗	void ConvexCell::clear() {
159		✗	nb_t_ = 0;
160		✗	nb_v_ = 0;
161		✗	first_free_ = END_OF_LIST;
162		✗	first_valid_ = END_OF_LIST;
163		✗	geometry_dirty_ = true;
164			#ifdef VBW_DEBUG
165			// Initialize all triangle flags with something
166			// different from VALID_TRIANGLE.
167			for(index_t t=0; t<max_t(); ++t) {
168			set_triangle_flags(t, END_OF_LIST);
169			}
170			#endif
171		✗	}
172
173			/***********************************************************************/
174
175		✗	void ConvexCell::init_with_box(
176			double xmin, double ymin, double zmin,
177			double xmax, double ymax, double zmax
178			) {
179		✗	clear();
180
181			// The vertex at infinity.
182		✗	plane_eqn_[0] = make_vec4(0,0,0,0);
183
184			// Offset for the 6 bounding box plane equations.
185			// Here they come first (offset is zero).
186			index_t boff = 1;
187
188			// The equations of the six faces of the bounding box.
189		✗	plane_eqn_[boff ] = make_vec4( 1.0, 0.0, 0.0, -xmin);
190		✗	plane_eqn_[boff+1] = make_vec4(-1.0, 0.0, 0.0, xmax);
191		✗	plane_eqn_[boff+2] = make_vec4( 0.0, 1.0, 0.0, -ymin);
192		✗	plane_eqn_[boff+3] = make_vec4( 0.0,-1.0, 0.0, ymax);
193		✗	plane_eqn_[boff+4] = make_vec4( 0.0, 0.0, 1.0, -zmin);
194		✗	plane_eqn_[boff+5] = make_vec4( 0.0, 0.0,-1.0, zmax);
195
196			// Create the 8 triangles that correspond to the
197			// 8 vertices of the bounding box.
198			// (Unused) adjacency info. ----------------.
199			// Triangle vertices -. |
200			// v v
201			new_triangle( boff+2,boff+5,boff+0, 1,4,2);
202			new_triangle( boff+5,boff+3,boff+0, 5,0,3);
203			new_triangle( boff+1,boff+5,boff+2, 0,6,3);
204			new_triangle( boff+5,boff+1,boff+3, 7,1,2);
205			new_triangle( boff+4,boff+2,boff+0, 0,5,6);
206			new_triangle( boff+4,boff+0,boff+3, 1,7,4);
207			new_triangle( boff+2,boff+4,boff+1, 7,2,4);
208			new_triangle( boff+4,boff+3,boff+1, 3,6,5);
209
210			// We already created 6 vertices (for the 6 bounding box
211			// plane equations) plus the vertex at infinity.
212		✗	nb_v_ = 7;
213
214		✗	geometry_dirty_ = true;
215		✗	}
216
217
218		✗	void ConvexCell::init_with_tet(
219			vec4 P0, vec4 P1, vec4 P2, vec4 P3
220			) {
221		✗	clear();
222
223			// The vertex at infinity.
224		✗	plane_eqn_[0] = make_vec4(0,0,0,0);
225
226			// Offset for the 4 plane equations.
227			// Plane 0 is vertex at infinity.
228			index_t boff = 1;
229
230		✗	plane_eqn_[boff ] = P0;
231		✗	plane_eqn_[boff+1] = P1;
232		✗	plane_eqn_[boff+2] = P2;
233		✗	plane_eqn_[boff+3] = P3;
234
235			// Create the 4 triangles (that correspond to
236			// the 4 vertices of the tetrahedron)
237			// (Unused) adjacency info. ----------.
238			// Triangle vertices -. |
239			// v v
240			new_triangle(boff+3, boff+2, boff+1, 3, 2, 1);
241			new_triangle(boff+3, boff+0, boff+2, 3, 0, 2);
242			new_triangle(boff+3, boff+1, boff+0, 3, 1, 0);
243			new_triangle(boff+2, boff+0, boff+1, 2, 0, 1);
244
245			// We already created 4 vertices (for the 4 facets
246			// plane equations) plus the vertex at infinity.
247		✗	nb_v_ = 5;
248
249		✗	geometry_dirty_ = true;
250		✗	}
251
252
253		✗	void ConvexCell::init_with_tet(
254			vec4 P0, vec4 P1, vec4 P2, vec4 P3,
255			global_index_t P0_global_index,
256			global_index_t P1_global_index,
257			global_index_t P2_global_index,
258			global_index_t P3_global_index
259			) {
260			geo_debug_assert(has_vglobal_);
261		✗	init_with_tet(P0, P1, P2, P3);
262
263			// Offset for the 4 plane equations.
264			// Plane 0 is vertex at infinity.
265			index_t boff = 1;
266
267		✗	vglobal_[boff ] = P0_global_index;
268		✗	vglobal_[boff+1] = P1_global_index;
269		✗	vglobal_[boff+2] = P2_global_index;
270		✗	vglobal_[boff+3] = P3_global_index;
271		✗	}
272
273
274			/***********************************************************************/
275
276		✗	void ConvexCell::save(const std::string& filename, double shrink) const {
277		✗	std::ofstream out(filename.c_str());
278		✗	save(out, 1, shrink);
279		✗	}
280
281
282		✗	index_t ConvexCell::save(
283			std::ostream& out, global_index_t v_offset,
284			double shrink, bool borders_only
285			) const {
286
287			vec3 g = make_vec3(0.0, 0.0, 0.0);
288		✗	if(shrink != 0.0) {
289		✗	const_cast<ConvexCell*>(this)->compute_geometry();
290		✗	g = barycenter();
291			}
292
293		✗	vector<index_t> v2t(nb_v(),index_t(-1));
294		✗	vector<index_t> t_index(nb_t(),index_t(-1));
295			index_t nt=0;
296
297			{
298		✗	index_t t = first_valid_;
299		✗	while(t != END_OF_LIST) {
300			TriangleWithFlags T = get_triangle_and_flags(t);
301			vec4 p;
302		✗	if(geometry_dirty_) {
303		✗	p = compute_triangle_point(t);
304		✗	p.x /= p.w;
305		✗	p.y /= p.w;
306		✗	p.z /= p.w;
307			p.w = 1.0;
308			} else {
309		✗	p.x = triangle_point_[t].x;
310		✗	p.y = triangle_point_[t].y;
311		✗	p.z = triangle_point_[t].z;
312			p.w = 1.0;
313			}
314
315		✗	if(shrink != 0.0) {
316		✗	p.x = shrink * g.x + (1.0 - shrink) * p.x;
317		✗	p.y = shrink * g.y + (1.0 - shrink) * p.y;
318		✗	p.z = shrink * g.z + (1.0 - shrink) * p.z;
319			}
320			out << "v " << p.x << " " << p.y << " " << p.z << std::endl;
321		✗	t_index[t] = nt;
322		✗	++nt;
323		✗	v2t[T.i] = t;
324		✗	v2t[T.j] = t;
325		✗	v2t[T.k] = t;
326		✗	t = index_t(T.flags);
327			}
328			}
329
330		✗	for(index_t v=1; v<nb_v(); ++v) {
331		✗	if(borders_only &&
332		✗	has_vglobal() &&
333		✗	v_global_index(v) != global_index_t(-1) &&
334			v_global_index(v) != global_index_t(-2)
335			// index_t(-2) is for fluid free bndry
336			) {
337		✗	continue;
338			}
339		✗	if(v2t[v] != index_t(-1)) {
340			index_t t = v2t[v];
341		✗	out << "f ";
342			do {
343		✗	out << (t_index[t]+v_offset) << " ";
344			index_t lv = triangle_find_vertex(t,v);
345		✗	t = triangle_adjacent(t, (lv + 1)%3);
346		✗	} while(t != v2t[v]);
347			out << std::endl;
348			}
349			}
350
351		✗	return nt;
352			}
353
354		✗	void ConvexCell::for_each_Voronoi_vertex(
355			index_t v,
356			std::function<void(index_t)> vertex
357			) {
358			geo_debug_assert(!geometry_dirty_);
359		✗	if(v2t_[v] != END_OF_LIST) {
360		✗	index_t t = index_t(v2t_[v]);
361			do {
362		✗	vertex(t);
363			index_t lv = triangle_find_vertex(t,v);
364		✗	t = triangle_adjacent(t, (lv + 1)%3);
365		✗	} while(t != v2t_[v]);
366			}
367		✗	}
368
369			#if !defined(STANDALONE_CONVEX_CELL) && !defined(GEOGRAM_PSM)
370
371		✗	void ConvexCell::append_to_mesh(
372			GEO::Mesh* mesh, double shrink, bool borders_only,
373			GEO::Attribute<GEO::index_t>* facet_attr
374			) const {
375
376			global_index_t v_offset = mesh->vertices.nb();
377
378			vec3 g = make_vec3(0.0, 0.0, 0.0);
379		✗	if(shrink != 0.0) {
380		✗	const_cast<ConvexCell*>(this)->compute_geometry();
381		✗	g = barycenter();
382			}
383
384		✗	vector<index_t> v2t(nb_v(),index_t(-1));
385		✗	vector<index_t> t_index(nb_t(),index_t(-1));
386			index_t nt=0;
387
388			{
389		✗	index_t t = first_valid_;
390		✗	while(t != END_OF_LIST) {
391			TriangleWithFlags T = get_triangle_and_flags(t);
392		✗	vec4 p = compute_triangle_point(t);
393		✗	p.x /= p.w;
394		✗	p.y /= p.w;
395		✗	p.z /= p.w;
396		✗	p.w = 1.0;
397		✗	if(shrink != 0.0) {
398		✗	p.x = shrink * g.x + (1.0 - shrink) * p.x;
399		✗	p.y = shrink * g.y + (1.0 - shrink) * p.y;
400		✗	p.z = shrink * g.z + (1.0 - shrink) * p.z;
401			}
402		✗	mesh->vertices.create_vertex(p.data());
403		✗	t_index[t] = nt;
404		✗	++nt;
405		✗	v2t[T.i] = t;
406		✗	v2t[T.j] = t;
407		✗	v2t[T.k] = t;
408		✗	t = index_t(T.flags);
409			}
410			}
411
412		✗	for(index_t v=1; v<nb_v(); ++v) {
413		✗	if(borders_only &&
414		✗	has_vglobal() &&
415		✗	v_global_index(v) != global_index_t(-1) &&
416			v_global_index(v) != global_index_t(-2) // This one for
417			// fluid free bndry
418			) {
419		✗	continue;
420			}
421			std::vector<global_index_t> facet_vertices;
422		✗	if(v2t[v] != index_t(-1)) {
423			index_t t = v2t[v];
424			do {
425		✗	facet_vertices.push_back(t_index[t]+v_offset);
426			index_t lv = triangle_find_vertex(t,v);
427		✗	t = triangle_adjacent(t, (lv + 1)%3);
428		✗	} while(t != v2t[v]);
429			}
430		✗	if(facet_vertices.size() < 3) {
431			continue;
432			}
433		✗	global_index_t f = mesh->facets.create_polygon(
434			GEO::index_t(facet_vertices.size())
435			);
436		✗	for(index_t i=0; i<facet_vertices.size(); ++i) {
437		✗	mesh->facets.set_vertex(f, i, facet_vertices[i]);
438			}
439		✗	if(facet_attr != nullptr && facet_attr->is_bound()) {
440		✗	(*facet_attr)[f] = v_global_index(v);
441			}
442			}
443
444		✗	}
445
446			#endif
447
448			/***********************************************************************/
449
450		✗	bool ConvexCell::has_v_global_index(global_index_t v) const {
451			vbw_assert(has_vglobal_);
452		✗	for(index_t i=0; i<nb_v(); ++i) {
453		✗	if(vglobal_[i] == v) {
454			return true;
455			}
456			}
457			return false;
458			}
459
460
461		✗	void ConvexCell::clip_by_plane(vec4 eqn, global_index_t j) {
462			vbw_assert(has_vglobal_);
463		✗	clip_by_plane(eqn);
464		✗	vglobal_[nb_v()-1] = j;
465		✗	}
466
467		✗	void ConvexCell::clip_by_plane(vec4 eqn) {
468		✗	geometry_dirty_ = true;
469
470		✗	index_t lv = nb_v_;
471		✗	if(lv == max_v()) {
472		✗	grow_v();
473			}
474		✗	plane_eqn_[lv] = eqn;
475			vbw_assert(lv < max_v());
476		✗	++nb_v_;
477
478			// Step 1: Find conflict zone and link conflicted triangles
479			// (to recycle them in free list).
480
481			index_t conflict_head = END_OF_LIST;
482			index_t conflict_tail = END_OF_LIST;
483
484			// Classify triangles, compute conflict list and valid list.
485			// Note: This could be done by climbing from a random triangle,
486			// but here we prefer complete linear scan for several reasons:
487			// - it is more robust to numerical errors (here we are not
488			// using exact predicates).
489			// - the code is simpler.
490			// - and more importantly, we got no more than a few tenths of
491			// vertices.
492			// The 'climbing from a random triangle' strategy is implemented
493			// in clip_by_plane_fast(). We keep both implementations for now,
494			// until we make sure than one is more efficient than the other one.
495
496		✗	index_t t = first_valid_;
497		✗	first_valid_ = END_OF_LIST;
498		✗	while(t != END_OF_LIST) {
499			TriangleWithFlags T = get_triangle_and_flags(t);
500		✗	if(triangle_is_in_conflict(T,eqn)) {
501			set_triangle_flags(
502		✗	t, ushort(conflict_head) | ushort(CONFLICT_MASK)
503			);
504			conflict_head = t;
505		✗	if(conflict_tail == END_OF_LIST) {
506			conflict_tail = t;
507			}
508			} else {
509		✗	set_triangle_flags(t, ushort(first_valid_));
510		✗	first_valid_ = t;
511			}
512		✗	t = index_t(T.flags);
513			}
514
515		✗	triangulate_conflict_zone(lv, conflict_head, conflict_tail);
516		✗	}
517
518			// This version of clip_by_plane(), with a user-defined predicate,
519			// is duplicated from the standard version above. May be fixed by
520			// moving it to the header file (to have faster invokation of the
521			// predicate) and make the default version call it with PCK, but
522			// I do not want to do that because:
523			// - will make the header more heavy, longer compilation time
524			// - not sure about the impact on performance
525			// - the function is short, not a big drama to duplicate it...
526		✗	void ConvexCell::clip_by_plane(
527			vec4 eqn, global_index_t global_index,
528			std::function<bool(ushort,ushort)> triangle_conflict_predicate
529			) {
530		✗	geometry_dirty_ = true;
531
532		✗	index_t lv = nb_v_;
533		✗	if(lv == max_v()) {
534		✗	grow_v();
535			}
536		✗	plane_eqn_[lv] = eqn;
537			vbw_assert(lv < max_v());
538		✗	++nb_v_;
539
540			// Note: it is unlikely that this function is used without
541			// global indices (because without global indices, it would
542			// mean we are only using geometry, then we should use the
543			// default predicate), so we could probably make it mandatory
544			// to have global indices here.
545		✗	if(has_vglobal_) {
546		✗	vglobal_[nb_v()-1] = global_index;
547			}
548
549
550			// Step 1: Find conflict zone and link conflicted triangles
551			// (to recycle them in free list).
552
553			index_t conflict_head = END_OF_LIST;
554			index_t conflict_tail = END_OF_LIST;
555
556			// Classify triangles, compute conflict list and valid list.
557			// Note: This could be done by climbing from a random triangle,
558			// but here we prefer complete linear scan for several reasons:
559			// - it is more robust to numerical errors (here we are not
560			// using exact predicates).
561			// - the code is simpler.
562			// - and more importantly, we got no more than a few tenths of
563			// vertices.
564			// The 'climbing from a random triangle' strategy is implemented
565			// in clip_by_plane_fast(). We keep both implementations for now,
566			// until we make sure than one is more efficient than the other one
567			// (and clip_by_plane_fast() does not have a version with the user-
568			// defined predicate for now).
569
570		✗	index_t t = first_valid_;
571		✗	first_valid_ = END_OF_LIST;
572		✗	while(t != END_OF_LIST) {
573			TriangleWithFlags T = get_triangle_and_flags(t);
574		✗	if(triangle_conflict_predicate(ushort(t), ushort(nb_v()-1))) {
575			set_triangle_flags(
576		✗	t, ushort(conflict_head) | ushort(CONFLICT_MASK)
577			);
578			conflict_head = t;
579		✗	if(conflict_tail == END_OF_LIST) {
580			conflict_tail = t;
581			}
582			} else {
583		✗	set_triangle_flags(t, ushort(first_valid_));
584		✗	first_valid_ = t;
585			}
586		✗	t = index_t(T.flags);
587			}
588
589		✗	triangulate_conflict_zone(lv, conflict_head, conflict_tail);
590		✗	}
591
592
593		✗	void ConvexCell::clip_by_plane_fast(vec4 P, global_index_t j) {
594			vbw_assert(has_vglobal_);
595		✗	clip_by_plane_fast(P);
596		✗	vglobal_[nb_v()-1] = j;
597		✗	}
598
599		✗	void ConvexCell::clip_by_plane_fast(vec4 P) {
600		✗	geometry_dirty_ = true;
601		✗	index_t lv = nb_v_;
602		✗	if(lv == max_v()) {
603		✗	grow_v();
604			}
605		✗	plane_eqn_[lv] = P;
606			vbw_assert(lv < max_v());
607		✗	++nb_v_;
608
609			// Step 1: Find a good seed triangle (likely to
610			// be in conflict). If it is not in conflict,
611			// it is not a big problem (it will be just a
612			// bit slower), so we use inexact predicates
613			// here.
614
615			index_t t_init = END_OF_LIST;
616
617			{
618			index_t t_pred = END_OF_LIST;
619		✗	index_t t = first_valid_;
620		✗	if(t == END_OF_LIST) {
621		✗	return;
622			}
623
624		✗	auto triangle_distance = [this](index_t t_in, vec4 P_in) {
625		✗	Triangle T = get_triangle(t_in);
626			vbw_assert(T.i != VERTEX_AT_INFINITY);
627			vbw_assert(T.j != VERTEX_AT_INFINITY);
628			vbw_assert(T.k != VERTEX_AT_INFINITY);
629			vec4 p1 = vertex_plane(T.i);
630			vec4 p2 = vertex_plane(T.j);
631			vec4 p3 = vertex_plane(T.k);
632		✗	return det4x4(
633			p1.x, p2.x, p3.x, P_in.x,
634			p1.y, p2.y, p3.y, P_in.y,
635			p1.z, p2.z, p3.z, P_in.z,
636			p1.w, p2.w, p3.w, P_in.w
637		✗	);
638		✗	};
639
640			index_t count = 100;
641		✗	double t_dist = triangle_distance(t,P);
642
643		✗	still_walking:
644		✗	for(index_t le=0; le<3; ++le) {
645			index_t t_next = triangle_adjacent(t, le);
646		✗	if(t_next == t_pred) {
647		✗	continue;
648			}
649		✗	double t_next_dist = triangle_distance(t_next,P);
650		✗	if(t_next_dist < t_dist) {
651		✗	continue;
652			}
653		✗	--count;
654			t_pred = t;
655			t = t_next;
656			t_dist = t_next_dist;
657		✗	if(count > 0 && t_dist < 0.0) {
658		✗	goto still_walking;
659			}
660			}
661			t_init = t;
662			}
663
664			// note: t_init is now one triangle, probably in conflict
665			// (but not always: there can be degenerate cases where count
666			// reach zero). When t_init is in conflict, it is in general
667			// not the "highest" triangle (as one expect in a Delaunay
668			// triangulation code, but here it is different).
669
670
671			// Step 2: mark all the triangles that have the same
672			// conflict status as t_init with CONFLICT_MASK
673			// (even if t_init was not in conflict, then we will
674			// swap confict mask)
675
676		✗	bool t_init_is_in_conflict = triangle_is_in_conflict(
677			get_triangle_and_flags(t_init), P
678			);
679
680			{
681		✗	ushort* buff = (ushort*)alloca(max_t_*sizeof(ushort));
682			SmallStack_ushort S(buff, int(max_t_));
683		✗	S.push(ushort(t_init));
684			set_triangle_flags(
685			t_init,
686			get_triangle_flags(t_init) |
687		✗	ushort(CONFLICT_MASK) |
688			ushort(MARKED_MASK)
689			);
690		✗	while(!S.empty()) {
691			index_t t = index_t(S.top());
692			S.pop();
693		✗	for(index_t le=0; le<3; ++le) {
694			index_t t_neigh = triangle_adjacent(t,le);
695		✗	if(
696			(get_triangle_flags(t_neigh) & ushort(MARKED_MASK))
697			== 0
698			) {
699		✗	if( triangle_is_in_conflict(
700			get_triangle_and_flags(t_neigh), P
701			) == t_init_is_in_conflict
702			) {
703			set_triangle_flags(
704			t_neigh,
705			get_triangle_flags(t_neigh) |
706		✗	ushort(CONFLICT_MASK) |
707			ushort(MARKED_MASK)
708			);
709			S.push(ushort(t_neigh));
710			} else {
711			set_triangle_flags(
712			t_neigh,
713		✗	get_triangle_flags(t_neigh) |
714			ushort(MARKED_MASK)
715			);
716			}
717			}
718			}
719			}
720			}
721
722			// Step 3: update conflict list and active triangles list
723			index_t conflict_head = END_OF_LIST;
724			index_t conflict_tail = END_OF_LIST;
725			{
726			index_t new_first_valid = END_OF_LIST;
727		✗	index_t t = first_valid_;
728		✗	while(t != END_OF_LIST) {
729			bool t_is_in_conflict = triangle_is_marked_as_conflict(t);
730			index_t t_next =
731		✗	index_t(
732			get_triangle_flags(t) &
733			~(CONFLICT_MASK | MARKED_MASK)
734		✗	);
735			// Flip conflict flag if t_init was not in conflict.
736		✗	if(!t_init_is_in_conflict) {
737		✗	t_is_in_conflict = !t_is_in_conflict;
738			}
739		✗	if(t_is_in_conflict) {
740			set_triangle_flags(
741		✗	t, ushort(conflict_head) | ushort(CONFLICT_MASK)
742			);
743			conflict_head = t;
744		✗	if(conflict_tail == END_OF_LIST) {
745			conflict_tail = t;
746			}
747			} else {
748		✗	set_triangle_flags(t, ushort(new_first_valid));
749			new_first_valid = t;
750			}
751			t = t_next;
752			}
753		✗	first_valid_ = new_first_valid;
754			}
755
756		✗	triangulate_conflict_zone(lv, conflict_head, conflict_tail);
757			}
758
759			/***********************************************************************/
760
761		✗	void ConvexCell::triangulate_conflict_zone(
762			index_t lv, index_t conflict_head, index_t conflict_tail
763			) {
764			// Special case: no triangle in conflict.
765		✗	if(conflict_head == END_OF_LIST) {
766			return;
767			}
768
769
770			// Link the vertices on the border of the conflict zone.
771			// Consider the edge e of triangle t such that:
772			// t is not marked as conflict
773			// triangle_adjacent(t,e) is marked as conflict
774			// Let v1 = triangle_vertex(t, (e+1)%3)
775			// v2 = triangle_vertex(t, (e+2)%3)
776			// We set:
777			// v2t_[v1] = t
778			// v2e_[v1] = e
779			// Once done for all the triangles, one can traverse the
780			// border of the conflict zone with:
781			//
782			// v = first_v_on_border]
783			// do {
784			// t = v2t_[v];
785			// e = v2t_[e];
786			// do something with t,e
787			// v = triangle_vertex(t, (e+2)%3);
788			// } while(v != first_v_on_border]);
789
790			geo_debug(index_t nb = 0); // for sanity check, nbr of vertices on border
791			// of conflict zone.
792
793			VBW::index_t first_v_on_border = END_OF_LIST;
794		✗	for(
795			VBW::ushort t = first_triangle();
796		✗	t != END_OF_LIST;
797			t = next_triangle(t)
798			) {
799			vbw_assert(!triangle_is_marked_as_conflict(t));
800
801		✗	if(triangle_is_marked_as_conflict(triangle_adjacent(t,0))) {
802			first_v_on_border = triangle_vertex(t,1);
803		✗	v2t_[first_v_on_border] = t;
804		✗	v2e_[first_v_on_border] = 0;
805			geo_debug(++nb);
806			}
807		✗	if(triangle_is_marked_as_conflict(triangle_adjacent(t,1))) {
808			first_v_on_border = triangle_vertex(t,2);
809		✗	v2t_[first_v_on_border] = t;
810		✗	v2e_[first_v_on_border] = 1;
811			geo_debug(++nb);
812			}
813		✗	if(triangle_is_marked_as_conflict(triangle_adjacent(t,2))) {
814			first_v_on_border = triangle_vertex(t,0);
815		✗	v2t_[first_v_on_border] = t;
816		✗	v2e_[first_v_on_border] = 2;
817			geo_debug(++nb);
818			}
819			}
820
821			geo_debug(index_t nb2 = 0); // for sanity check, number of vertices on border
822			// of conflict zone (should match nb).
823
824			// Traverse the list of edges on the border of the conflict zone
825			// (see previous comment block for explanations). For each edge
826			// on the border of the conflict zone, generate a new triangle.
827			// - Connect it to the triangle on the border of the conflict zone
828			// - Connect it to the previous new triangle
829			// - Special case: in the end, connect the first new triangle with
830			// the last one.
831
832			// check we are not in the special case with all triangles in conflict
833		✗	if(first_v_on_border != END_OF_LIST) {
834			VBW::index_t v = first_v_on_border;
835			VBW::ushort prev_new_t = VBW::ushort(-1);
836			VBW::ushort first_new_t = VBW::ushort(-1);
837			do {
838			geo_debug(++nb2);
839		✗	index_t t = v2t_[v];
840		✗	index_t e = v2e_[v];
841		✗	index_t v1 = triangle_vertex(t, (e+1)%3);
842		✗	index_t v2 = triangle_vertex(t, (e+2)%3);
843			vbw_assert(v1 == v);
844		✗	VBW::ushort new_t = VBW::ushort(new_triangle(lv, v2, v1));
845			set_triangle_adjacent(new_t, 0, t);
846			set_triangle_adjacent(t, e, new_t);
847		✗	if(prev_new_t == VBW::ushort(-1)) {
848			first_new_t = new_t;
849			} else {
850			set_triangle_adjacent(prev_new_t, 2, new_t);
851			set_triangle_adjacent(new_t, 1, prev_new_t);
852			}
853			prev_new_t = new_t;
854			v = v2;
855		✗	} while (v != first_v_on_border);
856			set_triangle_adjacent(prev_new_t, 2, first_new_t);
857			set_triangle_adjacent(first_new_t, 1, prev_new_t);
858			}
859
860			vbw_assert(nb2 == nb);
861
862			// Recycle triangles in conflict zone
863		✗	set_triangle_flags(conflict_tail, ushort(first_free_));
864		✗	first_free_ = conflict_head;
865			}
866
867
868			/***********************************************************************/
869
870		✗	bool ConvexCell::triangle_is_in_conflict(
871			TriangleWithFlags T, const vec4& eqn
872			) const {
873			#ifndef STANDALONE_CONVEX_CELL
874		✗	if(use_exact_predicates_) {
875		✗	if(T.i == VERTEX_AT_INFINITY) {
876		✗	vec3 E = make_vec3(eqn.x, eqn.y, eqn.z);
877			vec3 n2 = vertex_plane_normal(T.j);
878			vec3 n3 = vertex_plane_normal(T.k);
879		✗	return(GEO::PCK::det_3d(E.data(), n2.data(), n3.data()) <= 0);
880			}
881
882		✗	if(T.j == VERTEX_AT_INFINITY) {
883		✗	vec3 E = make_vec3(eqn.x, eqn.y, eqn.z);
884			vec3 n3 = vertex_plane_normal(T.k);
885			vec3 n1 = vertex_plane_normal(T.i);
886		✗	return(GEO::PCK::det_3d(n1.data(), E.data(), n3.data()) <= 0);
887			}
888
889		✗	if(T.k == VERTEX_AT_INFINITY) {
890		✗	vec3 E = make_vec3(eqn.x, eqn.y, eqn.z);
891			vec3 n1 = vertex_plane_normal(T.i);
892			vec3 n2 = vertex_plane_normal(T.j);
893		✗	return(GEO::PCK::det_3d(n1.data(), n2.data(), E.data()) <= 0);
894			}
895
896			// The triangle is in conflict with eqn if the
897			// result of compute_triangle_point(t) injected in eqn
898			// has a negative sign.
899			// Examining the formula in compute_triangle_point(),
900			// this corresponds to (minus) the 4x4 determinant of the
901			// 4 plane equations developed w.r.t. the 4th column.
902			// (see Edelsbrunner - Simulation of Simplicity for similar examples
903			// of computations).
904
905			vec4 p1 = vertex_plane(T.i);
906			vec4 p2 = vertex_plane(T.j);
907			vec4 p3 = vertex_plane(T.k);
908
909		✗	return(GEO::PCK::det_4d(
910			p1.data(), p2.data(), p3.data(), eqn.data()) >= 0
911		✗	);
912			}
913			#endif
914
915			double det = 0.0;
916
917			// If one of the vertices of the triangle is the
918			// vertex at infinity, then the triangle is in conflict
919			// with eqn if the oriented director vector of the intersection
920			// between the two planes of the two other vertices
921			// is opposite to the normal vector to eqn.
922
923		✗	if(T.i == VERTEX_AT_INFINITY) {
924			vec3 n2 = vertex_plane_normal(T.j);
925			vec3 n3 = vertex_plane_normal(T.k);
926		✗	det = -det3x3(
927		✗	eqn.x, n2.x, n3.x,
928		✗	eqn.y, n2.y, n3.y,
929		✗	eqn.z, n2.z, n3.z
930			);
931		✗	} else if(T.j == VERTEX_AT_INFINITY) {
932			vec3 n3 = vertex_plane_normal(T.k);
933			vec3 n1 = vertex_plane_normal(T.i);
934		✗	det = -det3x3(
935		✗	n1.x, eqn.x, n3.x,
936		✗	n1.y, eqn.y, n3.y,
937		✗	n1.z, eqn.z, n3.z
938			);
939		✗	} else if(T.k == VERTEX_AT_INFINITY) {
940			vec3 n1 = vertex_plane_normal(T.i);
941			vec3 n2 = vertex_plane_normal(T.j);
942		✗	det = -det3x3(
943		✗	n1.x, n2.x, eqn.x,
944		✗	n1.y, n2.y, eqn.y,
945		✗	n1.z, n2.z, eqn.z
946			);
947			} else {
948
949			// The triangle is in conflict with eqn if the
950			// result of compute_triangle_point(t) injected in eqn
951			// has a negative sign.
952			// Examining the formula in compute_triangle_point(),
953			// this corresponds to (minus) the 4x4 determinant of the 4 plane
954			// equations developed w.r.t. the 4th column.
955			// (see Edelsbrunner - Simulation of Simplicity for similar examples
956			// of computations).
957
958			vec4 p1 = vertex_plane(T.i);
959			vec4 p2 = vertex_plane(T.j);
960			vec4 p3 = vertex_plane(T.k);
961		✗	det = det4x4(
962		✗	p1.x, p2.x, p3.x, eqn.x,
963		✗	p1.y, p2.y, p3.y, eqn.y,
964		✗	p1.z, p2.z, p3.z, eqn.z,
965		✗	p1.w, p2.w, p3.w, eqn.w
966			);
967			}
968		✗	return (det > 0.0);
969			}
970
971		✗	vec4 ConvexCell::compute_triangle_point(index_t t) const {
972
973			double infinite_len = 16.0;
974			TriangleWithFlags T = get_triangle_and_flags(t);
975
976			// Special cases with one of the three vertices at infinity.
977		✗	if(T.i == VERTEX_AT_INFINITY) {
978			vec4 Pj = vertex_plane(T.j);
979			vec4 Pk = vertex_plane(T.k);
980		✗	vec3 Njk = normalize(cross(
981			make_vec3(Pj.x, Pj.y, Pj.z),
982			make_vec3(Pk.x, Pk.y, Pk.z)
983			));
984		✗	index_t t_adj = t_adj_[t].i; // vv2t(T.k, T.j);
985			vbw_assert(!triangle_is_infinite(t_adj));
986		✗	vec4 result = compute_triangle_point(t_adj);
987		✗	result.x += result.w * Njk.x * infinite_len;
988		✗	result.y += result.w * Njk.y * infinite_len;
989		✗	result.z += result.w * Njk.z * infinite_len;
990			return result;
991		✗	} else if(T.j == VERTEX_AT_INFINITY) {
992			vec4 Pk = vertex_plane(T.k);
993			vec4 Pi = vertex_plane(T.i);
994		✗	vec3 Nki = normalize(cross(
995			make_vec3(Pk.x, Pk.y, Pk.z),
996			make_vec3(Pi.x, Pi.y, Pi.z)
997			));
998		✗	index_t t_adj = t_adj_[t].j; // vv2t(T.i, T.k);
999			vbw_assert(!triangle_is_infinite(t_adj));
1000		✗	vec4 result = compute_triangle_point(t_adj);
1001		✗	result.x += result.w * Nki.x * infinite_len;
1002		✗	result.y += result.w * Nki.y * infinite_len;
1003		✗	result.z += result.w * Nki.z * infinite_len;
1004		✗	return result;
1005		✗	} else if(T.k == VERTEX_AT_INFINITY) {
1006			vec4 Pi = vertex_plane(T.i);
1007			vec4 Pj = vertex_plane(T.j);
1008		✗	vec3 Nij = normalize(cross(
1009			make_vec3(Pi.x, Pi.y, Pi.z),
1010			make_vec3(Pj.x, Pj.y, Pj.z)
1011			));
1012		✗	index_t t_adj = t_adj_[t].k; // vv2t(T.j, T.i);
1013			vbw_assert(!triangle_is_infinite(t_adj));
1014		✗	vec4 result = compute_triangle_point(t_adj);
1015		✗	result.x += result.w * Nij.x * infinite_len;
1016		✗	result.y += result.w * Nij.y * infinite_len;
1017		✗	result.z += result.w * Nij.z * infinite_len;
1018		✗	return result;
1019			}
1020
1021			// Get the plane equations associated with each vertex of t
1022
1023			vec4 pi1 = vertex_plane(T.i);
1024			vec4 pi2 = vertex_plane(T.j);
1025			vec4 pi3 = vertex_plane(T.k);
1026
1027			// Find the intersection of the three planes using Cramer's formula.
1028			// (see Edelsbrunner - Simulation of Simplicity for other examples).
1029			//
1030			// Cramer's formula: each component of the solution is obtained as
1031			// the ratio of two determinants:
1032			// - the determinant of the system where the ith column is replaced
1033			// with the rhs
1034			// divided by:
1035			// - the determinant of the system.
1036			//
1037			// System of equations to be solved:
1038			// pi1.x * x + pi1.y * y + pi1.z * z = -pi1.w
1039			// pi2.x * x + pi2.y * y + pi2.z * z = -pi2.w
1040			// pi3.x * x + pi3.y * y + pi3.z * z = -pi3.w
1041			//
1042			// Expression of the solution given by Cramer's formula:
1043			// | -pi1.w p1.y pi1.z | | pi1.x pi1.y pi1.z |
1044			// x = | -pi2.w p2.y pi2.z | / | pi2.x pi2.y pi2.z |
1045			// | -pi3.w p3.y pi3.z | | pi3.x pi3.y pi3.z |
1046			//
1047			// | pi1.x -p1.w pi1.z | | pi1.x pi1.y pi1.z |
1048			// y = | pi2.x -p2.w pi2.z | / | pi2.x pi2.y pi2.z |
1049			// | pi3.x -p3.w pi3.z | | pi3.x pi3.y pi3.z |
1050			//
1051			// | pi1.x p1.y -pi1.w | | pi1.x pi1.y pi1.z |
1052			// z = | pi2.x p2.y -pi2.w | / | pi2.x pi2.y pi2.z |
1053			// | pi3.x p3.y -pi3.w | | pi3.x pi3.y pi3.z |
1054
1055			vec4 result;
1056
1057		✗	result.x = -det3x3(
1058			pi1.w, pi1.y, pi1.z,
1059			pi2.w, pi2.y, pi2.z,
1060			pi3.w, pi3.y, pi3.z
1061			);
1062
1063		✗	result.y = -det3x3(
1064			pi1.x, pi1.w, pi1.z,
1065			pi2.x, pi2.w, pi2.z,
1066			pi3.x, pi3.w, pi3.z
1067			);
1068
1069		✗	result.z = -det3x3(
1070			pi1.x, pi1.y, pi1.w,
1071			pi2.x, pi2.y, pi2.w,
1072			pi3.x, pi3.y, pi3.w
1073			);
1074
1075			result.w = det3x3(
1076			pi1.x, pi1.y, pi1.z,
1077			pi2.x, pi2.y, pi2.z,
1078			pi3.x, pi3.y, pi3.z
1079			);
1080
1081		✗	return result;
1082			}
1083
1084			/***********************************************************************/
1085
1086		✗	void ConvexCell::grow_v() {
1087		✗	max_v_ *= 2;
1088		✗	plane_eqn_.resize(max_v_);
1089		✗	vglobal_.resize(max_v_, global_index_t(-1));
1090		✗	v2t_.resize(max_v_, ushort(-1));
1091		✗	v2e_.resize(max_v_, uchar(-1));
1092		✗	}
1093
1094		✗	void ConvexCell::grow_t() {
1095		✗	max_t_ *= 2;
1096		✗	t_.resize(max_t_);
1097		✗	t_adj_.resize(max_t_);
1098		✗	if(has_tflags_) {
1099		✗	tflags_.resize(max_t_,0);
1100			}
1101		✗	}
1102
1103			/***********************************************************************/
1104
1105		✗	void ConvexCell::kill_vertex(index_t v) {
1106		✗	for(index_t t=0; t<nb_t(); ++t) {
1107			Triangle T = get_triangle(t);
1108		✗	if(T.i == v) {
1109			T.i = VERTEX_AT_INFINITY;
1110			}
1111		✗	if(T.j == v) {
1112			T.j = VERTEX_AT_INFINITY;
1113			}
1114		✗	if(T.k == v) {
1115			T.k = VERTEX_AT_INFINITY;
1116			}
1117		✗	t_[t].i = T.i;
1118		✗	t_[t].j = T.j;
1119		✗	t_[t].k = T.k;
1120			}
1121		✗	}
1122
1123			/***********************************************************************/
1124
1125		✗	void ConvexCell::compute_geometry() {
1126		✗	if(!geometry_dirty_) {
1127			return;
1128			}
1129
1130		✗	triangle_point_.resize(nb_t());
1131		✗	v2t_.assign(max_v(),END_OF_LIST);
1132
1133		✗	index_t t = first_valid_;
1134		✗	while(t != END_OF_LIST) {
1135			TriangleWithFlags T = get_triangle_and_flags(t);
1136		✗	vec4 p = compute_triangle_point(t);
1137		✗	triangle_point_[t] = make_vec3(p.x/p.w, p.y/p.w, p.z/p.w);
1138		✗	v2t_[T.i] = ushort(t);
1139		✗	v2t_[T.j] = ushort(t);
1140		✗	v2t_[T.k] = ushort(t);
1141		✗	t = index_t(T.flags);
1142			}
1143
1144		✗	geometry_dirty_ = false;
1145			}
1146
1147		✗	inline double triangle_area(vec3 p1, vec3 p2, vec3 p3) {
1148		✗	double Ux = p2.x - p1.x;
1149		✗	double Uy = p2.y - p1.y;
1150		✗	double Uz = p2.z - p1.z;
1151		✗	double Vx = p3.x - p1.x;
1152		✗	double Vy = p3.y - p1.y;
1153		✗	double Vz = p3.z - p1.z;
1154		✗	double Wx = Uy * Vz - Uz * Vy;
1155		✗	double Wy = Uz * Vx - Ux * Vz;
1156		✗	double Wz = Ux * Vy - Uy * Vx;
1157		✗	return 0.5 * ::sqrt(Wx*Wx + Wy*Wy + Wz*Wz);
1158			}
1159
1160		✗	double ConvexCell::facet_area(index_t v) const {
1161			vbw_assert(v < nb_v());
1162			vbw_assert(!geometry_dirty_);
1163
1164			ushort t1t2[2];
1165			index_t cur=0;
1166			double result = 0.0;
1167
1168		✗	if(v2t_[v] != END_OF_LIST) {
1169		✗	index_t t = v2t_[v];
1170			index_t count = 0;
1171			do {
1172		✗	if(cur < 2) {
1173		✗	t1t2[cur] = ushort(t);
1174			} else {
1175		✗	result += triangle_area(
1176		✗	triangle_point_[t1t2[0]],
1177		✗	triangle_point_[t1t2[1]],
1178			triangle_point_[t]
1179			);
1180		✗	t1t2[1] = ushort(t);
1181			}
1182		✗	++cur;
1183			index_t lv = triangle_find_vertex(t,v);
1184		✗	t = triangle_adjacent(t, (lv + 1)%3);
1185			++count;
1186		✗	geo_assert(count < 100000);
1187		✗	} while(t != v2t_[v]);
1188			}
1189
1190		✗	return result;
1191			}
1192
1193		✗	inline double tet_volume(vec3 p1, vec3 p2, vec3 p3, vec3 p4) {
1194		✗	double Ux = p2.x - p1.x;
1195		✗	double Uy = p2.y - p1.y;
1196		✗	double Uz = p2.z - p1.z;
1197
1198		✗	double Vx = p3.x - p1.x;
1199		✗	double Vy = p3.y - p1.y;
1200		✗	double Vz = p3.z - p1.z;
1201
1202		✗	double Wx = p4.x - p1.x;
1203		✗	double Wy = p4.y - p1.y;
1204		✗	double Wz = p4.z - p1.z;
1205
1206		✗	double UVx = Uy * Vz - Uz * Vy;
1207		✗	double UVy = Uz * Vx - Ux * Vz;
1208		✗	double UVz = Ux * Vy - Uy * Vx;
1209
1210		✗	return ::fabs(
1211		✗	UVx * Wx + UVy * Wy + UVz * Wz
1212		✗	) / 6.0;
1213			}
1214
1215		✗	double ConvexCell::volume() const {
1216			vbw_assert(!geometry_dirty_);
1217			double result = 0.0;
1218
1219			ushort t_origin = END_OF_LIST;
1220		✗	for(index_t v=0; v<nb_v_; ++v) {
1221		✗	if(v2t_[v] == END_OF_LIST) {
1222		✗	continue;
1223			}
1224		✗	if(t_origin == END_OF_LIST) {
1225			t_origin = v2t_[v];
1226		✗	continue;
1227			}
1228			ushort t1t2[2];
1229			index_t cur=0;
1230		✗	index_t t = v2t_[v];
1231
1232			index_t count = 0;
1233			do {
1234		✗	if(cur < 2) {
1235		✗	t1t2[cur] = ushort(t);
1236			} else {
1237		✗	result += tet_volume(
1238			triangle_point_[t_origin],
1239		✗	triangle_point_[t1t2[0]],
1240		✗	triangle_point_[t1t2[1]],
1241			triangle_point_[t]
1242			);
1243		✗	t1t2[1] = ushort(t);
1244			}
1245		✗	++cur;
1246			index_t lv = triangle_find_vertex(t,v);
1247		✗	t = triangle_adjacent(t, (lv + 1)%3);
1248			++count;
1249		✗	geo_assert(count < 100000);
1250		✗	} while(t != v2t_[v]);
1251			}
1252		✗	return result;
1253			}
1254
1255		✗	vec3 ConvexCell::barycenter() const {
1256			vec3 result;
1257			double m;
1258		✗	compute_mg(m, result);
1259		✗	if(m != 0.0) {
1260		✗	result.x /= m;
1261		✗	result.y /= m;
1262		✗	result.z /= m;
1263			}
1264		✗	return result;
1265			}
1266
1267		✗	void ConvexCell::compute_mg(double& m, vec3& result) const {
1268			vbw_assert(!geometry_dirty_);
1269		✗	result = make_vec3(0.0, 0.0, 0.0);
1270		✗	m = 0.0;
1271
1272			ushort t_origin = END_OF_LIST;
1273		✗	for(index_t v=0; v<nb_v_; ++v) {
1274		✗	if(v2t_[v] == END_OF_LIST) {
1275		✗	continue;
1276			}
1277		✗	if(t_origin == END_OF_LIST) {
1278			t_origin = v2t_[v];
1279		✗	continue;
1280			}
1281			ushort t1t2[2];
1282			index_t cur=0;
1283		✗	index_t t = v2t_[v];
1284			index_t count = 0;
1285			do {
1286		✗	if(cur < 2) {
1287		✗	t1t2[cur] = ushort(t);
1288			} else {
1289		✗	vec3 p = triangle_point_[t_origin];
1290		✗	vec3 q = triangle_point_[t1t2[0]];
1291		✗	vec3 r = triangle_point_[t1t2[1]];
1292		✗	vec3 s = triangle_point_[t];
1293		✗	double cur_m = tet_volume(p,q,r,s);
1294		✗	m += cur_m;
1295		✗	result.x += cur_m*(p.x + q.x + r.x + s.x)/4.0;
1296		✗	result.y += cur_m*(p.y + q.y + r.y + s.y)/4.0;
1297		✗	result.z += cur_m*(p.z + q.z + r.z + s.z)/4.0;
1298		✗	t1t2[1] = ushort(t);
1299			}
1300		✗	++cur;
1301			index_t lv = triangle_find_vertex(t,v);
1302		✗	t = triangle_adjacent(t, (lv + 1)%3);
1303			++count;
1304		✗	geo_assert(count < 100000);
1305		✗	} while(t != v2t_[v]);
1306			}
1307		✗	}
1308
1309			/***********************************************************************/
1310
1311
1312		✗	double ConvexCell::squared_radius(vec3 center) const {
1313		✗	double result = 0.0;
1314		✗	index_t t = first_valid_;
1315		✗	while(t != END_OF_LIST) {
1316			TriangleWithFlags T = get_triangle_and_flags(t);
1317		✗	if(geometry_dirty_) {
1318		✗	vec4 p4 = compute_triangle_point(t);
1319		✗	vec3 p3 = make_vec3(
1320		✗	p4.x/p4.w, p4.y/p4.w, p4.z/p4.w
1321			);
1322		✗	result = std::max(result, squared_distance(center,p3));
1323			} else {
1324		✗	vec3 p = triangle_point_[t];
1325		✗	result = std::max(result, squared_distance(center,p));
1326			}
1327		✗	t = index_t(T.flags);
1328			}
1329		✗	return result;
1330			}
1331
1332		✗	double ConvexCell::squared_inner_radius(vec3 center) const {
1333		✗	double result = std::numeric_limits<double>::max();
1334		✗	for(index_t v=0; v<nb_v(); ++v) {
1335			vec4 P = vertex_plane(v);
1336			// Ignore vertex at infinity.
1337		✗	if(P.x == 0.0 && P.y == 0.0 && P.z == 0.0) {
1338			continue;
1339			}
1340		✗	result = std::min(
1341			result, squared_point_plane_distance(center, P)
1342			);
1343			}
1344		✗	return result;
1345			}
1346
1347
1348		✗	void ConvexCell::connect_triangles() {
1349
1350			// create array that maps vertices pairs to triangles.
1351			// size of the array is nb_v squared.
1352			// If nb_v is small, allocate it on the stack, else
1353			// allocate it on the heap (allocating on the stack is
1354			// interesting for multithreading).
1355
1356			const index_t MAX_NV_ON_STACK = 50;
1357			index_t NV = nb_v();
1358			ushort* vv2t =
1359		✗	(NV <= MAX_NV_ON_STACK) ? (ushort*)alloca(NV*NV*sizeof(ushort))
1360		✗	: new ushort[NV*NV]
1361			;
1362
1363			#ifdef GEO_DEBUG
1364			for(index_t i=0; i<NV*NV; ++i) {
1365			vv2t[i] = END_OF_LIST;
1366			}
1367			#endif
1368
1369			// For each triangle t, remember
1370			// that t is adjacent to the three
1371			// oriented edges (j,k), (k,i), (i,j)
1372		✗	for(
1373			ushort t = first_triangle();
1374		✗	t != END_OF_LIST;
1375			t = next_triangle(t)
1376			) {
1377		✗	Triangle T = t_[t];
1378			vbw_assert(T.i < nb_v());
1379			vbw_assert(T.j < nb_v());
1380			vbw_assert(T.k < nb_v());
1381		✗	vv2t[NV*T.j + T.k] = t;
1382		✗	vv2t[NV*T.k + T.i] = t;
1383		✗	vv2t[NV*T.i + T.j] = t;
1384			}
1385
1386			// For each triangle t, find the
1387			// three triangles adjacent to its
1388			// three edges (k,j), (i,k), (j,i)
1389			// (in reverse order because we
1390			// want to find the three triangles
1391			// on the other side of the edges)
1392			for(
1393			ushort t = first_triangle();
1394		✗	t != END_OF_LIST;
1395			t = next_triangle(t)
1396			) {
1397		✗	Triangle T = t_[t];
1398			vbw_assert(vv2t[NV*T.j + T.i] != END_OF_LIST);
1399			vbw_assert(vv2t[NV*T.k + T.j] != END_OF_LIST);
1400			vbw_assert(vv2t[NV*T.i + T.k] != END_OF_LIST);
1401		✗	t_adj_[t] = make_triangle(
1402		✗	vv2t[NV*T.k + T.j],
1403		✗	vv2t[NV*T.i + T.k],
1404		✗	vv2t[NV*T.j + T.i]
1405			);
1406			}
1407
1408		✗	if(NV > MAX_NV_ON_STACK) {
1409		✗	delete[] vv2t;
1410			}
1411		✗	}
1412
1413			/************************************************************************/
1414
1415
1416			}
1417