GCC Code Coverage Report

Directory: ./
File: lib/exploragram/optimal_transport/optimal_transport_2d.cpp
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 0 225 0.0%
Functions: 0 22 0.0%
Branches: 0 226 0.0%
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_2d.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
49 namespace {
50 using namespace GEO;
51
52 /**************************************************************************/
53
54 /**
55 * \brief Clips a polygon by a half-space.
56 * \param[in] P the polygon to be clipped
57 * \param[in] Pi the equation of the half-space,
58 * Pi.x*x + Pi.y*y + Pi.z >= 0
59 * \param[out] clipped the result
60 * \param[in,out] first index to be used for intersections
61 * \param[in] alloc the PointAllocator used to create the
62 * new vertices
63 */
64 ✗ void clip_polygon_by_halfplane(
65 const GEOGen::Polygon& P,
66 vec3 Pi,
67 GEOGen::Polygon& target,
68 index_t& n,
69 GEOGen::PointAllocator* alloc
70 ) {
71 target.clear();
72 ✗ if(P.nb_vertices() == 0) {
73 return;
74 }
75
76 // The predecessor of the first vertex is the last vertex
77 ✗ index_t prev_k = P.nb_vertices() - 1;
78 const GEOGen::Vertex* prev_vk = &(P.vertex(prev_k));
79 const double* geo_restrict prev_pk = prev_vk->point();
80
81 GEO::Sign prev_status = GEO::geo_sgn(
82 ✗ prev_pk[0]*Pi.x + prev_pk[1]*Pi.y + Pi.z
83 );
84
85 ✗ for(index_t k = 0; k < P.nb_vertices(); k++) {
86 const GEOGen::Vertex* vk = &(P.vertex(k));
87 const double* pk = vk->point();
88
89 ✗ GEO::Sign status = GEO::geo_sgn(pk[0]*Pi.x + pk[1]*Pi.y + Pi.z);
90
91 // If status of edge extremities differ,
92 // then there is an intersection.
93 if(status != prev_status && (prev_status != 0)) {
94 GEOGen::Vertex I;
95 ✗ double* Ipoint = alloc->new_item();
96 I.set_point(Ipoint);
97
98 // Compute lambda1 and lambda2, the
99 // barycentric coordinates of the intersection I
100 // in the segment [prev_vk vk]
101 // Note that d and l (used for the predicates)
102 // are reused here.
103
104
105 ✗ double denom = Pi.x * (
106 ✗ pk[0] - prev_pk[0]) + Pi.y * (pk[1] - prev_pk[1]
107 );
108 ✗ double lambda2 =
109 ✗ -(Pi.z + Pi.x * prev_pk[0] + Pi.y * prev_pk[1]);
110 ✗ double lambda1 = denom - lambda2;
111
112 // Shit happens ! [Forrest Gump]
113 ✗ if(::fabs(denom) < 1e-20) {
114 lambda1 = 0.5;
115 lambda2 = 0.5;
116 } else {
117 ✗ lambda1 = lambda1 / denom;
118 ✗ lambda2 = lambda2 / denom;
119 }
120
121 ✗ Ipoint[0] = lambda1 * prev_pk[0] + lambda2 * pk[0];
122 ✗ Ipoint[1] = lambda1 * prev_pk[1] + lambda2 * pk[1];
123
124 ✗ I.set_weight(
125 ✗ lambda1 * prev_vk->weight() + lambda2 * vk->weight()
126 );
127 ✗ if(status > 0) {
128 I.copy_edge_from(*prev_vk);
129 ✗ I.set_adjacent_seed(signed_index_t(n));
130 ✗ ++n;
131 } else {
132 I.set_flag(GEOGen::INTERSECT);
133 I.set_adjacent_seed(vk->adjacent_seed());
134 }
135 target.add_vertex(I);
136 }
137 ✗ if(status > 0) {
138 target.add_vertex(*vk);
139 }
140 prev_vk = vk;
141 prev_pk = pk;
142 prev_status = status;
143 prev_k = k;
144 }
145 }
146
147
148 /**
149 * \brief Clips a polygon by a ball.
150 * \param[in] P the polygon to be clipped
151 * \param[in] center the center of the ball
152 * \param[in] radius the radius of the ball
153 * \param[out] clipped the result
154 * \param[in] n first index to be used for intersections (number of
155 * vertices in Delaunay).
156 * \param[in] alloc the PointAllocator used to create the
157 * new vertices
158 * \param[out] work a temporary work zone
159 */
160 ✗ void clip_polygon_by_ball(
161 const GEOGen::Polygon& P,
162 vec2 center, double radius,
163 GEOGen::Polygon& clipped,
164 index_t n,
165 GEOGen::PointAllocator* alloc,
166 GEOGen::Polygon& work
167 ) {
168 const index_t N = 16;
169 const double dalpha = M_PI * 2.0 / double(N);
170 ✗ index_t first_new_index = n;
171 clipped.copy(P);
172 ✗ FOR(i,N) {
173 ✗ double s = ::sin(dalpha*i);
174 ✗ double c = ::cos(dalpha*i);
175 ✗ double off = (center.x+radius*c)*c + (center.y+radius*s)*s;
176 ✗ vec3 Pi(-c,-s,off);
177 ✗ clip_polygon_by_halfplane(
178 clipped, Pi, work, first_new_index, alloc
179 );
180 clipped.swap(work);
181 }
182 ✗ }
183
184 /**************************************************************************/
185
186 // For more details on the algorithm / structure of the objective function,
187 // see also implementation notes at the beginning of
188 // optimal_transport_3d.cpp.
189
190 /**
191 * \brief Computes the contribution of a polygon
192 * to the objective function minimized by a semi-discrete
193 * optimal transport map in 2D.
194 */
195 class OTMPolygonCallback :
196 public OptimalTransportMap::Callback,
197 public RVDPolygonCallback {
198 public:
199
200 /**
201 * \brief OTMPolygonCallback constructor.
202 * \param[in] OTM a pointer to the OptimalTransportMap2d
203 */
204 ✗ OTMPolygonCallback(OptimalTransportMap2d* OTM) :
205 ✗ OptimalTransportMap::Callback(OTM) {
206 ✗ }
207
208 /**
209 * \copydoc RVDPolygonCallback::operator()
210 */
211 ✗ void operator() (
212 index_t v,
213 index_t t,
214 const GEOGen::Polygon& P
215 ) const override {
216 geo_argused(t);
217 ✗ if(OTM_->air_fraction() != 0.0 && OTM_->nb_air_particles() == 0) {
218 ✗ if(v < OTM_->nb_points()) {
219 OptimalTransportMap2d* OTM =
220 static_cast<OptimalTransportMap2d*>(OTM_);
221 double R = OTM_->weight(v);
222 ✗ geo_assert(R > 0.0);
223 ✗ R = ::sqrt(R);
224 vec2 center(OTM_->point_ptr(v));
225 ✗ clip_polygon_by_ball(
226 ✗ P, center, R, OTM->clipped_,
227 ✗ OTM_->nb_points(),
228 ✗ OTM_->RVD()->point_allocator(),
229 ✗ OTM->work_
230 );
231 ✗ do_it(v,t,OTM->clipped_);
232 }
233 } else {
234 ✗ do_it(v,t,P);
235 }
236 ✗ }
237
238 ✗ void do_it(
239 index_t v,
240 index_t t,
241 const GEOGen::Polygon& P
242 ) const {
243
244 // v can be an air particle.
245 ✗ if(v >= n_) {
246 ✗ return;
247 }
248
249 ✗ if(P.nb_vertices() == 0) {
250 return;
251 }
252
253 geo_argused(t);
254 double m, mgx, mgy;
255 ✗ compute_m_and_mg(P, m, mgx, mgy);
256
257 ✗ if(spinlocks_ != nullptr) {
258 ✗ spinlocks_->acquire_spinlock(v);
259 }
260
261 // +m because we maximize F <=> minimize -F
262 ✗ g_[v] += m;
263
264 ✗ if(Newton_step_) {
265 // ... but here -m because Newton step =
266 // solve H p = -g (minus g in the RHS).
267 ✗ OTM_->add_i_right_hand_side(v,-m);
268 }
269
270 ✗ if(mg_ != nullptr) {
271 ✗ mg_[2*v] += mgx;
272 ✗ mg_[2*v+1] += mgy;
273 }
274
275 ✗ if(spinlocks_ != nullptr) {
276 ✗ spinlocks_->release_spinlock(v);
277 }
278
279 ✗ if(Newton_step_) {
280 // Spinlocks are managed internally by update_Hessian().
281 ✗ update_Hessian(P, v);
282 }
283
284
285 ✗ if(eval_F_) {
286 ✗ Thread* thread = Thread::current();
287 index_t current_thread_id =
288 (thread == nullptr) ? 0 : thread->id();
289 ✗ double F = weighted_ ? eval_F_weighted(P, v) : eval_F(P, v);
290 const_cast<OTMPolygonCallback*>(this)->
291 funcval_[current_thread_id] += F;
292 }
293 }
294
295 protected:
296
297 /**
298 * \brief Computes the mass and mass times centroid of the
299 * current intersection polygon.
300 * \details Weights are taken into account if present.
301 * \param[in] P a const reference to the current intersection polygon.
302 * \param[out] m , mgx , mgy the mass and the mass times the
303 * centroid of the ConvexCell. mgx and mgy are not computed
304 * if mg_ is nullptr.
305 */
306 ✗ void compute_m_and_mg(
307 const GEOGen::Polygon& P,
308 double& m, double& mgx, double& mgy
309 ) const {
310 ✗ m = 0.0;
311 ✗ mgx = 0.0;
312 ✗ mgy = 0.0;
313 const GEOGen::Vertex& V0 = P.vertex(0);
314 const double* p0 = V0.point();
315 ✗ for(index_t i=1; i+1<P.nb_vertices(); ++i) {
316 const GEOGen::Vertex& V1 = P.vertex(i);
317 const double* p1 = V1.point();
318 const GEOGen::Vertex& V2 = P.vertex(i+1);
319 const double* p2 = V2.point();
320 ✗ double cur_m = triangle_mass(V0,V1,V2);
321 ✗ m += cur_m;
322 ✗ if(mg_ != nullptr) {
323 ✗ if(weighted_) {
324 double w0 = V0.weight();
325 double w1 = V1.weight();
326 double w2 = V2.weight();
327 ✗ double s = cur_m/(w0+w1+w2);
328 ✗ mgx += s * (w0*p0[0] + w1*p1[0] + w2*p2[0]);
329 ✗ mgy += s * (w0*p0[1] + w1*p1[1] + w2*p2[1]);
330 } else {
331 ✗ mgx += cur_m * (p0[0] + p1[0] + p2[0]) / 3.0;
332 ✗ mgy += cur_m * (p0[1] + p1[1] + p2[1]) / 3.0;
333 }
334 }
335 }
336 ✗ }
337
338 /**
339 * \brief Updates the Hessian according to the current intersection
340 * polygon.
341 * \param[in] P a const reference to the current intersection polygon.
342 * \param[in] i the current seed
343 */
344 ✗ void update_Hessian(
345 const GEOGen::Polygon& P, index_t i
346 ) const {
347
348 // The coefficient of the Hessian associated to a pair of
349 // adjacent cells Lag(i),Lag(j) is :
350 // - mass(Lag(i) /\ Lag(j)) / (2*distance(pi,pj))
351
352 ✗ const double* pi = OTM_->point_ptr(i);
353
354 ✗ for(index_t k1=0; k1<P.nb_vertices(); ++k1) {
355 ✗ index_t k2 = k1+1;
356 ✗ if(k2 == P.nb_vertices()) {
357 k2 = 0;
358 }
359 // Note:
360 // It is P.vertex(k2).adjacent_seed(),
361 // not P.vertex(k1).adjacent_seed() !!!
362 ✗ index_t j = index_t(P.vertex(k2).adjacent_seed());
363 ✗ if(j != index_t(-1)) {
364 double hij = 0.0;
365 ✗ if(j < n_) {
366 ✗ const double* pj = OTM_->point_ptr(j);
367 ✗ hij =
368 ✗ edge_mass(P.vertex(k1), P.vertex(k2)) /
369 ✗ (2.0 * GEO::Geom::distance(pi,pj,2)) ;
370 ✗ } else if(OTM_->nb_air_particles() != 0) {
371 const double* pj = OTM_->point_ptr(j); // points and air
372 ✗ hij =
373 ✗ edge_mass(P.vertex(k1), P.vertex(k2)) /
374 ✗ (2.0 * GEO::Geom::distance(pi,pj,2)) ;
375 } else {
376 double R = OTM_->weight(i);
377 ✗ geo_assert(R >= 0.0);
378 ✗ R = ::sqrt(R);
379 ✗ hij = edge_mass(P.vertex(k1), P.vertex(k2)) / (2.0 * R);
380 }
381
382 // -hij because we maximize F <=> minimize -F
383 ✗ if(hij != 0.0) {
384 ✗ if(spinlocks_ != nullptr) {
385 ✗ spinlocks_->acquire_spinlock(i);
386 }
387 // Diagonal is positive, extra-diagonal
388 // coefficients are negative,
389 // this is a convex function.
390 ✗ if(j < n_) {
391 ✗ OTM_->add_ij_coefficient(i, j, -hij);
392 }
393 ✗ OTM_->add_ij_coefficient(i, i, hij);
394 ✗ if(spinlocks_ != nullptr) {
395 ✗ spinlocks_->release_spinlock(i);
396 }
397 }
398 }
399 }
400 ✗ }
401
402 /**
403 * \brief Computes the contribution of the current polygon
404 * to the objective function, in the uniform (non-weighted)
405 * case.
406 * \param[in] P a const reference to the current polygon.
407 * \param[in] i the current seed
408 */
409 ✗ double eval_F(const GEOGen::Polygon& P, index_t i) const {
410 geo_debug_assert(!weighted_);
411 geo_argused(P);
412 geo_argused(i);
413 // Not implemented yet.
414 ✗ geo_assert_not_reached;
415 }
416
417 /**
418 * \brief Computes the contribution of the current polygon
419 * to the objective function, in the weighted case.
420 * \param[in] P a const reference to the current polygon.
421 * \param[in] i the current seed
422 */
423 ✗ double eval_F_weighted(const GEOGen::Polygon& P, index_t i) const {
424 geo_argused(P);
425 geo_argused(i);
426 // Not implemented yet.
427 ✗ geo_assert_not_reached;
428 }
429
430 /**
431 * \brief Computes the mass of a triangle.
432 * \details Weights are taken into account in weighted_ mode.
433 * \param[in] V0 , V1 , V2 the three vertices of the triangle,
434 * given as RVD vertices.
435 * \return the area of the triangle in 3D times the average value of the
436 * three weights.
437 */
438 ✗ double triangle_mass(
439 const GEOGen::Vertex& V0,
440 const GEOGen::Vertex& V1,
441 const GEOGen::Vertex& V2
442 ) const {
443 double m =
444 Geom::triangle_area_2d(V0.point(), V1.point(), V2.point());
445 ✗ if(weighted_) {
446 ✗ m *= ((V0.weight() + V1.weight() + V2.weight())/3.0);
447 }
448 ✗ return m;
449 }
450
451 /**
452 * \brief Computes the mass of an edge.
453 * \details Weights are taken into account in weighted_ mode.
454 * \param[in] V0 , V1 the two vertices of the edge,
455 * given as RVD vertices.
456 * \return the length of the edge in 3D times the average value of the
457 * two weights.
458 */
459 ✗ double edge_mass(
460 const GEOGen::Vertex& V0,
461 const GEOGen::Vertex& V1
462 ) const {
463 ✗ double m = Geom::distance(V0.point(), V1.point(), 2);
464 ✗ if(weighted_) {
465 ✗ m *= ((V0.weight() + V1.weight())/2.0);
466 }
467 ✗ return m;
468 }
469 };
470
471 /********************************************************************/
472
473 /**
474 * \brief A RVDPolygonCallback that stores the Restricted Voronoi
475 * Diagram in a Mesh.
476 */
477 class ComputeRVDPolygonCallback : public RVDPolygonCallback {
478 public:
479 ✗ ComputeRVDPolygonCallback(OptimalTransportMap* OTM, Mesh* target) :
480 ✗ OTM_(OTM), target_(target) {
481 ✗ target_->clear();
482 ✗ target_->vertices.set_dimension(3);
483 ✗ chart_.bind(target_->facets.attributes(), "chart");
484 ✗ }
485
486 ✗ ~ComputeRVDPolygonCallback() override {
487 ✗ chart_.unbind();
488 ✗ }
489
490 ✗ void operator() (
491 index_t v,
492 index_t t,
493 const GEOGen::Polygon& P
494 ) const override {
495 geo_argused(t);
496 ✗ if(OTM_->air_fraction() != 0.0 && OTM_->nb_air_particles() == 0) {
497 ✗ if(v < OTM_->nb_points()) {
498 OptimalTransportMap2d* OTM =
499 static_cast<OptimalTransportMap2d*>(OTM_);
500 double R = OTM_->weight(v);
501 ✗ if(R < 0.0) {
502 ✗ std::cerr << '-' << std::flush;
503 }
504 ✗ R = R > 0.0 ? ::sqrt(R) : 0.0;
505 ✗ vec2 center(OTM_->point_ptr(v));
506 ✗ clip_polygon_by_ball(
507 ✗ P, center, R, OTM->clipped_,
508 ✗ OTM_->nb_points(),
509 ✗ OTM_->RVD()->point_allocator(),
510 ✗ OTM->work_
511 );
512 ✗ do_it(v,t,OTM->clipped_);
513 }
514 } else {
515 ✗ do_it(v,t,P);
516 }
517 ✗ }
518
519 ✗ void do_it(
520 index_t v,
521 index_t t,
522 const GEOGen::Polygon& P
523 ) const {
524 geo_argused(v);
525 geo_argused(t);
526
527 ✗ if(P.nb_vertices() == 0) {
528 return;
529 }
530
531 ✗ index_t voffset = target_->vertices.nb();
532 ✗ FOR(i,P.nb_vertices()) {
533 const double* p = P.vertex(i).point();
534 ✗ if(OTM_->dimension() == 2) {
535 ✗ target_->vertices.create_vertex(
536 ✗ vec3(p[0], p[1], 0.0).data()
537 );
538 } else {
539 ✗ target_->vertices.create_vertex(
540 ✗ vec3(p[0], p[1], p[2]).data()
541 );
542 }
543 }
544 ✗ index_t f = target_->facets.create_polygon(P.nb_vertices());
545 ✗ FOR(i,P.nb_vertices()) {
546 ✗ target_->facets.set_vertex(f,i,voffset+i);
547 }
548 ✗ const_cast<Attribute<index_t>&>(chart_)[f] = v;
549 }
550
551 private:
552 OptimalTransportMap* OTM_;
553 Mesh* target_;
554 Attribute<index_t> chart_;
555 };
556
557 /********************************************************************/
558 }
559
560 namespace GEO {
561
562 ✗ OptimalTransportMap2d::OptimalTransportMap2d(
563 Mesh* mesh, const std::string& delaunay, bool BRIO
564 ✗ ) :
565 OptimalTransportMap(
566 2,
567 mesh,
568 ✗ (delaunay == "default") ? "BPOW2d" : delaunay,
569 BRIO
570 ✗ ) {
571 ✗ callback_ = new OTMPolygonCallback(this);
572 ✗ total_mass_ = total_mesh_mass();
573 ✗ }
574
575 ✗ OptimalTransportMap2d::~OptimalTransportMap2d() {
576 ✗ }
577
578 ✗ void OptimalTransportMap2d::get_RVD(Mesh& RVD_mesh) {
579 ✗ ComputeRVDPolygonCallback callback(this, &RVD_mesh);
580 ✗ RVD()->for_each_polygon(callback, false, false, false);
581 /*
582 // NOTE: Does not work, TODO: determine why
583 Attribute<index_t> tet_region(RVD_mesh.cells.attributes(),"region");
584 RVD()->compute_RVD(
585 RVD_mesh,
586 0, // dim (0 means use default)
587 false, // borders_only
588 show_RVD_seed_ // integration_simplices
589 );
590 */
591 ✗ }
592
593 ✗ void OptimalTransportMap2d::compute_Laguerre_centroids(double* centroids) {
594 ✗ vector<double> g(nb_points(), 0.0);
595 ✗ Memory::clear(centroids, nb_points()*sizeof(double)*2);
596
597 ✗ callback_->set_Laguerre_centroids(centroids);
598 callback_->set_g(g.data());
599 {
600 Stopwatch* W = nullptr;
601 ✗ if(newton_ && verbose_) {
602 ✗ W = new Stopwatch("RVD");
603 ✗ Logger::out("OTM") << "In RVD (centroids)..." << std::endl;
604 }
605 ✗ RVD_->for_each_polygon(
606 ✗ *dynamic_cast<RVDPolygonCallback*>(callback_),
607 false, false, true
608 );
609 ✗ if(newton_ && verbose_) {
610 ✗ delete W;
611 }
612 }
613
614 ✗ callback_->set_Laguerre_centroids(nullptr);
615
616 ✗ for(index_t v=0; v<nb_points(); ++v) {
617 ✗ centroids[2*v ] /= g[v];
618 ✗ centroids[2*v+1] /= g[v];
619 }
620 ✗ }
621
622 ✗ double OptimalTransportMap2d::total_mesh_mass() const {
623 double result = 0.0;
624
625 // This is terribly confusing, the parameters for
626 // a power diagram are called "weights", and the
627 // standard attribute name for vertices density is
628 // also called "weight" (and is unrelated).
629 // In this program, what is called weight corresponds
630 // to the parameters of the power diagram (except the
631 // name of the attribute), and everything that corresponds
632 // to mass/density is called mass.
633 Attribute<double> vertex_mass;
634 ✗ vertex_mass.bind_if_is_defined(
635 ✗ mesh_->vertices.attributes(), "weight"
636 );
637
638 ✗ for(index_t t: mesh_->facets) {
639 double tri_mass = GEO::Geom::triangle_area(
640 vec2(mesh_->vertices.point_ptr(mesh_->facets.vertex(t, 0))),
641 vec2(mesh_->vertices.point_ptr(mesh_->facets.vertex(t, 1))),
642 vec2(mesh_->vertices.point_ptr(mesh_->facets.vertex(t, 2)))
643 );
644 if(vertex_mass.is_bound()) {
645 ✗ tri_mass *= (
646 ✗ vertex_mass[mesh_->facets.vertex(t, 0)] +
647 ✗ vertex_mass[mesh_->facets.vertex(t, 1)] +
648 vertex_mass[mesh_->facets.vertex(t, 2)]
649 ✗ ) / 3.0;
650 }
651 ✗ result += tri_mass;
652 }
653 ✗ return result;
654 }
655
656 ✗ void OptimalTransportMap2d::call_callback_on_RVD() {
657 ✗ RVD_->for_each_polygon(
658 ✗ *dynamic_cast<RVDPolygonCallback*>(callback_),
659 false, // symbolic
660 false, // connected components priority
661 ✗ !clip_by_balls_ // parallel
662 );
663 // clip_by_balls deactivates parallel mode, because it needs to access
664 // the PointAllocator of the current thread, and we do not have any
665 // access (for now).
666 ✗ }
667
668 /**********************************************************************/
669
670 ✗ void compute_Laguerre_centroids_2d(
671 Mesh* omega,
672 index_t nb_points,
673 const double* points,
674 double* centroids,
675 Mesh* RVD,
676 bool verbose,
677 index_t nb_air_particles,
678 const double* air_particles,
679 index_t air_particles_stride,
680 double air_fraction,
681 const double* weights_in,
682 double* weights_out,
683 index_t nb_iter
684 ) {
685 // Omega can be either 2d or 3d with third coordinate set to
686 // zero.
687 index_t omega_dim_backup = omega->vertices.dimension();
688 ✗ omega->vertices.set_dimension(3);
689
690 // false = no BRIO
691 // (OTM does not use multilevel and lets Delaunay
692 // reorder the vertices)
693 OptimalTransportMap2d OTM(
694 omega,
695 ✗ std::string("BPOW2d"),
696 false
697 ✗ );
698
699 static bool initialized = false;
700 static bool has_cholmod = false;
701 ✗ if(!initialized) {
702 ✗ initialized = true;
703 ✗ has_cholmod = (nlInitExtension("CHOLMOD") == NL_TRUE);
704 }
705
706 ✗ if(has_cholmod) {
707 OTM.set_regularization(1e-3);
708 OTM.set_linear_solver(OT_CHOLMOD);
709 }
710
711 OTM.set_Newton(true);
712 OTM.set_air_particles(
713 nb_air_particles, air_particles, air_particles_stride, air_fraction
714 );
715 ✗ OTM.set_points(nb_points, points);
716 ✗ if(weights_in != nullptr) {
717 ✗ FOR(i, nb_points) {
718 ✗ OTM.set_initial_weight(i, weights_in[i]);
719 }
720 }
721 OTM.set_epsilon(0.01);
722 OTM.set_Laguerre_centroids(centroids);
723 OTM.set_verbose(verbose);
724 ✗ OTM.optimize(nb_iter);
725
726 ✗ if(RVD != nullptr) {
727 ✗ OTM.get_RVD(*RVD);
728 }
729
730 ✗ omega->vertices.set_dimension(omega_dim_backup);
731
732 ✗ if(weights_out != nullptr) {
733 ✗ FOR(v, OTM.nb_points()) {
734 ✗ weights_out[v] = OTM.weight(v);
735 }
736 }
737 ✗ }
738 }
739