GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* Copyright (c) 2000-2022 Inria
3			* All rights reserved.
4			*
5			* Redistribution and use in source and binary forms, with or without
6			* modification, are permitted provided that the following conditions are met:
7			*
8			* * Redistributions of source code must retain the above copyright notice,
9			* this list of conditions and the following disclaimer.
10			* * Redistributions in binary form must reproduce the above copyright notice,
11			* this list of conditions and the following disclaimer in the documentation
12			* and/or other materials provided with the distribution.
13			* * Neither the name of the ALICE Project-Team nor the names of its
14			* contributors may be used to endorse or promote products derived from this
15			* software without specific prior written permission.
16			*
17			* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18			* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20			* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21			* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22			* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23			* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24			* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25			* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26			* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27			* POSSIBILITY OF SUCH DAMAGE.
28			*
29			* Contact: Bruno Levy
30			*
31			* https://www.inria.fr/fr/bruno-levy
32			*
33			* Inria,
34			* Domaine de Voluceau,
35			* 78150 Le Chesnay - Rocquencourt
36			* FRANCE
37			*
38			*/
39
40			#include <exploragram/optimal_transport/optimal_transport_on_surface.h>
41			#include <geogram/voronoi/generic_RVD_vertex.h>
42			#include <geogram/voronoi/generic_RVD_polygon.h>
43			#include <geogram/voronoi/RVD_callback.h>
44			#include <geogram/basic/geometry.h>
45			#include <geogram/basic/geometry_nd.h>
46			#include <geogram/basic/stopwatch.h>
47
48			#include <geogram/bibliography/bibliography.h>
49
50			namespace {
51			using namespace GEO;
52
53			/**********************************************************************/
54
55			// For more details on the algorithm / structure of the objective function,
56			//see also implementation notes at the beginning of optimal_transport_3d.cpp.
57
58			/**
59			* \brief Computes the contribution of a polygon
60			* to the objective function minimized by a semi-discrete
61			* optimal transport map in 2D.
62			*/
63			class SurfaceOTMPolygonCallback :
64			public OptimalTransportMap::Callback,
65			public RVDPolygonCallback {
66			public:
67
68			/**
69			* \brief OTMPolygonCallback constructor.
70			* \param[in] OTM a pointer to the OptimalTransportMapOnSurface
71			*/
72		✗	SurfaceOTMPolygonCallback(OptimalTransportMapOnSurface* OTM) :
73		✗	OptimalTransportMap::Callback(OTM) {
74		✗	}
75
76			/**
77			* \copydoc RVDPolygonCallback::operator()
78			*/
79		✗	void operator() (
80			index_t v,
81			index_t t,
82			const GEOGen::Polygon& P
83			) const override {
84			// v can be an air particle.
85		✗	if(v >= n_) {
86		✗	return;
87			}
88
89		✗	if(P.nb_vertices() == 0) {
90			return;
91			}
92
93			double m, mgx, mgy, mgz;
94		✗	compute_m_and_mg(P, m, mgx, mgy, mgz);
95
96		✗	if(spinlocks_ != nullptr) {
97		✗	spinlocks_->acquire_spinlock(v);
98			}
99
100			// +m because we maximize F <=> minimize -F
101		✗	g_[v] += m;
102
103		✗	if(Newton_step_) {
104			// ... but here -m because Newton step =
105			// solve H p = -g (minus g in the RHS).
106		✗	OTM_->add_i_right_hand_side(v,-m);
107			}
108
109		✗	if(mg_ != nullptr) {
110		✗	mg_[3*v] += mgx;
111		✗	mg_[3*v+1] += mgy;
112		✗	mg_[3*v+2] += mgz;
113			}
114
115		✗	if(spinlocks_ != nullptr) {
116		✗	spinlocks_->release_spinlock(v);
117			}
118
119		✗	if(Newton_step_) {
120			// Spinlocks are managed internally by update_Hessian().
121		✗	update_Hessian(P, v, t);
122			}
123
124
125		✗	if(eval_F_) {
126		✗	Thread* thread = Thread::current();
127			index_t current_thread_id = (thread==nullptr) ? 0 : thread->id();
128		✗	double F = weighted_ ? eval_F_weighted(P, v) : eval_F(P, v);
129			const_cast<SurfaceOTMPolygonCallback*>(this)->
130			funcval_[current_thread_id] += F;
131			}
132			}
133
134			protected:
135
136			/**
137			* \brief Computes the mass and mass times centroid of the
138			* current intersection polygon.
139			* \details Weights are taken into account if present.
140			* \param[in] P a const reference to the current intersection polygon.
141			* \param[out] m , mgx , mgy , mgz the mass and the mass times the
142			* centroid of the ConvexCell. mgx, mgy and mgy are not computed
143			* if mg_ is nullptr.
144			*/
145		✗	void compute_m_and_mg(
146			const GEOGen::Polygon& P,
147			double& m, double& mgx, double& mgy, double& mgz
148			) const {
149		✗	m = 0.0;
150		✗	mgx = 0.0;
151		✗	mgy = 0.0;
152		✗	mgz = 0.0;
153			const GEOGen::Vertex& V0 = P.vertex(0);
154			const double* p0 = V0.point();
155		✗	for(index_t i=1; i+1<P.nb_vertices(); ++i) {
156			const GEOGen::Vertex& V1 = P.vertex(i);
157			const double* p1 = V1.point();
158			const GEOGen::Vertex& V2 = P.vertex(i+1);
159			const double* p2 = V2.point();
160		✗	double cur_m = triangle_mass(V0,V1,V2);
161		✗	m += cur_m;
162		✗	if(mg_ != nullptr) {
163		✗	if(weighted_) {
164			double w0 = V0.weight();
165			double w1 = V1.weight();
166			double w2 = V2.weight();
167		✗	double s = cur_m/(w0+w1+w2);
168		✗	mgx += s * (w0*p0[0] + w1*p1[0] + w2*p2[0]);
169		✗	mgy += s * (w0*p0[1] + w1*p1[1] + w2*p2[1]);
170		✗	mgz += s * (w0*p0[2] + w1*p1[2] + w2*p2[2]);
171			} else {
172		✗	mgx += cur_m * (p0[0] + p1[0] + p2[0]) / 3.0;
173		✗	mgy += cur_m * (p0[1] + p1[1] + p2[1]) / 3.0;
174		✗	mgz += cur_m * (p0[2] + p1[2] + p2[2]) / 3.0;
175			}
176			}
177			}
178		✗	}
179
180			/**
181			* \brief Gets the unit normal vector to a facet.
182			*/
183		✗	vec3 facet_normal(index_t t) const {
184		✗	Mesh& M = OTM_->mesh();
185			index_t i = M.facets.vertex(t,0);
186			index_t j = M.facets.vertex(t,1);
187			index_t k = M.facets.vertex(t,2);
188			const double* p0 = M.vertices.point_ptr(i);
189			const double* p1 = M.vertices.point_ptr(j);
190			const double* p2 = M.vertices.point_ptr(k);
191		✗	double x1 = p1[0] - p0[0];
192		✗	double y1 = p1[1] - p0[1];
193		✗	double z1 = p1[2] - p0[2];
194		✗	double x2 = p2[0] - p0[0];
195		✗	double y2 = p2[1] - p0[1];
196		✗	double z2 = p2[2] - p0[2];
197		✗	double x = y1*z2 - z1*y2;
198		✗	double y = z1*x2 - x1*z2;
199		✗	double z = x1*y2 - y1*x2;
200		✗	return vec3(x,y,z);
201			}
202
203			/**
204			* \brief Updates the Hessian according to the current intersection
205			* polygon.
206			* \param[in] P a const reference to the current intersection polygon.
207			* \param[in] i the current seed
208			* \param[in] t the current mesh facet
209			*/
210		✗	void update_Hessian(
211			const GEOGen::Polygon& P, index_t i, index_t t
212			) const {
213
214		✗	vec3 N = normalize(facet_normal(t));
215
216			// The coefficient of the Hessian associated to a pair of
217			// adjacent cells Lag(i),Lag(j) is :
218			// - mass(Lag(i) /\ Lag(j)) / (2*distance(pi,pj))
219
220		✗	const double* pi = OTM_->point_ptr(i);
221
222		✗	for(index_t k1=0; k1<P.nb_vertices(); ++k1) {
223		✗	index_t k2 = k1+1;
224		✗	if(k2 == P.nb_vertices()) {
225			k2 = 0;
226			}
227			// Note:
228			// It is P.vertex(k2).adjacent_seed(),
229			// not P.vertex(k1).adjacent_seed() !!!
230		✗	index_t j = index_t(P.vertex(k2).adjacent_seed());
231		✗	if(j != index_t(-1)) {
232		✗	const double* pj = OTM_->point_ptr(j);
233
234		✗	vec3 pij(pj[0] - pi[0],pj[1] - pi[1],pj[2] - pi[2]);
235			// Remove from pij the part that is normal
236			// to the current polygon.
237			pij -= dot(pij,N)*N;
238			double lij = length(pij);
239
240			double hij =
241		✗	edge_mass(P.vertex(k1), P.vertex(k2)) / (2.0 * lij);
242
243			// -hij because we maximize F <=> minimize -F
244		✗	if(hij != 0.0) {
245		✗	if(spinlocks_ != nullptr) {
246		✗	spinlocks_->acquire_spinlock(i);
247			}
248			// Diagonal is positive, extra-diagonal
249			// coefficients are negative,
250			// this is a convex function.
251		✗	if(j < n_) {
252		✗	OTM_->add_ij_coefficient(i, j, -hij);
253			}
254		✗	OTM_->add_ij_coefficient(i, i, hij);
255		✗	if(spinlocks_ != nullptr) {
256		✗	spinlocks_->release_spinlock(i);
257			}
258			}
259			}
260			}
261		✗	}
262
263			/**
264			* \brief Computes the contribution of the current polygon
265			* to the objective function, in the uniform (non-weighted)
266			* case.
267			* \param[in] P a const reference to the current polygon.
268			* \param[in] i the current seed
269			*/
270		✗	double eval_F(const GEOGen::Polygon& P, index_t i) const {
271			geo_debug_assert(!weighted_);
272			geo_argused(P);
273			geo_argused(i);
274			// Not implemented yet.
275		✗	geo_assert_not_reached;
276			}
277
278			/**
279			* \brief Computes the contribution of the current polygon
280			* to the objective function, in the weighted case.
281			* \param[in] P a const reference to the current polygon.
282			* \param[in] i the current seed
283			*/
284		✗	double eval_F_weighted(const GEOGen::Polygon& P, index_t i) const {
285			geo_argused(P);
286			geo_argused(i);
287			// Not implemented yet.
288		✗	geo_assert_not_reached;
289			}
290
291			/**
292			* \brief Computes the mass of a triangle.
293			* \details Weights are taken into account in weighted_ mode.
294			* \param[in] V0 , V1 , V2 the three vertices of the triangle,
295			* given as RVD vertices.
296			* \return the area of the triangle in 3D times the average value of the
297			* three weights.
298			*/
299		✗	double triangle_mass(
300			const GEOGen::Vertex& V0,
301			const GEOGen::Vertex& V1,
302			const GEOGen::Vertex& V2
303			) const {
304			double m =
305		✗	Geom::triangle_area_3d(V0.point(), V1.point(), V2.point());
306		✗	if(weighted_) {
307		✗	m *= ((V0.weight() + V1.weight() + V2.weight())/3.0);
308			}
309		✗	return m;
310			}
311
312			/**
313			* \brief Computes the mass of an edge.
314			* \details Weights are taken into account in weighted_ mode.
315			* \param[in] V0 , V1 the two vertices of the edge,
316			* given as RVD vertices.
317			* \return the length of the edge in 3D times the average value of the
318			* two weights.
319			*/
320		✗	double edge_mass(
321			const GEOGen::Vertex& V0,
322			const GEOGen::Vertex& V1
323			) const {
324		✗	double m = Geom::distance(V0.point(), V1.point(), 3);
325		✗	if(weighted_) {
326		✗	m *= ((V0.weight() + V1.weight())/2.0);
327			}
328		✗	return m;
329			}
330			};
331
332			/**********************************************************************/
333
334			/**
335			* \brief A RVDPolygonCallback that stores the Restricted Voronoi
336			* Diagram in a Mesh.
337			*/
338			class ComputeRVDPolygonCallback : public RVDPolygonCallback {
339			public:
340		✗	ComputeRVDPolygonCallback(OptimalTransportMap* OTM, Mesh* target) :
341		✗	OTM_(OTM), target_(target) {
342		✗	target_->clear();
343		✗	target_->vertices.set_dimension(3);
344		✗	chart_.bind(target_->facets.attributes(), "chart");
345		✗	}
346
347		✗	~ComputeRVDPolygonCallback() override {
348		✗	chart_.unbind();
349		✗	}
350
351		✗	void operator() (
352			index_t v,
353			index_t t,
354			const GEOGen::Polygon& P
355			) const override {
356			geo_argused(v);
357			geo_argused(t);
358
359		✗	if(P.nb_vertices() == 0) {
360			return;
361			}
362
363		✗	index_t voffset = target_->vertices.nb();
364		✗	FOR(i,P.nb_vertices()) {
365			const double* p = P.vertex(i).point();
366		✗	if(OTM_->dimension() == 2) {
367		✗	target_->vertices.create_vertex(
368		✗	vec3(p[0], p[1], 0.0).data()
369			);
370			} else {
371		✗	target_->vertices.create_vertex(
372		✗	vec3(p[0], p[1], p[2]).data()
373			);
374			}
375			}
376		✗	index_t f = target_->facets.create_polygon(P.nb_vertices());
377		✗	FOR(i,P.nb_vertices()) {
378		✗	target_->facets.set_vertex(f,i,voffset+i);
379			}
380		✗	const_cast<Attribute<index_t>&>(chart_)[f] = v;
381			}
382
383			private:
384			OptimalTransportMap* OTM_;
385			Mesh* target_;
386			Attribute<index_t> chart_;
387			};
388
389
390
391			/**
392			* \brief Gets the Delaunay implementation that best fits
393			* user desire.
394			* \param[in] user_delaunay the desired implementation.
395			* \return the available implementation nearest to user desire.
396			*/
397		✗	std::string default_delaunay(const std::string& user_delaunay) {
398		✗	std::string result = "BPOW";
399		✗	if(user_delaunay != "default") {
400			result = user_delaunay;
401			}
402			#ifndef GEOGRAM_WITH_PDEL
403		✗	if(result == "PDEL") {
404			result = "BPOW";
405			}
406			#endif
407		✗	return result;
408			}
409
410
411			}
412
413			/**********************************************************************/
414
415			namespace GEO {
416
417		✗	OptimalTransportMapOnSurface::OptimalTransportMapOnSurface(
418			Mesh* mesh, const std::string& delaunay, bool BRIO
419		✗	) :
420			OptimalTransportMap(
421			3,
422			mesh,
423		✗	default_delaunay(delaunay),
424			BRIO
425		✗	) {
426		✗	callback_ = new SurfaceOTMPolygonCallback(this);
427		✗	total_mass_ = total_mesh_mass();
428		✗	geo_cite("DBLP:journals/corr/MerigotMT17");
429		✗	}
430
431		✗	OptimalTransportMapOnSurface::~OptimalTransportMapOnSurface() {
432		✗	}
433
434		✗	void OptimalTransportMapOnSurface::get_RVD(Mesh& RVD_mesh) {
435		✗	ComputeRVDPolygonCallback callback(this, &RVD_mesh);
436		✗	RVD()->for_each_polygon(callback, false, false, false);
437			/*
438			// NOTE: Does not work, TODO: determine why
439			Attribute<index_t> tet_region(RVD_mesh.cells.attributes(),"region");
440			RVD()->compute_RVD(
441			RVD_mesh,
442			0, // dim (0 means use default)
443			false, // borders_only
444			show_RVD_seed_ // integration_simplices
445			);
446			*/
447		✗	}
448
449		✗	void OptimalTransportMapOnSurface::compute_Laguerre_centroids(
450			double* centroids
451			) {
452		✗	vector<double> g(nb_points(), 0.0);
453		✗	Memory::clear(centroids, nb_points()*sizeof(double)*3);
454
455		✗	callback_->set_Laguerre_centroids(centroids);
456			callback_->set_g(g.data());
457			{
458			Stopwatch* W = nullptr;
459		✗	if(newton_ && verbose_) {
460		✗	W = new Stopwatch("RVD");
461		✗	Logger::out("OTM") << "In RVD (centroids)..." << std::endl;
462			}
463		✗	RVD_->for_each_polygon(
464		✗	*dynamic_cast<RVDPolygonCallback*>(callback_), false, false, true
465			);
466		✗	if(newton_ && verbose_) {
467		✗	delete W;
468			}
469			}
470
471		✗	callback_->set_Laguerre_centroids(nullptr);
472
473		✗	for(index_t v=0; v<nb_points(); ++v) {
474		✗	centroids[3*v ] /= g[v];
475		✗	centroids[3*v+1] /= g[v];
476		✗	centroids[3*v+2] /= g[v];
477			}
478		✗	}
479
480		✗	double OptimalTransportMapOnSurface::total_mesh_mass() const {
481			double result = 0.0;
482
483			// This is terribly confusing, the parameters for
484			// a power diagram are called "weights", and the
485			// standard attribute name for vertices density is
486			// also called "weight" (and is unrelated).
487			// In this program, what is called weight corresponds
488			// to the parameters of the power diagram (except the
489			// name of the attribute), and everything that corresponds
490			// to mass/density is called mass.
491			Attribute<double> vertex_mass;
492		✗	vertex_mass.bind_if_is_defined(
493		✗	mesh_->vertices.attributes(), "weight"
494			);
495
496		✗	for(index_t t: mesh_->facets) {
497			double tri_mass = GEO::Geom::triangle_area(
498		✗	vec3(mesh_->vertices.point_ptr(mesh_->facets.vertex(t, 0))),
499		✗	vec3(mesh_->vertices.point_ptr(mesh_->facets.vertex(t, 1))),
500		✗	vec3(mesh_->vertices.point_ptr(mesh_->facets.vertex(t, 2)))
501			);
502			if(vertex_mass.is_bound()) {
503		✗	tri_mass *= (
504		✗	vertex_mass[mesh_->facets.vertex(t, 0)] +
505		✗	vertex_mass[mesh_->facets.vertex(t, 1)] +
506			vertex_mass[mesh_->facets.vertex(t, 2)]
507		✗	) / 3.0;
508			}
509		✗	result += tri_mass;
510			}
511		✗	return result;
512			}
513
514		✗	void OptimalTransportMapOnSurface::call_callback_on_RVD() {
515		✗	RVD_->for_each_polygon(
516		✗	*dynamic_cast<RVDPolygonCallback*>(callback_),
517			false, // symbolic
518			false, // connected components priority
519		✗	!clip_by_balls_ // parallel
520			);
521			// clip_by_balls deactivates parallel mode, because it needs to access
522			// the PointAllocator of the current thread, and we do not have any
523			// access (for now).
524		✗	}
525
526			/**********************************************************************/
527
528		✗	void compute_Laguerre_centroids_on_surface(
529			Mesh* omega,
530			index_t nb_points,
531			const double* points,
532			double* centroids,
533			Mesh* RVD,
534			bool verbose
535			) {
536
537		✗	omega->vertices.set_dimension(4);
538
539			// false = no BRIO
540			// (OTM does not use multilevel and lets Delaunay
541			// reorder the vertices)
542			OptimalTransportMapOnSurface OTM(
543			omega,
544		✗	"BPOW", // "PDEL" will not be much faster because
545			// we are on a surface (lots of threads interactions).
546			false
547		✗	);
548
549			static bool initialized = false;
550			static bool has_cholmod = false;
551		✗	if(!initialized) {
552		✗	initialized = true;
553		✗	has_cholmod = (nlInitExtension("CHOLMOD") == NL_TRUE);
554			}
555		✗	if(has_cholmod) {
556			OTM.set_regularization(1e-3);
557			OTM.set_linear_solver(OT_CHOLMOD);
558			}
559
560			OTM.set_Newton(true);
561		✗	OTM.set_points(nb_points, points);
562			OTM.set_epsilon(0.01);
563			OTM.set_Laguerre_centroids(centroids);
564			OTM.set_verbose(verbose);
565		✗	OTM.optimize(1000);
566
567		✗	if(RVD != nullptr) {
568		✗	OTM.get_RVD(*RVD);
569			}
570
571		✗	omega->vertices.set_dimension(3);
572		✗	}
573
574			/**********************************************************************/
575
576			}
577