From: El Hadi Moussi Date: Thu, 8 Aug 2024 15:51:39 +0000 (+0200) Subject: Rename NodeCurvatureCalculator to NodesBuilder X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=bba7ee9b171ea00a424939ad3323d2b8f93cbc73;p=tools%2Fmedcoupling.git Rename NodeCurvatureCalculator to NodesBuilder --- diff --git a/src/ShapeRecogn/CMakeLists.txt b/src/ShapeRecogn/CMakeLists.txt index a4dc5717e..9c4b506d8 100644 --- a/src/ShapeRecogn/CMakeLists.txt +++ b/src/ShapeRecogn/CMakeLists.txt @@ -44,7 +44,7 @@ INCLUDE_DIRECTORIES( SET(shaperecogn_SOURCES MathOps.cxx Nodes.cxx - NodeCurvatureCalculator.cxx + NodesBuilder.cxx Areas.cxx AreasBuilder.cxx ShapeRecognMesh.cxx diff --git a/src/ShapeRecogn/NodeCurvatureCalculator.cxx b/src/ShapeRecogn/NodeCurvatureCalculator.cxx deleted file mode 100644 index e6e950d36..000000000 --- a/src/ShapeRecogn/NodeCurvatureCalculator.cxx +++ /dev/null @@ -1,285 +0,0 @@ - -#include "NodeCurvatureCalculator.hxx" -#include "Nodes.hxx" -#include "MathOps.hxx" -#include "ShapeRecongConstants.hxx" - -using namespace MEDCoupling; - -NodeCurvatureCalculator::NodeCurvatureCalculator( - const MEDCouplingUMesh *mesh) : mesh(mesh) -{ -} - -Nodes *NodeCurvatureCalculator::compute() -{ - DataArrayInt64 *neighbors; - DataArrayInt64 *neighborsIdx; - mesh->computeNeighborsOfNodes(neighbors, neighborsIdx); - nodes = new Nodes(mesh, neighbors, neighborsIdx); - computeNormals(); - computeCurvatures(); - return nodes; -} - -void NodeCurvatureCalculator::computeNormals() -{ - mcIdType nbNodes = mesh->getNumberOfNodes(); - mcIdType nbCells = mesh->getNumberOfCells(); - std::array cellNormal; - cellNormal.fill(0.0); - // product of normal AreaByCell - std::vector prodNormalAreaByCell(3 * nbCells, 0.0); - for (int cellId = 0; cellId < nbCells; cellId++) - { - std::vector nodeIds; - mesh->getNodeIdsOfCell(cellId, nodeIds); - computeCellNormal(nodeIds, cellNormal); - prodNormalAreaByCell[3 * cellId + 0] = cellNormal[0] / 2.0; - prodNormalAreaByCell[3 * cellId + 1] = cellNormal[1] / 2.0; - prodNormalAreaByCell[3 * cellId + 2] = cellNormal[2] / 2.0; - } - // - nodes->normals.resize(3 * nbNodes, 1.0); - DataArrayInt64 *revNodal = DataArrayInt64::New(); - DataArrayInt64 *revNodalIdx = DataArrayInt64::New(); - mesh->getReverseNodalConnectivity(revNodal, revNodalIdx); - for (mcIdType nodeId = 0; nodeId < nbNodes; nodeId++) - { - int nbCells = revNodalIdx->getIJ(nodeId + 1, 0) - revNodalIdx->getIJ(nodeId, 0); - std::vector cellIds(nbCells, 0); - int start = revNodalIdx->getIJ(nodeId, 0); - for (size_t i = 0; i < cellIds.size(); ++i) - cellIds[i] = revNodal->getIJ(start + i, 0); - double normal = 0.0; - for (size_t i = 0; i < 3; i++) - { - nodes->normals[3 * nodeId + i] = 0.0; - for (mcIdType j = 0; j < nbCells; j++) - { - nodes->normals[3 * nodeId + i] += - prodNormalAreaByCell[3 * cellIds[j] + i]; - } - nodes->normals[3 * nodeId + i] /= (double)nbCells; - normal += pow(nodes->normals[3 * nodeId + i], 2); - } - for (size_t i = 0; i < 3; i++) - nodes->normals[3 * nodeId + i] /= sqrt(normal); - } - revNodal->decrRef(); - revNodalIdx->decrRef(); -} - -void NodeCurvatureCalculator::computeCurvatures(double tol) -{ - mcIdType nbNodes = mesh->getNumberOfNodes(); - nodes->k1.resize(nbNodes); - nodes->k2.resize(nbNodes); - nodes->adimK1.resize(nbNodes); - nodes->adimK2.resize(nbNodes); - nodes->primitives.resize(nbNodes); - nodes->weakDirections.resize(3 * nbNodes); - nodes->mainDirections.resize(3 * nbNodes); - for (mcIdType nodeId = 0; nodeId < nbNodes; nodeId++) - computeCurvatures(nodeId, tol); -} - -void NodeCurvatureCalculator::computeCurvatures(mcIdType nodeId, double tol) -{ - std::array normal = nodes->getNormal(nodeId); - std::vector neighborIds = nodes->getNeighbors(nodeId); - double theta0 = 0.0; - double k1 = 0.0; - double k2 = 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}; - if (neighborIds.size() > 0) - { - std::vector discreteCurvatures = computeDiscreteCurvatures(nodeId, neighborIds); - std::vector tangents = computeTangentDirections(nodeId, neighborIds); - mcIdType maxCurvatureId = std::distance( - discreteCurvatures.begin(), - std::max_element(discreteCurvatures.begin(), discreteCurvatures.end())); - std::array e1{ - tangents[3 * maxCurvatureId], - tangents[3 * maxCurvatureId + 1], - tangents[3 * maxCurvatureId + 2]}; - std::array e2 = MathOps::normalize(MathOps::cross(e1, normal)); - std::vector coeffs = computeNormalCurvatureCoefficients( - discreteCurvatures, tangents, normal, e1); - double a = coeffs[0], b = coeffs[1], c = coeffs[2]; - double h = (a + c) / 2.0; - double kg = a * c - pow(b, 2) / 4.0; - k1 = h + sqrt(fabs(pow(h, 2) - kg)); - k2 = h - sqrt(fabs(pow(h, 2) - kg)); - if (fabs(k1 - k2) >= tol) - theta0 = 0.5 * asin(b / (k1 - k2)); - std::array direction1{0.0, 0.0, 0.0}; - std::array direction2{0.0, 0.0, 0.0}; - for (size_t i = 0; i < 3; ++i) - { - 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)) - { - adimKdiff0 = adimK1; - adimKis0 = adimK2; - mainDirection = direction1; - weakDirection = direction2; - } - else - { - adimKdiff0 = adimK2; - adimKis0 = adimK1; - mainDirection = direction2; - weakDirection = direction1; - } - primitive = findPrimitiveType(adimK1, adimK2, adimKdiff0, adimKis0); - } - // Populate nodes - nodes->k1[nodeId] = k1; - nodes->k2[nodeId] = k2; - nodes->adimK1[nodeId] = adimK1; - nodes->adimK2[nodeId] = adimK2; - for (size_t i = 0; i < 3; ++i) - { - nodes->weakDirections[3 * nodeId + i] = weakDirection[i]; - nodes->mainDirections[3 * nodeId + i] = mainDirection[i]; - } - nodes->primitives[nodeId] = primitive; -} - -PrimitiveType NodeCurvatureCalculator::findPrimitiveType(double k1, double k2, double kdiff0, double kis0) const -{ - if ((fabs(k1 - k2) < EPSILON_PRIMITIVE) && (fabs((k1 + k2) / 2) < EPSILON_PRIMITIVE)) - return PrimitiveType::Plane; - else if ((fabs(k1 - k2) < EPSILON_PRIMITIVE) && (fabs((k1 + k2) / 2) > EPSILON_PRIMITIVE)) - return PrimitiveType::Sphere; - else if ((fabs(kdiff0) > EPSILON_PRIMITIVE) && (fabs(kis0) < EPSILON_PRIMITIVE)) - return PrimitiveType::Cylinder; - else if ((fabs(k1) > EPSILON_PRIMITIVE) && (fabs(k2) > EPSILON_PRIMITIVE)) - return PrimitiveType::Torus; - else - return PrimitiveType::Unknown; -} - -std::vector NodeCurvatureCalculator::computeNormalCurvatureCoefficients( - const std::vector &discreteCurvatures, - const std::vector &tangents, - const std::array &normal, - const std::array &e1) const -{ - size_t nbNeighbors = discreteCurvatures.size(); - std::vector a(3 * nbNeighbors, 0.0); - for (size_t i = 0; i < nbNeighbors; ++i) - { - std::array tangent{tangents[3 * i], tangents[3 * i + 1], tangents[3 * i + 2]}; - double theta = MathOps::computeOrientedAngle(normal, tangent, e1); - double cos_theta = cos(theta); - double sin_theta = sin(theta); - a[i] = pow(cos_theta, 2); - a[nbNeighbors + i] = cos_theta * sin_theta; - a[2 * nbNeighbors + i] = pow(sin_theta, 2); - } - return MathOps::lstsq(a, discreteCurvatures); -} - -void NodeCurvatureCalculator::computeCellNormal( - const std::vector &nodeIds, - std::array &cellNormal) const -{ - std::vector point1; - std::vector point2; - std::vector point3; - mesh->getCoordinatesOfNode(nodeIds[0], point1); - mesh->getCoordinatesOfNode(nodeIds[1], point2); - mesh->getCoordinatesOfNode(nodeIds[2], point3); - std::array a; - a.fill(3); - std::array b; - b.fill(3); - for (int i = 0; i < 3; i++) - { - a[i] = point2[i] - point1[i]; - b[i] = point3[i] - point1[i]; - } - cellNormal[0] = a[1] * b[2] - a[2] * b[1]; - cellNormal[1] = a[2] * b[0] - a[0] * b[2]; - cellNormal[2] = a[0] * b[1] - a[1] * b[0]; -} - -double NodeCurvatureCalculator::computeAverageDistance(mcIdType nodeId, const std::vector &neighborIds) const -{ - double distance = 0.0; - std::array - nodeCoords = nodes->getCoordinates(nodeId); - for (size_t i = 0; i < neighborIds.size(); ++i) - { - std::array neighborCoords = nodes->getCoordinates(neighborIds[i]); - double distanceToNeighbor = 0.0; - for (size_t j = 0; j < 3; ++j) - distanceToNeighbor += pow(neighborCoords[j] - nodeCoords[j], 2); - distance += sqrt(distanceToNeighbor); - } - return distance / (double)neighborIds.size(); -} - -std::vector NodeCurvatureCalculator::computeDiscreteCurvatures( - mcIdType nodeId, - const std::vector &neighborIds) const -{ - std::vector discreteCurvatures(neighborIds.size(), 0.0); - for (size_t i = 0; i < neighborIds.size(); ++i) - discreteCurvatures[i] = computeDiscreteCurvature(nodeId, neighborIds[i]); - return discreteCurvatures; -} - -double NodeCurvatureCalculator::computeDiscreteCurvature( - mcIdType nodeId, - mcIdType neighborId) const -{ - double curvature = 0.0; - double n = 0.0; - for (size_t i = 0; i < 3; i++) - { - double ni = nodes->coords->getIJ(neighborId, i) - nodes->coords->getIJ(nodeId, i); - curvature += ni * (nodes->normals[3 * neighborId + i] - nodes->normals[3 * nodeId + i]); - n += ni * ni; - } - return curvature / n; -} - -std::vector NodeCurvatureCalculator::computeTangentDirections( - mcIdType nodeId, - const std::vector &neighborIds) const -{ - size_t nbNeighbors = neighborIds.size(); - std::vector tangent(3 * nbNeighbors, 0.0); - for (size_t i = 0; i < nbNeighbors; ++i) - { - mcIdType neighborId = neighborIds[i]; - double s = 0.0; - for (size_t j = 0; j < 3; j++) - { - tangent[3 * i + j] = nodes->coords->getIJ(neighborId, j) - nodes->coords->getIJ(nodeId, j); - s += tangent[3 * i + j] * nodes->normals[3 * nodeId + j]; - } - double n = 0.0; - for (size_t j = 0; j < 3; j++) - { - tangent[3 * i + j] -= s * nodes->normals[3 * nodeId + j]; - n += tangent[3 * i + j] * tangent[3 * i + j]; - } - for (size_t j = 0; j < 3; j++) - tangent[3 * i + j] /= sqrt(n); - } - return tangent; -} diff --git a/src/ShapeRecogn/NodeCurvatureCalculator.hxx b/src/ShapeRecogn/NodeCurvatureCalculator.hxx deleted file mode 100644 index ab15e2f87..000000000 --- a/src/ShapeRecogn/NodeCurvatureCalculator.hxx +++ /dev/null @@ -1,61 +0,0 @@ -// Copyright (C) 2007-2024 CEA, EDF -// -// This library is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com -// - -#ifndef __NODECURVATURECALCULATOR_HXX__ -#define __NODECURVATURECALCULATOR_HXX__ - -#include -#include -#include -#include "MEDCouplingUMesh.hxx" -#include "PrimitiveType.hxx" - -namespace MEDCoupling -{ - class Nodes; - - class NodeCurvatureCalculator - { - public: - NodeCurvatureCalculator(const MEDCouplingUMesh *); - - Nodes *compute(); - - private: - void computeNormals(); - void computeCurvatures(double tol = 0.000001); - void computeCurvatures(mcIdType nodeId, double tol); - PrimitiveType findPrimitiveType(double k1, double k2, double kdiff0, double kis0) const; - std::vector computeNormalCurvatureCoefficients( - const std::vector &discreteCurvatures, - const std::vector &tangents, - const std::array &normal, - const std::array &e1) const; - void computeCellNormal(const std::vector &nodeIds, std::array &cellNormal) const; - double computeAverageDistance(mcIdType nodeId, const std::vector &neighborIds) const; - std::vector computeDiscreteCurvatures(mcIdType nodeId, const std::vector &neighborIds) const; - double computeDiscreteCurvature(mcIdType nodeId, mcIdType neighborId) const; - std::vector computeTangentDirections(mcIdType nodeId, const std::vector &neighborIds) const; - - const MEDCouplingUMesh *mesh; - Nodes *nodes; - }; -}; - -#endif //__NODECURVATURECALCULATOR_HXX__ \ No newline at end of file diff --git a/src/ShapeRecogn/Nodes.hxx b/src/ShapeRecogn/Nodes.hxx index c22dc91d6..1d47d3060 100644 --- a/src/ShapeRecogn/Nodes.hxx +++ b/src/ShapeRecogn/Nodes.hxx @@ -29,7 +29,7 @@ namespace MEDCoupling class Nodes { public: - friend class NodeCurvatureCalculator; + friend class NodesBuilder; Nodes(const MEDCouplingUMesh *mesh, const DataArrayInt64 *neighbors, const DataArrayInt64 *neighborsIdx); diff --git a/src/ShapeRecogn/NodesBuilder.cxx b/src/ShapeRecogn/NodesBuilder.cxx new file mode 100644 index 000000000..1d261f27a --- /dev/null +++ b/src/ShapeRecogn/NodesBuilder.cxx @@ -0,0 +1,285 @@ + +#include "NodesBuilder.hxx" +#include "Nodes.hxx" +#include "MathOps.hxx" +#include "ShapeRecongConstants.hxx" + +using namespace MEDCoupling; + +NodesBuilder::NodesBuilder( + const MEDCouplingUMesh *mesh) : mesh(mesh) +{ +} + +Nodes *NodesBuilder::compute() +{ + DataArrayInt64 *neighbors; + DataArrayInt64 *neighborsIdx; + mesh->computeNeighborsOfNodes(neighbors, neighborsIdx); + nodes = new Nodes(mesh, neighbors, neighborsIdx); + computeNormals(); + computeCurvatures(); + return nodes; +} + +void NodesBuilder::computeNormals() +{ + mcIdType nbNodes = mesh->getNumberOfNodes(); + mcIdType nbCells = mesh->getNumberOfCells(); + std::array cellNormal; + cellNormal.fill(0.0); + // product of normal AreaByCell + std::vector prodNormalAreaByCell(3 * nbCells, 0.0); + for (int cellId = 0; cellId < nbCells; cellId++) + { + std::vector nodeIds; + mesh->getNodeIdsOfCell(cellId, nodeIds); + computeCellNormal(nodeIds, cellNormal); + prodNormalAreaByCell[3 * cellId + 0] = cellNormal[0] / 2.0; + prodNormalAreaByCell[3 * cellId + 1] = cellNormal[1] / 2.0; + prodNormalAreaByCell[3 * cellId + 2] = cellNormal[2] / 2.0; + } + // + nodes->normals.resize(3 * nbNodes, 1.0); + DataArrayInt64 *revNodal = DataArrayInt64::New(); + DataArrayInt64 *revNodalIdx = DataArrayInt64::New(); + mesh->getReverseNodalConnectivity(revNodal, revNodalIdx); + for (mcIdType nodeId = 0; nodeId < nbNodes; nodeId++) + { + int nbCells = revNodalIdx->getIJ(nodeId + 1, 0) - revNodalIdx->getIJ(nodeId, 0); + std::vector cellIds(nbCells, 0); + int start = revNodalIdx->getIJ(nodeId, 0); + for (size_t i = 0; i < cellIds.size(); ++i) + cellIds[i] = revNodal->getIJ(start + i, 0); + double normal = 0.0; + for (size_t i = 0; i < 3; i++) + { + nodes->normals[3 * nodeId + i] = 0.0; + for (mcIdType j = 0; j < nbCells; j++) + { + nodes->normals[3 * nodeId + i] += + prodNormalAreaByCell[3 * cellIds[j] + i]; + } + nodes->normals[3 * nodeId + i] /= (double)nbCells; + normal += pow(nodes->normals[3 * nodeId + i], 2); + } + for (size_t i = 0; i < 3; i++) + nodes->normals[3 * nodeId + i] /= sqrt(normal); + } + revNodal->decrRef(); + revNodalIdx->decrRef(); +} + +void NodesBuilder::computeCurvatures(double tol) +{ + mcIdType nbNodes = mesh->getNumberOfNodes(); + nodes->k1.resize(nbNodes); + nodes->k2.resize(nbNodes); + nodes->adimK1.resize(nbNodes); + nodes->adimK2.resize(nbNodes); + nodes->primitives.resize(nbNodes); + nodes->weakDirections.resize(3 * nbNodes); + nodes->mainDirections.resize(3 * nbNodes); + for (mcIdType nodeId = 0; nodeId < nbNodes; nodeId++) + computeCurvatures(nodeId, tol); +} + +void NodesBuilder::computeCurvatures(mcIdType nodeId, double tol) +{ + std::array normal = nodes->getNormal(nodeId); + std::vector neighborIds = nodes->getNeighbors(nodeId); + double theta0 = 0.0; + double k1 = 0.0; + double k2 = 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}; + if (neighborIds.size() > 0) + { + std::vector discreteCurvatures = computeDiscreteCurvatures(nodeId, neighborIds); + std::vector tangents = computeTangentDirections(nodeId, neighborIds); + mcIdType maxCurvatureId = std::distance( + discreteCurvatures.begin(), + std::max_element(discreteCurvatures.begin(), discreteCurvatures.end())); + std::array e1{ + tangents[3 * maxCurvatureId], + tangents[3 * maxCurvatureId + 1], + tangents[3 * maxCurvatureId + 2]}; + std::array e2 = MathOps::normalize(MathOps::cross(e1, normal)); + std::vector coeffs = computeNormalCurvatureCoefficients( + discreteCurvatures, tangents, normal, e1); + double a = coeffs[0], b = coeffs[1], c = coeffs[2]; + double h = (a + c) / 2.0; + double kg = a * c - pow(b, 2) / 4.0; + k1 = h + sqrt(fabs(pow(h, 2) - kg)); + k2 = h - sqrt(fabs(pow(h, 2) - kg)); + if (fabs(k1 - k2) >= tol) + theta0 = 0.5 * asin(b / (k1 - k2)); + std::array direction1{0.0, 0.0, 0.0}; + std::array direction2{0.0, 0.0, 0.0}; + for (size_t i = 0; i < 3; ++i) + { + 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)) + { + adimKdiff0 = adimK1; + adimKis0 = adimK2; + mainDirection = direction1; + weakDirection = direction2; + } + else + { + adimKdiff0 = adimK2; + adimKis0 = adimK1; + mainDirection = direction2; + weakDirection = direction1; + } + primitive = findPrimitiveType(adimK1, adimK2, adimKdiff0, adimKis0); + } + // Populate nodes + nodes->k1[nodeId] = k1; + nodes->k2[nodeId] = k2; + nodes->adimK1[nodeId] = adimK1; + nodes->adimK2[nodeId] = adimK2; + for (size_t i = 0; i < 3; ++i) + { + nodes->weakDirections[3 * nodeId + i] = weakDirection[i]; + nodes->mainDirections[3 * nodeId + i] = mainDirection[i]; + } + nodes->primitives[nodeId] = primitive; +} + +PrimitiveType NodesBuilder::findPrimitiveType(double k1, double k2, double kdiff0, double kis0) const +{ + if ((fabs(k1 - k2) < EPSILON_PRIMITIVE) && (fabs((k1 + k2) / 2) < EPSILON_PRIMITIVE)) + return PrimitiveType::Plane; + else if ((fabs(k1 - k2) < EPSILON_PRIMITIVE) && (fabs((k1 + k2) / 2) > EPSILON_PRIMITIVE)) + return PrimitiveType::Sphere; + else if ((fabs(kdiff0) > EPSILON_PRIMITIVE) && (fabs(kis0) < EPSILON_PRIMITIVE)) + return PrimitiveType::Cylinder; + else if ((fabs(k1) > EPSILON_PRIMITIVE) && (fabs(k2) > EPSILON_PRIMITIVE)) + return PrimitiveType::Torus; + else + return PrimitiveType::Unknown; +} + +std::vector NodesBuilder::computeNormalCurvatureCoefficients( + const std::vector &discreteCurvatures, + const std::vector &tangents, + const std::array &normal, + const std::array &e1) const +{ + size_t nbNeighbors = discreteCurvatures.size(); + std::vector a(3 * nbNeighbors, 0.0); + for (size_t i = 0; i < nbNeighbors; ++i) + { + std::array tangent{tangents[3 * i], tangents[3 * i + 1], tangents[3 * i + 2]}; + double theta = MathOps::computeOrientedAngle(normal, tangent, e1); + double cos_theta = cos(theta); + double sin_theta = sin(theta); + a[i] = pow(cos_theta, 2); + a[nbNeighbors + i] = cos_theta * sin_theta; + a[2 * nbNeighbors + i] = pow(sin_theta, 2); + } + return MathOps::lstsq(a, discreteCurvatures); +} + +void NodesBuilder::computeCellNormal( + const std::vector &nodeIds, + std::array &cellNormal) const +{ + std::vector point1; + std::vector point2; + std::vector point3; + mesh->getCoordinatesOfNode(nodeIds[0], point1); + mesh->getCoordinatesOfNode(nodeIds[1], point2); + mesh->getCoordinatesOfNode(nodeIds[2], point3); + std::array a; + a.fill(3); + std::array b; + b.fill(3); + for (int i = 0; i < 3; i++) + { + a[i] = point2[i] - point1[i]; + b[i] = point3[i] - point1[i]; + } + cellNormal[0] = a[1] * b[2] - a[2] * b[1]; + cellNormal[1] = a[2] * b[0] - a[0] * b[2]; + cellNormal[2] = a[0] * b[1] - a[1] * b[0]; +} + +double NodesBuilder::computeAverageDistance(mcIdType nodeId, const std::vector &neighborIds) const +{ + double distance = 0.0; + std::array + nodeCoords = nodes->getCoordinates(nodeId); + for (size_t i = 0; i < neighborIds.size(); ++i) + { + std::array neighborCoords = nodes->getCoordinates(neighborIds[i]); + double distanceToNeighbor = 0.0; + for (size_t j = 0; j < 3; ++j) + distanceToNeighbor += pow(neighborCoords[j] - nodeCoords[j], 2); + distance += sqrt(distanceToNeighbor); + } + return distance / (double)neighborIds.size(); +} + +std::vector NodesBuilder::computeDiscreteCurvatures( + mcIdType nodeId, + const std::vector &neighborIds) const +{ + std::vector discreteCurvatures(neighborIds.size(), 0.0); + for (size_t i = 0; i < neighborIds.size(); ++i) + discreteCurvatures[i] = computeDiscreteCurvature(nodeId, neighborIds[i]); + return discreteCurvatures; +} + +double NodesBuilder::computeDiscreteCurvature( + mcIdType nodeId, + mcIdType neighborId) const +{ + double curvature = 0.0; + double n = 0.0; + for (size_t i = 0; i < 3; i++) + { + double ni = nodes->coords->getIJ(neighborId, i) - nodes->coords->getIJ(nodeId, i); + curvature += ni * (nodes->normals[3 * neighborId + i] - nodes->normals[3 * nodeId + i]); + n += ni * ni; + } + return curvature / n; +} + +std::vector NodesBuilder::computeTangentDirections( + mcIdType nodeId, + const std::vector &neighborIds) const +{ + size_t nbNeighbors = neighborIds.size(); + std::vector tangent(3 * nbNeighbors, 0.0); + for (size_t i = 0; i < nbNeighbors; ++i) + { + mcIdType neighborId = neighborIds[i]; + double s = 0.0; + for (size_t j = 0; j < 3; j++) + { + tangent[3 * i + j] = nodes->coords->getIJ(neighborId, j) - nodes->coords->getIJ(nodeId, j); + s += tangent[3 * i + j] * nodes->normals[3 * nodeId + j]; + } + double n = 0.0; + for (size_t j = 0; j < 3; j++) + { + tangent[3 * i + j] -= s * nodes->normals[3 * nodeId + j]; + n += tangent[3 * i + j] * tangent[3 * i + j]; + } + for (size_t j = 0; j < 3; j++) + tangent[3 * i + j] /= sqrt(n); + } + return tangent; +} diff --git a/src/ShapeRecogn/NodesBuilder.hxx b/src/ShapeRecogn/NodesBuilder.hxx new file mode 100644 index 000000000..dde6312bd --- /dev/null +++ b/src/ShapeRecogn/NodesBuilder.hxx @@ -0,0 +1,61 @@ +// Copyright (C) 2007-2024 CEA, EDF +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#ifndef __NODECURVATURECALCULATOR_HXX__ +#define __NODECURVATURECALCULATOR_HXX__ + +#include +#include +#include +#include "MEDCouplingUMesh.hxx" +#include "PrimitiveType.hxx" + +namespace MEDCoupling +{ + class Nodes; + + class NodesBuilder + { + public: + NodesBuilder(const MEDCouplingUMesh *); + + Nodes *compute(); + + private: + void computeNormals(); + void computeCurvatures(double tol = 0.000001); + void computeCurvatures(mcIdType nodeId, double tol); + PrimitiveType findPrimitiveType(double k1, double k2, double kdiff0, double kis0) const; + std::vector computeNormalCurvatureCoefficients( + const std::vector &discreteCurvatures, + const std::vector &tangents, + const std::array &normal, + const std::array &e1) const; + void computeCellNormal(const std::vector &nodeIds, std::array &cellNormal) const; + double computeAverageDistance(mcIdType nodeId, const std::vector &neighborIds) const; + std::vector computeDiscreteCurvatures(mcIdType nodeId, const std::vector &neighborIds) const; + double computeDiscreteCurvature(mcIdType nodeId, mcIdType neighborId) const; + std::vector computeTangentDirections(mcIdType nodeId, const std::vector &neighborIds) const; + + const MEDCouplingUMesh *mesh; + Nodes *nodes; + }; +}; + +#endif //__NODECURVATURECALCULATOR_HXX__ \ No newline at end of file diff --git a/src/ShapeRecogn/ShapeRecognMesh.cxx b/src/ShapeRecogn/ShapeRecognMesh.cxx index d5910e949..b70ee9539 100644 --- a/src/ShapeRecogn/ShapeRecognMesh.cxx +++ b/src/ShapeRecogn/ShapeRecognMesh.cxx @@ -1,6 +1,6 @@ #include "ShapeRecognMesh.hxx" -#include "NodeCurvatureCalculator.hxx" +#include "NodesBuilder.hxx" #include "AreasBuilder.hxx" #include "MEDLoader.hxx" #include "MEDCouplingFieldDouble.hxx" @@ -25,8 +25,8 @@ ShapeRecognMesh::~ShapeRecognMesh() void ShapeRecognMesh::recognize() { mesh->incrRef(); - NodeCurvatureCalculator nodeCurvatureCalc(mesh); - nodes = nodeCurvatureCalc.compute(); + NodesBuilder nodesBuilder(mesh); + nodes = nodesBuilder.compute(); AreasBuilder areasBuilder(nodes); areasBuilder.explore(); areasBuilder.expand();