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			#ifndef GEOGRAM_VORONOI_CONVEX_CELL
41			#define GEOGRAM_VORONOI_CONVEX_CELL
42
43			#ifndef STANDALONE_CONVEX_CELL
44			#include <geogram/basic/common.h>
45			#include <geogram/basic/memory.h>
46			#include <geogram/basic/numeric.h>
47			#include <geogram/basic/geometry.h>
48			# ifndef GEOGRAM_PSM
49			# include <geogram/basic/attributes.h>
50			# endif
51			#endif
52
53			#include <string>
54			#include <vector>
55			#include <iostream>
56			#include <cmath>
57			#include <cassert>
58
59
60
61			/**
62			* \file geogram/voronoi/convex_cell.h
63			* \brief Class to compute the intersection of a set of half-spaces in 3D.
64			* \details Has its own types for points and vectors because it can be used
65			* independently from Geogram. In that case, define STANDALONE_CONVEX_CELL
66			*/
67
68			#ifndef STANDALONE_CONVEX_CELL
69			namespace GEO {
70			class Mesh;
71			class PeriodicDelaunay3d;
72			}
73			#endif
74
75
76			namespace VBW {
77
78			#ifdef STANDALONE_CONVEX_CELL
79			using std::vector;
80			typedef unsigned int index_t;
81			typedef unsigned int global_index_t;
82			# define vbw_assert(x) assert(x)
83			struct vec2 {
84			double x;
85			double y;
86			};
87			struct vec3 {
88			double x;
89			double y;
90			double z;
91			};
92			struct vec4 {
93			double x;
94			double y;
95			double z;
96			double w;
97			};
98			#else
99			using GEO::vector;
100			typedef unsigned int index_t; // Always 32 bits
101			typedef GEO::index_t global_index_t; // Possibly 64 bits in GARGANTUA mode
102			# define vbw_assert(x) geo_debug_assert(x)
103			using GEO::vec2;
104			using GEO::vec3;
105			using GEO::vec4;
106			#endif
107
108			/******************************************************************************/
109
110
111			/**
112			* \brief Creates a vec2 from its components.
113			* \param[in] x , y the components of the
114			* vector.
115			* \return the created vector.
116			*/
117			inline vec2 make_vec2(
118			double x, double y
119			) {
120			vec2 result;
121			result.x = x;
122			result.y = y;
123			return result;
124			}
125
126
127			/**
128			* \brief Creates a vec3 from its components.
129			* \param[in] x , y , z the components of the
130			* vector.
131			* \return the created vector.
132			*/
133			inline vec3 make_vec3(
134			double x, double y, double z
135			) {
136			vec3 result;
137		✗	result.x = x;
138		✗	result.y = y;
139		✗	result.z = z;
140			return result;
141			}
142
143
144			/**
145			* \brief Computes the cross product between
146			* two vectors.
147			* \param[in] v1 , v2 the two vectors.
148			* \return the cross product between \p v1 and
149			* \p v2.
150			*/
151			inline vec3 cross(vec3 v1, vec3 v2) {
152		✗	return make_vec3(
153		✗	v1.y*v2.z - v1.z*v2.y,
154		✗	v1.z*v2.x - v1.x*v2.z,
155		✗	v1.x*v2.y - v1.y*v2.x
156			);
157			}
158
159			/**
160			* \brief Computes the dot product between
161			* two vectors.
162			* \param[in] v1 , v2 the two vectors.
163			* \return the dot product between \p v1 and
164			* \p v2.
165			*/
166			inline double dot(vec3 v1, vec3 v2) {
167			return (
168			v1.x*v2.x + v1.y*v2.y + v1.z*v2.z
169			);
170			}
171
172			/**
173			* \brief Computes the squared length of a vector.
174			* \param[in] v the vector.
175			* \return the squared length of \p v.
176			*/
177			inline double squared_length(vec3 v) {
178		✗	return (v.x*v.x + v.y*v.y + v.z*v.z);
179			}
180
181			/**
182			* \brief Computes the squared distance between two points.
183			* \param[in] v , w the two points.
184			* \return the squared distance between \p v and \p w.
185			*/
186			inline double squared_distance(vec3 v, vec3 w) {
187		✗	double dx = w.x-v.x;
188		✗	double dy = w.y-v.y;
189		✗	double dz = w.z-v.z;
190		✗	return (dx*dx+dy*dy+dz*dz);
191			}
192
193			/**
194			* \brief Computes the length of a vector.
195			* \param[in] v the vector.
196			* \return the length of \p v.
197			*/
198			inline double length(vec3 v) {
199		✗	return ::sqrt(squared_length(v));
200			}
201
202			/**
203			* \brief Computes a normalized vector.
204			* \param[in] v the vector.
205			* \return a vector with the same direction
206			* as \v and unit length.
207			*/
208		✗	inline vec3 normalize(vec3 v) {
209		✗	double s = 1.0/length(v);
210		✗	return make_vec3(
211			s*v.x, s*v.y, s*v.z
212		✗	);
213			}
214
215			/**
216			* \brief Creates a vec4 from its components.
217			* \param[in] x , y , z , w the components of the
218			* vector.
219			* \return the created vector.
220			*/
221			inline vec4 make_vec4(
222			double x, double y, double z, double w
223			) {
224			vec4 result;
225			result.x = x;
226			result.y = y;
227			result.z = z;
228			result.w = w;
229			return result;
230			}
231
232			/**
233			* \brief Computes the dot product between
234			* two vectors.
235			* \param[in] v1 , v2 the two vectors.
236			* \return the dot product between \p v1 and
237			* \p v2.
238			*/
239			inline double dot(vec4 v1, vec4 v2) {
240			return (
241			v1.x*v2.x + v1.y*v2.y +
242			v1.z*v2.z + v1.w*v2.w
243			);
244			}
245
246			/**
247			* \brief Computes the squared length of a vector.
248			* \param[in] v the vector.
249			* \return the squared length of \p v.
250			*/
251			inline double squared_length(vec4 v) {
252			return (
253			v.x*v.x + v.y*v.y +
254			v.z*v.z + v.w*v.w
255			);
256			}
257
258			/**
259			* \brief Computes the length of a vector.
260			* \param[in] v the vector.
261			* \return the length of \p v.
262			*/
263			inline double length(vec4 v) {
264			return ::sqrt(squared_length(v));
265			}
266
267			/**
268			* \brief Computes the squared distance between a point and a plane
269			* \param[in] p the point
270			* \param[in] P the plane equation
271			* \return the squared distance between p and P
272			*/
273			inline double squared_point_plane_distance(VBW::vec3 p, VBW::vec4 P) {
274		✗	double result = P.x*p.x + P.y*p.y + P.z*p.z + P.w;
275		✗	result = (result*result) / (P.x*P.x + P.y*P.y + P.z*P.z);
276			return result;
277			}
278
279			/**
280			* \brief Some constants for the flags
281			* in TriangleWithFlags.
282			* \see TriangleWithFlags.
283			*/
284			enum {
285			CONFLICT_MASK = 32768, /**< \brief The mask for conflict triangles. */
286			MARKED_MASK = 16384, /**< \brief The mask for marked triangles. */
287			END_OF_LIST = 16383, /**< \brief Constant to indicate end of list.*/
288			VERTEX_AT_INFINITY = 0 /**< \brief Vertex at infinity. */
289			};
290
291
292			/**
293			* \brief Type for flags.
294			*/
295			typedef unsigned char uchar;
296
297			/**
298			* \brief Type for local indices.
299			* \details Valid values are between 0 and 32766.
300			* Full range is not used due to bookkeeping reasons,
301			* \see TriangleWithFlags.
302			*/
303			typedef unsigned short ushort;
304
305			/**
306			* \brief A triangle with the local indices of its three
307			* vertices.
308			*/
309			struct Triangle {
310			ushort i;
311			ushort j;
312			ushort k;
313			ushort operator[](unsigned int index) const {
314			vbw_assert(index < 3);
315		✗	return (&i)[index];
316			}
317			ushort& operator[](unsigned int index) {
318			vbw_assert(index < 3);
319		✗	return (&i)[index];
320			}
321			};
322
323			/**
324			* \brief Creates a triangle from its three vertices.
325			* \param[in] i , j , k the local indices of the three
326			* vertices.
327			* \return The created triangle.
328			*/
329			inline Triangle make_triangle(
330			ushort i, ushort j, ushort k
331			) {
332			Triangle result;
333			result.i = i;
334			result.j = j;
335			result.k = k;
336			return result;
337			}
338
339			/**
340			* \brief A triangle with flags.
341			* \details The flags are used for two purposes:
342			* - bits [0..15] are used to chain the triangles:
343			* there are two lists of triangles, the valid triangles
344			* and the free list. End of list is indicated by value 32767.
345			* - bit 16 (32768) is set if the triangle is in conflict.
346			*/
347			struct TriangleWithFlags : public Triangle {
348			ushort flags;
349			};
350
351			inline TriangleWithFlags make_triangle_with_flags(
352			ushort i, ushort j, ushort k, ushort f
353			) {
354			TriangleWithFlags result;
355			result.i = i;
356			result.j = j;
357			result.k = k;
358			result.flags = f;
359			return result;
360			}
361
362
363			/******************************************************************************/
364
365			inline double det2x2(
366			double a11, double a12,
367			double a21, double a22
368			) {
369		✗	return a11*a22 - a12*a21;
370			}
371
372			inline double det3x3(
373			double a11, double a12, double a13,
374			double a21, double a22, double a23,
375			double a31, double a32, double a33
376			) {
377			return
378		✗	a11*det2x2(a22,a23,a32,a33)
379		✗	-a21*det2x2(a12,a13,a32,a33)
380		✗	+a31*det2x2(a12,a13,a22,a23);
381			}
382
383		✗	inline double det4x4(
384			double a11, double a12, double a13, double a14,
385			double a21, double a22, double a23, double a24,
386			double a31, double a32, double a33, double a34,
387			double a41, double a42, double a43, double a44
388			) {
389		✗	double m12 = a21*a12 - a11*a22;
390		✗	double m13 = a31*a12 - a11*a32;
391		✗	double m14 = a41*a12 - a11*a42;
392		✗	double m23 = a31*a22 - a21*a32;
393		✗	double m24 = a41*a22 - a21*a42;
394		✗	double m34 = a41*a32 - a31*a42;
395
396		✗	double m123 = m23*a13 - m13*a23 + m12*a33;
397		✗	double m124 = m24*a13 - m14*a23 + m12*a43;
398		✗	double m134 = m34*a13 - m14*a33 + m13*a43;
399		✗	double m234 = m34*a23 - m24*a33 + m23*a43;
400
401		✗	return (m234*a14 - m134*a24 + m124*a34 - m123*a44);
402			}
403
404			/******************************************************************************/
405
406			enum ConvexCellFlag {
407			None = 0, /**< \brief default */
408			WithVGlobal = 1, /**< \brief store global vertex indices */
409			WithTFlags = 2 /**< \brief store user triange flags */
410			};
411
412			typedef index_t ConvexCellFlags;
413
414			/**
415			* \brief Computes the intersection between a set of halfplanes using
416			* Bowyer-Watson algorithm.
417			* \details Do not use with a large number of planes.
418			*/
419			class GEOGRAM_API ConvexCell {
420			public:
421
422			/**
423			* \brief ConvexCell constructor.
424			* \param[in] flags a combination of WithVGlobal, WithTFlags
425			*/
426			ConvexCell(ConvexCellFlags flags = None);
427
428			#ifndef STANDALONE_CONVEX_CELL
429			/**
430			* \brief Specifies whether exact predicates should be used.
431			* \param[in] x true if exact predicates should be used.
432			* \details Not supported if ConvexCell distributed
433			* as standalone file.
434			*/
435			void use_exact_predicates(bool x) {
436		✗	use_exact_predicates_ = x;
437			}
438			#endif
439
440			/**
441			* \brief Tests whether global vertex indices are stored.
442			* \retval true if global vertex indices are stored.
443			* \retval false otherwise.
444			*/
445			bool has_vglobal() const {
446		✗	return has_vglobal_;
447			}
448
449			/**
450			* \brief Tests whether triangle flags are stored.
451			* \retval true if triangle flags are stored.
452			* \retval false otherwise.
453			*/
454			bool has_tflags() const {
455			return has_tflags_;
456			}
457
458			/**
459			* \brief Creates vertex global indices if they are
460			* not present.
461			*/
462			void create_vglobal() {
463		✗	if(!has_vglobal()) {
464		✗	has_vglobal_ = true;
465		✗	vglobal_.assign(max_v(), global_index_t(-1));
466			}
467			}
468
469			/**
470			* \brief Removes all vertices and triangles from this
471			* ConvexCell.
472			* \details Keeps allocated memory for future use.
473			*/
474			void clear();
475
476			/**
477			* \brief Initializes this ConvexCell to an axis-aligned
478			* box.
479			* \details Previous contents of this ConvexCell are
480			* discarded. Vertex 0 is vertex at infinity.
481			* \param[in] xmin , ymin , zmin , xmax , ymax , zmax
482			* the coordinates of the box.
483			*/
484			void init_with_box(
485			double xmin, double ymin, double zmin,
486			double xmax, double ymax, double zmax
487			);
488
489			/**
490			* \brief Initializes this ConvexCell to a tetrahedron.
491			* \details Previous contents of this ConvexCell are
492			* discarded. Vertex 0 is vertex at infinity.
493			* \param[in] P0 , P1 , P2 , P3 the plane equations of
494			* the four faces of the tetrahedron.
495			*/
496			void init_with_tet(
497			vec4 P0, vec4 P1, vec4 P2, vec4 P3
498			);
499
500			/**
501			* \brief Initializes this ConvexCell to a tetrahedron.
502			* \details Previous contents of this ConvexCell are
503			* discarded. Vertex 0 is vertex at infinity.
504			* \param[in] P0 , P1 , P2 , P3 the plane equations of
505			* the four faces of the tetrahedron.
506			* \param[in] P0_global_index , P1_global_index ,
507			* P1_global_index , P2_global_index the global
508			* indices associated with the plane equations.
509			* \pre has_vglobal()
510			*/
511			void init_with_tet(
512			vec4 P0, vec4 P1, vec4 P2, vec4 P3,
513			global_index_t P0_global_index,
514			global_index_t P1_global_index,
515			global_index_t P2_global_index,
516			global_index_t P3_global_index
517			);
518
519			/**
520			* \brief Saves the computed cell in alias wavefront
521			* file format.
522			* \param[in] filename the name of the file where to
523			* save the cell.
524			* \param[in] shrink shrinking factor to ease visualization.
525			*/
526			void save(const std::string& filename, double shrink=0.0) const;
527
528
529			/**
530			* \brief Saves the computed cell in alias wavefront
531			* file format.
532			* \param[out] out a stream where to save the output.
533			* \param[in] v_offset offset applied to vertex indices.
534			* \param[in] shrink shrinking factor to ease visualization.
535			* \param[in] borders_only if set, only facets that correspond
536			* to vertex global index -1 are saved.
537			* \return the number of created vertices.
538			*/
539			index_t save(
540			std::ostream& out, global_index_t v_offset=1, double shrink=0.0,
541			bool borders_only=false
542			) const;
543
544			#if !defined(STANDALONE_CONVEX_CELL) && !defined(GEOGRAM_PSM)
545			/**
546			* \brief Appends the computed cell to a GEO::Mesh.
547			* \param[out] mesh a pointer to the mesh.
548			* \param[in] shrink shrinking factor to ease visualization.
549			* \param[in] borders_only if set, only facets that correspond
550			* to vertex global index -1 are saved.
551			* \param[in] facet_attr optional facet attribute that stores
552			* global facet (dual vertex) ids.
553			* \note One needs to call mesh->facets.connect() afterwards to
554			* have facets adjacencies. It is not called because one may
555			* want to append multiple cells to the same mesh.
556			*/
557			void append_to_mesh(
558			GEO::Mesh* mesh,
559			double shrink=0.0, bool borders_only=false,
560			GEO::Attribute<GEO::index_t>* facet_attr=nullptr
561			) const;
562
563			#endif
564
565			/**
566			* \brief Calls a user-defined function for each vertex of a Voronoi
567			* facet.
568			* \details One needs to call compute_geometry() before calling this
569			* function.
570			* \param[in] v the index of the (dual) Voronoi Facet, that is a
571			* (primal) vertex, in [0..nb_v()-1]
572			* \param[in] vertex a function that takes an index_t as an argument,
573			* with the index of the triangle that corresponds to the current
574			* Voronoi vertex.
575			*/
576			void for_each_Voronoi_vertex(
577			index_t v,
578			std::function<void(index_t)> vertex
579			);
580
581			/**
582			* \brief Clips this convex cell by a new plane.
583			* \details The positive side of the plane equation corresponds to
584			* what is kept. In other words, the normal vector P.x, P.y, P.z
585			* points towards the interior of this ConvexCell.
586			* \param[in] P the plane equation.
587			*/
588			void clip_by_plane(vec4 P);
589
590			/**
591			* \brief Clips this convex cell by a new plane and stores
592			* the corresponding global index in the newly created vertex.
593			* \details The positive side of the plane equation corresponds to
594			* what is kept. In other words, the normal vector P.x, P.y, P.z
595			* points towards the interior of this ConvexCell.
596			* This function can only be called if global indices are stored.
597			* \param[in] P the plane equation.
598			* \param[in] j the global index of the plane.
599			*/
600			void clip_by_plane(vec4 P, global_index_t j);
601
602
603			/**
604			* \brief Clips this convex cell by a new plane, using a user-defined
605			* geometric predicate.
606			* \details It is useful to be able to have a user-defined geometric
607			* predicates when the vertices have a symbolic representation, stored
608			* in the global indices associated with the plane. It is used by
609			* the robust mesh boolean operations.
610			* The positive side of the plane equation corresponds to
611			* what is kept. In other words, the normal vector P.x, P.y, P.z
612			* points towards the interior of this ConvexCell.
613			* If global indices are stored, then j is stored as the global index
614			* of the plane equation.
615			* \param[in] P the plane equation.
616			* \param[in] P_global_index the global index of the plane.
617			* \param[in] triangle_conflict_predicate a function that takes as
618			* arguments a local triangle index and local vertex (plane eqn)
619			* index, and that returns true if the triangle is in conflict with
620			* the vertex.
621			*/
622			void clip_by_plane(
623			vec4 P, global_index_t P_global_index,
624			std::function<bool(ushort,ushort)> triangle_conflict_predicate
625			);
626
627			/**
628			* \brief Clips this convex cell by a new plane and stores
629			* the corresponding global index in the newly created vertex.
630			* \details For a ConvexCell with a large number of facets, this
631			* version is faster than clip_by_plane(). However, it cannot be
632			* used with a ConvexCell that has infinite faces.
633			* \param[in] P the plane equation.
634			* \see clip_by_plane()
635			*/
636			void clip_by_plane_fast(vec4 P);
637
638			/**
639			* \brief Clips this convex cell by a new plane and stores
640			* the corresponding global index in the newly created vertex.
641			* \details For a ConvexCell with a large number of facets, this
642			* version is faster than clip_by_plane(). However, it cannot be
643			* used with a ConvexCell that has infinite faces.
644			* \param[in] P the plane equation.
645			* \param[in] j the global index of the plane.
646			* \see clip_by_plane()
647			*/
648			void clip_by_plane_fast(vec4 P, global_index_t j);
649
650			/**
651			* \brief Gets the number of triangles.
652			* \return the number of created triangles.
653			* \details The created triangles are not
654			* necessarily valid ones. To get the valid triangles,
655			* one needs to traverse the list from first_valid_.
656			*/
657			index_t nb_t() const {
658		✗	return nb_t_;
659			}
660
661			/**
662			* \brief Gets the number of vertices.
663			* \return the number of vertices.
664			* \details Some vertices can be incident to no triangle.
665			* The first six vertices correspond to the facets of the
666			* initial axis aligned box passed to the constructor.
667			*/
668			index_t nb_v() const {
669		✗	return nb_v_;
670			}
671
672			/**
673			* \brief Directly creates a new vertex.
674			* \param[in] P the plane equation attached to the vertex.
675			* \return the index of the newly created vertex.
676			*/
677		✗	index_t create_vertex(vec4 P) {
678		✗	if(nb_v_ == max_v_) {
679		✗	grow_v();
680			}
681		✗	plane_eqn_[nb_v_] = P;
682			index_t result = nb_v_;
683		✗	++nb_v_;
684		✗	return result;
685			}
686
687			/**
688			* \brief Directly creates a new vertex.
689			* \param[in] P the plane equation attached to the vertex.
690			* \param[in] v the global index associated with the vertex.
691			* \return the index of the newly created vertex.
692			* \pre global vertex indices are stored
693			*/
694			index_t create_vertex(vec4 P, global_index_t v) {
695		✗	index_t result = create_vertex(P);
696		✗	vglobal_[nb_v()-1] = v;
697			return result;
698			}
699
700			/**
701			* \brief Directly creates a new triangle.
702			* \param[in] i , j, k the three vertices of the
703			* triangle.
704			* \details The triangle is inserted into the list
705			* of valid triangles.
706			* \return the index of the newly created triangle.
707			*/
708			index_t create_triangle(index_t i, index_t j, index_t k) {
709			vbw_assert(i < nb_v());
710			vbw_assert(j < nb_v());
711			vbw_assert(k < nb_v());
712		✗	return new_triangle(i,j,k);
713			}
714
715			/**
716			* \brief Replaces a vertex with the vertex at infinity
717			* in all facets.
718			* \param[in] v the vertex to be killed.
719			*/
720			void kill_vertex(index_t v);
721
722			/**
723			* \brief Tests whether a vertex has a corresponding
724			* facet in the cell.
725			* \details Calling compute_geometry() before makes
726			* this function faster.
727			*/
728			bool vertex_is_contributing(index_t v) const {
729			if(!geometry_dirty_) {
730			return v2t_[v] != END_OF_LIST;
731			}
732			index_t t = first_valid_;
733			while(t != END_OF_LIST) {
734			TriangleWithFlags T = get_triangle_and_flags(t);
735			if(T.i == v || T.j == v || T.k == v) {
736			return true;
737			}
738			t = index_t(T.flags);
739			}
740			return false;
741			}
742
743			/**
744			* \brief Gets a triangle incident to a vertex.
745			* \param[in] v vertex index.
746			* \return a triangle incident to v.
747			*/
748			index_t vertex_triangle(index_t v) const {
749			geo_assert(!geometry_dirty_);
750			return v2t_[v];
751			}
752
753			/**
754			* \brief Computes the geometry and some cached information.
755			* \details Needs to be called before volume(),
756			* facet_area() and barycenter().
757			*/
758			void compute_geometry();
759
760			/**
761			* \brief Gets the dual facet area of a given vertex.
762			* \details compute_geometry() needs to be called before.
763			* \param[in] v the vertex.
764			* \return the dual facet area associated with v.
765			* \details terminate() needs to be called before
766			* calling this function.
767			*/
768			double facet_area(index_t v) const;
769
770			/**
771			* \brief Computes the volume of this convex cell.
772			* \details compute_geometry() needs to be called before.
773			* \return the volume.
774			*/
775			double volume() const;
776
777			/**
778			* \brief Computes the barycenter of this convex cell.
779			* \details compute_geometry() needs to be called before.
780			* \return the barycenter.
781			*/
782			vec3 barycenter() const;
783
784			/**
785			* \brief Computes volume and barycenter.
786			* \param[out] m the computed volume
787			* \param[out] mg the computed volume times the barycenter
788			* \details compute_geometry() needs to be called before.
789			*/
790			void compute_mg(double& m, vec3& mg) const ;
791
792
793			/**
794			* \brief Computes the squared radius of the smallest sphere
795			* containing the cell and centered on a point.
796			* \return the maximum squared distance between center and
797			* all the vertices of the cell.
798			*/
799			double squared_radius(vec3 center) const;
800
801			/**
802			* \brief Computes the squared radius of the largest sphere contained
803			* in the cell and centered on a point.
804			* \return the minimum squared distance between center and
805			* all facets of the cell.
806			*/
807			double squared_inner_radius(vec3 center) const;
808
809
810			/**
811			* \brief Tests whether this ConvexCell is empty.
812			* \details ConvexCell can be empty if everything was
813			* clipped out.
814			* \retval true if this ConvexCell is empty.
815			* \retval false otherwise.
816			*/
817			bool empty() const {
818			return first_valid_ == END_OF_LIST;
819			}
820
821			/**
822			* \brief Gets the global vertex index from a local
823			* vertex index.
824			* \details Vertex indices correspond to planes (remember,
825			* we are in dual form).
826			* \param[in] lv the local vertex index
827			* \return the global vertex index that corresponds to
828			* lv.
829			*/
830			global_index_t v_global_index(index_t lv) const {
831			vbw_assert(has_vglobal_);
832			vbw_assert(lv < nb_v());
833		✗	return vglobal_[lv];
834			}
835
836			/**
837			* \brief Sets the global vertex index associated with a local
838			* vertex index.
839			* \details Vertex indices correspond to planes (remember,
840			* we are in dual form).
841			* \param[in] lv the local vertex index
842			* \param[in] v the global vertex index that corresponds to
843			* lv.
844			*/
845			void set_v_global_index(index_t lv, global_index_t v) {
846			vbw_assert(has_vglobal_);
847			vbw_assert(lv < nb_v());
848			vglobal_[lv] = v;
849			}
850
851			/**
852			* \brief Tests whether a vertex with a given global index
853			* exists in this ConvexCell.
854			* \param[in] v the global index.
855			* \retval true if there exists in this ConvexCell a vertex with
856			* global index \p v.
857			* \retval false otherwise.
858			*/
859			bool has_v_global_index(global_index_t v) const;
860
861			/**
862			* \brief Gets the first triangle.
863			* \return the index of the first triangle, or END_OF_LIST
864			* if this ConvexCell is empty.
865			*/
866			ushort first_triangle() const {
867		✗	return ushort(first_valid_);
868			}
869
870			/**
871			* \brief Gets the successor of a triangle.
872			* \param[in] t the index of a valid triangle.
873			* \return the index of the successor of \p t, or END_OF_LIST
874			* if \p t is the last triangle.
875			*/
876			ushort next_triangle(ushort t) const {
877			return get_triangle_flags(t);
878			}
879
880			/**
881			* \brief Gets the point that corresponds to a triangle.
882			* \details If compute_geometry() was called, this gets
883			* the previously computed point, else it is computed
884			* and returned.
885			* \param[in] t the index of the triangle.
886			* \return the point that corresponds to triangle \p t.
887			*/
888			vec3 triangle_point(ushort t) const {
889			if(geometry_dirty_) {
890			vec4 result = compute_triangle_point(t);
891			vbw_assert(result.w != 0.0);
892			return make_vec3(
893			result.x/result.w, result.y/result.w, result.z/result.w
894			);
895			}
896			return triangle_point_[t];
897			}
898
899			/**
900			* \brief Gets the global index of a triangle vertex.
901			* \param[in] t the triangle.
902			* \param[in] llv one of 0,1,2.
903			* \return the global index of the vertex.
904			* \pre global indices are stored.
905			*/
906			global_index_t triangle_v_global_index(ushort t, index_t llv) const {
907			Triangle T = get_triangle(t);
908			ushort lv = ushort((llv==0)*T.i + (llv==1)*T.j + (llv==2)*T.k);
909			return v_global_index(lv);
910			}
911
912			/**
913			* \brief Gets the local index of a triangle vertex.
914			* \param[in] t the triangle.
915			* \param[in] llv one of 0,1,2.
916			* \return the local index of the vertex, in 0..nb_v()-1
917			*/
918			index_t triangle_v_local_index(ushort t, index_t llv) const {
919			Triangle T = get_triangle(t);
920			return index_t((llv==0)*T.i + (llv==1)*T.j + (llv==2)*T.k);
921			}
922
923			/**
924			* \brief Tests whether a triangle is marked by the user.
925			* \param[in] t the triangle.
926			* \retval true if the triangle is marked.
927			* \retval false otherwise.
928			* \pre triangle flags are stored.
929			*/
930			bool triangle_is_user_marked(ushort t) {
931			vbw_assert(has_tflags_);
932			vbw_assert(t < max_t_);
933			return (tflags_[t] != 0);
934			}
935
936			/**
937			* \brief Sets the user mark on a triangle.
938			* \param[in] t the triangle.
939			* \pre triangle flags are stored.
940			*/
941			void triangle_user_mark(ushort t) {
942			vbw_assert(has_tflags_);
943			vbw_assert(t < max_t_);
944			tflags_[t] = 1;
945			}
946
947			/**
948			* \brief Resets the user mark on a triangle.
949			* \param[in] t the triangle.
950			* \pre triangle flags are stored.
951			*/
952			void triangle_user_unmark(ushort t) {
953			vbw_assert(has_tflags_);
954			vbw_assert(t < max_t_);
955			tflags_[t] = 0;
956			}
957
958			/**
959			* \brief Tests whether a cell has at least one vertex in conflict with
960			* a halfspace.
961			* \param[in] P the equation of the halfspace.
962			* \retval true if there exists a triangle t such that
963			* triangle_is_in_conflict(P)
964			* \retval false otherwise.
965			*/
966			bool cell_has_conflict(const vec4& P) {
967			for(
968			ushort t = first_triangle();
969			t!=END_OF_LIST; t=next_triangle(t)
970			) {
971			TriangleWithFlags T = get_triangle_and_flags(t);
972			if(triangle_is_in_conflict(T,P)) {
973			return true;
974			}
975			}
976			return false;
977			}
978
979			/**
980			* \brief Tests whether a cell has all its vertices in conflict
981			* with a plane.
982			* \param[in] P the equation of the halfspace.
983			* \retval true if all the triangles are in conflict with P.
984			* \retval false otherwise.
985			*/
986			bool cell_is_totally_in_conflict(const vec4& P) {
987			for(
988			ushort t = first_triangle();
989			t!=END_OF_LIST; t=next_triangle(t)
990			) {
991			TriangleWithFlags T = get_triangle_and_flags(t);
992			if(!triangle_is_in_conflict(T,P)) {
993			return false;
994			}
995			}
996			return true;
997			}
998
999			/**
1000			* \brief Gets a triangle adjacent to another triangle by edge
1001			* local index.
1002			* \param[in] t a triangle.
1003			* \param[in] le local index of an edge of \p t (in 0..2).
1004			* \return the triangle adjacent to \p t along \ p e.
1005			*/
1006			index_t triangle_adjacent(index_t t, index_t le) const {
1007			vbw_assert(t < max_t());
1008			vbw_assert(le < 3);
1009		✗	return t_adj_[t][le];
1010			}
1011
1012
1013			/**
1014			* \brief Sets triangle to triangle adjacency.
1015			* \param[in] t1 a triangle.
1016			* \param[in] le local index of an edge of \p t (in 0..2).
1017			* \param[in] t2 triangle to be made adjacent to \p t1 along edge \p le.
1018			*/
1019			void set_triangle_adjacent(index_t t1, index_t le, index_t t2) {
1020			vbw_assert(t1 < max_t());
1021			vbw_assert(le < 3);
1022			vbw_assert(t2 < max_t());
1023		✗	t_adj_[t1][le] = VBW::ushort(t2);
1024		✗	}
1025
1026
1027
1028			/**
1029			* \brief Gets a triangle vertex.
1030			* \param[in] t a triangle.
1031			* \param[in] lv local index of a vertex of \p t (in 0..2).
1032			* \return the vertex
1033			*/
1034			index_t triangle_vertex(index_t t, index_t lv) const {
1035			vbw_assert(t < max_t());
1036			vbw_assert(lv < 3);
1037		✗	return t_[t][lv];
1038			}
1039
1040
1041			/**
1042			* \brief Gets the local index of a vertex in a triangle.
1043			* \param[in] t a triangle.
1044			* \param[in] v a vertex index.
1045			* \return the local index of \p v in \p t (in 0..2).
1046			*/
1047			index_t triangle_find_vertex(index_t t, index_t v) const {
1048			vbw_assert(t < max_t());
1049			Triangle T = get_triangle(t);
1050		✗	index_t result = index_t((T.j == v) + 2*(T.k == v));
1051			vbw_assert(triangle_vertex(t,result) == v);
1052			return result;
1053			}
1054
1055			/**
1056			* \brief Gets the edge on witch a triangle is adjacent to another one
1057			* \param[in] t1 a triangle.
1058			* \param[in] t2 a triangle adjacent to t1
1059			* \return the edge index e such that triangle_adjacent(t1,e)=t2
1060			*/
1061			index_t triangle_find_adjacent(index_t t1, index_t t2) const {
1062			vbw_assert(t1 < max_t());
1063			vbw_assert(t2 < max_t());
1064			Triangle T = t_adj_[t1];
1065			index_t result = index_t((T.j == t2) + 2*(T.k == t2));
1066			vbw_assert(triangle_adjacent(t1,result) == t2);
1067			return result;
1068			}
1069
1070			/**
1071			* \brief Tests whether a triangle is infinite.
1072			* \param[in] t the triangle
1073			* \retval true if t is incident to the vertex at
1074			* infinity.
1075			* \retval false otherwise.
1076			*/
1077			bool triangle_is_infinite(index_t t) const {
1078			vbw_assert(t < max_t());
1079			Triangle T = get_triangle(t);
1080			return (
1081			T.i == VERTEX_AT_INFINITY ||
1082			T.j == VERTEX_AT_INFINITY ||
1083			T.k == VERTEX_AT_INFINITY
1084			);
1085			}
1086
1087			/**
1088			* \brief Gets the equation of a plane associated with a vertex.
1089			* \details The first six equations correspond to the six
1090			* facets of a cube.
1091			* \param[in] v the local index of the vertex.
1092			*/
1093			vec4 vertex_plane(index_t v) const {
1094			vbw_assert(v < max_v());
1095		✗	return plane_eqn_[v];
1096			}
1097
1098			/**
1099			* \brief Gets the normal to the plane associated with a vertex.
1100			* \details The first six equations correspond to the six
1101			* facets of a cube.
1102			* \param[in] v the local index of the vertex.
1103			*/
1104			vec3 vertex_plane_normal(index_t v) const {
1105			vbw_assert(v != VERTEX_AT_INFINITY);
1106			vbw_assert(v < max_v());
1107			return make_vec3(
1108		✗	plane_eqn_[v].x,
1109		✗	plane_eqn_[v].y,
1110		✗	plane_eqn_[v].z
1111			);
1112			}
1113
1114			/**
1115			* \brief Tests whether a triangle is marked as conflict.
1116			* \param[in] t a triangle.
1117			* \retval true if \p t is marked as conflict.
1118			* \retval false otherwise.
1119			*/
1120			bool triangle_is_marked_as_conflict(index_t t) const {
1121			vbw_assert(t < max_t());
1122			return (get_triangle_flags(t) & ushort(CONFLICT_MASK)) != 0;
1123			}
1124
1125			/**
1126			* \brief Tests whether a triangle is in conflict with a plane.
1127			* \details A triangle is in conflict with a plane if feeding the point
1128			* associated with the triangle in the equation of the plane yields
1129			* a negative number.
1130			* \param[in] T a triangle.
1131			* \param[in] eqn the four coefficients of the equation of the plane.
1132			* \retval true if \p t is in conflict with \p eqn.
1133			* \retval false otherwise.
1134			*/
1135			bool triangle_is_in_conflict(
1136			TriangleWithFlags T, const vec4& eqn
1137			) const;
1138
1139			/**
1140			* \brief Creates a new triangle.
1141			* \param[in] i , j , k the three vertices of the triangle.
1142			*/
1143		✗	index_t new_triangle(index_t i, index_t j, index_t k) {
1144		✗	index_t result = first_free_;
1145		✗	if(result == END_OF_LIST) {
1146		✗	result = nb_t_;
1147		✗	++nb_t_;
1148		✗	if(nb_t_ > max_t()) {
1149		✗	grow_t();
1150			}
1151			} else {
1152		✗	first_free_ = index_t(
1153			get_triangle_flags(first_free_) & ~ushort(CONFLICT_MASK)
1154			);
1155			}
1156			vbw_assert(result < max_t());
1157		✗	t_[result] = make_triangle_with_flags(
1158		✗	ushort(i), ushort(j), ushort(k), ushort(first_valid_)
1159			);
1160		✗	first_valid_ = result;
1161		✗	if(has_tflags_) {
1162		✗	tflags_[result] = 0;
1163			}
1164		✗	return result;
1165			}
1166
1167			/**
1168			* \brief Creates a new triangle.
1169			* \details Adjacency information is not used (kept for reference).
1170			* \param[in] i , j , k the three vertices of the triangle.
1171			* \param[in] adj0 , adj1 , adj2 the three adjacent triangles
1172			* (unused in this version).
1173			* \return the index of the new triangle.
1174			*/
1175			index_t new_triangle(
1176			index_t i, index_t j, index_t k,
1177			index_t adj0, index_t adj1, index_t adj2
1178			) {
1179		✗	index_t result = new_triangle(i, j, k);
1180		✗	t_adj_[result] = make_triangle(
1181			ushort(adj0), ushort(adj1), ushort(adj2)
1182			);
1183			return result;
1184			}
1185
1186			/**
1187			* \brief Computes the coordinates of the point
1188			* associated with a triangle.
1189			* \param[in] t the triangle.
1190			* \return the intersection between the three planes
1191			* associated with the three vertices of the triangle,
1192			* in homogeneous coordinates.
1193			*/
1194			vec4 compute_triangle_point(index_t t) const;
1195
1196			/**
1197			* \brief Gets the three vertices of a triangle.
1198			* \param[in] t the triangle.
1199			* \return a Triangle with the indices of the three vertices
1200			* of the triangle.
1201			*/
1202			Triangle get_triangle(index_t t) const {
1203			vbw_assert(t < max_t());
1204		✗	return t_[t];
1205			}
1206
1207			/**
1208			* \brief Gets the flags associated with a triangle.
1209			* \details Contains both the conflict flag and the
1210			* chaining.
1211			* \param[in] t the triangle.
1212			* \return the flags associated with \p t.
1213			*/
1214			ushort get_triangle_flags(index_t t) const {
1215			vbw_assert(t < max_t());
1216		✗	return t_[t].flags;
1217			}
1218
1219			/**
1220			* \brief Sets the flags of a triangle.
1221			* \param[in] t the triangle.
1222			* \param[in] flags the flags to be set.
1223			*/
1224			void set_triangle_flags(index_t t, ushort flags) {
1225			vbw_assert(t < max_t());
1226		✗	t_[t].flags = flags;
1227		✗	}
1228
1229			/**
1230			* \brief Gets the three vertices of a triangle and its flags.
1231			* \param[in] t the triangle.
1232			* \return a TriangleWithFlags with the indices of the three vertices
1233			* of the triangle and the flags.
1234			*/
1235			TriangleWithFlags get_triangle_and_flags(index_t t) const {
1236			vbw_assert(t < max_t());
1237		✗	return t_[t];
1238			}
1239
1240			/**
1241			* \brief Tests whether a given triangle is in the conflict zone.
1242			*/
1243			bool triangle_is_marked_as_conflict(index_t t) {
1244			vbw_assert(t < max_t());
1245			ushort flg = get_triangle_flags(t);
1246		✗	return ((flg & ushort(CONFLICT_MASK)) != 0);
1247			}
1248
1249			/**
1250			* \brief Gets the maximum valid index for a triangle.
1251			* \return the maximum valid index of a triangle.
1252			*/
1253			index_t max_t() const {
1254		✗	return max_t_;
1255			}
1256
1257			/**
1258			* \brief Gets the maximum valid index for a vertex.
1259			* \return the maximum valid index of a vertex.
1260			*/
1261			index_t max_v() const {
1262		✗	return max_v_;
1263			}
1264
1265			/**
1266			* \brief Allocates more space for triangles.
1267			* \details Makes max_t_ twice bigger.
1268			*/
1269			void grow_t();
1270
1271			/**
1272			* \brief Allocates more space for vertices.
1273			* \details Makes max_v_ twice bigger.
1274			*/
1275			void grow_v();
1276
1277
1278			/**
1279			* \brief Swaps two ConvexCells.
1280			* \param[in] other the ConvexCell to be
1281			* exchanged with this ConvexCell.
1282			*/
1283			void swap(ConvexCell& other) {
1284			std::swap(max_t_,other.max_t_);
1285			std::swap(max_v_,other.max_v_);
1286			std::swap(t_,other.t_);
1287			std::swap(t_adj_,other.t_adj_);
1288			std::swap(plane_eqn_,other.plane_eqn_);
1289			std::swap(nb_t_,other.nb_t_);
1290			std::swap(nb_v_,other.nb_v_);
1291			std::swap(first_free_,other.first_free_);
1292			std::swap(first_valid_,other.first_valid_);
1293			std::swap(geometry_dirty_,other.geometry_dirty_);
1294			std::swap(triangle_point_,other.triangle_point_);
1295			std::swap(v2t_,other.v2t_);
1296			std::swap(v2e_,other.v2e_);
1297			std::swap(vglobal_,other.vglobal_);
1298			std::swap(has_vglobal_,other.has_vglobal_);
1299			std::swap(tflags_,other.tflags_);
1300			std::swap(has_tflags_,other.has_tflags_);
1301			#ifndef STANDALONE_CONVEX_CELL
1302			std::swap(use_exact_predicates_,other.use_exact_predicates_);
1303			#endif
1304			}
1305
1306			/**
1307			* \brief Gets a modifiable reference to a triangle point.
1308			* \param[in] t the index
1309			* \return a modifiable reference to the stored point
1310			*/
1311			vec3& stored_triangle_point(ushort t) {
1312			return triangle_point_[t];
1313			}
1314
1315			protected:
1316
1317			/**
1318			* \brief finds all triangle-triangle adjacency relations.
1319			* \details Client code should not need to call this function. It is used
1320			* by PeriodicDelaunay3d::copy_Laguerre_cell_from_Delaunay().
1321			*/
1322			void connect_triangles();
1323
1324
1325			/**
1326			* \brief Triangulates the conflict zone.
1327			* \param[in] lv the local index of the new vertex
1328			* \param[in] conflict_head , conflict tail the first
1329			* and last triangle of the conflict zone stored
1330			* as a linked list.
1331			*/
1332			void triangulate_conflict_zone(
1333			index_t lv, index_t conflict_head, index_t conflict_tail
1334			);
1335
1336			/**
1337			* \brief Changes a vertex plane equation.
1338			* \param[in] v the vertex.
1339			* \param[in] P the plane equation.
1340			* \details Does not update combinatorics.
1341			* \note Use with care, for experts only.
1342			*/
1343			void set_vertex_plane(index_t v, vec4 P) {
1344			vbw_assert(v < max_v());
1345			plane_eqn_[v] = P;
1346			geometry_dirty_ = true;
1347			}
1348
1349
1350			private:
1351
1352			/** \brief number of allocated triangles */
1353			index_t max_t_;
1354
1355			/** \brief number of allocated vertices */
1356			index_t max_v_;
1357
1358			/** \brief indices of triangle vertices and flags */
1359			vector<TriangleWithFlags> t_;
1360
1361			/** \brief adjacency of each triangle */
1362			vector<Triangle> t_adj_;
1363
1364			/**
1365			* \brief plane equation attached to each vertex,
1366			* as specified by clip_by_plane().
1367			*/
1368			vector<vec4> plane_eqn_;
1369
1370			/** \brief number of used triangles. */
1371			index_t nb_t_;
1372
1373			/** \brief number of used vertices. */
1374			index_t nb_v_;
1375
1376			/** \brief Head of the linked list of free triangles. */
1377			index_t first_free_;
1378
1379			/** \brief Head of the linked list of valid triangles. */
1380			index_t first_valid_;
1381
1382			/**
1383			* \brief true if triangle_point_ and t2v_ are
1384			* not up to date.
1385			*/
1386			bool geometry_dirty_;
1387
1388			/**
1389			* \brief dual vertex attached to each triangle.
1390			*/
1391			vector<vec3> triangle_point_;
1392
1393			/**
1394			* \brief One triangle incident to each vertex,
1395			* or END_OF_LIST if there is no such triangle.
1396			* Used also to store linked list of vertices
1397			* around conflict zone.
1398			*/
1399			vector<ushort> v2t_;
1400
1401			/**
1402			* \brief Used by linked list of vertices around
1403			* conflict zone. Indicates which edge of
1404			* v2t_[v] is incident to the conflict zone.
1405			*/
1406			vector<uchar> v2e_;
1407
1408			/**
1409			* \brief Optional vector of gloval vertex indices.
1410			*/
1411			vector<global_index_t> vglobal_;
1412
1413			/**
1414			* \brief True if global vertex indices are stored.
1415			*/
1416			bool has_vglobal_;
1417
1418			/**
1419			* \brief Optional flags attached to the triangles.
1420			*/
1421			vector<uchar> tflags_;
1422
1423			/**
1424			* \brief True if triangle flags are stored.
1425			*/
1426			bool has_tflags_;
1427
1428			#ifndef STANDALONE_CONVEX_CELL
1429			/**
1430			* \brief True if exact predicates should be used.
1431			*/
1432			bool use_exact_predicates_;
1433			#endif
1434
1435			friend class GEO::PeriodicDelaunay3d;
1436			};
1437			}
1438
1439			namespace GEO {
1440			using VBW::ConvexCell;
1441			}
1442
1443			#endif
1444