--- /dev/null
+#include "Areas.hxx"
+#include "MathOps.hxx"
+
+#include <random>
+
+using namespace MEDCoupling;
+
+Areas::Areas(const Nodes *nodes) : nodes(nodes), areaIdByNodes(nodes->getNbNodes(), -1)
+{
+}
+
+mcIdType Areas::addArea(PrimitiveType primitive)
+{
+ Area area;
+ area.primitive = primitive;
+ areas.push_back(area);
+ return areas.size() - 1;
+}
+
+void Areas::cleanArea(mcIdType areaId)
+{
+ cleanArea(areaId, -1);
+}
+
+void Areas::cancelArea(mcIdType areaId, PrimitiveType primitive)
+{
+ cleanArea(areaId, -2);
+ Area &area = areas[areaId];
+ area.primitive = primitive;
+}
+
+void Areas::removeArea(mcIdType areaId)
+{
+ for (size_t nodeId = 0; nodeId < areaIdByNodes.size(); ++nodeId)
+ {
+ if (areaIdByNodes[nodeId] == areaId)
+ areaIdByNodes[nodeId] = -1;
+ else if (areaIdByNodes[nodeId] > areaId)
+ areaIdByNodes[nodeId] -= 1;
+ }
+ areas.erase(areas.begin() + areaId);
+}
+
+void Areas::addNode(mcIdType areaId, mcIdType nodeId)
+{
+ removeNode(nodeId);
+ areaIdByNodes[nodeId] = areaId;
+ Area &area = areas[areaId];
+ area.nodeIds.push_back(nodeId);
+ size_t nbNodes = area.nodeIds.size();
+ area.k1 = ((nbNodes - 1) * area.k1 + nodes->getK1(nodeId)) / nbNodes;
+ area.k2 = ((nbNodes - 1) * area.k2 + nodes->getK2(nodeId)) / nbNodes;
+ area.kdiff0 = ((nbNodes - 1) * area.kdiff0 + nodes->getKdiff0(nodeId)) / nbNodes;
+}
+
+void Areas::cleanInvalidNodeAreas()
+{
+ for (mcIdType nodeId = 0; nodeId < (mcIdType)areaIdByNodes.size(); ++nodeId)
+ {
+ if (areaIdByNodes[nodeId] <= -2)
+ areaIdByNodes[nodeId] = -1;
+ }
+}
+
+mcIdType Areas::getAreaId(mcIdType nodeId) const
+{
+ return areaIdByNodes[nodeId];
+}
+
+bool Areas::isEmpty(mcIdType areaId) const
+{
+ return areas[areaId].nodeIds.empty();
+}
+
+size_t Areas::getNumberOfAreas() const
+{
+ return areas.size();
+}
+
+size_t Areas::getNumberOfNodes(mcIdType areaId) const
+{
+ return areas[areaId].nodeIds.size();
+}
+
+bool Areas::isNodeCompatible(mcIdType areaId, mcIdType nodeId) const
+{
+ PrimitiveType areaType = getPrimitiveType(areaId);
+ PrimitiveType nodeType = nodes->getPrimitiveType(nodeId);
+ return (
+ (areaType != PrimitiveType::Cylinder ||
+ nodeType != PrimitiveType::Plane) &&
+ (areaType != PrimitiveType::Sphere ||
+ nodeType != PrimitiveType::Plane));
+}
+
+PrimitiveType Areas::getPrimitiveType(mcIdType areaId) const
+{
+ const Area &area = areas[areaId];
+ if (area.primitive != PrimitiveType::Unknown)
+ return area.primitive;
+ else if (area.nodeIds.empty())
+ return PrimitiveType::Unknown;
+ else
+ return nodes->getPrimitiveType(area.nodeIds[0]);
+}
+
+std::string Areas::getPrimitiveTypeName(mcIdType areaId) const
+{
+ return convertPrimitiveToString(getPrimitiveType(areaId));
+}
+
+int Areas::getPrimitiveTypeInt(mcIdType areaId) const
+{
+ return convertPrimitiveToInt(getPrimitiveType(areaId));
+}
+
+const std::vector<mcIdType> &Areas::getNodeIds(mcIdType areaId) const
+{
+ return areas[areaId].nodeIds;
+}
+
+double Areas::getK1(mcIdType areaId) const
+{
+
+ return areas[areaId].k1;
+}
+
+double Areas::getK2(mcIdType areaId) const
+{
+ return areas[areaId].k2;
+}
+
+double Areas::getKdiff0(mcIdType areaId) const
+{
+ return areas[areaId].kdiff0;
+}
+
+void Areas::computeProperties(mcIdType areaId)
+{
+ switch (getPrimitiveType(areaId))
+ {
+ case PrimitiveType::Plane:
+ computePlaneProperties(areaId);
+ break;
+ case PrimitiveType::Sphere:
+ computeSphereProperties(areaId);
+ break;
+ case PrimitiveType::Cylinder:
+ computeCylinderProperties(areaId);
+ break;
+ case PrimitiveType::Cone:
+ computeConeProperties(areaId);
+ break;
+ case PrimitiveType::Torus:
+ computeTorusProperties(areaId);
+ break;
+ case PrimitiveType::Unknown:
+ default:
+ break;
+ }
+}
+
+double Areas::getRadius(mcIdType areaId) const
+{
+ const Area &area = areas[areaId];
+ return area.radius;
+}
+
+double Areas::getAngle(mcIdType areaId) const
+{
+ return areas[areaId].angle;
+}
+
+const std::array<double, 3> &Areas::getNormal(mcIdType areaId) const
+{
+ return areas[areaId].normal;
+}
+
+const std::array<double, 3> &Areas::getCenter(mcIdType areaId) const
+{
+ const Area &area = areas[areaId];
+ return area.center;
+}
+
+const std::array<double, 3> &Areas::getAxis(mcIdType areaId) const
+{
+ return areas[areaId].axis;
+}
+
+const std::array<double, 3> &Areas::getAxisPoint(mcIdType areaId) const
+{
+ return areas[areaId].axisPoint;
+}
+
+const std::array<double, 3> &Areas::getApex(mcIdType areaId) const
+{
+ return areas[areaId].apex;
+}
+
+std::array<double, 3> Areas::getAffinePoint(mcIdType areaId) const
+{
+ const Area &area = areas[areaId];
+ if (area.nodeIds.empty())
+ return {0.0, 0.0, 0.0};
+ else
+ return nodes->getCoordinates(area.nodeIds[0]);
+}
+
+void Areas::cleanArea(mcIdType areaId, mcIdType newAreaId = -1)
+{
+ Area &area = areas[areaId];
+ for (mcIdType nodeId : area.nodeIds)
+ areaIdByNodes[nodeId] = newAreaId;
+ area.primitive = PrimitiveType::Unknown;
+ area.k1 = 0.0;
+ area.k2 = 0.0;
+ area.kdiff0 = 0.0;
+ area.radius = 0.0;
+ area.angle = 0.0;
+ area.normal = {0.0, 0.0, 0.0};
+ area.center = {0.0, 0.0, 0.0};
+ area.axis = {0.0, 0.0, 0.0};
+ area.axisPoint = {0.0, 0.0, 0.0};
+ area.apex = {0.0, 0.0, 0.0};
+ area.nodeIds.clear();
+}
+
+void Areas::removeNode(mcIdType nodeId)
+{
+ mcIdType areaId = areaIdByNodes[nodeId];
+ if (areaId > -1)
+ {
+ Area &area = areas[areaId];
+ area.nodeIds.erase(
+ std::remove(
+ area.nodeIds.begin(),
+ area.nodeIds.end(),
+ nodeId),
+ area.nodeIds.end());
+ areaIdByNodes[nodeId] = -1;
+ // TODO: Update the parameters of the area ?
+ }
+}
+
+void Areas::computePlaneProperties(mcIdType areaId)
+{
+ Area &area = areas[areaId];
+ const std::vector<double> &normals = nodes->getNormals();
+ mcIdType nbNodes = area.nodeIds.size();
+ area.normal = {0.0, 0.0, 0.0};
+ for (mcIdType nodeId : area.nodeIds)
+ {
+ for (size_t i = 0; i < 3; ++i)
+ area.normal[i] += normals[3 * nodeId + i];
+ }
+ for (size_t i = 0; i < 3; ++i)
+ area.normal[i] /= nbNodes;
+}
+
+void Areas::computeSphereProperties(mcIdType areaId)
+{
+ Area &area = areas[areaId];
+ const std::vector<double> &normals = nodes->getNormals();
+ area.radius = (2 / (area.k1 + area.k2));
+ std::array<double, 3> center{0.0, 0.0, 0.0};
+ if (!area.nodeIds.empty())
+ {
+ size_t nbNodes = area.nodeIds.size();
+ for (mcIdType nodeId : area.nodeIds)
+ {
+ std::array<double, 3> nodeCoords = nodes->getCoordinates(nodeId);
+ for (size_t i = 0; i < 3; ++i)
+ center[i] += nodeCoords[i] - area.radius * normals[3 * nodeId + i];
+ }
+ for (size_t i = 0; i < 3; ++i)
+ center[i] /= nbNodes;
+ }
+ area.center = center;
+ area.radius = fabs(area.radius);
+}
+
+void Areas::computeCylinderProperties(mcIdType areaId)
+{
+ Area &area = areas[areaId];
+ size_t nbNodes = area.nodeIds.size();
+ // Project the nodes to the central axis of the cylinder
+ std::vector<double> projectedNodes(3 * nbNodes, 0.0);
+ area.radius = 0;
+ area.axisPoint.fill(0.0);
+ for (size_t i = 0; i < nbNodes; ++i)
+ {
+ mcIdType nodeId = area.nodeIds[i];
+ std::array<double, 3> nodeCoords = nodes->getCoordinates(nodeId);
+ double approxRadius = 1.0 / nodes->getKdiff0(nodeId);
+ for (size_t j = 0; j < 3; ++j)
+ {
+ projectedNodes[3 * i + j] =
+ nodeCoords[j] - approxRadius * nodes->getNormals()[3 * nodeId + j];
+ area.axisPoint[j] += projectedNodes[3 * i + j];
+ }
+ area.radius += approxRadius;
+ }
+ // Axis point is the mean of the projected nodes
+ for (size_t i = 0; i < 3; ++i)
+ area.axisPoint[i] /= nbNodes;
+ // Radius of the cylinder is the mean of the approximate radius of each node
+ area.radius = fabs(area.radius / nbNodes);
+ // Compute the axis of the cylinder
+ area.axis = MathOps::computePCAFirstAxis(projectedNodes);
+}
+
+void Areas::computeConeProperties(mcIdType areaId)
+{
+ Area &area = areas[areaId];
+ size_t nbNodes = area.nodeIds.size();
+ // Project the nodes to the central axis of the cone
+ std::vector<double> projectedNodes(3 * nbNodes, 0.0);
+ std::vector<double> radiusNodes(nbNodes, 0.0);
+ area.axisPoint.fill(0.0);
+ for (size_t i = 0; i < nbNodes; ++i)
+ {
+ mcIdType nodeId = area.nodeIds[i];
+ std::array<double, 3> nodeCoords = nodes->getCoordinates(nodeId);
+ radiusNodes[i] = 1.0 / nodes->getKdiff0(nodeId);
+ for (size_t j = 0; j < 3; ++j)
+ {
+ projectedNodes[3 * i + j] =
+ nodeCoords[j] - radiusNodes[i] * nodes->getNormals()[3 * nodeId + j];
+ area.axisPoint[j] += projectedNodes[3 * i + j];
+ }
+ }
+ // Axis point is the mean of the projected nodes
+ for (size_t i = 0; i < 3; ++i)
+ area.axisPoint[i] /= nbNodes;
+ // Compute the axis of the cone
+ area.axis = MathOps::computePCAFirstAxis(projectedNodes);
+ double normAxis = MathOps::computeNorm(area.axis);
+ for (size_t i = 0; i < 3; ++i)
+ area.axis[i] /= normAxis;
+ // Compute the angle of the cone
+ const std::vector<double> &weakDirections = nodes->getWeakDirections();
+ std::vector<double> weakDirectionNodes(3 * nbNodes, 0.0);
+ for (size_t i = 0; i < nbNodes; ++i)
+ {
+ mcIdType nodeId = area.nodeIds[i];
+ for (size_t j = 0; j < 3; ++j)
+ weakDirectionNodes[3 * i + j] = weakDirections[3 * nodeId + j];
+ }
+ std::vector<double> angles = MathOps::computeAngles(
+ weakDirectionNodes, area.axis);
+ // Correct the angles > pi/2 with the supplementary angle
+ for (size_t i = 0; i < angles.size(); ++i)
+ {
+ if (angles[i] > M_PI_2)
+ angles[i] = M_PI - angles[i];
+ }
+ area.angle = MathOps::mean(angles);
+ // Compute the radius of the cone which is the mean of the approximate radius of each node
+ area.radius = MathOps::mean(radiusNodes);
+ // Select extrem nodes
+ double q1 = MathOps::computeQuantile(radiusNodes, 0.1);
+ double q2 = MathOps::computeQuantile(radiusNodes, 0.9);
+ std::vector<mcIdType> q1_indices;
+ for (mcIdType idx = 0; idx < (mcIdType)radiusNodes.size(); ++idx)
+ {
+ if (radiusNodes[idx] < q1)
+ q1_indices.push_back(idx);
+ }
+ std::vector<mcIdType> q2_indices;
+ for (mcIdType idx = 0; idx < (mcIdType)radiusNodes.size(); ++idx)
+ {
+ if (radiusNodes[idx] > q2)
+ q2_indices.push_back(idx);
+ }
+ std::random_device rd;
+ std::mt19937 g(rd());
+ std::shuffle(q2_indices.begin(), q2_indices.end(), g);
+ // Compute the height of the cone
+ // std::vector<double> heights(q1_indices.size(), 0.0);
+ // std::vector<double> distancesToApex(q1_indices.size(), 0.0);
+ std::array<double, 3> p{0.0, 0.0, 0.0};
+ for (size_t i = 0; i < q1_indices.size(); ++i)
+ {
+ for (size_t j = 0; j < 3; ++j)
+ p[j] = (projectedNodes[3 * q1_indices[i] + j] - projectedNodes[3 * q2_indices[i] + j]);
+ double height = MathOps::computeNorm(p);
+ double distanceToApex = height * radiusNodes[q1_indices[i]] / (radiusNodes[q2_indices[i]] - radiusNodes[q1_indices[i]]);
+ double orientationAxis = MathOps::dot(p, area.axis);
+ for (size_t j = 0; j < 3; ++j)
+ {
+ area.apex[j] += projectedNodes[3 * q1_indices[i] + j];
+ if (orientationAxis >= 0)
+ area.apex[j] += distanceToApex * area.axis[j];
+ else
+ area.apex[j] -= distanceToApex * area.axis[j];
+ }
+ }
+ for (size_t j = 0; j < 3; ++j)
+ area.apex[j] /= q1_indices.size();
+}
+
+void Areas::computeTorusProperties(mcIdType areaId)
+{
+ Area &area = areas[areaId];
+ size_t n = area.nodeIds.size();
+ std::vector<double> areaNodesK1(n, 0.0);
+ std::vector<double> areaNodesK2(n, 0.0);
+ for (size_t i = 0; i < n; ++i)
+ {
+ 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 minorCurvature;
+ if (var1 > var2)
+ minorCurvature = MathOps::mean(areaNodesK2);
+ else
+ minorCurvature = MathOps::mean(areaNodesK1);
+ double 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;
+ }
+ std::array<double, 3> meanMajorRadiusNodes = MathOps::meanCoordinates(majorRadiusNodes);
+ for (size_t i = 0; i < n; ++i)
+ {
+ for (size_t j = 0; j < 3; ++j)
+ 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,
+ }
+ )
+ */
+}
+
+const std::vector<mcIdType> &Areas::getAreaIdByNodes() const
+{
+ return areaIdByNodes;
+}
--- /dev/null
+// 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 __AREAS_HXX__
+#define __AREAS_HXX__
+
+#include "PrimitiveType.hxx"
+#include "Nodes.hxx"
+
+#include <array>
+
+namespace MEDCoupling
+{
+
+ struct Area
+ {
+ PrimitiveType primitive = PrimitiveType::Unknown;
+ double k1 = 0.0;
+ double k2 = 0.0;
+ double kdiff0 = 0.0;
+ double radius = 0.0;
+ double angle = 0.0;
+ std::array<double, 3> normal{0.0, 0.0, 0.0};
+ std::array<double, 3> center{0.0, 0.0, 0.0};
+ std::array<double, 3> axis{0.0, 0.0, 0.0};
+ std::array<double, 3> axisPoint{0.0, 0.0, 0.0};
+ std::array<double, 3> apex{0.0, 0.0, 0.0};
+ std::vector<mcIdType> nodeIds;
+ };
+
+ class Areas
+ {
+ public:
+ Areas(const Nodes *nodes);
+
+ mcIdType addArea(PrimitiveType primitive = PrimitiveType::Unknown);
+ void cleanArea(mcIdType areaId);
+ void cancelArea(mcIdType areaId, PrimitiveType primitive = PrimitiveType::Unknown);
+ void removeArea(mcIdType areaId);
+ void addNode(mcIdType areaId, mcIdType nodeId);
+ void cleanInvalidNodeAreas();
+
+ mcIdType getAreaId(mcIdType nodeId) const;
+ const std::vector<mcIdType> &getAreaIdByNodes() const;
+
+ bool isEmpty(mcIdType areaId) const;
+ size_t getNumberOfAreas() const;
+ size_t getNumberOfNodes(mcIdType areaId) const;
+
+ bool isNodeCompatible(mcIdType areaId, mcIdType nodeId) const;
+
+ PrimitiveType getPrimitiveType(mcIdType areaId) const;
+ std::string getPrimitiveTypeName(mcIdType areaId) const;
+ int getPrimitiveTypeInt(mcIdType areaId) const;
+
+ const std::vector<mcIdType> &getNodeIds(mcIdType areaId) const;
+
+ double getK1(mcIdType areaId) const;
+ double getK2(mcIdType areaId) const;
+ double getKdiff0(mcIdType areaId) const;
+
+ void computeProperties(mcIdType areaId);
+ double getRadius(mcIdType areaId) const;
+ double getAngle(mcIdType areaId) const;
+ const std::array<double, 3> &getNormal(mcIdType areaId) const;
+ const std::array<double, 3> &getCenter(mcIdType areaId) const;
+ const std::array<double, 3> &getAxis(mcIdType areaId) const;
+ const std::array<double, 3> &getAxisPoint(mcIdType areaId) const;
+ const std::array<double, 3> &getApex(mcIdType areaId) const;
+ std::array<double, 3> getAffinePoint(mcIdType areaId) const;
+
+ private:
+ void cleanArea(mcIdType areaId, mcIdType newAreaId);
+ void removeNode(mcIdType nodeId);
+
+ void computePlaneProperties(mcIdType areaId);
+ void computeSphereProperties(mcIdType areaId);
+ void computeCylinderProperties(mcIdType areaId);
+ void computeConeProperties(mcIdType areaId);
+ void computeTorusProperties(mcIdType areaId);
+
+ std::vector<Area> areas;
+ const Nodes *nodes;
+ std::vector<mcIdType> areaIdByNodes;
+ };
+};
+
+#endif // __AREAS_HXX__
--- /dev/null
+#include "AreasBuilder.hxx"
+#include "MathOps.hxx"
+#include "ShapeRecongConstants.hxx"
+
+#include <unordered_set>
+
+using namespace MEDCoupling;
+
+AreasBuilder::AreasBuilder(const Nodes *nodes) : nodes(nodes), areas(new Areas(nodes))
+{
+ size_t nbNodes = nodes->getNbNodes();
+ threshold = std::max<size_t>(THRESHOLD_MIN_NB_NODES, nbNodes / THRESHOLD_MAX_NB_AREAS);
+}
+
+void AreasBuilder::explore()
+{
+ exploreAreas();
+ // std::cout << "exploreAreas ended\n";
+ areas->cleanInvalidNodeAreas();
+ // std::cout << "cleanInvalidNodeAreas ended\n";
+ filterHighPass();
+ // std::cout << "explore ended\n";
+ // std::cout << "-------\n";
+ // for (mcIdType areaId = 0; areaId < areas->getNumberOfAreas(); ++areaId)
+ // {
+ // std::cout << "area " << areaId << ": "
+ // << areas->getNumberOfNodes(areaId) << ", "
+ // << areas->getPrimitiveTypeName(areaId)
+ // << "\n";
+ // }
+}
+
+void AreasBuilder::expand()
+{
+ expandAreas();
+ // std::cout << "expandAreas ended\n";
+ filterHighPass();
+ // std::cout << "expand ended\n";
+ // std::cout << "-------\n";
+ // for (mcIdType areaId = 0; areaId < areas->getNumberOfAreas(); ++areaId)
+ // {
+ // std::cout << "area " << areaId << ": "
+ // << areas->getNumberOfNodes(areaId) << ", "
+ // << areas->getPrimitiveTypeName(areaId)
+ // << "\n";
+ // }
+}
+
+void AreasBuilder::rebuild()
+{
+ rebuildInvalidAreas();
+ // std::cout << "rebuildInvalidAreas ended\n";
+ filterHighPass();
+ // std::cout << "filterHighPass ended\n";
+ expandAreasByType(PrimitiveType::Cone);
+ // std::cout << "expandAreas Cone ended\n";
+ filterHighPass();
+ // std::cout << "rebuild ended\n";
+}
+
+Areas *AreasBuilder::getAreas() const
+{
+ return areas;
+}
+
+void AreasBuilder::exploreAreas()
+{
+ // int nbNodesExplored = 0;
+ std::vector<bool> exploredNodeIds(nodes->getNbNodes(), false);
+ std::unordered_set<mcIdType> nodesToExplore;
+ // Reserve a set with the size of nodes to avoid reallocation for each insert/erase
+ // TODO: Improve the size ? Nb of Nodes is too much ?
+ nodesToExplore.reserve(nodes->getNbNodes());
+ mcIdType areaId = -1;
+ for (mcIdType nodeId = 0; nodeId < nodes->getNbNodes(); ++nodeId)
+ {
+ if (!exploredNodeIds[nodeId] &&
+ nodes->getPrimitiveType(nodeId) != PrimitiveType::Unknown)
+ {
+ exploredNodeIds[nodeId] = true;
+ if (areaId != -1 && areas->getNumberOfNodes(areaId) < threshold)
+ areas->cancelArea(areaId, nodes->getPrimitiveType(nodeId));
+ else
+ areaId = areas->addArea(nodes->getPrimitiveType(nodeId));
+ areas->addNode(areaId, nodeId);
+ // if (nbNodesExplored % 1 == 0)
+ // {
+ // std::cout << "nodeId: " << nodeId << " " << convertPrimitiveToString(nodes->getPrimitiveType(nodeId)) << "\n";
+ // std::cout << "nbNodesExplored (current area: "
+ // << areaId
+ // << ", "
+ // << areas->getPrimitiveTypeName(areaId)
+ // << "): "
+ // << nbNodesExplored
+ // << "\n";
+ // }
+ const std::vector<mcIdType> neighbors = nodes->getNeighbors(nodeId);
+ for (mcIdType neighborId : neighbors)
+ {
+ if (nodes->getPrimitiveType(neighborId) == areas->getPrimitiveType(areaId) &&
+ areas->getAreaId(neighborId) <= -1)
+ nodesToExplore.insert(neighborId);
+ }
+ // Explore all the neighbors matching the area
+ while (!nodesToExplore.empty())
+ {
+ mcIdType neighborId = *nodesToExplore.begin();
+ nodesToExplore.erase(neighborId);
+ if (doesItMatch(areaId, neighborId))
+ {
+ // if (!exploredNodeIds[neighborId])
+ // nbNodesExplored += 1;
+ // if (areas->getNumberOfNodes(areaId) % 500 == 0)
+ // {
+ // std::cout << "nbNodesExplored (current area: "
+ // << areaId
+ // << ", " << areas->getPrimitiveTypeName(areaId)
+ // << ", " << areas->getNumberOfNodes(areaId)
+ // << "): "
+ // << nbNodesExplored
+ // << "\n";
+ // std::cout << "neighborId: " << neighborId << " " << convertPrimitiveToString(nodes->getPrimitiveType(neighborId)) << "\n";
+ // }
+ exploredNodeIds[neighborId] = true;
+ areas->addNode(areaId, neighborId);
+ const std::vector<mcIdType> neighborsOfNeighbor = nodes->getNeighbors(neighborId);
+ for (mcIdType neighborIdOfNeighbor : neighborsOfNeighbor)
+ {
+ if (!exploredNodeIds[neighborIdOfNeighbor] &&
+ areas->getAreaId(neighborIdOfNeighbor) <= -1 &&
+ // Already in doesItMatch but avoid useless insertion
+ nodes->getPrimitiveType(neighborIdOfNeighbor) == areas->getPrimitiveType(areaId))
+ nodesToExplore.insert(neighborIdOfNeighbor);
+ }
+ }
+ }
+ }
+ // if (!exploredNodeIds[nodeId])
+ // nbNodesExplored += 1;
+ exploredNodeIds[nodeId] = true;
+ }
+}
+
+void AreasBuilder::expandAreas()
+{
+ // Expand by topological order
+ expandAreasByType(PrimitiveType::Plane);
+ expandAreasByType(PrimitiveType::Sphere);
+ expandAreasByType(PrimitiveType::Cylinder);
+ expandAreasByType(PrimitiveType::Cone);
+ expandAreasByType(PrimitiveType::Torus);
+}
+
+void AreasBuilder::expandAreasByType(PrimitiveType primitive)
+{
+ std::unordered_set<mcIdType> nodesToExplore;
+ // Reserve a set with the size of nodes to avoid reallocation for each insert/erase
+ // TODO: Improve the size ? Nb of Nodes is too much ?
+ nodesToExplore.reserve(nodes->getNbNodes());
+ for (mcIdType areaId = 0; areaId < (mcIdType)areas->getNumberOfAreas(); ++areaId)
+ {
+ if (areas->getPrimitiveType(areaId) == primitive)
+ {
+ std::vector<bool> exploredNodeIds(nodes->getNbNodes(), false);
+ areas->computeProperties(areaId);
+ const std::vector<mcIdType> &nodeIds = areas->getNodeIds(areaId);
+ for (mcIdType nodeId : nodeIds)
+ {
+ exploredNodeIds[nodeId] = true;
+ nodesToExplore.insert(nodeId);
+ }
+ while (!nodesToExplore.empty())
+ {
+ mcIdType nodeId = *nodesToExplore.begin();
+ nodesToExplore.erase(nodeId);
+ if (doesItBelong(areaId, nodeId))
+ {
+ // TODO: Is the properties need to be updated after adding a node ?
+ // It gives bad results for the cone and the cylinder
+ // mcIdType oldAreaId = areas->getAreaId(nodeId);
+ areas->addNode(areaId, nodeId);
+ // areas->computeProperties(areaId);
+ // areas->computeProperties(oldAreaId);
+ const std::vector<mcIdType> neighborIds = nodes->getNeighbors(nodeId);
+ for (mcIdType neighborId : neighborIds)
+ {
+ if (!exploredNodeIds[neighborId])
+ nodesToExplore.insert(neighborId);
+ }
+ }
+ exploredNodeIds[nodeId] = true;
+ }
+ }
+ }
+}
+
+void AreasBuilder::rebuildInvalidAreas()
+{
+ std::vector<mcIdType> exploredNodeIds(nodes->getNbNodes(), false);
+ std::vector<bool> isIinvalidNodes(nodes->getNbNodes(), false);
+ std::unordered_set<mcIdType> nodesToExplore;
+ // Reserve a set with the size of nodes to avoid reallocation for each insert/erase
+ // TODO: Improve the size ? Nb of Nodes is too much ?
+ nodesToExplore.reserve(nodes->getNbNodes());
+ for (mcIdType nodeId = 0; nodeId < nodes->getNbNodes(); ++nodeId)
+ isIinvalidNodes[nodeId] = isInvalidCylinderNode(nodeId);
+ for (mcIdType nodeId = 0; nodeId < nodes->getNbNodes(); ++nodeId)
+ {
+ if (isIinvalidNodes[nodeId] && !exploredNodeIds[nodeId])
+ {
+ mcIdType areaId = areas->addArea(PrimitiveType::Cone);
+ areas->addNode(areaId, nodeId);
+ const std::vector<mcIdType> neighbors = nodes->getNeighbors(nodeId);
+ for (mcIdType neighborId : neighbors)
+ nodesToExplore.insert(neighborId);
+ while (!nodesToExplore.empty())
+ {
+ mcIdType neighborId = *nodesToExplore.begin();
+ nodesToExplore.erase(neighborId);
+ if (
+ (areas->getAreaId(neighborId) == -1 ||
+ isIinvalidNodes[neighborId]))
+ {
+ exploredNodeIds[neighborId] = true;
+ areas->addNode(areaId, neighborId);
+ const std::vector<mcIdType> neighborsOfNeighbor = nodes->getNeighbors(neighborId);
+ for (mcIdType neighborIdOfNeighbor : neighborsOfNeighbor)
+ {
+ if (!exploredNodeIds[neighborIdOfNeighbor])
+ nodesToExplore.insert(neighborIdOfNeighbor);
+ }
+ }
+ }
+ }
+ }
+}
+
+void AreasBuilder::filterHighPass()
+{
+ mcIdType nbAreas = areas->getNumberOfAreas();
+ for (mcIdType areaId = (nbAreas - 1); areaId >= 0; --areaId)
+ {
+ if (areas->getNumberOfNodes(areaId) < threshold)
+ areas->removeArea(areaId);
+ }
+}
+
+bool AreasBuilder::doesItMatch(mcIdType areaId, mcIdType nodeId) const
+{
+ PrimitiveType areaType = areas->getPrimitiveType(areaId);
+ PrimitiveType neighborType = nodes->getPrimitiveType(nodeId);
+ bool isMatching = false;
+ if (areaType == neighborType)
+ {
+ switch (areaType)
+ {
+ case PrimitiveType::Plane:
+ case PrimitiveType::Torus:
+ isMatching = true;
+ break;
+ case PrimitiveType::Sphere:
+ {
+ double kmoy = fabs(
+ (areas->getK1(areaId) + areas->getK2(areaId)) / 2.0);
+ double nodeKmoy = fabs(
+ (nodes->getK1(nodeId) + nodes->getK2(nodeId)) / 2.0);
+ // TODO: Replace by relative tolerance
+ isMatching = (kmoy + TOL_MATCH_SPHERE > nodeKmoy) &&
+ (nodeKmoy > kmoy - TOL_MATCH_SPHERE);
+ }
+ break;
+ case PrimitiveType::Cylinder:
+ isMatching =
+ fabs((areas->getKdiff0(areaId) - nodes->getKdiff0(nodeId)) /
+ areas->getKdiff0(areaId)) < EPSILON_PROP;
+ break;
+ case PrimitiveType::Cone:
+ case PrimitiveType::Unknown:
+ default:
+ break;
+ }
+ }
+ return isMatching;
+}
+
+bool AreasBuilder::doesItBelong(mcIdType areaId, mcIdType nodeId) const
+{
+ bool isClose = false;
+ if (areas->isNodeCompatible(areaId, nodeId))
+ {
+ PrimitiveType areaType = areas->getPrimitiveType(areaId);
+ switch (areaType)
+ {
+ case PrimitiveType::Plane:
+ {
+ isClose = distanceToPlane(
+ nodes->getCoordinates(nodeId),
+ areas->getAffinePoint(areaId),
+ areas->getNormal(areaId)) < DELTA_PLANE;
+ }
+ break;
+ case PrimitiveType::Sphere:
+ {
+ double distanceToCenter = distanceToSphere(
+ nodes->getCoordinates(nodeId),
+ areas->getCenter(areaId));
+ isClose = fabs(distanceToCenter - areas->getRadius(areaId)) < DELTA_SPHERE;
+ }
+ break;
+ case PrimitiveType::Cylinder:
+ {
+ double distance = distanceToCylinder(
+ nodes->getCoordinates(nodeId),
+ areas->getAxis(areaId),
+ areas->getAxisPoint(areaId));
+ double radius = areas->getRadius(areaId);
+ isClose = fabs((distance - radius) / radius) < DELTA_CYLINDER;
+ }
+ break;
+ case PrimitiveType::Cone:
+ {
+ isClose = distanceToCone(
+ nodes->getCoordinates(nodeId),
+ areas->getAxis(areaId),
+ areas->getApex(areaId),
+ areas->getAngle(areaId),
+ areas->getRadius(areaId)) < DELTA_CONE;
+ }
+ break;
+ case PrimitiveType::Torus:
+ case PrimitiveType::Unknown:
+ default:
+ break;
+ }
+ }
+ return isClose;
+}
+
+bool AreasBuilder::isInvalidCylinderNode(mcIdType nodeId) const
+{
+ mcIdType areaId = areas->getAreaId(nodeId);
+ if (areaId != -1 &&
+ nodes->getPrimitiveType(nodeId) == PrimitiveType::Cylinder &&
+ areas->getPrimitiveType(areaId) == PrimitiveType::Cylinder)
+ {
+ double angle = MathOps::computeAngle(
+ nodes->getWeakDirection(nodeId),
+ areas->getAxis(areaId));
+ return angle >= THETA_MAX_CYLINDER && angle <= (M_PI - THETA_MAX_CYLINDER);
+ }
+ return false;
+}
+
+double AreasBuilder::distanceToPlane(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &point,
+ const std::array<double, 3> &normal)
+{
+ std::array<double, 3> vec{0.0, 0.0, 0.0};
+ for (size_t i = 0; i < nodeCoords.size(); ++i)
+ vec[i] = nodeCoords[i] - point[i];
+ return fabs(MathOps::dot(vec, normal)) / MathOps::computeNorm(normal);
+}
+
+double AreasBuilder::distanceToSphere(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> ¢er)
+{
+ return MathOps::computeNorm(
+ std::array<double, 3>{
+ nodeCoords[0] - center[0],
+ nodeCoords[1] - center[1],
+ nodeCoords[2] - center[2]});
+}
+
+double AreasBuilder::distanceToCylinder(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &axis,
+ const std::array<double, 3> &axisPoint)
+{
+
+ std::array<double, 3> pa = {
+ axisPoint[0] - nodeCoords[0],
+ axisPoint[1] - nodeCoords[1],
+ axisPoint[2] - nodeCoords[2]};
+ double innerProduct = MathOps::dot(pa, axis);
+ return MathOps::computeNorm(
+ std::array<double, 3>({pa[0] - innerProduct * axis[0],
+ pa[1] - innerProduct * axis[1],
+ pa[2] - innerProduct * axis[2]}));
+}
+
+double AreasBuilder::distanceToCone(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &axis,
+ const std::array<double, 3> &apex,
+ double angle,
+ double radius)
+{
+ std::array<double, 3> ps{
+ apex[0] - nodeCoords[0],
+ apex[1] - nodeCoords[1],
+ apex[2] - nodeCoords[2]};
+ std::array<double, 3> v(axis);
+ if (MathOps::dot(axis, ps) <= 0)
+ {
+ for (size_t i = 0; i < 3; ++i)
+ v[i] *= -1;
+ }
+ double a = MathOps::computeNorm(
+ MathOps::cross(ps, v));
+ double b = MathOps::dot(ps, v);
+ return fabs(a * cos(angle) - b * sin(angle)) / fabs(radius);
+}
--- /dev/null
+// 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 __AREASBUILDER_HXX__
+#define __AREASBUILDER_HXX__
+
+#include "Nodes.hxx"
+#include "Areas.hxx"
+
+namespace MEDCoupling
+{
+ class AreasBuilder
+ {
+ public:
+ AreasBuilder(const Nodes *nodes);
+
+ void explore();
+ void expand();
+ void rebuild();
+
+ Areas *getAreas() const;
+
+ private:
+ void exploreAreas();
+ void expandAreas();
+ void expandAreasByType(PrimitiveType primitive);
+ void rebuildInvalidAreas();
+ void filterHighPass();
+ bool doesItMatch(mcIdType areaId, mcIdType nodeId) const;
+ bool doesItBelong(mcIdType areaId, mcIdType nodeId) const;
+ bool isInvalidCylinderNode(mcIdType nodeId) const;
+
+ static double distanceToPlane(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &point,
+ const std::array<double, 3> &normal);
+ static double distanceToSphere(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> ¢er);
+ static double distanceToCylinder(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &axis,
+ const std::array<double, 3> &axisPoint);
+ static double distanceToCone(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &axis,
+ const std::array<double, 3> &apex,
+ double angle,
+ double radius);
+
+ const Nodes *nodes;
+ Areas *areas;
+
+ size_t threshold = 5;
+ };
+};
+
+#endif // __AREASBUILDER_HXX__
\ No newline at end of file
${LAPACKE_INCLUDE_DIRS}
${CMAKE_CURRENT_SOURCE_DIR}/../MEDCoupling
${CMAKE_CURRENT_SOURCE_DIR}/../MEDLoader
+ ${CMAKE_CURRENT_SOURCE_DIR}/../INTERP_KERNEL
+ ${CMAKE_CURRENT_SOURCE_DIR}/../INTERP_KERNEL/Bases
)
SET(shaperecogn_SOURCES
- )
+ MathOps.cxx
+ Nodes.cxx
+ NodeCurvatureCalculator.cxx
+ Areas.cxx
+ AreasBuilder.cxx
+ ShapeRecognMesh.cxx
+)
ADD_LIBRARY(shaperecogn ${shaperecogn_SOURCES})
SET_TARGET_PROPERTIES(shaperecogn PROPERTIES COMPILE_FLAGS "")
--- /dev/null
+// 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
+//
+
+#include "MathOps.hxx"
+
+#include <algorithm>
+#include <lapacke.h>
+#include <cblas.h>
+#include <iostream>
+#include <cfloat>
+
+using namespace MEDCoupling;
+
+std::vector<double> MathOps::lstsq(
+ std::vector<double> &a, const std::vector<double> &b)
+{
+ int m = (int)b.size();
+ int n = 3;
+ int nrhs = 1;
+ return lstsq(a, b, m, n, nrhs);
+}
+
+std::vector<double> MathOps::lstsq(
+ 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> s(lds, 0.0);
+ double rcond = DBL_EPSILON * (double)ldb; // same value as numpy.linalg.lstsq
+ int rank = 0;
+ int info = LAPACKE_dgelsd(
+ LAPACK_COL_MAJOR,
+ m, n, nrhs,
+ a.data(), m,
+ x.data(), ldb,
+ 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};
+ c[0] = a[1] * b[2] - a[2] * b[1];
+ c[1] = a[2] * b[0] - a[0] * b[2];
+ c[2] = a[0] * b[1] - a[1] * b[0];
+ return c;
+}
+
+std::array<double, 3> MathOps::normalize(const std::array<double, 3> &a)
+{
+ std::array<double, 3> an{0.0, 0.0, 0.0};
+ double n = computeNorm(a);
+ for (size_t i = 0; i < 3; i++)
+ an[i] = a[i] / n;
+ return an;
+}
+
+double MathOps::computeNorm(const std::array<double, 3> &a)
+{
+ double n = 0;
+ for (size_t i = 0; i < 3; i++)
+ n += a[i] * a[i];
+ return sqrt(n);
+}
+
+std::vector<double> MathOps::computeNorm(
+ const std::vector<double> &a)
+{
+ size_t s = a.size() / 3;
+ std::vector<double> n(s, 0.0);
+ for (size_t i = 0; i < s; i++)
+ {
+ for (size_t j = 0; j < 3; j++)
+ n[i] += a[3 * i + j] * a[3 * i + j];
+ n[i] = sqrt(n[i]);
+ }
+ return n;
+}
+
+double MathOps::dot(const std::array<double, 3> &a, const std::array<double, 3> &b)
+{
+ double d = 0.0;
+ for (size_t i = 0; i < 3; i++)
+ d += a[i] * b[i];
+ return d;
+}
+
+std::vector<double> MathOps::dot(const std::vector<double> &a, const std::array<double, 3> &b)
+{
+ size_t nbNodes = a.size() / 3;
+ std::vector<double> d(nbNodes, 0.0);
+ cblas_dgemv(
+ CBLAS_LAYOUT::CblasRowMajor,
+ CBLAS_TRANSPOSE::CblasNoTrans,
+ (int)nbNodes, 3, 1.0,
+ a.data(), 3,
+ b.data(), 1,
+ 0.0, d.data(), 1);
+ return d;
+}
+
+double MathOps::mean(const std::vector<double> &values)
+{
+ double mean = 0.0;
+ for (double value : values)
+ mean += value;
+ return mean / values.size();
+}
+
+std::array<double, 3> MathOps::meanCoordinates(const std::vector<double> &coordinates)
+{
+ std::array<double, 3> coordsMean{0.0, 0.0, 0.0};
+ size_t nbNodes = coordinates.size() / 3;
+ for (size_t nodeId = 0; nodeId < nbNodes; ++nodeId)
+ {
+ coordsMean[0] += coordinates[3 * nodeId];
+ coordsMean[1] += coordinates[3 * nodeId + 1];
+ coordsMean[2] += coordinates[3 * nodeId + 2];
+ }
+ coordsMean[0] /= nbNodes;
+ coordsMean[1] /= nbNodes;
+ coordsMean[2] /= nbNodes;
+ return coordsMean;
+}
+
+std::array<double, 9> MathOps::computeCov(
+ const std::vector<double> &coordinates)
+{
+ std::array<double, 9> covMatrix;
+ covMatrix.fill(0);
+ size_t nbNodes = coordinates.size() / 3;
+ // Center the coordinates
+ std::vector<double> coordsCentered(coordinates);
+ std::array<double, 3> coordsMean = meanCoordinates(coordinates);
+ for (size_t nodeId = 0; nodeId < nbNodes; ++nodeId)
+ {
+ coordsCentered[3 * nodeId] -= coordsMean[0];
+ coordsCentered[3 * nodeId + 1] -= coordsMean[1];
+ coordsCentered[3 * nodeId + 2] -= coordsMean[2];
+ }
+ cblas_dgemm(
+ CBLAS_LAYOUT::CblasColMajor,
+ CBLAS_TRANSPOSE::CblasNoTrans,
+ CBLAS_TRANSPOSE::CblasTrans,
+ 3, 3, 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;
+ }
+ return covMatrix;
+}
+
+std::array<double, 9> MathOps::computePCA(const std::vector<double> &coordinates)
+{
+ std::array<double, 9> covMatrix = computeCov(coordinates);
+ std::array<double, 3> eigenValues{0.0, 0.0, 0.0};
+ LAPACKE_dsyevd(
+ LAPACK_COL_MAJOR,
+ 'V',
+ 'U',
+ 3,
+ covMatrix.data(),
+ 3,
+ eigenValues.data());
+ return covMatrix;
+}
+
+std::array<double, 3> MathOps::computePCAFirstAxis(const std::vector<double> &coordinates)
+{
+ std::array<double, 9> pca = computePCA(coordinates);
+ // The eignvalues are in ascending order so the first axis correspond to the last vector
+ return {pca[6], pca[7], pca[8]};
+}
+
+std::array<double, 3> MathOps::computePCAThirdAxis(const std::vector<double> &coordinates)
+{
+ std::array<double, 9> pca = computePCA(coordinates);
+ // The eignvalues are in ascending order so the third axis correspond to the first vector
+ return {pca[0], pca[1], pca[2]};
+}
+
+double MathOps::computeQuantile(const std::vector<double> &values, double q)
+{
+ std::vector<double> sortedValues(values);
+ std::sort(sortedValues.begin(), sortedValues.end());
+ double pos = q * (sortedValues.size() - 1);
+ size_t index = static_cast<size_t>(pos);
+ if (pos == index)
+ return sortedValues[index];
+ else
+ {
+ double frac = pos - index;
+ return sortedValues[index] * (1 - frac) + sortedValues[index + 1] * frac;
+ }
+}
+
+double MathOps::computeAngle(const std::array<double, 3> &direction, std::array<double, 3> axis)
+{
+ double angle = dot(direction, axis);
+ double normAxis = computeNorm(axis);
+ double normDirection = computeNorm(direction);
+ angle /= normDirection * normAxis;
+ if (fabs(angle) >= 1)
+ return 0;
+ else
+ return acos(angle);
+}
+
+std::vector<double> MathOps::computeAngles(
+ const std::vector<double> &directions, std::array<double, 3> axis)
+{
+ size_t nbDirections = directions.size() / 3;
+ std::vector<double> angles(nbDirections, 0.0);
+ cblas_dgemv(
+ CBLAS_LAYOUT::CblasRowMajor,
+ CBLAS_TRANSPOSE::CblasNoTrans,
+ nbDirections, 3, 1.0,
+ directions.data(), 3,
+ axis.data(), 1,
+ 0.0, angles.data(), 1);
+ double normAxis = computeNorm(axis);
+ std::vector<double> normDirections = computeNorm(directions);
+ for (size_t i = 0; i < nbDirections; ++i)
+ {
+ angles[i] /= normAxis * normDirections[i];
+ if (fabs(angles[i]) >= 1.0)
+ angles[i] = 0.0;
+ angles[i] = acos(angles[i]);
+ }
+ return angles;
+}
+
+double MathOps::computeOrientedAngle(const std::array<double, 3> &normal, const std::array<double, 3> &vector1, const std::array<double, 3> &vector2)
+{
+ double angle = computeAngle(vector1, vector2);
+ if (dot(cross(vector1, vector2), normal) >= 0.0)
+ return angle;
+ else
+ return -angle;
+}
+
+double MathOps::computeVar(std::vector<double> values)
+{
+ int 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;
+}
--- /dev/null
+// 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 __MATHOPS_HXX__
+#define __MATHOPS_HXX__
+
+#include <vector>
+#include <array>
+
+namespace MEDCoupling
+{
+ class MathOps
+ {
+ public:
+ static std::vector<double> lstsq(std::vector<double> &a, const std::vector<double> &b);
+ static std::vector<double> lstsq(
+ std::vector<double> &a,
+ const std::vector<double> &b,
+ int m,
+ int n = 3,
+ int nrhs = 1);
+ 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);
+ static std::vector<double> computeNorm(const std::vector<double> &a);
+ static double dot(const std::array<double, 3> &a, const std::array<double, 3> &b);
+ static std::vector<double> dot(const std::vector<double> &a, const std::array<double, 3> &b);
+ static double mean(const std::vector<double> &values);
+ static std::array<double, 3> meanCoordinates(const std::vector<double> &coordinates);
+ static std::array<double, 9> computeCov(const std::vector<double> &coordinates);
+ static std::array<double, 9> computePCA(const std::vector<double> &coordinates);
+ static std::array<double, 3> computePCAFirstAxis(const std::vector<double> &coordinates);
+ static std::array<double, 3> computePCAThirdAxis(const std::vector<double> &coordinates);
+ static double computeQuantile(
+ const std::vector<double> &values,
+ double q);
+ static double computeAngle(
+ const std::array<double, 3> &direction,
+ std::array<double, 3> axis);
+ static std::vector<double> computeAngles(
+ const std::vector<double> &directions,
+ std::array<double, 3> axis);
+ static double computeOrientedAngle(
+ 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);
+ };
+};
+
+#endif //__MATHOPS_HXX__
--- /dev/null
+
+#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();
+ // std::cout << "computeNormals ended\n";
+ computeCurvatures();
+ // std::cout << "computeCurvatures ended\n";
+ return nodes;
+}
+
+void NodeCurvatureCalculator::computeNormals()
+{
+ mcIdType nbNodes = mesh->getNumberOfNodes();
+ mcIdType nbCells = mesh->getNumberOfCells();
+ std::array<double, 3> cellNormal;
+ cellNormal.fill(0.0);
+ // product of normal AreaByCell
+ std::vector<double> prodNormalAreaByCell(3 * nbCells, 0.0);
+ for (int cellId = 0; cellId < nbCells; cellId++)
+ {
+ std::vector<mcIdType> 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++)
+ {
+ // nodeIds[0] = nodeId;
+ // auto cellIds = mesh->getCellIdsLyingOnNodes(nodeIds.data(), nodeIds.data() + 1, false);
+ int nbCells = revNodalIdx->getIJ(nodeId + 1, 0) - revNodalIdx->getIJ(nodeId, 0);
+ std::vector<mcIdType> 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->theta0.resize(nbNodes);
+ nodes->isConvex.resize(nbNodes);
+ nodes->kdiff0.resize(nbNodes);
+ nodes->primitives.resize(nbNodes);
+ nodes->weakDirections.resize(3 * nbNodes);
+ for (mcIdType nodeId = 0; nodeId < nbNodes; nodeId++)
+ computeCurvatures(nodeId, tol);
+}
+
+void NodeCurvatureCalculator::computeCurvatures(mcIdType nodeId, double tol)
+{
+ std::array<double, 3> normal = nodes->getNormal(nodeId);
+ std::vector<mcIdType> neighborIds = nodes->getNeighbors(nodeId);
+ std::vector<double> discreteCurvatures = computeDiscreteCurvatures(nodeId, neighborIds);
+ std::vector<double> tangents = computeTangentDirections(nodeId, neighborIds);
+ mcIdType maxCurvatureId = std::distance(
+ discreteCurvatures.begin(),
+ std::max_element(discreteCurvatures.begin(), discreteCurvatures.end()));
+ std::array<double, 3> e1{
+ tangents[3 * maxCurvatureId],
+ tangents[3 * maxCurvatureId + 1],
+ tangents[3 * maxCurvatureId + 2]};
+ std::array<double, 3> e2 = MathOps::normalize(MathOps::cross(e1, normal));
+ std::vector<double> 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;
+ double k1 = h + sqrt(fabs(pow(h, 2) - kg));
+ double k2 = h - sqrt(fabs(pow(h, 2) - kg));
+ double theta0 = 0.0;
+ if (fabs(k1 - k2) >= tol)
+ theta0 = 0.5 * asin(b / (k1 - k2));
+ std::array<double, 3> direction1{0.0, 0.0, 0.0};
+ std::array<double, 3> 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];
+ }
+ bool isConvex;
+ double kdiff0;
+ double kis0;
+ std::array<double, 3> mainDirection{0.0, 0.0, 0.0};
+ std::array<double, 3> weakDirection{0.0, 0.0, 0.0};
+ if (fabs(k1) > fabs(k2))
+ {
+ isConvex = true;
+ kdiff0 = k1;
+ kis0 = k2;
+ mainDirection = direction1;
+ weakDirection = direction2;
+ }
+ else
+ {
+ isConvex = false;
+ kdiff0 = k2;
+ kis0 = k1;
+ mainDirection = direction2;
+ weakDirection = direction1;
+ }
+ PrimitiveType primitive = findPrimitiveType(k1, k2, kdiff0, kis0);
+ // Populate nodes
+ nodes->isConvex[nodeId] = isConvex;
+ nodes->theta0[nodeId] = theta0;
+ nodes->k1[nodeId] = k1;
+ nodes->k2[nodeId] = k2;
+ nodes->kdiff0[nodeId] = kdiff0;
+ for (size_t i = 0; i < 3; ++i)
+ nodes->weakDirections[3 * nodeId + i] = weakDirection[i];
+ nodes->primitives[nodeId] = primitive;
+}
+
+PrimitiveType NodeCurvatureCalculator::findPrimitiveType(double k1, double k2, double kdiff0, double kis0) const
+{
+ if ((fabs(k1 - k2) < EPSILON_SC) && (fabs((k1 + k2) / 2) < EPSILON_PL))
+ return PrimitiveType::Plane;
+ else if ((fabs(k1 - k2) < EPSILON_SC) && (fabs((k1 + k2) / 2) > EPSILON_PL))
+ return PrimitiveType::Sphere;
+ else if ((fabs(kdiff0) > EPSILON_CY) && (fabs(kis0) < EPSILON_CY))
+ return PrimitiveType::Cylinder;
+ else if ((fabs(k1) > EPSILON_CY) && (fabs(k2) > EPSILON_CY))
+ return PrimitiveType::Torus;
+ else
+ return PrimitiveType::Unknown;
+}
+
+std::vector<double> NodeCurvatureCalculator::computeNormalCurvatureCoefficients(
+ const std::vector<double> &discreteCurvatures,
+ const std::vector<double> &tangents,
+ const std::array<double, 3> &normal,
+ const std::array<double, 3> &e1) const
+{
+ size_t nbNeighbors = discreteCurvatures.size();
+ std::vector<double> a(3 * nbNeighbors, 0.0);
+ for (size_t i = 0; i < nbNeighbors; ++i)
+ {
+ std::array<double, 3> 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<mcIdType> &nodeIds,
+ std::array<double, 3> &cellNormal) const
+{
+ std::vector<double> point1;
+ std::vector<double> point2;
+ std::vector<double> point3;
+ mesh->getCoordinatesOfNode(nodeIds[0], point1);
+ mesh->getCoordinatesOfNode(nodeIds[1], point2);
+ mesh->getCoordinatesOfNode(nodeIds[2], point3);
+ std::array<double, 3> a;
+ a.fill(3);
+ std::array<double, 3> 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];
+}
+
+std::vector<double> NodeCurvatureCalculator::computeDiscreteCurvatures(
+ mcIdType nodeId,
+ const std::vector<mcIdType> &neighborIds) const
+{
+ std::vector<double> 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<double> NodeCurvatureCalculator::computeTangentDirections(
+ mcIdType nodeId,
+ const std::vector<mcIdType> &neighborIds) const
+{
+ size_t nbNeighbors = neighborIds.size();
+ std::vector<double> 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;
+}
--- /dev/null
+// 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 <string>
+#include <vector>
+#include <array>
+#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<double> computeNormalCurvatureCoefficients(
+ const std::vector<double> &discreteCurvatures,
+ const std::vector<double> &tangents,
+ const std::array<double, 3> &normal,
+ const std::array<double, 3> &e1) const;
+ void computeCellNormal(const std::vector<mcIdType> &nodeIds, std::array<double, 3> &cellNormal) const;
+ std::vector<double> computeDiscreteCurvatures(mcIdType nodeId, const std::vector<mcIdType> &neighborIds) const;
+ double computeDiscreteCurvature(mcIdType nodeId, mcIdType neighborId) const;
+ std::vector<double> computeTangentDirections(mcIdType nodeId, const std::vector<mcIdType> &neighborIds) const;
+
+ const MEDCouplingUMesh *mesh;
+ Nodes *nodes;
+ };
+};
+
+#endif //__NODECURVATURECALCULATOR_HXX__
\ No newline at end of file
--- /dev/null
+#include "Nodes.hxx"
+
+using namespace MEDCoupling;
+
+Nodes::Nodes(
+ const MEDCouplingUMesh *mesh,
+ const DataArrayInt64 *neighbors,
+ const DataArrayInt64 *neighborsIdx)
+ : mesh(mesh), coords(mesh->getCoords()), neighbors(neighbors), neighborsIdx(neighborsIdx)
+{
+}
+
+Nodes::~Nodes()
+{
+ mesh->decrRef();
+ neighbors->decrRef();
+ neighborsIdx->decrRef();
+}
+
+mcIdType Nodes::getNbNodes() const
+{
+ return coords->getNumberOfTuples();
+}
+
+const std::vector<double> &Nodes::getNormals() const
+{
+ return normals;
+}
+
+const std::vector<double> &Nodes::getWeakDirections() const
+{
+ return weakDirections;
+}
+
+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]};
+}
+
+double Nodes::getK1(mcIdType nodeId) const
+{
+ return k1[nodeId];
+}
+
+double Nodes::getK2(mcIdType nodeId) const
+{
+ return k2[nodeId];
+}
+
+double Nodes::getKdiff0(mcIdType nodeId) const
+{
+ return kdiff0[nodeId];
+}
+
+const std::vector<mcIdType> Nodes::getNeighbors(mcIdType nodeId) const
+{
+ mcIdType start = neighborsIdx->getIJ(nodeId, 0);
+ mcIdType nbNeighbors = neighborsIdx->getIJ(nodeId + 1, 0) - start;
+
+ std::vector<mcIdType> neighborNodes(nbNeighbors, 0);
+ for (mcIdType i = 0; i < nbNeighbors; i++)
+ neighborNodes[i] = neighbors->getIJ(start + i, 0);
+ return neighborNodes;
+}
+
+PrimitiveType Nodes::getPrimitiveType(mcIdType nodeId) const
+{
+ return primitives[nodeId];
+}
+
+std::array<double, 3> Nodes::getWeakDirection(mcIdType nodeId) const
+{
+ return {weakDirections[3 * nodeId],
+ weakDirections[3 * nodeId + 1],
+ weakDirections[3 * nodeId + 2]};
+}
+
+std::array<double, 3> Nodes::getCoordinates(mcIdType nodeId) const
+{
+ std::array<double, 3> nodeCoords;
+ coords->getTuple(nodeId, nodeCoords.data());
+ return nodeCoords;
+}
--- /dev/null
+// 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 __NODES_HXX__
+#define __NODES_HXX__
+
+#include <vector>
+#include "MEDCouplingUMesh.hxx"
+#include "PrimitiveType.hxx"
+
+namespace MEDCoupling
+{
+ class Nodes
+ {
+ public:
+ friend class NodeCurvatureCalculator;
+ Nodes(const MEDCouplingUMesh *mesh,
+ const DataArrayInt64 *neighbors,
+ const DataArrayInt64 *neighborsIdx);
+ ~Nodes();
+
+ mcIdType getNbNodes() const;
+ const std::vector<double> &getNormals() const;
+ const std::vector<double> &getWeakDirections() 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;
+ const std::vector<mcIdType> getNeighbors(mcIdType nodeId) const;
+ PrimitiveType getPrimitiveType(mcIdType nodeId) const;
+ std::array<double, 3> getWeakDirection(mcIdType nodeId) const;
+ std::array<double, 3> getCoordinates(mcIdType nodeId) const;
+
+ private:
+ const MEDCouplingUMesh *mesh;
+ const DataArrayDouble *coords;
+ // normals 3 * nbNodes
+ std::vector<double> normals;
+ // neighbors
+ 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> weakDirections;
+ std::vector<PrimitiveType> primitives;
+ };
+};
+
+#endif //__NODES_HXX__
--- /dev/null
+// 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 __PRIMITIVETYPE_HXX__
+#define __PRIMITIVETYPE_HXX__
+
+#include <string>
+namespace MEDCoupling
+{
+ enum PrimitiveType
+ {
+ Plane = 0,
+ Sphere = 1,
+ Cylinder = 2,
+ Cone = 3,
+ Torus = 4,
+ Unknown = 5
+ };
+
+ inline std::string convertPrimitiveToString(PrimitiveType type)
+ {
+ std::string typeName = "";
+ switch (type)
+ {
+ case PrimitiveType::Plane:
+ typeName = "Plane";
+ break;
+ case PrimitiveType::Sphere:
+ typeName = "Sphere";
+ break;
+ case PrimitiveType::Cylinder:
+ typeName = "Cylinder";
+ break;
+ case PrimitiveType::Cone:
+ typeName = "Cone";
+ break;
+ case PrimitiveType::Torus:
+ typeName = "Torus";
+ break;
+ case PrimitiveType::Unknown:
+ typeName = "Unknown";
+ break;
+ default:
+ break;
+ }
+ return typeName;
+ };
+
+ inline int convertPrimitiveToInt(PrimitiveType type)
+ {
+ int typeInt = 5;
+ switch (type)
+ {
+ case PrimitiveType::Plane:
+ typeInt = 0;
+ break;
+ case PrimitiveType::Sphere:
+ typeInt = 1;
+ break;
+ case PrimitiveType::Cylinder:
+ typeInt = 2;
+ break;
+ case PrimitiveType::Cone:
+ typeInt = 3;
+ break;
+ case PrimitiveType::Torus:
+ typeInt = 4;
+ break;
+ case PrimitiveType::Unknown:
+ default:
+ break;
+ }
+ return typeInt;
+ };
+};
+
+#endif // __PRIMITIVETYPE_HXX__
\ No newline at end of file
--- /dev/null
+#include "ShapeRecognMesh.hxx"
+
+#include "NodeCurvatureCalculator.hxx"
+#include "AreasBuilder.hxx"
+#include "MEDLoader.hxx"
+#include "MEDCouplingFieldDouble.hxx"
+
+using namespace MEDCoupling;
+
+ShapeRecognMesh::ShapeRecognMesh(const std::string &fileName)
+{
+ mesh = ReadUMeshFromFile(fileName);
+ // TODO: check if the elements of the mesh are triangle ?
+}
+
+ShapeRecognMesh::~ShapeRecognMesh()
+{
+ if (areas != 0)
+ delete areas;
+ if (nodes != 0)
+ delete nodes;
+ mesh->decrRef();
+}
+
+const Nodes *ShapeRecognMesh::getNodes() const
+{
+ return nodes;
+}
+
+const Areas *ShapeRecognMesh::getAreas() const
+{
+ return areas;
+}
+
+void ShapeRecognMesh::recognize()
+{
+ mesh->incrRef();
+ NodeCurvatureCalculator nodeCurvatureCalc(mesh);
+ nodes = nodeCurvatureCalc.compute();
+ AreasBuilder areasBuilder(nodes);
+ areasBuilder.explore();
+ areasBuilder.expand();
+ areasBuilder.rebuild();
+ areas = areasBuilder.getAreas();
+ // for (mcIdType areaId = 0; areaId < areas->getNumberOfAreas(); ++areaId)
+ // {
+ // std::cout << "area ("
+ // << areaId
+ // << ") : "
+ // << areas->getPrimitiveTypeName(areaId)
+ // << ", " << areas->getNumberOfNodes(areaId) << " "
+ // << "[";
+ // // for (mcIdType nodeId : areas->getNodeIds(areaId))
+ // // {
+ // // std::cout << nodeId << ",";
+ // // }
+ // std::cout << "]\n";
+ // }
+}
+
+void ShapeRecognMesh::recognize(const std::string &outputFile)
+{
+ recognize();
+ // Nodes
+ // - k1
+ MEDCouplingFieldDouble *fieldNodeK1 = MEDCouplingFieldDouble::New(
+ ON_NODES);
+ fieldNodeK1->setName("node_k1");
+ fieldNodeK1->setMesh(mesh);
+ DataArrayDouble *nodeK1 = DataArrayDouble::New();
+ nodeK1->setName("node_k1");
+ nodeK1->alloc(nodes->getNbNodes(), 1);
+ for (mcIdType nodeId = 0; nodeId < nodes->getNbNodes(); ++nodeId)
+ nodeK1->setIJ(nodeId, 0, nodes->getK1(nodeId));
+ fieldNodeK1->setArray(nodeK1);
+ nodeK1->decrRef();
+ WriteField(outputFile, fieldNodeK1, true);
+ fieldNodeK1->decrRef();
+ // - k2
+ MEDCouplingFieldDouble *fieldNodeK2 = MEDCouplingFieldDouble::New(
+ ON_NODES);
+ fieldNodeK2->setName("node_k2");
+ fieldNodeK2->setMesh(mesh);
+ DataArrayDouble *nodeK2 = DataArrayDouble::New();
+ nodeK2->setName("node_k2");
+ nodeK2->alloc(nodes->getNbNodes(), 1);
+ for (mcIdType nodeId = 0; nodeId < nodes->getNbNodes(); ++nodeId)
+ nodeK2->setIJ(nodeId, 0, nodes->getK2(nodeId));
+ fieldNodeK2->setArray(nodeK2);
+ nodeK2->decrRef();
+ WriteField(outputFile, fieldNodeK2, false);
+ fieldNodeK2->decrRef();
+ // - primitive types
+ MEDCouplingFieldDouble *fieldNodePrimitiveTypes = MEDCouplingFieldDouble::New(
+ ON_NODES);
+ fieldNodePrimitiveTypes->setName("node_primitive_types");
+ fieldNodePrimitiveTypes->setMesh(mesh);
+ DataArrayDouble *nodePrimitiveTypes = DataArrayDouble::New();
+ nodePrimitiveTypes->setName("node_primitive_types");
+ nodePrimitiveTypes->alloc(nodes->getNbNodes(), 1);
+ for (mcIdType nodeId = 0; nodeId < nodes->getNbNodes(); ++nodeId)
+ nodePrimitiveTypes->setIJ(nodeId, 0, (double)convertPrimitiveToInt(nodes->getPrimitiveType(nodeId)));
+ fieldNodePrimitiveTypes->setArray(nodePrimitiveTypes);
+ nodePrimitiveTypes->decrRef();
+ WriteField(outputFile, fieldNodePrimitiveTypes, false);
+ fieldNodePrimitiveTypes->decrRef();
+ // Areas
+ const std::vector<mcIdType> &areaIdByNodes = areas->getAreaIdByNodes();
+ // - area id
+ MEDCouplingFieldDouble *fieldAreaIds = MEDCouplingFieldDouble::New(
+ ON_NODES);
+ fieldAreaIds->setName("area_ids");
+ fieldAreaIds->setMesh(mesh);
+ DataArrayDouble *areaIds = DataArrayDouble::New();
+ areaIds->setName("area_ids");
+ areaIds->alloc(nodes->getNbNodes(), 1);
+ for (mcIdType nodeId = 0; nodeId < areaIdByNodes.size(); ++nodeId)
+ areaIds->setIJ(nodeId, 0, (double)areaIdByNodes[nodeId]);
+ fieldAreaIds->setArray(areaIds);
+ areaIds->decrRef();
+ WriteField(outputFile, fieldAreaIds, false);
+ fieldAreaIds->decrRef();
+ // - primitive types
+ MEDCouplingFieldDouble *fieldAreaPrimitiveTypes = MEDCouplingFieldDouble::New(
+ ON_NODES);
+ fieldAreaPrimitiveTypes->setName("area_primitive_types");
+ fieldAreaPrimitiveTypes->setMesh(mesh);
+ DataArrayDouble *areaPrimitiveTypes = DataArrayDouble::New();
+ areaPrimitiveTypes->setName("area_primitive_types");
+ areaPrimitiveTypes->alloc(nodes->getNbNodes(), 1);
+ for (mcIdType nodeId = 0; nodeId < areaIdByNodes.size(); ++nodeId)
+ areaPrimitiveTypes->setIJ(nodeId, 0, (double)areas->getPrimitiveTypeInt(areaIdByNodes[nodeId]));
+ fieldAreaPrimitiveTypes->setArray(areaPrimitiveTypes);
+ areaPrimitiveTypes->decrRef();
+ WriteField(outputFile, fieldAreaPrimitiveTypes, false);
+ fieldAreaPrimitiveTypes->decrRef();
+}
--- /dev/null
+// 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 __SHAPERECOGNMESH_HXX__
+#define __SHAPERECOGNMESH_HXX__
+
+#include <string>
+
+#include "MEDCouplingUMesh.hxx"
+
+namespace MEDCoupling
+{
+ class Nodes;
+ class Areas;
+
+ class ShapeRecognMesh
+ {
+ public:
+ ShapeRecognMesh(const std::string &fileName);
+ ~ShapeRecognMesh();
+
+ const Areas *getAreas() const;
+ const Nodes *getNodes() const;
+
+ void recognize();
+ void recognize(const std::string &outputFile);
+
+ private:
+ const MEDCouplingUMesh *mesh;
+ Nodes *nodes = 0;
+ Areas *areas = 0;
+ };
+};
+
+#endif // __SHAPERECOGNMESH_HXX__
--- /dev/null
+// 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 __SHAPERECOGNCONSTANTS_HXX__
+#define __SHAPERECOGNCONSTANTS_HXX__
+
+namespace MEDCoupling
+{
+ // Nodes
+ constexpr double EPSILON_SC = 0.005; // For |k1 - k2|
+ constexpr double EPSILON_PL = 0.001; // For |k1 + k2| / 2 = 1 / R
+ constexpr double EPSILON_CY = 0.005; // For k1 or k2
+ // Areas
+ constexpr double EPSILON_PROP = 0.05;
+ constexpr double DELTA_PLANE = 0.01;
+ constexpr double DELTA_SPHERE = 0.1;
+ constexpr double DELTA_CYLINDER = 0.02;
+ constexpr double DELTA_CONE = 0.05;
+ constexpr double THETA_MAX_CYLINDER = 0.05;
+ constexpr double TOL_MATCH_SPHERE = 0.1;
+ constexpr int THRESHOLD_MIN_NB_NODES = 5;
+ constexpr int THRESHOLD_MAX_NB_AREAS = 500;
+};
+
+#endif //__SHAPERECOGNCONSTANTS_HXX__
\ No newline at end of file
