Finalize the computation of the properties of a torus

diff --git a/src/ShapeRecogn/Areas.cxx b/src/ShapeRecogn/Areas.cxx
index a9b1e0dff7bfd1532bc0b028b47769eb5c7f48d9..f8441cb48df6e663626f647e763fcc056beb9904 100644 (file)
--- a/src/ShapeRecogn/Areas.cxx
+++ b/src/ShapeRecogn/Areas.cxx
@@ -2,6 +2,7 @@
 #include "MathOps.hxx"
 
 #include <random>
+#include <cblas.h>
 
 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<double> areaNodesK1(n, 0.0);
     std::vector<double> 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<double> majorRadiusNodes(3 * n, 0.0);
     for (size_t i = 0; i < n; ++i)
     {
         std::array<double, 3> coords = nodes->getCoordinates(area.nodeIds[i]);
         std::array<double, 3> 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<double, 3> 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<double, 3> 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<double, 6> base2d = MathOps::computeBaseFromNormal(normal);
+    std::vector<double> 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<double> A(3 * n, 1.0);
+    std::vector<double> 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<double> 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<mcIdType> &Areas::getAreaIdByNodes() const

diff --git a/src/ShapeRecogn/Areas.hxx b/src/ShapeRecogn/Areas.hxx
index b8c765594843c6b85b97fbf1e14cde99706a3f59..b1b6f6bca438b5fc81e4aed5e985ba599be53da2 100644 (file)
--- 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<double, 3> 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<double, 3> &getNormal(mcIdType areaId) const;

diff --git a/src/ShapeRecogn/MathOps.cxx b/src/ShapeRecogn/MathOps.cxx
index 910b8cbe664f7ffbc7662d3bb73cc2ddea548b6d..8f1705ace4d49e69732e4929858a1c876f4b5c04 100644 (file)
--- a/src/ShapeRecogn/MathOps.cxx
+++ b/src/ShapeRecogn/MathOps.cxx
@@ -37,11 +37,15 @@ std::vector<double> MathOps::lstsq(
 }
 
 std::vector<double> MathOps::lstsq(
-    std::vector<double> &a, const std::vector<double> &b, int m, int n, int nrhs)
+    std::vector<double> &a,
+    const std::vector<double> &b,
+    int m, int n, int nrhs)
 {
     int ldb = std::max<int>(m, n);
     int lds = std::min<int>(m, n);
-    std::vector<double> x(b.begin(), b.end());
+    std::vector<double> x(ldb, 0.0);
+    for (size_t i = 0; i < b.size(); ++i)
+        x[i] = b[i];
     std::vector<double> 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<double> MathOps::lstsq(
     return x;
 }
 
+std::vector<double> MathOps::lstsqRow(
+    std::vector<double> &a,
+    const std::vector<double> &b)
+{
+    int m = b.size();
+    int n = 3;
+    int nrhs = 1;
+    int ldb = std::max<int>(m, n);
+    int lds = std::min<int>(m, n);
+    std::vector<double> x(ldb, 0.0);
+    for (size_t i = 0; i < b.size(); ++i)
+        x[i] = b[i];
+    std::vector<double> 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<double, 3> MathOps::cross(const std::array<double, 3> &a, const std::array<double, 3> &b)
 {
     std::array<double, 3> c{0.0, 0.0, 0.0};
@@ -161,13 +190,13 @@ std::array<double, 9> 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<double, 3> &normal, const
         return -angle;
 }
 
-double MathOps::computeVar(std::vector<double> values)
+double MathOps::computeVariance(std::vector<double> 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<double, 6> MathOps::computeBaseFromNormal(std::array<double, 3> normal)
+{
+    std::array<double, 3> n_normal = normalize(normal);
+    std::array<double, 3> s;
+    std::array<double, 1> u;
+    std::array<double, 9> v;
+    LAPACKE_dgesdd(
+        LAPACK_COL_MAJOR,
+        'A', 1, 3, n_normal.data(), 1,
+        s.data(),
+        u.data(), 1,
+        v.data(), 3);
+    std::array<double, 3> v1 = {v[3], v[4], v[5]};
+    std::array<double, 3> v2 = {v[6], v[7], v[8]};
+    std::array<double, 3> u1 = normalize(v1);
+    std::array<double, 3> 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 267f9170d40a46d197e427685febbf794157898b..e511845da16639aa48ac30f69ef68cc3541b429a 100644 (file)
--- 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<double> lstsqRow(
+            std::vector<double> &a,
+            const std::vector<double> &b);
         static std::array<double, 3> cross(const std::array<double, 3> &a, const std::array<double, 3> &b);
         static std::array<double, 3> normalize(const std::array<double, 3> &a);
         static double computeNorm(const std::array<double, 3> &a);
@@ -60,7 +63,8 @@ namespace MEDCoupling
             const std::array<double, 3> &normal,
             const std::array<double, 3> &vector1,
             const std::array<double, 3> &vector2);
-        static double computeVar(std::vector<double> values);
+        static double computeVariance(std::vector<double> values);
+        static std::array<double, 6> computeBaseFromNormal(std::array<double, 3> normal);
     };
 };
 

diff --git a/src/ShapeRecogn/Test/MathOpsTest.cxx b/src/ShapeRecogn/Test/MathOpsTest.cxx
index 17bef13fc989099cae6664ef8f3e698e42034bcf..033f9921b763c0327eee131a85f046020da0cf7e 100644 (file)
--- 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<double, 3> normal = {1.0, 2.0, 3.0};
+    std::array<double, 6> base = MathOps::computeBaseFromNormal(normal);
+    std::array<double, 6> 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 4ff787f543029db65c6d806f129accb3bbb05f5d..2988a950be40784d4b976052596b664ce75f046e 100644 (file)
--- 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();
     };
 };