Store adimensional principal curvatures and use them in AreasBuilder

diff --git a/src/ShapeRecogn/Areas.cxx b/src/ShapeRecogn/Areas.cxx
index cd4709e31693b8058ea6dccdcb40bd1f1e6dabc1..a9b1e0dff7bfd1532bc0b028b47769eb5c7f48d9 100644 (file)
--- a/src/ShapeRecogn/Areas.cxx
+++ b/src/ShapeRecogn/Areas.cxx
@@ -50,7 +50,9 @@ void Areas::addNode(mcIdType areaId, mcIdType nodeId)
     size_t nbNodes = area.nodeIds.size();
     area.k1 = ((double)(nbNodes - 1) * area.k1 + nodes->getK1(nodeId)) / (double)nbNodes;
     area.k2 = ((double)(nbNodes - 1) * area.k2 + nodes->getK2(nodeId)) / (double)nbNodes;
-    area.kdiff0 = ((double)(nbNodes - 1) * area.kdiff0 + nodes->getKdiff0(nodeId)) / (double)nbNodes;
+    area.adimK1 = ((double)(nbNodes - 1) * area.adimK1 + nodes->getAdimK1(nodeId)) / (double)nbNodes;
+    area.adimK2 = ((double)(nbNodes - 1) * area.adimK2 + nodes->getAdimK2(nodeId)) / (double)nbNodes;
+    area.adimKdiff0 = ((double)(nbNodes - 1) * area.adimKdiff0 + nodes->getAdimKdiff0(nodeId)) / (double)nbNodes;
 }
 
 void Areas::cleanInvalidNodeAreas()
@@ -119,20 +121,20 @@ const std::vector<mcIdType> &Areas::getNodeIds(mcIdType areaId) const
     return areas[areaId].nodeIds;
 }
 
-double Areas::getK1(mcIdType areaId) const
+double Areas::getAdimK1(mcIdType areaId) const
 {
 
-    return areas[areaId].k1;
+    return areas[areaId].adimK1;
 }
 
-double Areas::getK2(mcIdType areaId) const
+double Areas::getAdimK2(mcIdType areaId) const
 {
-    return areas[areaId].k2;
+    return areas[areaId].adimK2;
 }
 
-double Areas::getKdiff0(mcIdType areaId) const
+double Areas::getAdimKdiff0(mcIdType areaId) const
 {
-    return areas[areaId].kdiff0;
+    return areas[areaId].adimKdiff0;
 }
 
 void Areas::computeProperties(mcIdType areaId)
@@ -214,7 +216,9 @@ void Areas::cleanArea(mcIdType areaId, mcIdType newAreaId = -1)
     area.primitive = PrimitiveType::Unknown;
     area.k1 = 0.0;
     area.k2 = 0.0;
-    area.kdiff0 = 0.0;
+    area.adimK1 = 0.0;
+    area.adimK2 = 0.0;
+    area.adimKdiff0 = 0.0;
     area.radius = 0.0;
     area.angle = 0.0;
     area.normal = {0.0, 0.0, 0.0};

diff --git a/src/ShapeRecogn/Areas.hxx b/src/ShapeRecogn/Areas.hxx
index e0d8ae440705ecfc63b0488441e8f0b929e0f41a..b8c765594843c6b85b97fbf1e14cde99706a3f59 100644 (file)
--- a/src/ShapeRecogn/Areas.hxx
+++ b/src/ShapeRecogn/Areas.hxx
@@ -33,7 +33,9 @@ namespace MEDCoupling
         PrimitiveType primitive = PrimitiveType::Unknown;
         double k1 = 0.0;
         double k2 = 0.0;
-        double kdiff0 = 0.0;
+        double adimK1 = 0.0;
+        double adimK2 = 0.0;
+        double adimKdiff0 = 0.0;
         double radius = 0.0;
         double angle = 0.0;
         std::array<double, 3> normal{0.0, 0.0, 0.0};
@@ -71,9 +73,9 @@ namespace MEDCoupling
 
         const std::vector<mcIdType> &getNodeIds(mcIdType areaId) const;
 
-        double getK1(mcIdType areaId) const;
-        double getK2(mcIdType areaId) const;
-        double getKdiff0(mcIdType areaId) const;
+        double getAdimK1(mcIdType areaId) const;
+        double getAdimK2(mcIdType areaId) const;
+        double getAdimKdiff0(mcIdType areaId) const;
 
         void computeProperties(mcIdType areaId);
         double getRadius(mcIdType areaId) const;

