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/VSDM.h>
41			#include <geogram/mesh/mesh_geometry.h>
42			#include <geogram/mesh/mesh_subdivision.h>
43			#include <geogram/basic/progress.h>
44			#include <geogram/bibliography/bibliography.h>
45
46			namespace {
47			using namespace GEO;
48
49			/**
50			* \brief A callback class for mesh subdivision that interpolates
51			* attributes and constructs the subdivision matrix.
52			*/
53			class SymbolicMeshSplitCallbacks : public MeshSplitCallbacks {
54			public:
55
56			/**
57			* \brief SymbolicMeshSplitCallbacks constructor.
58			* \param[in] mesh a pointer to the mesh.
59			* \param[in] matrix a pointer to the matrix. Should contain
60			* the identity matrix with nv*dim columns, where nv denotes the
61			* number of vertices in the mesh and dim the dimension of the
62			* vertices (mostly 3 in the present case).
63			*/
64		✗	SymbolicMeshSplitCallbacks(Mesh* mesh, NLSparseMatrix* matrix) :
65		✗	MeshSplitCallbacks(mesh), matrix_(matrix) {
66		✗	}
67
68			/**
69			* \copydoc MeshSplitCallbacks::create_vertex()
70			*/
71		✗	index_t create_vertex() override {
72		✗	FOR(i,mesh_->vertices.dimension()) {
73		✗	nlSparseMatrixAddRow(matrix_);
74			}
75		✗	return MeshSplitCallbacks::create_vertex();
76			}
77
78			/**
79			* \copydoc MeshSplitCallbacks::scale_vertex()
80			*/
81		✗	void scale_vertex(index_t v, double s) override {
82		✗	index_t dim = mesh_->vertices.dimension();
83		✗	FOR(i,dim) {
84		✗	nlSparseMatrixScaleRow(matrix_, v*dim+i, s);
85			}
86		✗	MeshSplitCallbacks::scale_vertex(v,s);
87		✗	}
88
89			/**
90			* \copydoc MeshSplitCallbacks::zero_vertex()
91			*/
92		✗	void zero_vertex(index_t v) override {
93		✗	index_t dim = mesh_->vertices.dimension();
94		✗	FOR(i,dim) {
95		✗	nlSparseMatrixZeroRow(matrix_, v*dim+i);
96			}
97		✗	MeshSplitCallbacks::zero_vertex(v);
98		✗	}
99
100			/**
101			* \copydoc MeshSplitCallbacks::madd_vertex()
102			*/
103		✗	void madd_vertex(
104			index_t v1, double s, index_t v2
105			) override {
106		✗	index_t dim = mesh_->vertices.dimension();
107		✗	FOR(i,dim) {
108		✗	nlSparseMatrixMAddRow(matrix_, v1*dim+i, s, v2*dim+i);
109			}
110		✗	MeshSplitCallbacks::madd_vertex(v1,s,v2);
111		✗	}
112
113			private:
114			NLSparseMatrix* matrix_;
115			};
116			}
117
118			namespace GEO {
119
120			VSDM* VSDM::instance_ = nullptr;
121
122		✗	VSDM::VSDM(Mesh* S, Mesh* T):
123		✗	S_(S),
124		✗	T_(T),
125		✗	affinity_(1.0),
126		✗	progress_(nullptr),
127		✗	nb_iter_(0),
128		✗	cur_iter_(0),
129		✗	subd_(nullptr) {
130		✗	compute_graph_Laplacian(S_,&L_);
131		✗	temp_V1_.resize(S_->vertices.nb());
132		✗	temp_V2_.resize(S_->vertices.nb());
133		✗	delaunay_ = Delaunay::create(3);
134		✗	RVD_ = RestrictedVoronoiDiagram::create(delaunay_, T_);
135		✗	optimizer_ = Optimizer::create("HLBFGS");
136		✗	subd_matrix_ = nullptr;
137		✗	geo_cite("DBLP:conf/imr/NivoliersYL11");
138		✗	geo_cite("DBLP:journals/cgf/LuLW12");
139		✗	geo_cite("DBLP:journals/ewc/NivoliersYL14");
140		✗	}
141
142		✗	VSDM::~VSDM() {
143		✗	nlSparseMatrixDestroy(&L_);
144		✗	nlDeleteMatrix(subd_matrix_);
145		✗	}
146
147		✗	void VSDM::optimize(index_t nb_iter) {
148		✗	if(nb_iter == 0) {
149		✗	return;
150			}
151
152		✗	if(progress_ != nullptr) {
153		✗	progress_->reset(nb_iter);
154			}
155
156		✗	double S = Geom::mesh_area(*T_);
157		✗	double R = ::pow(S, 1.0 / 2.0);
158		✗	double CVT_normalization = 1.0 / pow(R, 4.0);
159		✗	double affinity_normalization = 0.001 / pow(R, 2.0);
160		✗	affinity_scaling_ =
161		✗	affinity_ * affinity_normalization / CVT_normalization;
162
163		✗	nb_iter_ = nb_iter;
164		✗	cur_iter_ = 0;
165		✗	instance_ = this;
166		✗	index_t n = S_->vertices.nb() * 3;
167		✗	index_t m = 7;
168		✗	double* x = S_->vertices.point_ptr(0);
169		✗	optimizer_->set_epsg(0.0);
170		✗	optimizer_->set_epsf(0.0);
171		✗	optimizer_->set_epsx(0.0);
172		✗	optimizer_->set_newiteration_callback(VSDM::newiteration_CB);
173		✗	optimizer_->set_funcgrad_callback(VSDM::funcgrad_CB);
174		✗	optimizer_->set_N(n);
175		✗	optimizer_->set_M(m);
176		✗	optimizer_->set_max_iter(nb_iter);
177		✗	optimizer_->optimize(x);
178		✗	instance_ = nullptr;
179			}
180
181		✗	void VSDM::funcgrad(index_t n, double* x, double& f, double* g) {
182		✗	f = 0.0;
183		✗	Memory::clear(g, n * sizeof(double));
184		✗	if(subd_ == nullptr) {
185		✗	delaunay_->set_vertices(n/3, x);
186		✗	RVD_->compute_CVT_func_grad(f, g);
187			} else {
188		✗	subd_g_.assign(subd_->vertices.nb()*3, 0.0);
189		✗	nlMultMatrixVector(subd_matrix_, x, subd_->vertices.point_ptr(0));
190		✗	delaunay_->set_vertices(
191		✗	subd_->vertices.nb(), subd_->vertices.point_ptr(0)
192			);
193		✗	RVD_->compute_CVT_func_grad(f, subd_g_.data());
194
195			// g = transpose(subd_matrix_) * subd_g_
196			{
197		✗	NLCRSMatrix* CRS = (NLCRSMatrix*)(subd_matrix_);
198		✗	FOR(i,CRS->m) {
199		✗	for(index_t jj=CRS->rowptr[i]; jj<CRS->rowptr[i+1]; ++jj) {
200		✗	index_t j = CRS->colind[jj];
201		✗	double a = CRS->val[jj];
202		✗	g[j] += a * subd_g_[i];
203			}
204			}
205			}
206			}
207		✗	if(affinity_ != 0.0) {
208		✗	add_funcgrad_affinity(n,x,f,g);
209			}
210		✗	}
211
212		✗	void VSDM::add_funcgrad_affinity(
213			index_t n, double* x, double& f, double* g
214			) {
215		✗	geo_assert(L_.n*3 == n);
216		✗	FOR(coord,3) {
217		✗	FOR(i,L_.n) {
218		✗	temp_V1_[i] = x[3*i+coord];
219			}
220		✗	nlSparseMatrixMult(&L_, temp_V1_.data(), temp_V2_.data());
221		✗	double F = 0.0;
222		✗	FOR(i,L_.n) {
223		✗	F += temp_V1_[i] * temp_V2_[i];
224			}
225		✗	f += affinity_scaling_ * F;
226		✗	FOR(i,L_.n) {
227		✗	g[3*i + coord] += 2.0 * affinity_scaling_ * temp_V2_[i];
228			}
229			}
230		✗	}
231
232		✗	void VSDM::newiteration() {
233		✗	cur_iter_++;
234		✗	if(cur_iter_ <= nb_iter_) {
235		✗	Logger::out("VSDM")
236		✗	<< "Iter: " << cur_iter_ << "/" << nb_iter_ << std::endl;
237			}
238		✗	if(progress_ != nullptr) {
239		✗	progress_->next();
240			}
241		✗	}
242
243		✗	void VSDM::funcgrad_CB(index_t n, double* x, double& f, double* g) {
244		✗	geo_assert(instance_ != nullptr);
245		✗	instance_->funcgrad(n, x, f, g);
246		✗	}
247
248		✗	void VSDM::newiteration_CB(
249			index_t n, const double* x, double f, const double* g, double gnorm
250			) {
251		✗	geo_argused(n);
252		✗	geo_argused(x);
253		✗	geo_argused(f);
254		✗	geo_argused(g);
255		✗	geo_argused(gnorm);
256		✗	geo_assert(instance_ != nullptr);
257		✗	instance_->newiteration();
258		✗	}
259
260		✗	void VSDM::compute_graph_Laplacian(Mesh* S, NLSparseMatrix* L) {
261		✗	index_t n = S->vertices.nb();
262		✗	nlSparseMatrixConstruct(L, n, n, NL_MATRIX_STORE_ROWS);
263		✗	vector<index_t> v_degree(S->vertices.nb(),0);
264		✗	FOR(c,S->facet_corners.nb()) {
265		✗	++v_degree[S->facet_corners.vertex(c)];
266			}
267		✗	FOR(f,S->facets.nb()) {
268		✗	for(
269		✗	index_t c1 = S->facets.corners_begin(f);
270		✗	c1 < S->facets.corners_end(f); ++c1) {
271		✗	index_t c2 = S->facets.next_corner_around_facet(f,c1);
272		✗	index_t v1 = S->facet_corners.vertex(c1);
273		✗	index_t v2 = S->facet_corners.vertex(c2);
274		✗	double a = 2.0 / (double(v_degree[v1]) + double(v_degree[v2]));
275		✗	nlSparseMatrixAdd(L, v1, v2, -a);
276		✗	nlSparseMatrixAdd(L, v1, v1, a);
277			}
278			}
279		✗	}
280
281		✗	void VSDM::set_subdivision_surface(Mesh* mesh, index_t nb_subdiv) {
282		✗	subd_ = mesh;
283		✗	index_t n = S_->vertices.nb()*3;
284		✗	subd_matrix_ = nlSparseMatrixNew(n, n, NL_MATRIX_STORE_ROWS);
285		✗	FOR(i,n) {
286		✗	nlSparseMatrixAdd((NLSparseMatrix*)subd_matrix_, i, i, 1.0);
287			}
288		✗	subd_->clear();
289		✗	subd_->copy(*S_);
290		✗	SymbolicMeshSplitCallbacks cb(subd_, (NLSparseMatrix*)subd_matrix_);
291		✗	FOR(i,nb_subdiv) {
292		✗	mesh_split_catmull_clark(*subd_, &cb);
293			}
294		✗	nlMatrixCompress(&subd_matrix_);
295		✗	subd_g_.assign(subd_->vertices.nb()*3, 0);
296		✗	}
297			}
298