From f9db763cf21ae8aafccbba3fbc1ed65fce4cd017 Mon Sep 17 00:00:00 2001 From: El Hadi Moussi Date: Wed, 7 Aug 2024 18:41:18 +0200 Subject: [PATCH] Finalize the computation of the properties of a torus --- src/ShapeRecogn/Areas.cxx | 61 ++++++++++++++++------- src/ShapeRecogn/Areas.hxx | 2 + src/ShapeRecogn/MathOps.cxx | 73 +++++++++++++++++++++++++--- src/ShapeRecogn/MathOps.hxx | 6 ++- src/ShapeRecogn/Test/MathOpsTest.cxx | 15 ++++++ src/ShapeRecogn/Test/MathOpsTest.hxx | 2 + 6 files changed, 133 insertions(+), 26 deletions(-) diff --git a/src/ShapeRecogn/Areas.cxx b/src/ShapeRecogn/Areas.cxx index a9b1e0dff..f8441cb48 100644 --- a/src/ShapeRecogn/Areas.cxx +++ b/src/ShapeRecogn/Areas.cxx @@ -2,6 +2,7 @@ #include "MathOps.hxx" #include +#include using namespace MEDCoupling; @@ -162,6 +163,12 @@ void Areas::computeProperties(mcIdType areaId) } } +double Areas::getMinorRadius(mcIdType areaId) const +{ + const Area &area = areas[areaId]; + return area.minorRadius; +} + double Areas::getRadius(mcIdType areaId) const { const Area &area = areas[areaId]; @@ -219,6 +226,7 @@ void Areas::cleanArea(mcIdType areaId, mcIdType newAreaId = -1) area.adimK1 = 0.0; area.adimK2 = 0.0; area.adimKdiff0 = 0.0; + area.minorRadius = 0.0; area.radius = 0.0; area.angle = 0.0; area.normal = {0.0, 0.0, 0.0}; @@ -408,6 +416,8 @@ void Areas::computeTorusProperties(mcIdType areaId) { Area &area = areas[areaId]; size_t n = area.nodeIds.size(); + if (n == 0) + return; std::vector areaNodesK1(n, 0.0); std::vector areaNodesK2(n, 0.0); for (size_t i = 0; i < n; ++i) @@ -415,21 +425,21 @@ void Areas::computeTorusProperties(mcIdType areaId) areaNodesK1[i] = nodes->getK1(area.nodeIds[i]); areaNodesK2[i] = nodes->getK2(area.nodeIds[i]); } - double var1 = MathOps::computeVar(areaNodesK1); - double var2 = MathOps::computeVar(areaNodesK2); + double var1 = MathOps::computeVariance(areaNodesK1); + double var2 = MathOps::computeVariance(areaNodesK2); double minorCurvature; if (var1 > var2) minorCurvature = MathOps::mean(areaNodesK2); else minorCurvature = MathOps::mean(areaNodesK1); - double minorRadius = 1.0 / minorCurvature; + area.minorRadius = 1.0 / minorCurvature; std::vector majorRadiusNodes(3 * n, 0.0); for (size_t i = 0; i < n; ++i) { std::array coords = nodes->getCoordinates(area.nodeIds[i]); std::array normal = nodes->getNormal(area.nodeIds[i]); for (size_t j = 0; j < 3; ++j) - majorRadiusNodes[3 * i + j] = coords[j] - normal[j] * minorRadius; + majorRadiusNodes[3 * i + j] = coords[j] - normal[j] * area.minorRadius; } std::array meanMajorRadiusNodes = MathOps::meanCoordinates(majorRadiusNodes); for (size_t i = 0; i < n; ++i) @@ -438,20 +448,35 @@ void Areas::computeTorusProperties(mcIdType areaId) majorRadiusNodes[3 * i + j] -= meanMajorRadiusNodes[j]; } std::array normal = MathOps::computePCAThirdAxis(majorRadiusNodes); - /* - u, v = get_normal_plane_vectors(normal) - p_matrix = np.array([u, v]).T - p_major_radius_points = major_radius_points @ p_matrix - centre_2D, self.major_radius = fitCircle2D(p_major_radius_points) - self.centre = centre_2D[0] * u + centre_2D[1] * v + P_mean - self.parameters.update( - { - "minor_radius": self.minor_radius, - "major_radius": self.major_radius, - "centre": self.centre, - } - ) - */ + std::array base2d = MathOps::computeBaseFromNormal(normal); + std::vector projectedMajorRadiusNodes(2 * n, 0.0); + cblas_dgemm( + CBLAS_LAYOUT::CblasRowMajor, + CBLAS_TRANSPOSE::CblasNoTrans, + CBLAS_TRANSPOSE::CblasTrans, + (int)n, 2, 3, + 1.0, + majorRadiusNodes.data(), 3, + base2d.data(), 3, + 0.0, projectedMajorRadiusNodes.data(), 2); + std::vector A(3 * n, 1.0); + std::vector B(n, 0.0); + for (size_t i = 0; i < n; ++i) + { + for (size_t j = 0; j < 2; ++j) + A[3 * i + j] = projectedMajorRadiusNodes[2 * i + j]; + B[i] = pow(projectedMajorRadiusNodes[2 * i], 2) + + pow(projectedMajorRadiusNodes[2 * i + 1], 2); + } + std::vector fit = MathOps::lstsqRow(A, B); + double a = fit[0]; + double b = fit[1]; + double c = fit[2]; + double xc = a / 2.0; + double yc = b / 2.0; + area.radius = sqrt(4.0 * c + pow(a, 2) + pow(b, 2)) / 2.0; + for (size_t i = 0; i < 3; ++i) + area.center[i] = xc * base2d[i] + yc * base2d[3 + i] + meanMajorRadiusNodes[i]; } const std::vector &Areas::getAreaIdByNodes() const diff --git a/src/ShapeRecogn/Areas.hxx b/src/ShapeRecogn/Areas.hxx index b8c765594..b1b6f6bca 100644 --- a/src/ShapeRecogn/Areas.hxx +++ b/src/ShapeRecogn/Areas.hxx @@ -36,6 +36,7 @@ namespace MEDCoupling double adimK1 = 0.0; double adimK2 = 0.0; double adimKdiff0 = 0.0; + double minorRadius = 0.0; double radius = 0.0; double angle = 0.0; std::array normal{0.0, 0.0, 0.0}; @@ -78,6 +79,7 @@ namespace MEDCoupling double getAdimKdiff0(mcIdType areaId) const; void computeProperties(mcIdType areaId); + double getMinorRadius(mcIdType areaId) const; double getRadius(mcIdType areaId) const; double getAngle(mcIdType areaId) const; const std::array &getNormal(mcIdType areaId) const; diff --git a/src/ShapeRecogn/MathOps.cxx b/src/ShapeRecogn/MathOps.cxx index 910b8cbe6..8f1705ace 100644 --- a/src/ShapeRecogn/MathOps.cxx +++ b/src/ShapeRecogn/MathOps.cxx @@ -37,11 +37,15 @@ std::vector MathOps::lstsq( } std::vector MathOps::lstsq( - std::vector &a, const std::vector &b, int m, int n, int nrhs) + std::vector &a, + const std::vector &b, + int m, int n, int nrhs) { int ldb = std::max(m, n); int lds = std::min(m, n); - std::vector x(b.begin(), b.end()); + std::vector x(ldb, 0.0); + for (size_t i = 0; i < b.size(); ++i) + x[i] = b[i]; std::vector s(lds, 0.0); double rcond = DBL_EPSILON * (double)ldb; // same value as numpy.linalg.lstsq int rank = 0; @@ -56,6 +60,31 @@ std::vector MathOps::lstsq( return x; } +std::vector MathOps::lstsqRow( + std::vector &a, + const std::vector &b) +{ + int m = b.size(); + int n = 3; + int nrhs = 1; + int ldb = std::max(m, n); + int lds = std::min(m, n); + std::vector x(ldb, 0.0); + for (size_t i = 0; i < b.size(); ++i) + x[i] = b[i]; + std::vector s(lds, 0.0); + double rcond = DBL_EPSILON * (double)ldb; // same value as numpy.linalg.lstsq + int rank = 0; + int info = LAPACKE_dgelsd( + LAPACK_ROW_MAJOR, + m, n, nrhs, + a.data(), n, + x.data(), nrhs, + s.data(), + rcond, + &rank); + return x; +} std::array MathOps::cross(const std::array &a, const std::array &b) { std::array c{0.0, 0.0, 0.0}; @@ -161,13 +190,13 @@ std::array MathOps::computeCov( CBLAS_LAYOUT::CblasColMajor, CBLAS_TRANSPOSE::CblasNoTrans, CBLAS_TRANSPOSE::CblasTrans, - 3, 3, nbNodes, 1.0, + 3, 3, (int)nbNodes, 1.0, coordsCentered.data(), 3, coordsCentered.data(), 3, 0.0, covMatrix.data(), 3); for (size_t i = 0; i < 9; ++i) { - covMatrix[i] /= nbNodes - 1; + covMatrix[i] /= (double)nbNodes - 1; } return covMatrix; } @@ -261,12 +290,42 @@ double MathOps::computeOrientedAngle(const std::array &normal, const return -angle; } -double MathOps::computeVar(std::vector values) +double MathOps::computeVariance(std::vector values) { - int n = values.size(); + size_t n = values.size(); double m = mean(values); double d2 = 0; for (size_t i = 0; i < n; ++i) d2 += pow(values[i] - m, 2); - return d2 / n; + return d2 / (double)n; +} + +std::array MathOps::computeBaseFromNormal(std::array normal) +{ + std::array n_normal = normalize(normal); + std::array s; + std::array u; + std::array v; + LAPACKE_dgesdd( + LAPACK_COL_MAJOR, + 'A', 1, 3, n_normal.data(), 1, + s.data(), + u.data(), 1, + v.data(), 3); + std::array v1 = {v[3], v[4], v[5]}; + std::array v2 = {v[6], v[7], v[8]}; + std::array u1 = normalize(v1); + std::array u2; + u2.fill(0.0); + double innerProd = dot(u1, v2); + for (size_t i = 0; i < 3; ++i) + u2[i] = v2[i] - innerProd * u1[i]; + u2 = normalize(u2); + double sign = dot(cross(u1, u2), normal); + if (sign < 0) + { + for (size_t i = 0; i < 3; ++i) + u1[i] *= -1; + } + return {u1[0], u1[1], u1[2], u2[0], u2[1], u2[2]}; } diff --git a/src/ShapeRecogn/MathOps.hxx b/src/ShapeRecogn/MathOps.hxx index 267f9170d..e511845da 100644 --- a/src/ShapeRecogn/MathOps.hxx +++ b/src/ShapeRecogn/MathOps.hxx @@ -35,6 +35,9 @@ namespace MEDCoupling int m, int n = 3, int nrhs = 1); + static std::vector lstsqRow( + std::vector &a, + const std::vector &b); static std::array cross(const std::array &a, const std::array &b); static std::array normalize(const std::array &a); static double computeNorm(const std::array &a); @@ -60,7 +63,8 @@ namespace MEDCoupling const std::array &normal, const std::array &vector1, const std::array &vector2); - static double computeVar(std::vector values); + static double computeVariance(std::vector values); + static std::array computeBaseFromNormal(std::array normal); }; }; diff --git a/src/ShapeRecogn/Test/MathOpsTest.cxx b/src/ShapeRecogn/Test/MathOpsTest.cxx index 17bef13fc..033f9921b 100644 --- a/src/ShapeRecogn/Test/MathOpsTest.cxx +++ b/src/ShapeRecogn/Test/MathOpsTest.cxx @@ -80,3 +80,18 @@ void MathOpsTest::testComputeAngles() CPPUNIT_ASSERT_DOUBLES_EQUAL(1.27845478, angles[2], 1E-6); CPPUNIT_ASSERT_DOUBLES_EQUAL(2.61880376, angles[3], 1E-6); } + +void MathOpsTest::testComputeBaseFromNormal() +{ + std::array normal = {1.0, 2.0, 3.0}; + std::array base = MathOps::computeBaseFromNormal(normal); + std::array baseRef = { + -0.53452248, 0.77454192, -0.33818712, + -0.80178373, -0.33818712, 0.49271932}; + for (size_t i = 0; i < 6; ++i) + { + std::ostringstream message; + message << "Mismatch at index " << i; + CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(message.str().c_str(), baseRef[i], base[i], 1E-6); + } +} diff --git a/src/ShapeRecogn/Test/MathOpsTest.hxx b/src/ShapeRecogn/Test/MathOpsTest.hxx index 4ff787f54..2988a950b 100644 --- a/src/ShapeRecogn/Test/MathOpsTest.hxx +++ b/src/ShapeRecogn/Test/MathOpsTest.hxx @@ -34,6 +34,7 @@ namespace MEDCoupling CPPUNIT_TEST(testComputeCov); CPPUNIT_TEST(testComputePCAFirstAxis); CPPUNIT_TEST(testComputeAngles); + CPPUNIT_TEST(testComputeBaseFromNormal); CPPUNIT_TEST_SUITE_END(); public: @@ -43,6 +44,7 @@ namespace MEDCoupling static void testComputeCov(); static void testComputePCAFirstAxis(); static void testComputeAngles(); + static void testComputeBaseFromNormal(); }; }; -- 2.39.2