diff --git a/src/ShapeRecogn/AreasBuilder.cxx b/src/ShapeRecogn/AreasBuilder.cxx
index e435c5d1e3baafa0160a45c76edc0b1ce75a2963..ef5d6ea20e9b7a403857d6623fcecbd1317562b1 100644 (file)
--- a/src/ShapeRecogn/AreasBuilder.cxx
+++ b/src/ShapeRecogn/AreasBuilder.cxx
@@ -261,16 +261,15 @@ bool AreasBuilder::doesItMatch(mcIdType areaId, mcIdType nodeId) const
         case PrimitiveType::Sphere:
         {
             double kmoy = fabs(
-                (areas->getK1(areaId) + areas->getK2(areaId)) / 2.0);
+                (areas->getAdimK1(areaId) + areas->getAdimK2(areaId)) / 2.0);
             double nodeKmoy = fabs(
-                (nodes->getK1(nodeId) + nodes->getK2(nodeId)) / 2.0);
-            isMatching = fabs((nodeKmoy - kmoy) / kmoy) < TOL_MATCH_SPHERE;
+                (nodes->getAdimK1(nodeId) + nodes->getAdimK2(nodeId)) / 2.0);
+            isMatching = fabs((nodeKmoy - kmoy)) < TOL_MATCH_SPHERE;
         }
         break;
         case PrimitiveType::Cylinder:
-            isMatching =
-                fabs((areas->getKdiff0(areaId) - nodes->getKdiff0(nodeId)) /
-                     areas->getKdiff0(areaId)) < EPSILON_PROP;
+            isMatching = fabs(areas->getAdimKdiff0(areaId) -
+                              nodes->getAdimKdiff0(nodeId)) < EPSILON_PROP;
             break;
         case PrimitiveType::Cone:
         case PrimitiveType::Unknown:

diff --git a/src/ShapeRecogn/NodeCurvatureCalculator.cxx b/src/ShapeRecogn/NodeCurvatureCalculator.cxx
index 7787240862a0dbaa652035179c9e021c8a7cae57..6a2f6528ad7ba967e73e7ff028bd858e6eb48834 100644 (file)
--- a/src/ShapeRecogn/NodeCurvatureCalculator.cxx
+++ b/src/ShapeRecogn/NodeCurvatureCalculator.cxx
@@ -77,9 +77,8 @@ void NodeCurvatureCalculator::computeCurvatures(double tol)
     mcIdType nbNodes = mesh->getNumberOfNodes();
     nodes->k1.resize(nbNodes);
     nodes->k2.resize(nbNodes);
-    nodes->theta0.resize(nbNodes);
-    nodes->isConvex.resize(nbNodes);
-    nodes->kdiff0.resize(nbNodes);
+    nodes->adimK1.resize(nbNodes);
+    nodes->adimK2.resize(nbNodes);
     nodes->primitives.resize(nbNodes);
     nodes->weakDirections.resize(3 * nbNodes);
     nodes->mainDirections.resize(3 * nbNodes);
@@ -92,11 +91,10 @@ void NodeCurvatureCalculator::computeCurvatures(mcIdType nodeId, double tol)
     std::array<double, 3> normal = nodes->getNormal(nodeId);
     std::vector<mcIdType> neighborIds = nodes->getNeighbors(nodeId);
     double theta0 = 0.0;
-    bool isConvex = false;
     double k1 = 0.0;
     double k2 = 0.0;
-    double kdiff0 = 0.0;
-    double kis0 = 0.0;
+    double adimK1 = 0.0;
+    double adimK2 = 0.0;
     PrimitiveType primitive = PrimitiveType::Unknown;
     std::array<double, 3> mainDirection{0.0, 0.0, 0.0};
     std::array<double, 3> weakDirection{0.0, 0.0, 0.0};
@@ -128,31 +126,31 @@ void NodeCurvatureCalculator::computeCurvatures(mcIdType nodeId, double tol)
             direction1[i] = cos(theta0) * e1[i] + sin(theta0) * e2[i];
             direction2[i] = -sin(theta0) * e1[i] + cos(theta0) * e2[i];
         }
+        double averageDistance = computeAverageDistance(nodeId, neighborIds);
+        double adimK1 = k1 * averageDistance;
+        double adimK2 = k2 * averageDistance;
+        double adimKdiff0, adimKis0;
         if (fabs(k1) > fabs(k2))
         {
-            isConvex = true;
-            kdiff0 = k1;
-            kis0 = k2;
+            adimKdiff0 = adimK1;
+            adimKis0 = adimK2;
             mainDirection = direction1;
             weakDirection = direction2;
         }
         else
         {
-            isConvex = false;
-            kdiff0 = k2;
-            kis0 = k1;
+            adimKdiff0 = adimK2;
+            adimKis0 = adimK1;
             mainDirection = direction2;
             weakDirection = direction1;
         }
