Visual Computing Library
Loading...
Searching...
No Matches
curvature.h
1/*****************************************************************************
2 * VCLib *
3 * Visual Computing Library *
4 * *
5 * Copyright(C) 2021-2025 *
6 * Visual Computing Lab *
7 * ISTI - Italian National Research Council *
8 * *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the Mozilla Public License Version 2.0 as published *
13 * by the Mozilla Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * Mozilla Public License Version 2.0 *
20 * (https://www.mozilla.org/en-US/MPL/2.0/) for more details. *
21 ****************************************************************************/
22
23#ifndef VCL_ALGORITHMS_MESH_UPDATE_CURVATURE_H
24#define VCL_ALGORITHMS_MESH_UPDATE_CURVATURE_H
25
26#include <vclib/algorithms/core/polygon.h>
27#include <vclib/algorithms/core/stat.h>
28#include <vclib/algorithms/mesh/intersection.h>
29#include <vclib/algorithms/mesh/point_sampling.h>
30#include <vclib/algorithms/mesh/stat.h>
31#include <vclib/algorithms/mesh/update/normal.h>
32#include <vclib/math/transform.h>
33#include <vclib/mesh/requirements.h>
34#include <vclib/misc/logger.h>
35#include <vclib/misc/parallel.h>
36#include <vclib/space/complex/grid.h>
37#include <vclib/space/complex/mesh_pos.h>
38#include <vclib/space/core/principal_curvature.h>
39#include <vclib/views/pointers.h>
40
41#include <mutex>
42
43namespace vcl {
44
45typedef enum {
46 VCL_PRINCIPAL_CURVATURE_TAUBIN95,
47 VCL_PRINCIPAL_CURVATURE_PCA
48} VCLibPrincipalCurvatureAlgorithm;
49
50template<FaceMeshConcept MeshType, LoggerConcept LogType = NullLogger>
51void updatePrincipalCurvatureTaubin95(MeshType& m, LogType& log)
52{
53 requirePerVertexPrincipalCurvature(m);
54 requirePerVertexAdjacentFaces(m);
56
57 using VertexType = MeshType::VertexType;
58 using CoordType = VertexType::CoordType;
59 using NormalType = VertexType::NormalType;
60 using ScalarType = CoordType::ScalarType;
61 using FaceType = MeshType::FaceType;
62
63 // Auxiliary data structure
64 struct AdjVertex
65 {
66 const VertexType* vert;
67 double doubleArea;
68 bool isBorder;
69 };
70
71 log.log(0, "Updating per vertex normals...");
72
73 updatePerVertexNormalsAngleWeighted(m);
74 normalizePerVertexNormals(m);
75
76 log.log(5, "Computing per vertex curvature...");
77 // log every 5%, starting from 5% to 100%
78 log.startProgress("", m.vertexNumber(), 5, 5, 100);
79
80 for (VertexType& v : m.vertices()) {
81 std::vector<ScalarType> weights;
82 std::vector<AdjVertex> vertices;
83
84 MeshPos<FaceType> pos(v.adjFace(0), &v);
85 const VertexType* firstVertex = pos.adjVertex();
86 const VertexType* tmpVertex;
87 float totalDoubleAreaSize = 0;
88
89 // compute the area of each triangle around the central vertex as well
90 // as their total area
91 do {
92 pos.nextEdgeAdjacentToV();
93 tmpVertex = pos.adjVertex();
94 AdjVertex adjV;
95 adjV.isBorder = pos.isEdgeOnBorder();
96 adjV.vert = tmpVertex;
97 adjV.doubleArea = faceArea(*pos.face()) * 2;
98 totalDoubleAreaSize += adjV.doubleArea;
99 vertices.push_back(adjV);
100 } while (tmpVertex != firstVertex);
101
102 weights.reserve(vertices.size());
103
104 for (int i = 0; i < vertices.size(); ++i) {
105 if (vertices[i].isBorder) {
106 weights.push_back(vertices[i].doubleArea / totalDoubleAreaSize);
107 }
108 else {
109 weights.push_back(
110 0.5f *
111 (vertices[i].doubleArea +
112 vertices[(i - 1) % vertices.size()].doubleArea) /
113 totalDoubleAreaSize);
114 }
115 assert(weights.back() < 1.0f);
116 }
117
118 // compute I-NN^t to be used for computing the T_i's
119 Matrix33<ScalarType> Tp;
120
121 NormalType n = v.normal();
122 for (int i = 0; i < 3; ++i)
123 Tp(i, i) = 1.0f - std::pow(n[i], 2);
124 Tp(0, 1) = Tp(1, 0) = -1.0f * (n[0] * n[1]);
125 Tp(1, 2) = Tp(2, 1) = -1.0f * (n[1] * n[2]);
126 Tp(0, 2) = Tp(2, 0) = -1.0f * (n[0] * n[2]);
127
128 // for all neighbors vi compute the directional curvatures k_i and the
129 // T_i compute M by summing all w_i k_i T_i T_i^t
130 Matrix33<ScalarType> tempMatrix;
131 Matrix33<ScalarType> M = Matrix33<ScalarType>::Zero();
132 for (size_t i = 0; i < vertices.size(); ++i) {
133 CoordType edge = (v.coord() - vertices[i].vert->coord());
134 float curvature =
135 (2.0f * (v.normal().dot(edge))) / edge.squaredNorm();
136 CoordType t = Tp * edge;
137 t.normalize();
138 tempMatrix = t.outerProduct(t);
139 M += tempMatrix * weights[i] * curvature;
140 }
141
142 // compute vector W for the Householder matrix
143 CoordType w;
144 CoordType e1(1.0f, 0.0f, 0.0f);
145 if ((e1 - v.normal()).squaredNorm() > (e1 + v.normal()).squaredNorm())
146 w = e1 - v.normal();
147 else
148 w = e1 + v.normal();
149 w.normalize();
150
151 // compute the Householder matrix I - 2WW^t
152 Matrix33<ScalarType> Q = Matrix33<ScalarType>::Identity();
153 tempMatrix = w.outerProduct(w);
154 Q -= tempMatrix * 2.0f;
155
156 // compute matrix Q^t M Q
157 Matrix33<ScalarType> QtMQ = (Q.transpose() * M * Q);
158
159 Eigen::Matrix<ScalarType, 1, 3> T1 = Q.col(1);
160 Eigen::Matrix<ScalarType, 1, 3> T2 = Q.col(2);
161
162 // find sin and cos for the Givens rotation
163 ScalarType s, c;
164 // Gabriel Taubin hint and Valentino Fiorin impementation
165 ScalarType alpha = QtMQ(1, 1) - QtMQ(2, 2);
166 ScalarType beta = QtMQ(2, 1);
167
168 ScalarType h[2];
169 ScalarType delta =
170 std::sqrt(4.0f * std::pow(alpha, 2) + 16.0f * std::pow(beta, 2));
171 h[0] = (2.0f * alpha + delta) / (2.0f * beta);
172 h[1] = (2.0f * alpha - delta) / (2.0f * beta);
173
174 ScalarType t[2];
175 ScalarType bestC, bestS;
176 ScalarType minError = std::numeric_limits<ScalarType>::infinity();
177 for (uint i = 0; i < 2; i++) {
178 delta = std::sqrt(std::pow(h[i], 2) + 4.0f);
179 t[0] = (h[i] + delta) / 2.0f;
180 t[1] = (h[i] - delta) / 2.0f;
181
182 for (uint j = 0; j < 2; j++) {
183 ScalarType squaredT = std::pow(t[j], 2);
184 ScalarType denominator = 1.0f + squaredT;
185
186 s = (2.0f * t[j]) / denominator;
187 c = (1 - squaredT) / denominator;
188
189 ScalarType approximation =
190 c * s * alpha + (std::pow(c, 2) - std::pow(s, 2)) * beta;
191 ScalarType angleSimilarity =
192 std::abs(std::acos(c) / std::asin(s));
193 ScalarType error =
194 std::abs(1.0f - angleSimilarity) + std::abs(approximation);
195 if (error < minError) {
196 minError = error;
197 bestC = c;
198 bestS = s;
199 }
200 }
201 }
202 c = bestC;
203 s = bestS;
204
205 Eigen::Matrix2f minor2x2;
206 Eigen::Matrix2f S;
207
208 // diagonalize M
209 minor2x2(0, 0) = QtMQ(1, 1);
210 minor2x2(0, 1) = QtMQ(1, 2);
211 minor2x2(1, 0) = QtMQ(2, 1);
212 minor2x2(1, 1) = QtMQ(2, 2);
213
214 S(0, 0) = S(1, 1) = c;
215 S(0, 1) = s;
216 S(1, 0) = -1.0f * s;
217
218 Eigen::Matrix2f StMS = S.transpose() * minor2x2 * S;
219
220 // compute curvatures and curvature directions
221 ScalarType principalCurv1 = (3.0f * StMS(0, 0)) - StMS(1, 1);
222 ScalarType principalCurv2 = (3.0f * StMS(1, 1)) - StMS(0, 0);
223
224 Eigen::Matrix<ScalarType, 1, 3> principalDir1 = T1 * c - T2 * s;
225 Eigen::Matrix<ScalarType, 1, 3> principalDir2 = T1 * s + T2 * c;
226
227 v.principalCurvature().maxDir() = principalDir1;
228 v.principalCurvature().minDir() = principalDir2;
229 v.principalCurvature().maxValue() = principalCurv1;
230 v.principalCurvature().minValue() = principalCurv2;
231
232 log.progress(m.index(v));
233 }
234
235 log.endProgress();
236 log.log(100, "Per vertex curvature computed.");
237}
238
251template<FaceMeshConcept MeshType, LoggerConcept LogType = NullLogger>
252void updatePrincipalCurvaturePCA(
253 MeshType& m,
254 typename MeshType::VertexType::CoordType::ScalarType radius,
255 bool montecarloSampling = true,
256 LogType& log = nullLogger)
257{
258 using VertexType = MeshType::VertexType;
259 using CoordType = VertexType::CoordType;
260 using ScalarType = CoordType::ScalarType;
261 using NormalType = VertexType::NormalType;
262 using FaceType = MeshType::FaceType;
263
264 using VGrid = StaticGrid3<VertexType*, ScalarType>;
265 using VGridIterator = VGrid::ConstIterator;
266
267 VGrid pGrid;
268 ScalarType area;
269
270 log.log(0, "Updating per vertex normals...");
271
272 updatePerVertexNormalsAngleWeighted(m);
273 normalizePerVertexNormals(m);
274
275 log.log(0, "Computing per vertex curvature...");
276 log.startProgress("", m.vertexNumber());
277
278 if (montecarloSampling) {
279 area = surfaceArea(m);
280 pGrid = VGrid(m.vertices() | views::addrOf);
281 pGrid.build();
282 }
283
284 parallelFor(m.vertices(), [&](VertexType& v) {
285 // for (VertexType& v : m.vertices()) {
286 Matrix33<ScalarType> A, eigenvectors;
287 CoordType bp, eigenvalues;
288 if (montecarloSampling) {
289 Sphere s(v.coord(), radius);
290 std::vector<VGridIterator> vec = pGrid.valuesInSphere(s);
291 std::vector<CoordType> points;
292 points.reserve(vec.size());
293 for (const auto& it : vec) {
294 points.push_back(it->second->coord());
295 }
296 A = covarianceMatrixOfPointCloud(points);
297 A *= area * area / 1000;
298 }
299 else {
300 Sphere<ScalarType> sph(v.coord(), radius);
301 MeshType tmpMesh = intersection(m, sph);
302
303 A = covarianceMatrixOfMesh(tmpMesh);
304 }
305
306 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<ScalarType, 3, 3>> eig(A);
307 eigenvalues = CoordType(eig.eigenvalues());
308 eigenvectors = eig.eigenvectors(); // eigenvector are stored as columns.
309 // get the estimate of curvatures from eigenvalues and eigenvectors
310 // find the 2 most tangent eigenvectors (by finding the one closest to
311 // the normal)
312 uint best = 0;
313 ScalarType bestv = std::abs(
314 v.normal().dot(CoordType(eigenvectors.col(0).normalized())));
315 for (uint i = 1; i < 3; ++i) {
316 ScalarType prod = std::abs(
317 v.normal().dot(CoordType(eigenvectors.col(i).normalized())));
318 if (prod > bestv) {
319 bestv = prod;
320 best = i;
321 }
322 }
323 v.principalCurvature().maxDir() =
324 (eigenvectors.col((best + 1) % 3).normalized());
325 v.principalCurvature().minDir() =
326 (eigenvectors.col((best + 2) % 3).normalized());
327
328 Matrix33<ScalarType> rot;
329 ScalarType angle;
330 angle = acos(v.principalCurvature().maxDir().dot(v.normal()));
331
332 rot = rotationMatrix<Matrix33<ScalarType>>(
333 CoordType(v.principalCurvature().maxDir().cross(v.normal())),
334 -(M_PI * 0.5 - angle));
335
336 v.principalCurvature().maxDir() = rot * v.principalCurvature().maxDir();
337
338 angle = acos(v.principalCurvature().minDir().dot(v.normal()));
339
340 rot = rotationMatrix<Matrix33<ScalarType>>(
341 CoordType(v.principalCurvature().minDir().cross(v.normal())),
342 -(M_PI * 0.5 - angle));
343
344 v.principalCurvature().minDir() = rot * v.principalCurvature().minDir();
345
346 // computes the curvature values
347 const ScalarType r5 = std::pow(radius, 5);
348 const ScalarType r6 = r5 * radius;
349 v.principalCurvature().maxValue() =
350 (2.0 / 5.0) *
351 (4.0 * M_PI * r5 + 15 * eigenvalues[(best + 2) % 3] -
352 45.0 * eigenvalues[(best + 1) % 3]) /
353 (M_PI * r6);
354 v.principalCurvature().minValue() =
355 (2.0 / 5.0) *
356 (4.0 * M_PI * r5 + 15 * eigenvalues[(best + 1) % 3] -
357 45.0 * eigenvalues[(best + 2) % 3]) /
358 (M_PI * r6);
359 if (v.principalCurvature().maxValue() <
360 v.principalCurvature().minValue()) {
361 std::swap(
362 v.principalCurvature().minValue(),
363 v.principalCurvature().maxValue());
364 std::swap(
365 v.principalCurvature().minDir(),
366 v.principalCurvature().maxDir());
367 }
368
369 log.progress(m.index(v));
370 //}
371 });
372
373 log.endProgress();
374 log.log(100, "Per vertex curvature computed.");
375}
376
377template<FaceMeshConcept MeshType, LoggerConcept LogType = NullLogger>
378void updatePrincipalCurvature(MeshType& m, LogType& log = nullLogger)
379{
380 requirePerVertexPrincipalCurvature(m);
381
382 updatePrincipalCurvatureTaubin95(m, log);
383}
384
385template<FaceMeshConcept MeshType, LoggerConcept LogType = NullLogger>
386void updatePrincipalCurvature(
387 MeshType& m,
388 VCLibPrincipalCurvatureAlgorithm alg = VCL_PRINCIPAL_CURVATURE_TAUBIN95,
389 LogType& log = nullLogger)
390{
391 requirePerVertexPrincipalCurvature(m);
392
393 double radius;
394 switch (alg) {
395 case VCL_PRINCIPAL_CURVATURE_TAUBIN95:
396 updatePrincipalCurvatureTaubin95(m, log);
397 break;
398 case VCL_PRINCIPAL_CURVATURE_PCA:
399 radius = boundingBox(m).diagonal() * 0.1;
400 updatePrincipalCurvaturePCA(m, radius, true, log);
401 }
402}
403
404} // namespace vcl
405
406#endif // VCL_ALGORITHMS_MESH_UPDATE_CURVATURE_H
Segment()
Default constructor. Creates a segment with endpoints at the origin.
Definition segment.h:68
auto boundingBox(const PointType &p)
Compute the bounding box of a single point.
Definition bounding_box.h:65
auto faceArea(const FaceType &f)
Computes the area of a face. Works both for triangle and polygonal faces, and it is optimized in case...
Definition geometry.h:101
void requirePerFaceAdjacentFaces(const MeshType &m)
This function asserts that a Mesh has a FaceContainer, the Face has a AdjacentFaces Component,...
Definition face_requirements.h:698
NullLogger nullLogger
The nullLogger object is an object of type NullLogger that is used as default argument in the functio...
Definition null_logger.h:125
constexpr detail::VerticesView vertices
A view that allows to iterate over the Vertex elements of an object.
Definition vertex.h:60
constexpr detail::AddressOfView addrOf
The addrOf view applies the address-of operator & on the input view.
Definition pointers.h:120