From 5d7f47ea4e7e7b3da0cf3f329b18c85bc765b40d Mon Sep 17 00:00:00 2001 From: El Hadi Moussi Date: Mon, 5 Aug 2024 11:34:42 +0200 Subject: [PATCH] Add ShapeRecogn sources --- src/ShapeRecogn/Areas.cxx | 455 ++++++++++++++++++++ src/ShapeRecogn/Areas.hxx | 104 +++++ src/ShapeRecogn/AreasBuilder.cxx | 414 ++++++++++++++++++ src/ShapeRecogn/AreasBuilder.hxx | 74 ++++ src/ShapeRecogn/CMakeLists.txt | 10 +- src/ShapeRecogn/MathOps.cxx | 272 ++++++++++++ src/ShapeRecogn/MathOps.hxx | 67 +++ src/ShapeRecogn/NodeCurvatureCalculator.cxx | 264 ++++++++++++ src/ShapeRecogn/NodeCurvatureCalculator.hxx | 60 +++ src/ShapeRecogn/Nodes.cxx | 88 ++++ src/ShapeRecogn/Nodes.hxx | 71 +++ src/ShapeRecogn/PrimitiveType.hxx | 93 ++++ src/ShapeRecogn/ShapeRecognMesh.cxx | 137 ++++++ src/ShapeRecogn/ShapeRecognMesh.hxx | 51 +++ src/ShapeRecogn/ShapeRecongConstants.hxx | 41 ++ 15 files changed, 2200 insertions(+), 1 deletion(-) create mode 100644 src/ShapeRecogn/Areas.cxx create mode 100644 src/ShapeRecogn/Areas.hxx create mode 100644 src/ShapeRecogn/AreasBuilder.cxx create mode 100644 src/ShapeRecogn/AreasBuilder.hxx create mode 100644 src/ShapeRecogn/MathOps.cxx create mode 100644 src/ShapeRecogn/MathOps.hxx create mode 100644 src/ShapeRecogn/NodeCurvatureCalculator.cxx create mode 100644 src/ShapeRecogn/NodeCurvatureCalculator.hxx create mode 100644 src/ShapeRecogn/Nodes.cxx create mode 100644 src/ShapeRecogn/Nodes.hxx create mode 100644 src/ShapeRecogn/PrimitiveType.hxx create mode 100644 src/ShapeRecogn/ShapeRecognMesh.cxx create mode 100644 src/ShapeRecogn/ShapeRecognMesh.hxx create mode 100644 src/ShapeRecogn/ShapeRecongConstants.hxx diff --git a/src/ShapeRecogn/Areas.cxx b/src/ShapeRecogn/Areas.cxx new file mode 100644 index 000000000..5c2559f8c --- /dev/null +++ b/src/ShapeRecogn/Areas.cxx @@ -0,0 +1,455 @@ +#include "Areas.hxx" +#include "MathOps.hxx" + +#include + +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 &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 &Areas::getNormal(mcIdType areaId) const +{ + return areas[areaId].normal; +} + +const std::array &Areas::getCenter(mcIdType areaId) const +{ + const Area &area = areas[areaId]; + return area.center; +} + +const std::array &Areas::getAxis(mcIdType areaId) const +{ + return areas[areaId].axis; +} + +const std::array &Areas::getAxisPoint(mcIdType areaId) const +{ + return areas[areaId].axisPoint; +} + +const std::array &Areas::getApex(mcIdType areaId) const +{ + return areas[areaId].apex; +} + +std::array 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 &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 &normals = nodes->getNormals(); + area.radius = (2 / (area.k1 + area.k2)); + std::array center{0.0, 0.0, 0.0}; + if (!area.nodeIds.empty()) + { + size_t nbNodes = area.nodeIds.size(); + for (mcIdType nodeId : area.nodeIds) + { + std::array 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 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 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 projectedNodes(3 * nbNodes, 0.0); + std::vector radiusNodes(nbNodes, 0.0); + area.axisPoint.fill(0.0); + for (size_t i = 0; i < nbNodes; ++i) + { + mcIdType nodeId = area.nodeIds[i]; + std::array 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 &weakDirections = nodes->getWeakDirections(); + std::vector 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 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 q1_indices; + for (mcIdType idx = 0; idx < (mcIdType)radiusNodes.size(); ++idx) + { + if (radiusNodes[idx] < q1) + q1_indices.push_back(idx); + } + std::vector 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 heights(q1_indices.size(), 0.0); + // std::vector distancesToApex(q1_indices.size(), 0.0); + std::array 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 areaNodesK1(n, 0.0); + std::vector 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 majorRadiusNodes(3 * n, 0.0); + for (size_t i = 0; i < n; ++i) + { + std::array coords = nodes->getCoordinates(area.nodeIds[i]); + std::array 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 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 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 &Areas::getAreaIdByNodes() const +{ + return areaIdByNodes; +} diff --git a/src/ShapeRecogn/Areas.hxx b/src/ShapeRecogn/Areas.hxx new file mode 100644 index 000000000..e0d8ae440 --- /dev/null +++ b/src/ShapeRecogn/Areas.hxx @@ -0,0 +1,104 @@ +// 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 + +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 normal{0.0, 0.0, 0.0}; + std::array center{0.0, 0.0, 0.0}; + std::array axis{0.0, 0.0, 0.0}; + std::array axisPoint{0.0, 0.0, 0.0}; + std::array apex{0.0, 0.0, 0.0}; + std::vector 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 &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 &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 &getNormal(mcIdType areaId) const; + const std::array &getCenter(mcIdType areaId) const; + const std::array &getAxis(mcIdType areaId) const; + const std::array &getAxisPoint(mcIdType areaId) const; + const std::array &getApex(mcIdType areaId) const; + std::array 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 areas; + const Nodes *nodes; + std::vector areaIdByNodes; + }; +}; + +#endif // __AREAS_HXX__ diff --git a/src/ShapeRecogn/AreasBuilder.cxx b/src/ShapeRecogn/AreasBuilder.cxx new file mode 100644 index 000000000..125c45bdb --- /dev/null +++ b/src/ShapeRecogn/AreasBuilder.cxx @@ -0,0 +1,414 @@ +#include "AreasBuilder.hxx" +#include "MathOps.hxx" +#include "ShapeRecongConstants.hxx" + +#include + +using namespace MEDCoupling; + +AreasBuilder::AreasBuilder(const Nodes *nodes) : nodes(nodes), areas(new Areas(nodes)) +{ + size_t nbNodes = nodes->getNbNodes(); + threshold = std::max(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 exploredNodeIds(nodes->getNbNodes(), false); + std::unordered_set 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 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 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 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 exploredNodeIds(nodes->getNbNodes(), false); + areas->computeProperties(areaId); + const std::vector &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 neighborIds = nodes->getNeighbors(nodeId); + for (mcIdType neighborId : neighborIds) + { + if (!exploredNodeIds[neighborId]) + nodesToExplore.insert(neighborId); + } + } + exploredNodeIds[nodeId] = true; + } + } + } +} + +void AreasBuilder::rebuildInvalidAreas() +{ + std::vector exploredNodeIds(nodes->getNbNodes(), false); + std::vector isIinvalidNodes(nodes->getNbNodes(), false); + std::unordered_set 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 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 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 &nodeCoords, + const std::array &point, + const std::array &normal) +{ + std::array 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 &nodeCoords, + const std::array ¢er) +{ + return MathOps::computeNorm( + std::array{ + nodeCoords[0] - center[0], + nodeCoords[1] - center[1], + nodeCoords[2] - center[2]}); +} + +double AreasBuilder::distanceToCylinder( + const std::array &nodeCoords, + const std::array &axis, + const std::array &axisPoint) +{ + + std::array pa = { + axisPoint[0] - nodeCoords[0], + axisPoint[1] - nodeCoords[1], + axisPoint[2] - nodeCoords[2]}; + double innerProduct = MathOps::dot(pa, axis); + return MathOps::computeNorm( + std::array({pa[0] - innerProduct * axis[0], + pa[1] - innerProduct * axis[1], + pa[2] - innerProduct * axis[2]})); +} + +double AreasBuilder::distanceToCone( + const std::array &nodeCoords, + const std::array &axis, + const std::array &apex, + double angle, + double radius) +{ + std::array ps{ + apex[0] - nodeCoords[0], + apex[1] - nodeCoords[1], + apex[2] - nodeCoords[2]}; + std::array 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); +} diff --git a/src/ShapeRecogn/AreasBuilder.hxx b/src/ShapeRecogn/AreasBuilder.hxx new file mode 100644 index 000000000..eb97c5714 --- /dev/null +++ b/src/ShapeRecogn/AreasBuilder.hxx @@ -0,0 +1,74 @@ +// 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 &nodeCoords, + const std::array &point, + const std::array &normal); + static double distanceToSphere( + const std::array &nodeCoords, + const std::array ¢er); + static double distanceToCylinder( + const std::array &nodeCoords, + const std::array &axis, + const std::array &axisPoint); + static double distanceToCone( + const std::array &nodeCoords, + const std::array &axis, + const std::array &apex, + double angle, + double radius); + + const Nodes *nodes; + Areas *areas; + + size_t threshold = 5; + }; +}; + +#endif // __AREASBUILDER_HXX__ \ No newline at end of file diff --git a/src/ShapeRecogn/CMakeLists.txt b/src/ShapeRecogn/CMakeLists.txt index c5717d293..5dac83f40 100644 --- a/src/ShapeRecogn/CMakeLists.txt +++ b/src/ShapeRecogn/CMakeLists.txt @@ -37,10 +37,18 @@ INCLUDE_DIRECTORIES( ${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 "") diff --git a/src/ShapeRecogn/MathOps.cxx b/src/ShapeRecogn/MathOps.cxx new file mode 100644 index 000000000..910b8cbe6 --- /dev/null +++ b/src/ShapeRecogn/MathOps.cxx @@ -0,0 +1,272 @@ +// 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 +#include +#include +#include +#include + +using namespace MEDCoupling; + +std::vector MathOps::lstsq( + std::vector &a, const std::vector &b) +{ + int m = (int)b.size(); + int n = 3; + int nrhs = 1; + return lstsq(a, b, m, n, nrhs); +} + +std::vector MathOps::lstsq( + std::vector &a, const std::vector &b, int m, int n, int nrhs) +{ + int ldb = std::max(m, n); + int lds = std::min(m, n); + std::vector x(b.begin(), b.end()); + std::vector 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 MathOps::cross(const std::array &a, const std::array &b) +{ + std::array 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 MathOps::normalize(const std::array &a) +{ + std::array 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 &a) +{ + double n = 0; + for (size_t i = 0; i < 3; i++) + n += a[i] * a[i]; + return sqrt(n); +} + +std::vector MathOps::computeNorm( + const std::vector &a) +{ + size_t s = a.size() / 3; + std::vector 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 &a, const std::array &b) +{ + double d = 0.0; + for (size_t i = 0; i < 3; i++) + d += a[i] * b[i]; + return d; +} + +std::vector MathOps::dot(const std::vector &a, const std::array &b) +{ + size_t nbNodes = a.size() / 3; + std::vector 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 &values) +{ + double mean = 0.0; + for (double value : values) + mean += value; + return mean / values.size(); +} + +std::array MathOps::meanCoordinates(const std::vector &coordinates) +{ + std::array 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 MathOps::computeCov( + const std::vector &coordinates) +{ + std::array covMatrix; + covMatrix.fill(0); + size_t nbNodes = coordinates.size() / 3; + // Center the coordinates + std::vector coordsCentered(coordinates); + std::array 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 MathOps::computePCA(const std::vector &coordinates) +{ + std::array covMatrix = computeCov(coordinates); + std::array eigenValues{0.0, 0.0, 0.0}; + LAPACKE_dsyevd( + LAPACK_COL_MAJOR, + 'V', + 'U', + 3, + covMatrix.data(), + 3, + eigenValues.data()); + return covMatrix; +} + +std::array MathOps::computePCAFirstAxis(const std::vector &coordinates) +{ + std::array 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 MathOps::computePCAThirdAxis(const std::vector &coordinates) +{ + std::array 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 &values, double q) +{ + std::vector sortedValues(values); + std::sort(sortedValues.begin(), sortedValues.end()); + double pos = q * (sortedValues.size() - 1); + size_t index = static_cast(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 &direction, std::array 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 MathOps::computeAngles( + const std::vector &directions, std::array axis) +{ + size_t nbDirections = directions.size() / 3; + std::vector 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 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 &normal, const std::array &vector1, const std::array &vector2) +{ + double angle = computeAngle(vector1, vector2); + if (dot(cross(vector1, vector2), normal) >= 0.0) + return angle; + else + return -angle; +} + +double MathOps::computeVar(std::vector 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; +} diff --git a/src/ShapeRecogn/MathOps.hxx b/src/ShapeRecogn/MathOps.hxx new file mode 100644 index 000000000..267f9170d --- /dev/null +++ b/src/ShapeRecogn/MathOps.hxx @@ -0,0 +1,67 @@ +// 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 +#include + +namespace MEDCoupling +{ + class MathOps + { + public: + static std::vector lstsq(std::vector &a, const std::vector &b); + static std::vector lstsq( + std::vector &a, + const std::vector &b, + int m, + int n = 3, + int nrhs = 1); + static std::array cross(const std::array &a, const std::array &b); + static std::array normalize(const std::array &a); + static double computeNorm(const std::array &a); + static std::vector computeNorm(const std::vector &a); + static double dot(const std::array &a, const std::array &b); + static std::vector dot(const std::vector &a, const std::array &b); + static double mean(const std::vector &values); + static std::array meanCoordinates(const std::vector &coordinates); + static std::array computeCov(const std::vector &coordinates); + static std::array computePCA(const std::vector &coordinates); + static std::array computePCAFirstAxis(const std::vector &coordinates); + static std::array computePCAThirdAxis(const std::vector &coordinates); + static double computeQuantile( + const std::vector &values, + double q); + static double computeAngle( + const std::array &direction, + std::array axis); + static std::vector computeAngles( + const std::vector &directions, + std::array axis); + static double computeOrientedAngle( + const std::array &normal, + const std::array &vector1, + const std::array &vector2); + static double computeVar(std::vector values); + }; +}; + +#endif //__MATHOPS_HXX__ diff --git a/src/ShapeRecogn/NodeCurvatureCalculator.cxx b/src/ShapeRecogn/NodeCurvatureCalculator.cxx new file mode 100644 index 000000000..46b545415 --- /dev/null +++ b/src/ShapeRecogn/NodeCurvatureCalculator.cxx @@ -0,0 +1,264 @@ + +#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 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++) + { + // 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 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 normal = nodes->getNormal(nodeId); + std::vector neighborIds = nodes->getNeighbors(nodeId); + 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; + 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 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]; + } + bool isConvex; + double kdiff0; + double kis0; + std::array mainDirection{0.0, 0.0, 0.0}; + std::array 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 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]; +} + +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 new file mode 100644 index 000000000..8b3a58515 --- /dev/null +++ b/src/ShapeRecogn/NodeCurvatureCalculator.hxx @@ -0,0 +1,60 @@ +// 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; + 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.cxx b/src/ShapeRecogn/Nodes.cxx new file mode 100644 index 000000000..845e48b5e --- /dev/null +++ b/src/ShapeRecogn/Nodes.cxx @@ -0,0 +1,88 @@ +#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 &Nodes::getNormals() const +{ + return normals; +} + +const std::vector &Nodes::getWeakDirections() const +{ + return weakDirections; +} + +const std::vector &Nodes::getK1() const +{ + return k1; +} + +std::array Nodes::getNormal(mcIdType nodeId) const +{ + return {normals[3 * nodeId], normals[3 * nodeId + 1], normals[3 * nodeId + 2]}; +} + +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 Nodes::getNeighbors(mcIdType nodeId) const +{ + mcIdType start = neighborsIdx->getIJ(nodeId, 0); + mcIdType nbNeighbors = neighborsIdx->getIJ(nodeId + 1, 0) - start; + + std::vector 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 Nodes::getWeakDirection(mcIdType nodeId) const +{ + return {weakDirections[3 * nodeId], + weakDirections[3 * nodeId + 1], + weakDirections[3 * nodeId + 2]}; +} + +std::array Nodes::getCoordinates(mcIdType nodeId) const +{ + std::array nodeCoords; + coords->getTuple(nodeId, nodeCoords.data()); + return nodeCoords; +} diff --git a/src/ShapeRecogn/Nodes.hxx b/src/ShapeRecogn/Nodes.hxx new file mode 100644 index 000000000..ed8d267c2 --- /dev/null +++ b/src/ShapeRecogn/Nodes.hxx @@ -0,0 +1,71 @@ +// 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 +#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 &getNormals() const; + const std::vector &getWeakDirections() const; + const std::vector &getK1() const; + + std::array getNormal(mcIdType nodeId) const; + double getK1(mcIdType nodeId) const; + double getK2(mcIdType nodeId) const; + double getKdiff0(mcIdType nodeId) const; + const std::vector getNeighbors(mcIdType nodeId) const; + PrimitiveType getPrimitiveType(mcIdType nodeId) const; + std::array getWeakDirection(mcIdType nodeId) const; + std::array getCoordinates(mcIdType nodeId) const; + + private: + const MEDCouplingUMesh *mesh; + const DataArrayDouble *coords; + // normals 3 * nbNodes + std::vector normals; + // neighbors + const DataArrayInt64 *neighbors; + const DataArrayInt64 *neighborsIdx; + // curvature + std::vector isConvex; + std::vector theta0; + std::vector k1; + std::vector k2; + std::vector kdiff0; + std::vector weakDirections; + std::vector primitives; + }; +}; + +#endif //__NODES_HXX__ diff --git a/src/ShapeRecogn/PrimitiveType.hxx b/src/ShapeRecogn/PrimitiveType.hxx new file mode 100644 index 000000000..b131da84a --- /dev/null +++ b/src/ShapeRecogn/PrimitiveType.hxx @@ -0,0 +1,93 @@ +// 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 +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 diff --git a/src/ShapeRecogn/ShapeRecognMesh.cxx b/src/ShapeRecogn/ShapeRecognMesh.cxx new file mode 100644 index 000000000..8dc73c44d --- /dev/null +++ b/src/ShapeRecogn/ShapeRecognMesh.cxx @@ -0,0 +1,137 @@ +#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 &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(); +} diff --git a/src/ShapeRecogn/ShapeRecognMesh.hxx b/src/ShapeRecogn/ShapeRecognMesh.hxx new file mode 100644 index 000000000..d4532e752 --- /dev/null +++ b/src/ShapeRecogn/ShapeRecognMesh.hxx @@ -0,0 +1,51 @@ +// 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 + +#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__ diff --git a/src/ShapeRecogn/ShapeRecongConstants.hxx b/src/ShapeRecogn/ShapeRecongConstants.hxx new file mode 100644 index 000000000..2c36f19fb --- /dev/null +++ b/src/ShapeRecogn/ShapeRecongConstants.hxx @@ -0,0 +1,41 @@ +// 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 -- 2.39.2