-        double averageDistance = computeAverageDistance(nodeId, neighborIds);
-        primitive = findPrimitiveType(k1 * averageDistance, k2 * averageDistance, kdiff0 * averageDistance, kis0 * averageDistance);
+        primitive = findPrimitiveType(adimK1, adimK2, adimKdiff0, adimKis0);
     }
     // Populate nodes
-    nodes->isConvex[nodeId] = isConvex;
-    nodes->theta0[nodeId] = theta0;
     nodes->k1[nodeId] = k1;
     nodes->k2[nodeId] = k2;
-    nodes->kdiff0[nodeId] = kdiff0;
+    nodes->adimK1[nodeId] = adimK1;
+    nodes->adimK2[nodeId] = adimK2;
     for (size_t i = 0; i < 3; ++i)
     {
         nodes->weakDirections[3 * nodeId + i] = weakDirection[i];

diff --git a/src/ShapeRecogn/Nodes.cxx b/src/ShapeRecogn/Nodes.cxx
index 31610b3d80c82bbb2b7a550b7d207753f0ba0de3..d1dd810c45bc8aa99ebd266f906f1da70e1abc06 100644 (file)
--- a/src/ShapeRecogn/Nodes.cxx
+++ b/src/ShapeRecogn/Nodes.cxx
@@ -37,11 +37,6 @@ const std::vector<double> &Nodes::getMainDirections() const
     return mainDirections;
 }
 
-const std::vector<double> &Nodes::getK1() const
-{
-    return k1;
-}
-
 std::array<double, 3> Nodes::getNormal(mcIdType nodeId) const
 {
     return {normals[3 * nodeId], normals[3 * nodeId + 1], normals[3 * nodeId + 2]};
@@ -59,7 +54,22 @@ double Nodes::getK2(mcIdType nodeId) const
 
 double Nodes::getKdiff0(mcIdType nodeId) const
 {
-    return kdiff0[nodeId];
+    return fabs(k1[nodeId]) > fabs(k2[nodeId]) ? k1[nodeId] : k2[nodeId];
+}
+
+double Nodes::getAdimK1(mcIdType nodeId) const
+{
+    return adimK1[nodeId];
+}
+
+double Nodes::getAdimK2(mcIdType nodeId) const
+{
+    return adimK2[nodeId];
+}
+
+double Nodes::getAdimKdiff0(mcIdType nodeId) const
+{
+    return fabs(adimK1[nodeId]) > fabs(adimK2[nodeId]) ? adimK1[nodeId] : adimK2[nodeId];
 }
 
 const std::vector<mcIdType> Nodes::getNeighbors(mcIdType nodeId) const

diff --git a/src/ShapeRecogn/Nodes.hxx b/src/ShapeRecogn/Nodes.hxx
index ef06d1540d8d9fadcf81ed87c33dd74fade8134b..273d757f7b7ab8c218e6c62c5a1ce1a0804f45be 100644 (file)
--- a/src/ShapeRecogn/Nodes.hxx
+++ b/src/ShapeRecogn/Nodes.hxx
@@ -39,12 +39,14 @@ namespace MEDCoupling
         const std::vector<double> &getNormals() const;
         const std::vector<double> &getWeakDirections() const;
         const std::vector<double> &getMainDirections() const;
-        const std::vector<double> &getK1() const;
 
         std::array<double, 3> getNormal(mcIdType nodeId) const;
         double getK1(mcIdType nodeId) const;
         double getK2(mcIdType nodeId) const;
         double getKdiff0(mcIdType nodeId) const;
+        double getAdimK1(mcIdType nodeId) const;
+        double getAdimK2(mcIdType nodeId) const;
+        double getAdimKdiff0(mcIdType nodeId) const;
         const std::vector<mcIdType> getNeighbors(mcIdType nodeId) const;
         PrimitiveType getPrimitiveType(mcIdType nodeId) const;
         std::array<double, 3> getWeakDirection(mcIdType nodeId) const;
@@ -60,11 +62,10 @@ namespace MEDCoupling
         const DataArrayInt64 *neighbors;
         const DataArrayInt64 *neighborsIdx;
         // curvature
-        std::vector<bool> isConvex;
-        std::vector<double> theta0;
         std::vector<double> k1;
         std::vector<double> k2;
-        std::vector<double> kdiff0;
+        std::vector<double> adimK1;
+        std::vector<double> adimK2;
         std::vector<double> weakDirections;
         std::vector<double> mainDirections;
         std::vector<PrimitiveType> primitives;