From 5d5b861db47d870d9321f267e1c7a4ddef433595 Mon Sep 17 00:00:00 2001 From: El Hadi Moussi Date: Wed, 7 Aug 2024 11:12:45 +0200 Subject: [PATCH] Store adimensional principal curvatures and use them in AreasBuilder --- src/ShapeRecogn/Areas.cxx | 20 +++++++------ src/ShapeRecogn/Areas.hxx | 10 ++++--- src/ShapeRecogn/AreasBuilder.cxx | 11 ++++--- src/ShapeRecogn/NodeCurvatureCalculator.cxx | 32 ++++++++++----------- src/ShapeRecogn/Nodes.cxx | 22 ++++++++++---- src/ShapeRecogn/Nodes.hxx | 9 +++--- 6 files changed, 59 insertions(+), 45 deletions(-) diff --git a/src/ShapeRecogn/Areas.cxx b/src/ShapeRecogn/Areas.cxx index cd4709e31..a9b1e0dff 100644 --- 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 &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 e0d8ae440..b8c765594 100644 --- 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 normal{0.0, 0.0, 0.0}; @@ -71,9 +73,9 @@ namespace MEDCoupling const std::vector &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 e435c5d1e..ef5d6ea20 100644 --- 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 778724086..6a2f6528a 100644 --- 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 normal = nodes->getNormal(nodeId); std::vector 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 mainDirection{0.0, 0.0, 0.0}; std::array 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 31610b3d8..d1dd810c4 100644 --- a/src/ShapeRecogn/Nodes.cxx +++ b/src/ShapeRecogn/Nodes.cxx @@ -37,11 +37,6 @@ const std::vector &Nodes::getMainDirections() const return mainDirections; } -const std::vector &Nodes::getK1() const -{ - return k1; -} - std::array 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 Nodes::getNeighbors(mcIdType nodeId) const diff --git a/src/ShapeRecogn/Nodes.hxx b/src/ShapeRecogn/Nodes.hxx index ef06d1540..273d757f7 100644 --- a/src/ShapeRecogn/Nodes.hxx +++ b/src/ShapeRecogn/Nodes.hxx @@ -39,12 +39,14 @@ namespace MEDCoupling const std::vector &getNormals() const; const std::vector &getWeakDirections() const; const std::vector &getMainDirections() const; - const std::vector &getK1() const; std::array 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 getNeighbors(mcIdType nodeId) const; PrimitiveType getPrimitiveType(mcIdType nodeId) const; std::array getWeakDirection(mcIdType nodeId) const; @@ -60,11 +62,10 @@ namespace MEDCoupling const DataArrayInt64 *neighbors; const DataArrayInt64 *neighborsIdx; // curvature - std::vector isConvex; - std::vector theta0; std::vector k1; std::vector k2; - std::vector kdiff0; + std::vector adimK1; + std::vector adimK2; std::vector weakDirections; std::vector mainDirections; std::vector primitives; -- 2.39.2