23#ifndef VCL_ALGORITHMS_MESH_UPDATE_CURVATURE_H
24#define VCL_ALGORITHMS_MESH_UPDATE_CURVATURE_H
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>
46 VCL_PRINCIPAL_CURVATURE_TAUBIN95,
47 VCL_PRINCIPAL_CURVATURE_PCA
48} VCLibPrincipalCurvatureAlgorithm;
50template<FaceMeshConcept MeshType, LoggerConcept LogType = NullLogger>
51void updatePrincipalCurvatureTaubin95(MeshType& m, LogType& log)
53 requirePerVertexPrincipalCurvature(m);
54 requirePerVertexAdjacentFaces(m);
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;
66 const VertexType* vert;
71 log.log(0,
"Updating per vertex normals...");
73 updatePerVertexNormalsAngleWeighted(m);
74 normalizePerVertexNormals(m);
76 log.log(5,
"Computing per vertex curvature...");
78 log.startProgress(
"", m.vertexNumber(), 5, 5, 100);
81 std::vector<ScalarType> weights;
84 MeshPos<FaceType> pos(v.adjFace(0), &v);
85 const VertexType* firstVertex = pos.adjVertex();
86 const VertexType* tmpVertex;
87 float totalDoubleAreaSize = 0;
92 pos.nextEdgeAdjacentToV();
93 tmpVertex = pos.adjVertex();
95 adjV.isBorder = pos.isEdgeOnBorder();
96 adjV.vert = tmpVertex;
97 adjV.doubleArea =
faceArea(*pos.face()) * 2;
98 totalDoubleAreaSize += adjV.doubleArea;
100 }
while (tmpVertex != firstVertex);
104 for (
int i = 0; i <
vertices.size(); ++i) {
105 if (vertices[i].isBorder) {
106 weights.push_back(vertices[i].doubleArea / totalDoubleAreaSize);
111 (vertices[i].doubleArea +
112 vertices[(i - 1) %
vertices.size()].doubleArea) /
113 totalDoubleAreaSize);
115 assert(weights.back() < 1.0f);
119 Matrix33<ScalarType> Tp;
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]);
130 Matrix33<ScalarType> tempMatrix;
132 for (
size_t i = 0; i <
vertices.size(); ++i) {
133 CoordType edge = (v.coord() -
vertices[i].vert->coord());
135 (2.0f * (v.normal().dot(edge))) / edge.squaredNorm();
136 CoordType t = Tp * edge;
138 tempMatrix = t.outerProduct(t);
139 M += tempMatrix * weights[i] * curvature;
144 CoordType e1(1.0f, 0.0f, 0.0f);
145 if ((e1 - v.normal()).squaredNorm() > (e1 + v.normal()).squaredNorm())
153 tempMatrix = w.outerProduct(w);
154 Q -= tempMatrix * 2.0f;
157 Matrix33<ScalarType> QtMQ = (Q.transpose() * M * Q);
159 Eigen::Matrix<ScalarType, 1, 3> T1 = Q.col(1);
160 Eigen::Matrix<ScalarType, 1, 3> T2 = Q.col(2);
165 ScalarType alpha = QtMQ(1, 1) - QtMQ(2, 2);
166 ScalarType beta = QtMQ(2, 1);
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);
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;
182 for (uint j = 0; j < 2; j++) {
183 ScalarType squaredT = std::pow(t[j], 2);
184 ScalarType denominator = 1.0f + squaredT;
186 s = (2.0f * t[j]) / denominator;
187 c = (1 - squaredT) / denominator;
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));
194 std::abs(1.0f - angleSimilarity) + std::abs(approximation);
195 if (error < minError) {
205 Eigen::Matrix2f minor2x2;
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);
214 S(0, 0) = S(1, 1) = c;
218 Eigen::Matrix2f StMS = S.transpose() * minor2x2 * S;
221 ScalarType principalCurv1 = (3.0f * StMS(0, 0)) - StMS(1, 1);
222 ScalarType principalCurv2 = (3.0f * StMS(1, 1)) - StMS(0, 0);
224 Eigen::Matrix<ScalarType, 1, 3> principalDir1 = T1 * c - T2 * s;
225 Eigen::Matrix<ScalarType, 1, 3> principalDir2 = T1 * s + T2 * c;
227 v.principalCurvature().maxDir() = principalDir1;
228 v.principalCurvature().minDir() = principalDir2;
229 v.principalCurvature().maxValue() = principalCurv1;
230 v.principalCurvature().minValue() = principalCurv2;
232 log.progress(m.index(v));
236 log.log(100,
"Per vertex curvature computed.");
251template<FaceMeshConcept MeshType, LoggerConcept LogType = NullLogger>
252void updatePrincipalCurvaturePCA(
254 typename MeshType::VertexType::CoordType::ScalarType radius,
255 bool montecarloSampling =
true,
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;
264 using VGrid = StaticGrid3<VertexType*, ScalarType>;
265 using VGridIterator = VGrid::ConstIterator;
270 log.log(0,
"Updating per vertex normals...");
272 updatePerVertexNormalsAngleWeighted(m);
273 normalizePerVertexNormals(m);
275 log.log(0,
"Computing per vertex curvature...");
276 log.startProgress(
"", m.vertexNumber());
278 if (montecarloSampling) {
279 area = surfaceArea(m);
284 parallelFor(m.vertices(), [&](VertexType& v) {
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());
296 A = covarianceMatrixOfPointCloud(points);
297 A *= area * area / 1000;
300 Sphere<ScalarType> sph(v.coord(), radius);
301 MeshType tmpMesh = intersection(m, sph);
303 A = covarianceMatrixOfMesh(tmpMesh);
306 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<ScalarType, 3, 3>> eig(A);
307 eigenvalues = CoordType(eig.eigenvalues());
308 eigenvectors = eig.eigenvectors();
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())));
323 v.principalCurvature().maxDir() =
324 (eigenvectors.col((best + 1) % 3).normalized());
325 v.principalCurvature().minDir() =
326 (eigenvectors.col((best + 2) % 3).normalized());
328 Matrix33<ScalarType> rot;
330 angle = acos(v.principalCurvature().maxDir().dot(v.normal()));
332 rot = rotationMatrix<Matrix33<ScalarType>>(
333 CoordType(v.principalCurvature().maxDir().cross(v.normal())),
334 -(M_PI * 0.5 - angle));
336 v.principalCurvature().maxDir() = rot * v.principalCurvature().maxDir();
338 angle = acos(v.principalCurvature().minDir().dot(v.normal()));
340 rot = rotationMatrix<Matrix33<ScalarType>>(
341 CoordType(v.principalCurvature().minDir().cross(v.normal())),
342 -(M_PI * 0.5 - angle));
344 v.principalCurvature().minDir() = rot * v.principalCurvature().minDir();
347 const ScalarType r5 = std::pow(radius, 5);
348 const ScalarType r6 = r5 * radius;
349 v.principalCurvature().maxValue() =
351 (4.0 * M_PI * r5 + 15 * eigenvalues[(best + 2) % 3] -
352 45.0 * eigenvalues[(best + 1) % 3]) /
354 v.principalCurvature().minValue() =
356 (4.0 * M_PI * r5 + 15 * eigenvalues[(best + 1) % 3] -
357 45.0 * eigenvalues[(best + 2) % 3]) /
359 if (v.principalCurvature().maxValue() <
360 v.principalCurvature().minValue()) {
362 v.principalCurvature().minValue(),
363 v.principalCurvature().maxValue());
365 v.principalCurvature().minDir(),
366 v.principalCurvature().maxDir());
369 log.progress(m.index(v));
374 log.log(100,
"Per vertex curvature computed.");
377template<FaceMeshConcept MeshType, LoggerConcept LogType = NullLogger>
378void updatePrincipalCurvature(MeshType& m, LogType& log = nullLogger)
380 requirePerVertexPrincipalCurvature(m);
382 updatePrincipalCurvatureTaubin95(m, log);
385template<FaceMeshConcept MeshType, LoggerConcept LogType = NullLogger>
386void updatePrincipalCurvature(
388 VCLibPrincipalCurvatureAlgorithm alg = VCL_PRINCIPAL_CURVATURE_TAUBIN95,
389 LogType& log = nullLogger)
391 requirePerVertexPrincipalCurvature(m);
395 case VCL_PRINCIPAL_CURVATURE_TAUBIN95:
396 updatePrincipalCurvatureTaubin95(m, log);
398 case VCL_PRINCIPAL_CURVATURE_PCA:
400 updatePrincipalCurvaturePCA(m, radius,
true, log);
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