From: El Hadi Moussi Date: Mon, 5 Aug 2024 09:34:26 +0000 (+0200) Subject: Reference version of Shape Recognition X-Git-Tag: V9_14_0a1~23 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=2ef39de72fbd6833f71756498a26c2891500776f;p=tools%2Fmedcoupling.git Reference version of Shape Recognition --- diff --git a/CMakeLists.txt b/CMakeLists.txt index 840bac26a..c080f174a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -88,6 +88,7 @@ OPTION(MEDCOUPLING_MICROMED "Build MED without MED file dependency." OFF) OPTION(MEDCOUPLING_ENABLE_PYTHON "Build PYTHON bindings." ON) OPTION(MEDCOUPLING_ENABLE_PARTITIONER "Build MEDPartitioner." ON) OPTION(MEDCOUPLING_ENABLE_RENUMBER "Build Renumber." ON) +OPTION(MEDCOUPLING_ENABLE_SHAPERECOGN "Build ShapeRecogn" OFF) OPTION(MEDCOUPLING_WITH_FILE_EXAMPLES "Install examples of files containing meshes and fields of different formats." ON) OPTION(MEDCOUPLING_USE_MPI "(Use MPI containers) - For MED this triggers the build of ParaMEDMEM." OFF) OPTION(MEDCOUPLING_BUILD_TESTS "Build MEDCoupling C++ tests." ON) @@ -201,6 +202,21 @@ IF(MEDCOUPLING_ENABLE_RENUMBER) SALOME_LOG_OPTIONAL_PACKAGE(Boost MEDCOUPLING_ENABLE_RENUMBER) ENDIF(MEDCOUPLING_ENABLE_RENUMBER) +IF(MEDCOUPLING_ENABLE_SHAPERECOGN) + FIND_PACKAGE(BLAS REQUIRED) + FIND_PACKAGE(LAPACK REQUIRED) + FIND_LIBRARY(LAPACKE_LIB NAMES lapacke REQUIRED) + SET(LAPACK_LIBRARIES ${LAPACKE_LIB} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) + FIND_PATH(LAPACKE_INCLUDE_DIRS NAMES lapacke.h HINTS ${LAPACK_LIBRARIES}) + IF(LAPACK_FOUND) + MESSAGE(STATUS "Lapacke libraries: ${LAPACK_LIBRARIES}") + MESSAGE(STATUS "Lapacke include dirs: ${LAPACKE_INCLUDE_DIRS}") + ELSE() + MESSAGE(FATAL_ERROR "Error in Lapacke detection ! lapacke not found !") + ENDIF(LAPACK_FOUND) + SALOME_LOG_OPTIONAL_PACKAGE(Lapack MEDCOUPLING_ENABLE_SHAPERECOGN) +ENDIF(MEDCOUPLING_ENABLE_SHAPERECOGN) + IF(MEDCOUPLING_ENABLE_PYTHON) FIND_PACKAGE(SalomePythonInterp) FIND_PACKAGE(SalomePythonLibs) @@ -339,6 +355,10 @@ IF(MEDCOUPLING_USE_MPI) ENDIF() ENDIF() ENDIF() + +IF(MEDCOUPLING_ENABLE_SHAPERECOGN) + LIST(APPEND _${PROJECT_NAME}_exposed_targets shaperecogn) +ENDIF(MEDCOUPLING_ENABLE_SHAPERECOGN) # Add all targets to the build-tree export set diff --git a/resources/CMakeLists.txt b/resources/CMakeLists.txt index 704221ca6..cf6ec2481 100644 --- a/resources/CMakeLists.txt +++ b/resources/CMakeLists.txt @@ -94,6 +94,11 @@ SET(MED_other_FILES test_MED_MAIL.sauv castem17_result_xdr.sauv castem17_result_ascii.sauv + ShapeRecognPlane.med + ShapeRecognCylinder.med + ShapeRecognCone.med + ShapeRecognSphere.med + ShapeRecognTorus.med ) SET(MED_RESOURCES_FILES ${MED_test_fig_files}) diff --git a/resources/ShapeRecognCone.med b/resources/ShapeRecognCone.med new file mode 100644 index 000000000..da32d3e8b Binary files /dev/null and b/resources/ShapeRecognCone.med differ diff --git a/resources/ShapeRecognCylinder.med b/resources/ShapeRecognCylinder.med new file mode 100755 index 000000000..f028da1ff Binary files /dev/null and b/resources/ShapeRecognCylinder.med differ diff --git a/resources/ShapeRecognPlane.med b/resources/ShapeRecognPlane.med new file mode 100755 index 000000000..c06ecfbf5 Binary files /dev/null and b/resources/ShapeRecognPlane.med differ diff --git a/resources/ShapeRecognSphere.med b/resources/ShapeRecognSphere.med new file mode 100755 index 000000000..4d7e7008f Binary files /dev/null and b/resources/ShapeRecognSphere.med differ diff --git a/resources/ShapeRecognTorus.med b/resources/ShapeRecognTorus.med new file mode 100755 index 000000000..7623db578 Binary files /dev/null and b/resources/ShapeRecognTorus.med differ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 14229c770..e4d70409c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -82,6 +82,10 @@ IF(MEDCOUPLING_USE_MPI) ENDIF(MEDCOUPLING_BUILD_TESTS) ENDIF(MEDCOUPLING_USE_MPI) +IF(MEDCOUPLING_ENABLE_SHAPERECOGN) + ADD_SUBDIRECTORY(ShapeRecogn) +ENDIF(MEDCOUPLING_ENABLE_SHAPERECOGN) + # Application tests configure_file(CTestTestfileInstall.cmake.in "CTestTestfileST.cmake" @ONLY) install(FILES ${CMAKE_CURRENT_BINARY_DIR}/CTestTestfileST.cmake DESTINATION ${MEDCOUPLING_INSTALL_TESTS} RENAME CTestTestfile.cmake) diff --git a/src/ShapeRecogn/Areas.cxx b/src/ShapeRecogn/Areas.cxx new file mode 100644 index 000000000..f8441cb48 --- /dev/null +++ b/src/ShapeRecogn/Areas.cxx @@ -0,0 +1,485 @@ +#include "Areas.hxx" +#include "MathOps.hxx" + +#include +#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 = ((double)(nbNodes - 1) * area.k1 + nodes->getK1(nodeId)) / (double)nbNodes; + area.k2 = ((double)(nbNodes - 1) * area.k2 + nodes->getK2(nodeId)) / (double)nbNodes; + area.adimK1 = ((double)(nbNodes - 1) * area.adimK1 + nodes->getAdimK1(nodeId)) / (double)nbNodes; + area.adimK2 = ((double)(nbNodes - 1) * area.adimK2 + nodes->getAdimK2(nodeId)) / (double)nbNodes; + area.adimKdiff0 = ((double)(nbNodes - 1) * area.adimKdiff0 + nodes->getAdimKdiff0(nodeId)) / (double)nbNodes; +} + +void Areas::cleanInvalidNodeAreas() +{ + 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::getAdimK1(mcIdType areaId) const +{ + + return areas[areaId].adimK1; +} + +double Areas::getAdimK2(mcIdType areaId) const +{ + return areas[areaId].adimK2; +} + +double Areas::getAdimKdiff0(mcIdType areaId) const +{ + return areas[areaId].adimKdiff0; +} + +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::getMinorRadius(mcIdType areaId) const +{ + const Area &area = areas[areaId]; + return area.minorRadius; +} + +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.adimK1 = 0.0; + area.adimK2 = 0.0; + area.adimKdiff0 = 0.0; + area.minorRadius = 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] /= (double)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] /= (double)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] /= (double)nbNodes; + // Radius of the cylinder is the mean of the approximate radius of each node + area.radius = fabs(area.radius / (double)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] /= (double)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}; + size_t min_q_size = std::min(q1_indices.size(), q2_indices.size()); + for (size_t i = 0; i < min_q_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] /= (double)min_q_size; +} + +void Areas::computeTorusProperties(mcIdType areaId) +{ + Area &area = areas[areaId]; + size_t n = area.nodeIds.size(); + if (n == 0) + return; + 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::computeVariance(areaNodesK1); + double var2 = MathOps::computeVariance(areaNodesK2); + double minorCurvature; + if (var1 > var2) + minorCurvature = MathOps::mean(areaNodesK2); + else + minorCurvature = MathOps::mean(areaNodesK1); + area.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] * area.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); + std::array base2d = MathOps::computeBaseFromNormal(normal); + std::vector projectedMajorRadiusNodes(2 * n, 0.0); + cblas_dgemm( + CBLAS_LAYOUT::CblasRowMajor, + CBLAS_TRANSPOSE::CblasNoTrans, + CBLAS_TRANSPOSE::CblasTrans, + (int)n, 2, 3, + 1.0, + majorRadiusNodes.data(), 3, + base2d.data(), 3, + 0.0, projectedMajorRadiusNodes.data(), 2); + std::vector A(3 * n, 1.0); + std::vector B(n, 0.0); + for (size_t i = 0; i < n; ++i) + { + for (size_t j = 0; j < 2; ++j) + A[3 * i + j] = projectedMajorRadiusNodes[2 * i + j]; + B[i] = pow(projectedMajorRadiusNodes[2 * i], 2) + + pow(projectedMajorRadiusNodes[2 * i + 1], 2); + } + std::vector fit = MathOps::lstsqRow(A, B); + double a = fit[0]; + double b = fit[1]; + double c = fit[2]; + double xc = a / 2.0; + double yc = b / 2.0; + area.radius = sqrt(4.0 * c + pow(a, 2) + pow(b, 2)) / 2.0; + for (size_t i = 0; i < 3; ++i) + area.center[i] = xc * base2d[i] + yc * base2d[3 + i] + meanMajorRadiusNodes[i]; +} + +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..b1b6f6bca --- /dev/null +++ b/src/ShapeRecogn/Areas.hxx @@ -0,0 +1,108 @@ +// 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 adimK1 = 0.0; + double adimK2 = 0.0; + double adimKdiff0 = 0.0; + double minorRadius = 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 getAdimK1(mcIdType areaId) const; + double getAdimK2(mcIdType areaId) const; + double getAdimKdiff0(mcIdType areaId) const; + + void computeProperties(mcIdType areaId); + double getMinorRadius(mcIdType areaId) const; + 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..947a54ed0 --- /dev/null +++ b/src/ShapeRecogn/AreasBuilder.cxx @@ -0,0 +1,371 @@ +#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::build() +{ + explore(); + expand(); + rebuild(); +} + +void AreasBuilder::explore() +{ + exploreAreas(); + areas->cleanInvalidNodeAreas(); + filterHighPass(); +} + +void AreasBuilder::expand() +{ + expandAreas(); + filterHighPass(); +} + +void AreasBuilder::rebuild() +{ + rebuildInvalidAreas(); + filterHighPass(); + expandAreasByType(PrimitiveType::Cone); + filterHighPass(); +} + +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); + 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)) + { + 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->getAdimK1(areaId) + areas->getAdimK2(areaId)) / 2.0); + double nodeKmoy = fabs( + (nodes->getAdimK1(nodeId) + nodes->getAdimK2(nodeId)) / 2.0); + isMatching = fabs((nodeKmoy - kmoy)) < TOL_MATCH_SPHERE; + } + break; + case PrimitiveType::Cylinder: + isMatching = fabs(areas->getAdimKdiff0(areaId) - + nodes->getAdimKdiff0(nodeId)) < TOL_MATCH_CYLINDER; + 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)); + double radius = areas->getRadius(areaId); + isClose = fabs((distanceToCenter - radius) / radius) < 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: + { + double radius = areas->getRadius(areaId); + isClose = distanceToCone( + nodes->getCoordinates(nodeId), + areas->getAxis(areaId), + areas->getApex(areaId), + areas->getAngle(areaId)) / + fabs(radius) < + 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) +{ + 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)); +} diff --git a/src/ShapeRecogn/AreasBuilder.hxx b/src/ShapeRecogn/AreasBuilder.hxx new file mode 100644 index 000000000..9742af648 --- /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 build(); + + Areas *getAreas() const; + + private: + void explore(); + void expand(); + void rebuild(); + 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); + + 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 new file mode 100644 index 000000000..1f63ae163 --- /dev/null +++ b/src/ShapeRecogn/CMakeLists.txt @@ -0,0 +1,63 @@ +# Copyright (C) 2012-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 +# + +ADD_DEFINITIONS(${HDF5_DEFINITIONS} ${MEDFILE_DEFINITIONS}) + +IF (NOT DEFINED MSVC) + ADD_DEFINITIONS(-Wsign-compare -Wconversion) +ENDIF() + +IF(MEDCOUPLING_ENABLE_PYTHON) + ADD_SUBDIRECTORY(Swig) +ENDIF(MEDCOUPLING_ENABLE_PYTHON) + +IF(MEDCOUPLING_BUILD_TESTS) + ADD_SUBDIRECTORY(Test) +ENDIF(MEDCOUPLING_BUILD_TESTS) + +INCLUDE_DIRECTORIES( + ${MEDFILE_INCLUDE_DIRS} + ${HDF5_INCLUDE_DIRS} + ${LAPACKE_INCLUDE_DIRS} + ${CMAKE_CURRENT_SOURCE_DIR}/../MEDLoader + ${CMAKE_CURRENT_SOURCE_DIR}/../MEDCoupling + ${CMAKE_CURRENT_SOURCE_DIR}/../INTERP_KERNEL + ${CMAKE_CURRENT_SOURCE_DIR}/../INTERP_KERNEL/Bases +) + +SET(shaperecogn_SOURCES + MathOps.cxx + Nodes.cxx + NodesBuilder.cxx + Areas.cxx + AreasBuilder.cxx + ShapeRecognMeshBuilder.cxx + ShapeRecognMesh.cxx +) + +ADD_LIBRARY(shaperecogn ${shaperecogn_SOURCES}) +SET_TARGET_PROPERTIES(shaperecogn PROPERTIES COMPILE_FLAGS "") +TARGET_LINK_LIBRARIES(shaperecogn medcouplingcpp medloader ${MEDFILE_C_LIBRARIES} ${HDF5_LIBRARIES} ${LAPACK_LIBRARIES}) +INSTALL(TARGETS shaperecogn EXPORT ${PROJECT_NAME}TargetGroup DESTINATION ${MEDCOUPLING_INSTALL_LIBS}) + +FILE(GLOB shaperecogn_HEADERS_HXX "${CMAKE_CURRENT_SOURCE_DIR}/*.hxx") +INSTALL(FILES ${shaperecogn_HEADERS_HXX} DESTINATION ${MEDCOUPLING_INSTALL_HEADERS}) + +# To allow usage as SWIG dependencies: +SET(shaperecogn_HEADERS_HXX PARENT_SCOPE) diff --git a/src/ShapeRecogn/MathOps.cxx b/src/ShapeRecogn/MathOps.cxx new file mode 100644 index 000000000..8f1705ace --- /dev/null +++ b/src/ShapeRecogn/MathOps.cxx @@ -0,0 +1,331 @@ +// 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(ldb, 0.0); + for (size_t i = 0; i < b.size(); ++i) + x[i] = b[i]; + 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::vector MathOps::lstsqRow( + std::vector &a, + const std::vector &b) +{ + int m = b.size(); + int n = 3; + int nrhs = 1; + int ldb = std::max(m, n); + int lds = std::min(m, n); + std::vector x(ldb, 0.0); + for (size_t i = 0; i < b.size(); ++i) + x[i] = b[i]; + 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_ROW_MAJOR, + m, n, nrhs, + a.data(), n, + x.data(), nrhs, + 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, (int)nbNodes, 1.0, + coordsCentered.data(), 3, + coordsCentered.data(), 3, + 0.0, covMatrix.data(), 3); + for (size_t i = 0; i < 9; ++i) + { + covMatrix[i] /= (double)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::computeVariance(std::vector values) +{ + size_t 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 / (double)n; +} + +std::array MathOps::computeBaseFromNormal(std::array normal) +{ + std::array n_normal = normalize(normal); + std::array s; + std::array u; + std::array v; + LAPACKE_dgesdd( + LAPACK_COL_MAJOR, + 'A', 1, 3, n_normal.data(), 1, + s.data(), + u.data(), 1, + v.data(), 3); + std::array v1 = {v[3], v[4], v[5]}; + std::array v2 = {v[6], v[7], v[8]}; + std::array u1 = normalize(v1); + std::array u2; + u2.fill(0.0); + double innerProd = dot(u1, v2); + for (size_t i = 0; i < 3; ++i) + u2[i] = v2[i] - innerProd * u1[i]; + u2 = normalize(u2); + double sign = dot(cross(u1, u2), normal); + if (sign < 0) + { + for (size_t i = 0; i < 3; ++i) + u1[i] *= -1; + } + return {u1[0], u1[1], u1[2], u2[0], u2[1], u2[2]}; +} diff --git a/src/ShapeRecogn/MathOps.hxx b/src/ShapeRecogn/MathOps.hxx new file mode 100644 index 000000000..e511845da --- /dev/null +++ b/src/ShapeRecogn/MathOps.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 __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::vector lstsqRow( + std::vector &a, + const std::vector &b); + 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 computeVariance(std::vector values); + static std::array computeBaseFromNormal(std::array normal); + }; +}; + +#endif //__MATHOPS_HXX__ diff --git a/src/ShapeRecogn/Nodes.cxx b/src/ShapeRecogn/Nodes.cxx new file mode 100644 index 000000000..615971ad4 --- /dev/null +++ b/src/ShapeRecogn/Nodes.cxx @@ -0,0 +1,125 @@ +#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::getK1() const +{ + return k1; +} + +const std::vector &Nodes::getK2() const +{ + return k2; +} + +const std::vector &Nodes::getPrimitiveType() const +{ + return primitives; +} + +const std::vector &Nodes::getNormals() const +{ + return normals; +} + +const std::vector &Nodes::getWeakDirections() const +{ + return weakDirections; +} + +const std::vector &Nodes::getMainDirections() const +{ + return mainDirections; +} + +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 fabs(k1[nodeId]) > fabs(k2[nodeId]) ? k1[nodeId] : k2[nodeId]; +} + +double Nodes::getAdimK1(mcIdType nodeId) const +{ + return adimK1[nodeId]; +} + +double Nodes::getAdimK2(mcIdType nodeId) const +{ + return adimK2[nodeId]; +} + +double Nodes::getAdimKdiff0(mcIdType nodeId) const +{ + return fabs(adimK1[nodeId]) > fabs(adimK2[nodeId]) ? adimK1[nodeId] : adimK2[nodeId]; +} + +const std::vector Nodes::getNeighbors(mcIdType nodeId) const +{ + 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::getMainDirection(mcIdType nodeId) const +{ + return {mainDirections[3 * nodeId], + mainDirections[3 * nodeId + 1], + mainDirections[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..1d47d3060 --- /dev/null +++ b/src/ShapeRecogn/Nodes.hxx @@ -0,0 +1,78 @@ +// 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 NodesBuilder; + Nodes(const MEDCouplingUMesh *mesh, + const DataArrayInt64 *neighbors, + const DataArrayInt64 *neighborsIdx); + ~Nodes(); + + mcIdType getNbNodes() const; + const std::vector &getK1() const; + const std::vector &getK2() const; + const std::vector &getPrimitiveType() const; + const std::vector &getNormals() const; + const std::vector &getWeakDirections() const; + const std::vector &getMainDirections() const; + + std::array getNormal(mcIdType nodeId) const; + double getK1(mcIdType nodeId) const; + double getK2(mcIdType nodeId) const; + double getKdiff0(mcIdType nodeId) const; + double getAdimK1(mcIdType nodeId) const; + double getAdimK2(mcIdType nodeId) const; + double getAdimKdiff0(mcIdType nodeId) const; + const std::vector getNeighbors(mcIdType nodeId) const; + PrimitiveType getPrimitiveType(mcIdType nodeId) const; + std::array getWeakDirection(mcIdType nodeId) const; + std::array getMainDirection(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 k1; + std::vector k2; + std::vector adimK1; + std::vector adimK2; + std::vector weakDirections; + std::vector mainDirections; + std::vector primitives; + }; +}; + +#endif //__NODES_HXX__ diff --git a/src/ShapeRecogn/NodesBuilder.cxx b/src/ShapeRecogn/NodesBuilder.cxx new file mode 100644 index 000000000..1302d3203 --- /dev/null +++ b/src/ShapeRecogn/NodesBuilder.cxx @@ -0,0 +1,311 @@ + +#include "NodesBuilder.hxx" +#include "Nodes.hxx" +#include "MathOps.hxx" +#include "ShapeRecongConstants.hxx" + +using namespace MEDCoupling; + +NodesBuilder::NodesBuilder( + const MEDCouplingUMesh *mesh) : mesh(mesh) +{ +} + +Nodes *NodesBuilder::build() +{ + DataArrayInt64 *neighbors; + DataArrayInt64 *neighborsIdx; + mesh->computeNeighborsOfNodes(neighbors, neighborsIdx); + nodes = new Nodes(mesh, neighbors, neighborsIdx); + computeNormals(); + computeCurvatures(); + return nodes; +} + +void NodesBuilder::computeNormals() +{ + mcIdType nbNodes = mesh->getNumberOfNodes(); + mcIdType nbCells = mesh->getNumberOfCells(); + std::array cellNormal; + cellNormal.fill(0.0); + // product of normal AreaByCell + std::vector prodNormalAreaByCell(3 * nbCells, 0.0); + for (int cellId = 0; cellId < nbCells; cellId++) + { + std::vector nodeIds; + mesh->getNodeIdsOfCell(cellId, nodeIds); + computeCellNormal(nodeIds, cellNormal); + prodNormalAreaByCell[3 * cellId + 0] = cellNormal[0] / 2.0; + prodNormalAreaByCell[3 * cellId + 1] = cellNormal[1] / 2.0; + prodNormalAreaByCell[3 * cellId + 2] = cellNormal[2] / 2.0; + } + // + nodes->normals.resize(3 * nbNodes, 1.0); + DataArrayInt64 *revNodal = DataArrayInt64::New(); + DataArrayInt64 *revNodalIdx = DataArrayInt64::New(); + mesh->getReverseNodalConnectivity(revNodal, revNodalIdx); + for (size_t nodeId = 0; nodeId < (size_t)nbNodes; nodeId++) + { + int nbCells = revNodalIdx->getIJ(nodeId + 1, 0) - + revNodalIdx->getIJ(nodeId, 0); + std::vector cellIds(nbCells, 0); + int start = revNodalIdx->getIJ(nodeId, 0); + for (size_t i = 0; i < cellIds.size(); ++i) + cellIds[i] = revNodal->getIJ(start + i, 0); + double normal = 0.0; + for (size_t i = 0; i < 3; i++) + { + nodes->normals[3 * nodeId + i] = 0.0; + for (mcIdType j = 0; j < nbCells; j++) + { + nodes->normals[3 * nodeId + i] += + prodNormalAreaByCell[3 * cellIds[j] + i]; + } + nodes->normals[3 * nodeId + i] /= (double)nbCells; + normal += pow(nodes->normals[3 * nodeId + i], 2); + } + for (size_t i = 0; i < 3; i++) + nodes->normals[3 * nodeId + i] /= sqrt(normal); + } + revNodal->decrRef(); + revNodalIdx->decrRef(); +} + +void NodesBuilder::computeCurvatures(double tol) +{ + mcIdType nbNodes = mesh->getNumberOfNodes(); + nodes->k1.resize(nbNodes); + nodes->k2.resize(nbNodes); + nodes->adimK1.resize(nbNodes); + nodes->adimK2.resize(nbNodes); + nodes->primitives.resize(nbNodes); + nodes->weakDirections.resize(3 * nbNodes); + nodes->mainDirections.resize(3 * nbNodes); + for (mcIdType nodeId = 0; nodeId < nbNodes; nodeId++) + computeCurvatures(nodeId, tol); +} + +void NodesBuilder::computeCurvatures(mcIdType nodeId, double tol) +{ + std::array normal = nodes->getNormal(nodeId); + std::vector neighborIds = nodes->getNeighbors(nodeId); + double theta0 = 0.0; + double k1 = 0.0; + double k2 = 0.0; + double adimK1 = 0.0; + double adimK2 = 0.0; + PrimitiveType primitive = PrimitiveType::Unknown; + std::array mainDirection{0.0, 0.0, 0.0}; + std::array weakDirection{0.0, 0.0, 0.0}; + if (neighborIds.size() > 0) + { + std::vector discreteCurvatures = computeDiscreteCurvatures(nodeId, neighborIds); + std::vector tangents = computeTangentDirections(nodeId, neighborIds); + mcIdType maxCurvatureId = std::distance( + discreteCurvatures.begin(), + std::max_element(discreteCurvatures.begin(), discreteCurvatures.end())); + std::array e1{ + tangents[3 * maxCurvatureId], + tangents[3 * maxCurvatureId + 1], + tangents[3 * maxCurvatureId + 2]}; + std::array e2 = MathOps::normalize(MathOps::cross(e1, normal)); + std::vector coeffs = computeNormalCurvatureCoefficients( + discreteCurvatures, tangents, normal, e1); + double a = coeffs[0], b = coeffs[1], c = coeffs[2]; + double h = (a + c) / 2.0; + double kg = a * c - pow(b, 2) / 4.0; + k1 = h + sqrt(fabs(pow(h, 2) - kg)); + k2 = h - sqrt(fabs(pow(h, 2) - kg)); + if (fabs(k1 - k2) >= tol) + theta0 = 0.5 * asin(b / (k1 - k2)); + std::array direction1{0.0, 0.0, 0.0}; + std::array direction2{0.0, 0.0, 0.0}; + for (size_t i = 0; i < 3; ++i) + { + direction1[i] = cos(theta0) * e1[i] + sin(theta0) * e2[i]; + direction2[i] = -sin(theta0) * e1[i] + cos(theta0) * e2[i]; + } + double averageDistance = computeAverageDistance(nodeId, neighborIds); + double adimK1 = k1 * averageDistance; + double adimK2 = k2 * averageDistance; + double adimKdiff0, adimKis0; + if (fabs(k1) > fabs(k2)) + { + adimKdiff0 = adimK1; + adimKis0 = adimK2; + mainDirection = direction1; + weakDirection = direction2; + } + else + { + adimKdiff0 = adimK2; + adimKis0 = adimK1; + mainDirection = direction2; + weakDirection = direction1; + } + primitive = findPrimitiveType(adimK1, adimK2, adimKdiff0, adimKis0); + } + // Populate nodes + nodes->k1[nodeId] = k1; + nodes->k2[nodeId] = k2; + nodes->adimK1[nodeId] = adimK1; + nodes->adimK2[nodeId] = adimK2; + for (size_t i = 0; i < 3; ++i) + { + nodes->weakDirections[3 * nodeId + i] = weakDirection[i]; + nodes->mainDirections[3 * nodeId + i] = mainDirection[i]; + } + nodes->primitives[nodeId] = primitive; +} + +PrimitiveType NodesBuilder::findPrimitiveType(double k1, double k2, double kdiff0, double kis0) const +{ + if ((fabs(k1 - k2) < EPSILON_PRIMITIVE) && (fabs((k1 + k2) / 2) < EPSILON_PRIMITIVE)) + return PrimitiveType::Plane; + else if ((fabs(k1 - k2) < EPSILON_PRIMITIVE) && (fabs((k1 + k2) / 2) > EPSILON_PRIMITIVE)) + return PrimitiveType::Sphere; + else if ((fabs(kdiff0) > EPSILON_PRIMITIVE) && (fabs(kis0) < EPSILON_PRIMITIVE)) + return PrimitiveType::Cylinder; + else if ((fabs(k1) > EPSILON_PRIMITIVE) && (fabs(k2) > EPSILON_PRIMITIVE)) + return PrimitiveType::Torus; + else + return PrimitiveType::Unknown; +} + +PrimitiveType NodesBuilder::findPrimitiveType2(double k1, double k2, double kdiff0, double kis0) const +{ + double epsilon2 = pow(EPSILON_PRIMITIVE, 2); + double diffCurvature = fabs(k1 - k2); + double gaussianCurvature = k1 * k2; + double meanCurvature = (k1 + k2) / 2.0; + if (fabs(k1) < EPSILON_PRIMITIVE && + fabs(k2) < EPSILON_PRIMITIVE && + gaussianCurvature < epsilon2 && + meanCurvature < EPSILON_PRIMITIVE) + return PrimitiveType::Plane; + else if (diffCurvature < EPSILON_PRIMITIVE && k1 > EPSILON_PRIMITIVE && k2 > EPSILON_PRIMITIVE) + return PrimitiveType::Sphere; + else if ( + (fabs(k1) > EPSILON_PRIMITIVE && fabs(k2) < EPSILON_PRIMITIVE) || + (fabs(k1) < EPSILON_PRIMITIVE && fabs(k2) > EPSILON_PRIMITIVE)) + return PrimitiveType::Cylinder; + else if ( + std::signbit(k1) != std::signbit(k2) || + (fabs(k1) < EPSILON_PRIMITIVE && fabs(k2) < EPSILON_PRIMITIVE)) + return PrimitiveType::Torus; + else + return PrimitiveType::Unknown; +} + +std::vector NodesBuilder::computeNormalCurvatureCoefficients( + const std::vector &discreteCurvatures, + const std::vector &tangents, + const std::array &normal, + const std::array &e1) const +{ + size_t nbNeighbors = discreteCurvatures.size(); + std::vector a(3 * nbNeighbors, 0.0); + for (size_t i = 0; i < nbNeighbors; ++i) + { + std::array tangent{tangents[3 * i], tangents[3 * i + 1], tangents[3 * i + 2]}; + double theta = MathOps::computeOrientedAngle(normal, tangent, e1); + double cos_theta = cos(theta); + double sin_theta = sin(theta); + a[i] = pow(cos_theta, 2); + a[nbNeighbors + i] = cos_theta * sin_theta; + a[2 * nbNeighbors + i] = pow(sin_theta, 2); + } + return MathOps::lstsq(a, discreteCurvatures); +} + +void NodesBuilder::computeCellNormal( + const std::vector &nodeIds, + std::array &cellNormal) const +{ + std::vector point1; + std::vector point2; + std::vector point3; + mesh->getCoordinatesOfNode(nodeIds[0], point1); + mesh->getCoordinatesOfNode(nodeIds[1], point2); + mesh->getCoordinatesOfNode(nodeIds[2], point3); + std::array a; + a.fill(3); + std::array b; + b.fill(3); + for (int i = 0; i < 3; i++) + { + a[i] = point2[i] - point1[i]; + b[i] = point3[i] - point1[i]; + } + cellNormal[0] = a[1] * b[2] - a[2] * b[1]; + cellNormal[1] = a[2] * b[0] - a[0] * b[2]; + cellNormal[2] = a[0] * b[1] - a[1] * b[0]; +} + +double NodesBuilder::computeAverageDistance(mcIdType nodeId, const std::vector &neighborIds) const +{ + double distance = 0.0; + std::array + nodeCoords = nodes->getCoordinates(nodeId); + for (size_t i = 0; i < neighborIds.size(); ++i) + { + std::array neighborCoords = nodes->getCoordinates(neighborIds[i]); + double distanceToNeighbor = 0.0; + for (size_t j = 0; j < 3; ++j) + distanceToNeighbor += pow(neighborCoords[j] - nodeCoords[j], 2); + distance += sqrt(distanceToNeighbor); + } + return distance / (double)neighborIds.size(); +} + +std::vector NodesBuilder::computeDiscreteCurvatures( + mcIdType nodeId, + const std::vector &neighborIds) const +{ + std::vector discreteCurvatures(neighborIds.size(), 0.0); + for (size_t i = 0; i < neighborIds.size(); ++i) + discreteCurvatures[i] = computeDiscreteCurvature(nodeId, neighborIds[i]); + return discreteCurvatures; +} + +double NodesBuilder::computeDiscreteCurvature( + mcIdType nodeId, + mcIdType neighborId) const +{ + double curvature = 0.0; + double n = 0.0; + for (size_t i = 0; i < 3; i++) + { + double ni = nodes->coords->getIJ(neighborId, i) - nodes->coords->getIJ(nodeId, i); + curvature += ni * (nodes->normals[3 * neighborId + i] - nodes->normals[3 * nodeId + i]); + n += ni * ni; + } + return curvature / n; +} + +std::vector NodesBuilder::computeTangentDirections( + mcIdType nodeId, + const std::vector &neighborIds) const +{ + size_t nbNeighbors = neighborIds.size(); + std::vector tangent(3 * nbNeighbors, 0.0); + for (size_t i = 0; i < nbNeighbors; ++i) + { + mcIdType neighborId = neighborIds[i]; + double s = 0.0; + for (size_t j = 0; j < 3; j++) + { + tangent[3 * i + j] = nodes->coords->getIJ(neighborId, j) - nodes->coords->getIJ(nodeId, j); + s += tangent[3 * i + j] * nodes->normals[3 * nodeId + j]; + } + double n = 0.0; + for (size_t j = 0; j < 3; j++) + { + tangent[3 * i + j] -= s * nodes->normals[3 * nodeId + j]; + n += tangent[3 * i + j] * tangent[3 * i + j]; + } + for (size_t j = 0; j < 3; j++) + tangent[3 * i + j] /= sqrt(n); + } + return tangent; +} diff --git a/src/ShapeRecogn/NodesBuilder.hxx b/src/ShapeRecogn/NodesBuilder.hxx new file mode 100644 index 000000000..190118064 --- /dev/null +++ b/src/ShapeRecogn/NodesBuilder.hxx @@ -0,0 +1,62 @@ +// Copyright (C) 2007-2024 CEA, EDF +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#ifndef __NODECURVATURECALCULATOR_HXX__ +#define __NODECURVATURECALCULATOR_HXX__ + +#include +#include +#include +#include "MEDCouplingUMesh.hxx" +#include "PrimitiveType.hxx" + +namespace MEDCoupling +{ + class Nodes; + + class NodesBuilder + { + public: + NodesBuilder(const MEDCouplingUMesh *); + + Nodes *build(); + + 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; + PrimitiveType findPrimitiveType2(double k1, double k2, double kdiff0, double kis0) const; + std::vector computeNormalCurvatureCoefficients( + const std::vector &discreteCurvatures, + const std::vector &tangents, + const std::array &normal, + const std::array &e1) const; + void computeCellNormal(const std::vector &nodeIds, std::array &cellNormal) const; + double computeAverageDistance(mcIdType nodeId, const std::vector &neighborIds) const; + std::vector computeDiscreteCurvatures(mcIdType nodeId, const std::vector &neighborIds) const; + double computeDiscreteCurvature(mcIdType nodeId, mcIdType neighborId) const; + std::vector computeTangentDirections(mcIdType nodeId, const std::vector &neighborIds) const; + + const MEDCouplingUMesh *mesh; + Nodes *nodes; + }; +}; + +#endif //__NODECURVATURECALCULATOR_HXX__ \ No newline at end of file diff --git a/src/ShapeRecogn/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/README.md b/src/ShapeRecogn/README.md new file mode 100644 index 000000000..61b9eb4d0 --- /dev/null +++ b/src/ShapeRecogn/README.md @@ -0,0 +1,106 @@ +# ShapeRecogn +A tool leveraging the MEDCoupling library for recognizing canonical shapes in 3D mesh files. + +## Table Of Contents + +1. [Quickstart](#quickstart) +2. [Installation](#installation) +3. [Description](#description) + +## Quickstart + +### With output file + +```python +>> import ShapeRecogn as sr +>> +>> shape_recogn_builder = sr.ShapeRecognMeshBuilder("resources/ShapeRecognCone.med") +>> shape_recogn = shape_recogn_builder.recognize() +>> shape_recogn.save("ShapeRecognCone_areas.med") +``` + +### Without output file + +```python +>> import ShapeRecogn as sr +>> +>> shape_recogn_builder = sr.ShapeRecognMeshBuilder("resources/ShapeRecognCone.med") +>> shape_recogn = shape_recogn_builder.recognize() +>> radius_field = shape_recogn.getRadius() +>> radius_field +MEDCouplingFieldDouble C++ instance at 0x55f3418c6700. Name : "Radius (Area)". +Nature of field : NoNature. +P1 spatial discretization. + +Mesh info : MEDCouplingUMesh C++ instance at 0x55f34171a2d0. Name : "Mesh_1". Mesh dimension : 2. Space dimension : 3. + +Array info : DataArrayDouble C++ instance at 0x55f3418b40e0. Number of tuples : 957. Number of components : 1. +[9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, ... ] +``` + +## Installation + +### Prerequisites + +**LAPACKE Library**: ShapeRecogn requires the [LAPACKE](https://www.netlib.org/lapack/lapacke.html) library for numerical computations. Ensure it is installed on your system. + +### CMAKE configuration + +Run CMake to configure the build. Make sure to enable ShapeRecogn by setting the `-DMEDCOUPLING_ENABLE_SHAPERECOGN=ON` option. + +## Description + +ShapeRecogn is a tool that leverages the MEDCoupling library to recognize +canonical shapes from 3D meshes using triangle cells. + +The tool is based on the thesis work of Roseline Bénière *[Reconstruction d’un modèle B-Rep à partir d’un maillage 3D](https://theses.fr/2012MON20034)*, and it recognizes five canonical shapes: + - Plane (0) + - Sphere (1) + - Cylinder (2) + - Cone (3) + - Torus (4) + - Unknown (5) (When the algorithm failed the recognition) + +The recognized shapes are divided into areas within the mesh. +The tool also generates an output file with detailed fields describing the areas and their properties. This makes it easier to analyze or further process the mesh data. + - Area Id: To distinguish the areas + - Primitive Type (Area): One of the canonical shape with the id describe above + - Normal (Area): + * Normal of a plane area + - Minor Radius (Area) + * Minor radius of a torus area + - Radius (Area) + * Radius of a sphere area + * Radius of a cylinder area + * Radius of the base of a cone area + * Major radius of a torus area + - Angle (Area) + * Angle between the central axis and the slant of a cone area + - Center (Area) + * Center of a sphere area + * Center of a torus area + - Axis (Area) + * Central axis of a cylinder area + * Central axis of a cone area + - Apex (Area) + * Apex of a cone area + +Some intermediate computation values concerning the nodes are also available: + - K1 (Node) and K2 (Node): the *[principal curvatures](https://en.wikipedia.org/wiki/Principal_curvature)* + - Primitive Type (Node) : One of the canonical shape with the id describe above. The primitive type is deduced usint K1 and K2 values. + - Normal (Node): Normal of the nodes using neighbor nodes + +Each field can be retrieved as a `MEDCouplingDoubleField` using the following methods of the `ShapeRecognMesh` class: + - `getAreaId()` + - `getAreaPrimitiveType()` + - `getAreaNormal()` + - `getMinorRadius()` + - `getRadius()` + - `getAngle()` + - `getCenter()` + - `getAxis()` + - `getApex()` + - `getNodeK1()` + - `getNodeK2()` + - `getNodePrimitiveType()` + - `getNodeNormal()` diff --git a/src/ShapeRecogn/ShapeRecognMesh.cxx b/src/ShapeRecogn/ShapeRecognMesh.cxx new file mode 100644 index 000000000..4c8ee5f9f --- /dev/null +++ b/src/ShapeRecogn/ShapeRecognMesh.cxx @@ -0,0 +1,155 @@ +#include "ShapeRecognMesh.hxx" + +#include "MEDLoader.hxx" + +using namespace MEDCoupling; + +ShapeRecognMesh::ShapeRecognMesh() + : nodeK1(0), nodeK2(0), nodePrimitiveType(0), + nodeNormal(0), areaId(0), areaPrimitiveType(0), + areaNormal(0), minorRadius(0), radius(0), + angle(0), center(0), axis(0), apex(0) +{ +} +ShapeRecognMesh::~ShapeRecognMesh() +{ + nodeK1->decrRef(); + nodeK2->decrRef(); + nodePrimitiveType->decrRef(); + nodeNormal->decrRef(); + areaId->decrRef(); + areaPrimitiveType->decrRef(); + areaNormal->decrRef(); + minorRadius->decrRef(); + radius->decrRef(); + angle->decrRef(); + center->decrRef(); + axis->decrRef(); + apex->decrRef(); +} + +std::size_t ShapeRecognMesh::getHeapMemorySizeWithoutChildren() const +{ + return 0; +} + +std::vector ShapeRecognMesh::getDirectChildrenWithNull() const +{ + std::vector ret; + ret.push_back(nodeK1); + ret.push_back(nodeK2); + ret.push_back(nodePrimitiveType); + ret.push_back(nodeNormal); + ret.push_back(areaId); + ret.push_back(areaPrimitiveType); + ret.push_back(areaNormal); + ret.push_back(minorRadius); + ret.push_back(radius); + ret.push_back(angle); + ret.push_back(center); + ret.push_back(axis); + ret.push_back(apex); + return ret; +} + +ShapeRecognMesh *ShapeRecognMesh::New() +{ + return new ShapeRecognMesh; +} + +void ShapeRecognMesh::save(const std::string &outputFile, bool writeFromScratch) const +{ + // Nodes + // - k1 + WriteField(outputFile, nodeK1, writeFromScratch); + // - k2 + WriteField(outputFile, nodeK2, false); + // - primitive types + WriteField(outputFile, nodePrimitiveType, false); + // - Normal + WriteField(outputFile, nodeNormal, false); + // Areas + // - Area Id + WriteField(outputFile, areaId, false); + // - Primitive Types + WriteField(outputFile, areaPrimitiveType, false); + // - Normal + WriteField(outputFile, areaNormal, false); + // - Minor Radius + WriteField(outputFile, minorRadius, false); + // - Radius + WriteField(outputFile, radius, false); + // - Angle + WriteField(outputFile, angle, false); + // - Center + WriteField(outputFile, center, false); + // - Axis + WriteField(outputFile, axis, false); + // - Apex + WriteField(outputFile, apex, false); +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getNodeK1() const +{ + return nodeK1; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getNodeK2() const +{ + return nodeK2; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getNodePrimitiveType() const +{ + return nodePrimitiveType; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getNodeNormal() const +{ + return nodeNormal; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getAreaId() const +{ + return areaId; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getAreaPrimitiveType() const +{ + return areaPrimitiveType; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getAreaNormal() const +{ + return areaNormal; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getMinorRadius() const +{ + return minorRadius; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getRadius() const +{ + return radius; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getAngle() const +{ + return angle; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getCenter() const +{ + return center; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getAxis() const +{ + return axis; +} + +MEDCouplingFieldDouble *ShapeRecognMesh::getApex() const +{ + return apex; +} diff --git a/src/ShapeRecogn/ShapeRecognMesh.hxx b/src/ShapeRecogn/ShapeRecognMesh.hxx new file mode 100644 index 000000000..02983e34a --- /dev/null +++ b/src/ShapeRecogn/ShapeRecognMesh.hxx @@ -0,0 +1,80 @@ +// 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" +#include "MEDCouplingFieldDouble.hxx" +#include "MEDCouplingRefCountObject.hxx" + +namespace MEDCoupling +{ + class ShapeRecognMesh : public RefCountObject + { + friend class ShapeRecognMeshBuilder; + + public: + static ShapeRecognMesh *New(); + std::size_t getHeapMemorySizeWithoutChildren() const; + std::vector getDirectChildrenWithNull() const; + + void save(const std::string &outputFile, bool writeFromScratch = true) const; + + // Node properties + MEDCoupling::MEDCouplingFieldDouble *getNodeK1() const; + MEDCoupling::MEDCouplingFieldDouble *getNodeK2() const; + MEDCoupling::MEDCouplingFieldDouble *getNodePrimitiveType() const; + MEDCoupling::MEDCouplingFieldDouble *getNodeNormal() const; + + // Area properties + MEDCoupling::MEDCouplingFieldDouble *getAreaId() const; + MEDCoupling::MEDCouplingFieldDouble *getAreaPrimitiveType() const; + MEDCoupling::MEDCouplingFieldDouble *getAreaNormal() const; + MEDCoupling::MEDCouplingFieldDouble *getMinorRadius() const; + MEDCoupling::MEDCouplingFieldDouble *getRadius() const; + MEDCoupling::MEDCouplingFieldDouble *getAngle() const; + MEDCoupling::MEDCouplingFieldDouble *getCenter() const; + MEDCoupling::MEDCouplingFieldDouble *getAxis() const; + MEDCoupling::MEDCouplingFieldDouble *getApex() const; + + protected: + ShapeRecognMesh(); + ~ShapeRecognMesh(); + + private: + MEDCoupling::MEDCouplingFieldDouble *nodeK1; + MEDCoupling::MEDCouplingFieldDouble *nodeK2; + MEDCoupling::MEDCouplingFieldDouble *nodePrimitiveType; + MEDCoupling::MEDCouplingFieldDouble *nodeNormal; + MEDCoupling::MEDCouplingFieldDouble *areaId; + MEDCoupling::MEDCouplingFieldDouble *areaPrimitiveType; + MEDCoupling::MEDCouplingFieldDouble *areaNormal; + MEDCoupling::MEDCouplingFieldDouble *minorRadius; + MEDCoupling::MEDCouplingFieldDouble *radius; + MEDCoupling::MEDCouplingFieldDouble *angle; + MEDCoupling::MEDCouplingFieldDouble *center; + MEDCoupling::MEDCouplingFieldDouble *axis; + MEDCoupling::MEDCouplingFieldDouble *apex; + }; +}; + +#endif // __SHAPERECOGNMESH_HXX__ diff --git a/src/ShapeRecogn/ShapeRecognMeshBuilder.cxx b/src/ShapeRecogn/ShapeRecognMeshBuilder.cxx new file mode 100644 index 000000000..f0a2b0581 --- /dev/null +++ b/src/ShapeRecogn/ShapeRecognMeshBuilder.cxx @@ -0,0 +1,250 @@ +#include "ShapeRecognMeshBuilder.hxx" + +#include "NodesBuilder.hxx" +#include "AreasBuilder.hxx" +#include "MEDLoader.hxx" +#include "ShapeRecognMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" + +using namespace MEDCoupling; + +ShapeRecognMeshBuilder::ShapeRecognMeshBuilder(const std::string &fileName, int meshDimRelToMax) +{ + mesh = ReadUMeshFromFile(fileName, meshDimRelToMax); + if (mesh->getMeshDimension() != 2) + throw INTERP_KERNEL::Exception("Expect a mesh with a dimension equal to 2"); + if (mesh->getNumberOfCellsWithType(INTERP_KERNEL::NORM_TRI3) != mesh->getNumberOfCells()) + throw INTERP_KERNEL::Exception("Expect a mesh containing exclusively triangular cells"); +} + +ShapeRecognMeshBuilder::~ShapeRecognMeshBuilder() +{ + if (areas != nullptr) + delete areas; + if (nodes != nullptr) + delete nodes; + mesh->decrRef(); +} + +ShapeRecognMesh *ShapeRecognMeshBuilder::recognize() +{ + mesh->incrRef(); + NodesBuilder nodesBuilder(mesh); + nodes = nodesBuilder.build(); + AreasBuilder areasBuilder(nodes); + areasBuilder.build(); + areas = areasBuilder.getAreas(); + MCAuto recognMesh = ShapeRecognMesh::New(); + recognMesh->nodeK1 = buildNodeK1(); + recognMesh->nodeK2 = buildNodeK2(); + recognMesh->nodePrimitiveType = buildNodePrimitiveType(); + recognMesh->nodeNormal = buildNodeNormal(); + recognMesh->areaId = buildAreaId(); + recognMesh->areaPrimitiveType = buildAreaPrimitiveType(); + recognMesh->areaNormal = buildAreaNormal(); + recognMesh->minorRadius = buildMinorRadius(); + recognMesh->radius = buildRadius(); + recognMesh->angle = buildAngle(); + recognMesh->center = buildCenter(); + recognMesh->axis = buildAxis(); + recognMesh->apex = buildApex(); + return recognMesh.retn(); +} + +const Nodes *ShapeRecognMeshBuilder::getNodes() const +{ + return nodes; +} + +const Areas *ShapeRecognMeshBuilder::getAreas() const +{ + return areas; +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildNodeK1() const +{ + if (nodes == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + return buildField("K1 (Node)", 1, nodes->getK1()); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildNodeK2() const +{ + if (nodes == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + return buildField("K2 (Node)", 1, nodes->getK2()); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildNodePrimitiveType() const +{ + if (nodes == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + return buildField("Primitive Type (Node)", 1, nodes->getPrimitiveType()); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildNodeNormal() const +{ + if (nodes == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + return buildField("Normal (Node)", 3, nodes->getNormals()); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildAreaId() const +{ + if (areas == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + return buildField("Area Id", 1, areas->getAreaIdByNodes()); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildAreaPrimitiveType() const +{ + if (areas == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + double *values = buildAreaArray([](Areas *areas, mcIdType areaId) -> double + { return (double)areas->getPrimitiveType(areaId); }); + return buildField("Primitive Type (Area)", 1, values); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildAreaNormal() const +{ + if (areas == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + double *values = buildArea3DArray([](Areas *areas, mcIdType areaId) -> const std::array & + { return areas->getNormal(areaId); }); + return buildField("Normal (Area)", 3, values); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildMinorRadius() const +{ + if (areas == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + double *values = buildAreaArray([](Areas *areas, mcIdType areaId) -> double + { return areas->getMinorRadius(areaId); }); + return buildField("Minor Radius (Area)", 1, values); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildRadius() const +{ + if (areas == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + double *values = buildAreaArray([](Areas *areas, mcIdType areaId) -> double + { return areas->getRadius(areaId); }); + return buildField("Radius (Area)", 1, values); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildAngle() const +{ + if (areas == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + double *values = buildAreaArray([](Areas *areas, mcIdType areaId) -> double + { return areas->getAngle(areaId); }); + return buildField("Angle (Area)", 1, values); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildCenter() const +{ + if (areas == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + double *values = buildArea3DArray([](Areas *areas, mcIdType areaId) -> const std::array & + { return areas->getCenter(areaId); }); + return buildField("Center (Area)", 3, values); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildAxis() const +{ + if (areas == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + double *values = buildArea3DArray([](Areas *areas, mcIdType areaId) -> const std::array & + { return areas->getAxis(areaId); }); + return buildField("Axis (Area)", 3, values); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildApex() const +{ + if (areas == nullptr) + throw INTERP_KERNEL::Exception("recognize must be called before building any fields"); + double *values = buildArea3DArray([](Areas *areas, mcIdType areaId) -> const std::array & + { return areas->getApex(areaId); }); + return buildField("Apex (Area)", 3, values); +} + +template +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildField( + const std::string &name, + size_t nbOfCompo, + const std::vector &values) const +{ + DataArrayDouble *data = DataArrayDouble::New(); + data->setName(name); + data->alloc(nodes->getNbNodes(), nbOfCompo); + std::copy(values.begin(), values.end(), data->getPointer()); + data->declareAsNew(); + return buildField(name, nbOfCompo, data); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildField( + const std::string &name, + size_t nbOfCompo, + double *values) const +{ + DataArrayDouble *data = DataArrayDouble::New(); + data->setName(name); + data->useArray( + values, + true, + MEDCoupling::DeallocType::CPP_DEALLOC, + nodes->getNbNodes(), + nbOfCompo); + return buildField(name, nbOfCompo, data); +} + +MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildField( + const std::string &name, + size_t nbOfCompo, + DataArrayDouble *data) const +{ + MEDCouplingFieldDouble *field = MEDCouplingFieldDouble::New(ON_NODES); + field->setName(name); + field->setMesh(mesh); + if (nbOfCompo == 3) + { + data->setInfoOnComponent(0, "X"); + data->setInfoOnComponent(1, "Y"); + data->setInfoOnComponent(2, "Z"); + } + field->setArray(data); + data->decrRef(); + return field; +} + +double *ShapeRecognMeshBuilder::buildArea3DArray( + const std::array &(*areaFunc)(Areas *, mcIdType)) const +{ + double *values = new double[3 * nodes->getNbNodes()]; + const std::vector &areaIdByNodes = areas->getAreaIdByNodes(); + for (size_t nodeId = 0; nodeId < areaIdByNodes.size(); ++nodeId) + { + mcIdType areaId = areaIdByNodes[nodeId]; + if (areaId != -1) + { + const std::array &areaValues = areaFunc(areas, areaId); + values[3 * nodeId] = areaValues[0]; + values[3 * nodeId + 1] = areaValues[1]; + values[3 * nodeId + 2] = areaValues[2]; + } + } + return values; +} + +double *ShapeRecognMeshBuilder::buildAreaArray(double (*areaFunc)(Areas *, mcIdType)) const +{ + const std::vector &areaIdByNodes = areas->getAreaIdByNodes(); + double *values = new double[nodes->getNbNodes()]; + for (size_t nodeId = 0; nodeId < areaIdByNodes.size(); ++nodeId) + { + mcIdType areaId = areaIdByNodes[nodeId]; + if (areaId != -1) + values[nodeId] = areaFunc(areas, areaId); + } + return values; +} diff --git a/src/ShapeRecogn/ShapeRecognMeshBuilder.hxx b/src/ShapeRecogn/ShapeRecognMeshBuilder.hxx new file mode 100644 index 000000000..eee7c1507 --- /dev/null +++ b/src/ShapeRecogn/ShapeRecognMeshBuilder.hxx @@ -0,0 +1,85 @@ +// 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 __SHAPERECOGNMESHBUILDER_HXX__ +#define __SHAPERECOGNMESHBUILDER_HXX__ + +#include + +#include "MEDCouplingUMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" + +namespace MEDCoupling +{ + class Nodes; + class Areas; + class ShapeRecognMesh; + + class ShapeRecognMeshBuilder + { + public: + ShapeRecognMeshBuilder(const std::string &fileName, int meshDimRelToMax = 0); + ~ShapeRecognMeshBuilder(); + + const Nodes *getNodes() const; + const Areas *getAreas() const; + + ShapeRecognMesh *recognize(); + + private: + // Node properties + MEDCoupling::MEDCouplingFieldDouble *buildNodeK1() const; + MEDCoupling::MEDCouplingFieldDouble *buildNodeK2() const; + MEDCoupling::MEDCouplingFieldDouble *buildNodePrimitiveType() const; + MEDCoupling::MEDCouplingFieldDouble *buildNodeNormal() const; + + // Area properties + MEDCoupling::MEDCouplingFieldDouble *buildAreaId() const; + MEDCoupling::MEDCouplingFieldDouble *buildAreaPrimitiveType() const; + MEDCoupling::MEDCouplingFieldDouble *buildAreaNormal() const; + MEDCoupling::MEDCouplingFieldDouble *buildMinorRadius() const; + MEDCoupling::MEDCouplingFieldDouble *buildRadius() const; + MEDCoupling::MEDCouplingFieldDouble *buildAngle() const; + MEDCoupling::MEDCouplingFieldDouble *buildCenter() const; + MEDCoupling::MEDCouplingFieldDouble *buildAxis() const; + MEDCoupling::MEDCouplingFieldDouble *buildApex() const; + + template + MEDCouplingFieldDouble *buildField( + const std::string &name, + size_t nbOfCompo, + const std::vector &values) const; + MEDCouplingFieldDouble *buildField( + const std::string &name, + size_t nbOfCompo, + double *values) const; + MEDCouplingFieldDouble *buildField( + const std::string &name, + size_t nbOfCompo, + DataArrayDouble *values) const; + double *buildArea3DArray(const std::array &(*areaFunc)(Areas *, mcIdType)) const; + double *buildAreaArray(double (*areaFunc)(Areas *, mcIdType)) const; + + const MEDCouplingUMesh *mesh; + Nodes *nodes = nullptr; + Areas *areas = nullptr; + }; +}; + +#endif // __SHAPERECOGNMESHBUILDER_HXX__ diff --git a/src/ShapeRecogn/ShapeRecongConstants.hxx b/src/ShapeRecogn/ShapeRecongConstants.hxx new file mode 100644 index 000000000..801df4dfa --- /dev/null +++ b/src/ShapeRecogn/ShapeRecongConstants.hxx @@ -0,0 +1,43 @@ +// 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_PRIMITIVE = 0.005; + // Areas + // - Match + constexpr double TOL_MATCH_CYLINDER = 0.05; + constexpr double TOL_MATCH_SPHERE = 0.05; + // - Relative distance + constexpr double DELTA_PLANE = 0.05; + constexpr double DELTA_SPHERE = 0.05; + constexpr double DELTA_CYLINDER = 0.05; + constexpr double DELTA_CONE = 0.05; + // - Invalid Zones + constexpr double THETA_MAX_CYLINDER = 0.05; + // - Thresholds + constexpr int THRESHOLD_MIN_NB_NODES = 5; + constexpr int THRESHOLD_MAX_NB_AREAS = 500; +}; + +#endif //__SHAPERECOGNCONSTANTS_HXX__ \ No newline at end of file diff --git a/src/ShapeRecogn/Swig/CMakeLists.txt b/src/ShapeRecogn/Swig/CMakeLists.txt new file mode 100644 index 000000000..8e839d156 --- /dev/null +++ b/src/ShapeRecogn/Swig/CMakeLists.txt @@ -0,0 +1,83 @@ +# Copyright (C) 2012-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 +# + +FIND_PACKAGE(SWIG REQUIRED) +INCLUDE(${SWIG_USE_FILE}) + +ADD_DEFINITIONS(${PYTHON_DEFINITIONS}) + +SET_SOURCE_FILES_PROPERTIES(ShapeRecogn.i PROPERTIES CPLUSPLUS ON) +IF ("${PYTHON_VERSION_MAJOR}" STREQUAL "3") + SET_SOURCE_FILES_PROPERTIES(ShapeRecogn.i PROPERTIES SWIG_FLAGS "-py3") +ELSE() + SET_SOURCE_FILES_PROPERTIES(ShapeRecogn.i PROPERTIES SWIG_DEFINITIONS "-shadow") +ENDIF() + +SET(SWIG_MODULE_ShapeRecogn_EXTRA_FLAGS "") +IF(MEDCOUPLING_USE_64BIT_IDS) + STRING(APPEND SWIG_MODULE_ShapeRecogn_EXTRA_FLAGS "-DMEDCOUPLING_USE_64BIT_IDS") +ENDIF(MEDCOUPLING_USE_64BIT_IDS) + +SET (ShapeRecogn_SWIG_DPYS_FILES + ShapeRecogn.i +) + +INCLUDE_DIRECTORIES( + ${PYTHON_INCLUDE_DIRS} + ${MEDFILE_INCLUDE_DIRS} + ${HDF5_INCLUDE_DIRS} + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_SOURCE_DIR}/.. + ${CMAKE_CURRENT_SOURCE_DIR}/../../MEDLoader + ${CMAKE_CURRENT_SOURCE_DIR}/../../MEDCoupling_Swig + ${CMAKE_CURRENT_SOURCE_DIR}/../../MEDCoupling + ${CMAKE_CURRENT_SOURCE_DIR}/../../INTERP_KERNEL + ${CMAKE_CURRENT_SOURCE_DIR}/../../INTERP_KERNEL/Bases +) + +SET (SWIG_MODULE_ShapeRecogn_EXTRA_DEPS ${ShapeRecogn_SWIG_DPYS_FILES} + ${shaperecogn_HEADERS_HXX} + ${medcoupling_HEADERS_HXX} ${medcoupling_HEADERS_TXX} + ${interpkernel_HEADERS_HXX} ${interpkernel_HEADERS_TXX}) + +IF(WIN32) + SET_PROPERTY(SOURCE ShapeRecogn.i PROPERTY COMPILE_DEFINITIONS WIN32) +ENDIF() + +IF(${CMAKE_VERSION} VERSION_LESS "3.8.0") + SWIG_ADD_MODULE(ShapeRecogn python ShapeRecogn.i) +ELSE() + SWIG_ADD_LIBRARY(ShapeRecogn LANGUAGE python SOURCES ShapeRecogn.i) +ENDIF() + +SWIG_LINK_LIBRARIES(ShapeRecogn ${PYTHON_LIBRARIES} ${PLATFORM_LIBS} shaperecogn medloader medcouplingcpp) +SWIG_CHECK_GENERATION(ShapeRecogn) +IF(WIN32) + SET_TARGET_PROPERTIES(_ShapeRecogn PROPERTIES DEBUG_OUTPUT_NAME _ShapeRecogn_d) + # To increase the size of the .obj file on Windows because ShapeRecognPYTHON_wrap.cxx, generated by SWIG, is too big + TARGET_COMPILE_OPTIONS(_ShapeRecogn PRIVATE /bigobj) +ENDIF(WIN32) + +INSTALL(TARGETS _ShapeRecogn DESTINATION ${MEDCOUPLING_INSTALL_PYTHON}) +INSTALL(FILES ShapeRecogn.i DESTINATION ${MEDCOUPLING_INSTALL_HEADERS}) + +SALOME_INSTALL_SCRIPTS( + ${CMAKE_CURRENT_BINARY_DIR}/ShapeRecogn.py + ${MEDCOUPLING_INSTALL_PYTHON} + EXTRA_DPYS "${SWIG_MODULE_ShapeRecogn_REAL_NAME}") diff --git a/src/ShapeRecogn/Swig/ShapeRecogn.i b/src/ShapeRecogn/Swig/ShapeRecogn.i new file mode 100644 index 000000000..58893b0c6 --- /dev/null +++ b/src/ShapeRecogn/Swig/ShapeRecogn.i @@ -0,0 +1,15 @@ +%module ShapeRecogn + +%include "std_string.i" +%include "MEDCouplingCommon.i" + +%{ +#include "ShapeRecognMesh.hxx" +#include "ShapeRecognMeshBuilder.hxx" +using namespace MEDCoupling; +%} + +%ignore getAreas() const; +%ignore getNodes() const; +%include "ShapeRecognMesh.hxx" +%include "ShapeRecognMeshBuilder.hxx" diff --git a/src/ShapeRecogn/Test/CMakeLists.txt b/src/ShapeRecogn/Test/CMakeLists.txt new file mode 100644 index 000000000..2dc29d57c --- /dev/null +++ b/src/ShapeRecogn/Test/CMakeLists.txt @@ -0,0 +1,65 @@ +# Copyright (C) 2012-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_DIRECTORIES( + ${CPPUNIT_INCLUDE_DIRS} + ${HDF5_INCLUDE_DIRS} + ${MEDFILE_INCLUDE_DIRS} + ${CMAKE_CURRENT_SOURCE_DIR}/.. + ${CMAKE_CURRENT_SOURCE_DIR}/../../MEDLoader + ${CMAKE_CURRENT_SOURCE_DIR}/../../MEDCoupling + ${CMAKE_CURRENT_SOURCE_DIR}/../../INTERP_KERNEL + ${CMAKE_CURRENT_SOURCE_DIR}/../../INTERP_KERNEL/Bases + ${CMAKE_CURRENT_SOURCE_DIR}/../../INTERP_KERNELTest # For common CppUnitTest.hxx file and TestIKUtils.hxx + ) + +SET(TestShapeRecogn_SOURCES + TestShapeRecogn.cxx + MathOpsTest.cxx + PlaneTest.cxx + CylinderTest.cxx + ConeTest.cxx + SphereTest.cxx + TorusTest.cxx +) + +SALOME_ACCUMULATE_ENVIRONMENT(MEDCOUPLING_RESOURCE_DIR "${CMAKE_BINARY_DIR}/resources") +SALOME_GENERATE_TESTS_ENVIRONMENT(tests_env) + +ADD_EXECUTABLE(TestShapeRecogn ${TestShapeRecogn_SOURCES}) +TARGET_LINK_LIBRARIES(TestShapeRecogn shaperecogn InterpKernelTestUtils ${CPPUNIT_LIBRARIES} ${PLATFORM_LIBS}) + +INSTALL(TARGETS TestShapeRecogn DESTINATION ${MEDCOUPLING_INSTALL_BINS}) + +SET(BASE_TESTS TestShapeRecogn) + +FOREACH(test ${BASE_TESTS}) + ADD_TEST(NAME ${test} + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MCTestLauncher.py ${CMAKE_CURRENT_BINARY_DIR}/${test}) + SET_TESTS_PROPERTIES(${test} PROPERTIES ENVIRONMENT "${tests_env}") +ENDFOREACH() + +# Application tests + +SET(TEST_INSTALL_DIRECTORY ${MEDCOUPLING_INSTALL_TESTS}/ShapeRecogn) +INSTALL(TARGETS TestShapeRecogn DESTINATION ${TEST_INSTALL_DIRECTORY}) + +INSTALL(FILES CTestTestfileInstall.cmake + DESTINATION ${TEST_INSTALL_DIRECTORY} + RENAME CTestTestfile.cmake) diff --git a/src/ShapeRecogn/Test/CTestTestfileInstall.cmake b/src/ShapeRecogn/Test/CTestTestfileInstall.cmake new file mode 100644 index 000000000..ab32ede46 --- /dev/null +++ b/src/ShapeRecogn/Test/CTestTestfileInstall.cmake @@ -0,0 +1,31 @@ +# Copyright (C) 2015-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 +# + +SET(TEST_NAMES + TestShapeRecogn +) + +FOREACH(tfile ${TEST_NAMES}) + SET(TEST_NAME ${COMPONENT_NAME}_${tfile}) + ADD_TEST(${TEST_NAME} python3 MCTestLauncher.py ${tfile}) + SET_TESTS_PROPERTIES(${TEST_NAME} PROPERTIES + LABELS "${COMPONENT_NAME}" + TIMEOUT ${TIMEOUT} + ) +ENDFOREACH() diff --git a/src/ShapeRecogn/Test/ConeTest.cxx b/src/ShapeRecogn/Test/ConeTest.cxx new file mode 100644 index 000000000..d25f4a792 --- /dev/null +++ b/src/ShapeRecogn/Test/ConeTest.cxx @@ -0,0 +1,324 @@ +#include "ConeTest.hxx" + +#include "ShapeRecognMeshBuilder.hxx" +#include "Areas.hxx" +#include "MathOps.hxx" +#include "TestInterpKernelUtils.hxx" // getResourceFile() + +using namespace MEDCoupling; + +void ConeTest::setUp() +{ + std::string file = INTERP_TEST::getResourceFile("ShapeRecognCone.med", 3); + srMesh = new ShapeRecognMeshBuilder(file); + srMesh->recognize(); + areas = srMesh->getAreas(); +} + +void ConeTest::tearDown() +{ + if (srMesh != 0) + delete srMesh; + areas = 0; +} + +void ConeTest::testNumberOfAreas() +{ + CPPUNIT_ASSERT_EQUAL(3, (int)areas->getNumberOfAreas()); +} + +void ConeTest::testComputePlaneProperties() +{ + const Nodes *nodes = srMesh->getNodes(); + Areas areas(nodes); + mcIdType areaId = areas.addArea(); + std::vector nodeIds{560, 561, 562, 563, 564}; + // Check the coordinates + std::vector coordsRef{ + -0.13274028, + -1.19035036, + 7.0, + 0.62579096, + -0.97231524, + 7.0, + -0.75942852, + -0.99273623, + 7.0, + 0.76454326, + 0.85588873, + 7.0, + 1.12546475, + -0.04924272, + 7.0, + }; + for (size_t i = 0; i < nodeIds.size(); ++i) + { + mcIdType nodeId = nodeIds[i]; + areas.addNode(areaId, nodeId); + std::array coords = nodes->getCoordinates(nodeId); + CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i], coords[0], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 1], coords[1], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 2], coords[2], 1E-6); + } + areas.computeProperties(areaId); + // Check normal + std::array normal = areas.getNormal(areaId); + std::array normalRef{0.0, 0.0, 1.0}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(normalRef[i], normal[i], 1E-6); +} + +void ConeTest::testComputeCylinderProperties() +{ + const Nodes *nodes = srMesh->getNodes(); + Areas areas(nodes); + mcIdType areaId = areas.addArea(); + // check coordinates + std::vector coordsRef{ + 8.58402821e+00, + -2.10248053e-15, + 1.41597179e+00, + 7.87604232e+00, + -1.92907400e-15, + 2.12395768e+00, + 7.16805642e+00, + -1.75566747e-15, + 2.83194358e+00, + }; + // Add nodes + std::vector nodeIds{26, 27, 28}; + for (size_t i = 0; i < nodeIds.size(); ++i) + { + mcIdType nodeId = nodeIds[i]; + areas.addNode(areaId, nodeId); + std::array coords = nodes->getCoordinates(nodeId); + CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i], coords[0], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 1], coords[1], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 2], coords[2], 1E-6); + } + areas.computeProperties(areaId); + // Check radius + double radius = areas.getRadius(areaId); + double radiusRef = 11.016809459292505; + CPPUNIT_ASSERT_DOUBLES_EQUAL(radiusRef, radius, 1E-6); + // Check axis + std::array axis = areas.getAxis(areaId); + std::array axisRef{-2.85185021e-02, 2.47000552e-05, 9.99593264e-01}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(axisRef[i], axis[i], 1E-6); + // Check axis point + std::array axisPoint = areas.getAxisPoint(areaId); + std::array axisPointRef{9.55163389e-02, -2.27195412e-05, -5.67562579e+00}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(axisPointRef[i], axisPoint[i], 1E-6); +} + +void ConeTest::testComputeConeProperties() +{ + std::srand(0); + const Nodes *nodes = srMesh->getNodes(); + Areas areas(nodes); + mcIdType areaId = areas.addArea(PrimitiveType::Cone); + std::vector coordsRef{ + 8.58402821e+00, + -2.10248053e-15, + 1.41597179e+00, + 7.87604232e+00, + -1.92907400e-15, + 2.12395768e+00, + 7.16805642e+00, + -1.75566747e-15, + 2.83194358e+00, + 6.46007053e+00, + -1.58226094e-15, + 3.53992947e+00}; + // Add nodes + std::vector nodeIds{26, 27, 28, 29}; + for (size_t i = 0; i < nodeIds.size(); ++i) + { + mcIdType nodeId = nodeIds[i]; + areas.addNode(areaId, nodeId); + std::array coords = nodes->getCoordinates(nodeId); + CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i], coords[0], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 1], coords[1], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 2], coords[2], 1E-6); + } + areas.computeProperties(areaId); + // Radius + double radiusRef = 10.546213172989718; + CPPUNIT_ASSERT_DOUBLES_EQUAL(radiusRef, areas.getRadius(areaId), 1E-6); + // Angle + double angleRef = 0.7567719849399294; + CPPUNIT_ASSERT_DOUBLES_EQUAL(angleRef, areas.getAngle(areaId), 1E-6); + // Axis + std::array axisRef{-2.99093478e-02, 1.48776702e-05, 9.99552615e-01}; + std::array axis = areas.getAxis(areaId); + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(axisRef[i], axis[i], 1E-6); + // Axis Point + std::array axisPointRef{ + 7.43242419e-02, -1.70396559e-05, -4.98890939e+00}; + std::array axisPoint = areas.getAxisPoint(areaId); + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(axisPointRef[i], axisPoint[i], 1E-6); + // Apex + std::array apexRef{ + -3.87870047e-01, 1.98282729e-04, 1.03928162e+01}; + std::array apex = areas.getApex(areaId); + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(apexRef[i], apex[i], 1E-6); +} + +void ConeTest::testFirstArea() +{ + // primitive type + CPPUNIT_ASSERT_EQUAL(PrimitiveType::Plane, areas->getPrimitiveType(0)); + // node ids + std::vector nodeIdsRef{ + 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559, + 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, + 571}; + std::vector nodeIds = areas->getNodeIds(0); + std::sort(nodeIds.begin(), nodeIds.end()); + CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size()); + for (size_t i = 0; i < nodeIds.size(); ++i) + CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]); + // normal + std::array normal = areas->getNormal(0); + std::array normalRef{0.0, 0.0, 1.0}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(normalRef[i], normal[i], 1E-6); +} + +void ConeTest::testSecondArea() +{ + // primitive type + CPPUNIT_ASSERT_EQUAL(PrimitiveType::Plane, areas->getPrimitiveType(0)); + // node ids + std::vector nodeIdsRef = { + 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, + 584, 585, 586, 587, 588, 589, 590, 591, 592, 593, 594, 595, + 596, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606, 607, + 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619, + 620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, + 632, 633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, + 644, 645, 646, 647, 648, 649, 650, 651, 652, 653, 654, 655, + 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667, + 668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, + 680, 681, 682, 683, 684, 685, 686, 687, 688, 689, 690, 691, + 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, + 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715, + 716, 717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727, + 728, 729, 730, 731, 732, 733, 734, 735, 736, 737, 738, 739, + 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751, + 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762, 763, + 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775, + 776, 777, 778, 779, 780, 781, 782, 783, 784, 785, 786, 787, + 788, 789, 790, 791, 792, 793, 794, 795, 796, 797, 798, 799, + 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 810, 811, + 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823, + 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835, + 836, 837, 838, 839, 840, 841, 842, 843, 844, 845, 846, 847, + 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 858, 859, + 860, 861, 862, 863, 864, 865, 866, 867, 868, 869, 870, 871, + 872, 873, 874, 875, 876, 877, 878, 879, 880, 881, 882, 883, + 884, 885, 886, 887, 888, 889, 890, 891, 892, 893, 894, 895, + 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907, + 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919, + 920, 921, 922, 923, 924, 925, 926, 927, 928, 929, 930, 931, + 932, 933, 934, 935, 936, 937, 938, 939, 940, 941, 942, 943, + 944, 945, 946, 947, 948, 949, 950, 951, 952, 953, 954, 955, + 956}; + std::vector nodeIds = areas->getNodeIds(1); + std::sort(nodeIds.begin(), nodeIds.end()); + CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size()); + for (size_t i = 0; i < nodeIds.size(); ++i) + CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]); + // normal + std::array normal = areas->getNormal(1); + std::array normalRef = {0.0, 0.0, -1.0}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(normalRef[i], normal[i], 1E-6); +} + +void ConeTest::testThirdArea() +{ + // primitive type + CPPUNIT_ASSERT_EQUAL(PrimitiveType::Cone, areas->getPrimitiveType(2)); + // node ids + std::vector nodeIdsRef{ + 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, + 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, + 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, + 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, + 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, + 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, + 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, + 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, + 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, + 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, + 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, + 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, + 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, + 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, + 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, + 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, + 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, + 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, + 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, + 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, + 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, + 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, + 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, + 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, + 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, + 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, + 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, + 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, + 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, + 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, + 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, + 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, + 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, + 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, + 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, + 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, + 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, + 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, + 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, + 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, + 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, + 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, + 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, + 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, + 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, + 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, + 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, + 540, 541, 542, 543, 544, 545, 546, 547, 548}; + std::vector nodeIds = areas->getNodeIds(2); + std::sort(nodeIds.begin(), nodeIds.end()); + CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size()); + for (size_t i = 0; i < nodeIds.size(); ++i) + CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]); + // radius + CPPUNIT_ASSERT_DOUBLES_EQUAL(10.0, areas->getRadius(2), 5E-1); + // angle + CPPUNIT_ASSERT_DOUBLES_EQUAL(M_PI_4, areas->getAngle(2), 1E-2); + // axis + std::array axis = areas->getAxis(2); + std::array axisRef{0.0, 0.0, 1.0}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(axisRef[i], axis[i], 1E-2); + // // axis point + // std::array axisPoint = areas->getAxisPoint(2); + // std::array axisPointRef{ + // -6.50039462e-03, -5.76665576e-05, -3.99198765e+00}; + // for (size_t i = 0; i < 3; ++i) + // CPPUNIT_ASSERT_DOUBLES_EQUAL( + // axisPointRef[i], axisPoint[i], 1E-2); + // apex + std::array apex = areas->getApex(2); + std::array apexRef{0.0, 0.0, 10.0}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(apexRef[i], apex[i], 1E-1); +} diff --git a/src/ShapeRecogn/Test/ConeTest.hxx b/src/ShapeRecogn/Test/ConeTest.hxx new file mode 100644 index 000000000..20b1bef3f --- /dev/null +++ b/src/ShapeRecogn/Test/ConeTest.hxx @@ -0,0 +1,62 @@ +// 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 __CONETEST_HXX__ +#define __CONETEST_HXX__ + +#include +#include + +namespace MEDCoupling +{ + class ShapeRecognMeshBuilder; + class Areas; + + class ConeTest : public CppUnit::TestFixture + { + CPPUNIT_TEST_SUITE(ConeTest); + CPPUNIT_TEST(testNumberOfAreas); + CPPUNIT_TEST(testComputePlaneProperties); + CPPUNIT_TEST(testComputeCylinderProperties); + CPPUNIT_TEST(testComputeConeProperties); + CPPUNIT_TEST(testFirstArea); + CPPUNIT_TEST(testSecondArea); + CPPUNIT_TEST(testThirdArea); + CPPUNIT_TEST_SUITE_END(); + + public: + void setUp() override; + void tearDown() override; + + void testComputePlaneProperties(); + void testComputeCylinderProperties(); + void testComputeConeProperties(); + + void testNumberOfAreas(); + void testFirstArea(); + void testSecondArea(); + void testThirdArea(); + + private: + ShapeRecognMeshBuilder *srMesh = 0; + const Areas *areas; + }; +}; + +#endif // __CONETEST_HXX__ diff --git a/src/ShapeRecogn/Test/CylinderTest.cxx b/src/ShapeRecogn/Test/CylinderTest.cxx new file mode 100644 index 000000000..47f6ac3d3 --- /dev/null +++ b/src/ShapeRecogn/Test/CylinderTest.cxx @@ -0,0 +1,136 @@ +#include "CylinderTest.hxx" + +#include "ShapeRecognMeshBuilder.hxx" +#include "Areas.hxx" +#include "MathOps.hxx" +#include "TestInterpKernelUtils.hxx" // getResourceFile() + +using namespace MEDCoupling; + +void CylinderTest::setUp() +{ + std::string file = INTERP_TEST::getResourceFile("ShapeRecognCylinder.med", 3); + srMesh = new ShapeRecognMeshBuilder(file); + srMesh->recognize(); + areas = srMesh->getAreas(); +} + +void CylinderTest::tearDown() +{ + if (srMesh != 0) + delete srMesh; + areas = 0; +} + +void CylinderTest::testNumberOfAreas() +{ + CPPUNIT_ASSERT_EQUAL(3, (int)areas->getNumberOfAreas()); +} + +void CylinderTest::testFirstArea() +{ + // primitive type + CPPUNIT_ASSERT_EQUAL(PrimitiveType::Cylinder, areas->getPrimitiveType(0)); + // node ids + std::vector nodeIdsRef{ + 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, + 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, + 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, + 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, + 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, + 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, + 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, + 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, + 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, + 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, + 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, + 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, + 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, + 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, + 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, + 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, + 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, + 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, + 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, + 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, + 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, + 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, + 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, + 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, + 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, + 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, + 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, + 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, + 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, + 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, + 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, + 369}; + std::vector nodeIds = areas->getNodeIds(0); + std::sort(nodeIds.begin(), nodeIds.end()); + CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size()); + for (size_t i = 0; i < nodeIds.size(); ++i) + CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]); + // radius + CPPUNIT_ASSERT_DOUBLES_EQUAL(4.993494657779537, areas->getRadius(0), 1E-2); + // axis + std::array axis = areas->getAxis(0); + std::array axisRef{7.66631075e-04, -1.59966800e-04, 9.99999693e-01}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(axisRef[i], axis[i], 1E-4); + // axis point + std::array axisPoint = areas->getAxisPoint(0); + std::array axisPointRef{3.78927537e-03, -2.03409415e-03, 5.03355802e+00}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL( + axisPointRef[i], axisPoint[i], 1E-2); +} + +void CylinderTest::testSecondArea() +{ + // primitive type + CPPUNIT_ASSERT_EQUAL(PrimitiveType::Plane, areas->getPrimitiveType(1)); + // node ids + std::vector nodeIdsRef{ + 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, + 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, + 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, + 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, + 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, + 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, + 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, + 440, 441, 442}; + std::vector nodeIds = areas->getNodeIds(1); + std::sort(nodeIds.begin(), nodeIds.end()); + CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size()); + for (size_t i = 0; i < nodeIds.size(); ++i) + CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]); + // normal + std::array normal = areas->getNormal(1); + std::array normalRef{0.0, 0.0, 1.0}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(normalRef[i], normal[i], 1E-6); +} + +void CylinderTest::testThirdArea() +{ + // primitive type + CPPUNIT_ASSERT_EQUAL(PrimitiveType::Plane, areas->getPrimitiveType(2)); + // node ids + std::vector nodeIdsRef{ + 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, + 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, + 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, + 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, + 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, + 493, 494, 495, 496, 497, 498, 499}; + std::vector nodeIds = areas->getNodeIds(2); + std::sort(nodeIds.begin(), nodeIds.end()); + CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size()); + for (size_t i = 0; i < nodeIds.size(); ++i) + CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]); + // normal + std::array normal = areas->getNormal(2); + std::array normalRef{0.0, 0.0, -1.0}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(normalRef[i], normal[i], 1E-6); +} diff --git a/src/ShapeRecogn/Test/CylinderTest.hxx b/src/ShapeRecogn/Test/CylinderTest.hxx new file mode 100644 index 000000000..a35ac1481 --- /dev/null +++ b/src/ShapeRecogn/Test/CylinderTest.hxx @@ -0,0 +1,55 @@ +// 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 __CYLINDERTEST_HXX__ +#define __CYLINDERTEST_HXX__ + +#include +#include + +namespace MEDCoupling +{ + class ShapeRecognMeshBuilder; + class Areas; + + class CylinderTest : public CppUnit::TestFixture + { + CPPUNIT_TEST_SUITE(CylinderTest); + CPPUNIT_TEST(testNumberOfAreas); + CPPUNIT_TEST(testFirstArea); + CPPUNIT_TEST(testSecondArea); + CPPUNIT_TEST(testThirdArea); + CPPUNIT_TEST_SUITE_END(); + + public: + void setUp() override; + void tearDown() override; + + void testNumberOfAreas(); + void testFirstArea(); + void testSecondArea(); + void testThirdArea(); + + private: + ShapeRecognMeshBuilder *srMesh = 0; + const Areas *areas; + }; +}; + +#endif // __CYLINDERTEST_HXX__ diff --git a/src/ShapeRecogn/Test/MathOpsTest.cxx b/src/ShapeRecogn/Test/MathOpsTest.cxx new file mode 100644 index 000000000..033f9921b --- /dev/null +++ b/src/ShapeRecogn/Test/MathOpsTest.cxx @@ -0,0 +1,97 @@ +#include "MathOpsTest.hxx" +#include "MathOps.hxx" + +using namespace MEDCoupling; + +void MathOpsTest::testLstsq() +{ + std::vector a = { + 0.69473263, 0.83318004, 0.60822673, 0.59243878, 0.82872553, + 0.84048546, 0.95698819, 0.02171218, 0.27683381, 0.20628928, + 0.80895323, 0.4207767, 0.37468575, 0.86258204, 0.42571846}; + std::vector b = { + 0.91167508, 0.95878824, 0.7234827, 0.51753917, 0.18603306}; + std::vector x = MathOps::lstsq(a, b); + std::array xRef = {0.35719095, 0.5134345, 0.26039343}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(xRef[i], x[i], 1E-6); +} + +void MathOpsTest::testLstsq2() +{ + std::vector a = { + 0.4564562, 0.3517006, + 0.28928215, 0.72309086, + 0.05944836, 0.56024464}; + std::vector b = {0.98902712, 0.46791812}; + std::vector x = MathOps::lstsq(a, b); + std::array xRef = {2.10752524, 0.2636243, -0.82807416}; + for (size_t i = 0; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(xRef[i], x[i], 1E-6); +} + +void MathOpsTest::testLstsqBig() +{ + std::vector a(5000000 * 3, 1.0); + std::vector b(5000000, 1.0); + std::vector x = MathOps::lstsq(a, b); + std::array xRef = {1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}; + for (size_t i = 2; i < 3; ++i) + CPPUNIT_ASSERT_DOUBLES_EQUAL(xRef[i], x[i], 1E-6); +} + +void MathOpsTest::testComputeCov() +{ + std::vector coordinates{ + 1, 2, -3, -8, 6, -3, 9, 2, -1, -10, -11, -120}; + std::array + covMatrix = MathOps::computeCov(coordinates); + CPPUNIT_ASSERT_DOUBLES_EQUAL(76.666666, covMatrix[0], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(26.666666, covMatrix[1], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(319.333333, covMatrix[2], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(26.666666, covMatrix[3], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(54.916666, covMatrix[4], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(420.75, covMatrix[5], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(319.333333, covMatrix[6], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(420.75, covMatrix[7], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(3462.25, covMatrix[8], 1E-6); +} + +void MathOpsTest::testComputePCAFirstAxis() +{ + std::vector coordinates{ + 1, 2, -3, -8, 6, -3, 9, 2, -1, -10, -11, -120}; + std::array + axis = MathOps::computePCAFirstAxis(coordinates); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09198798, axis[0], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.11994164, axis[1], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.9885101, axis[2], 1E-6); +} + +void MathOpsTest::testComputeAngles() +{ + std::vector directions{ + 1, 2, -3, -8, 6, -3, 9, 2, -1, -10, -11, -120}; + std::array axis{1, 2, 3}; + std::vector angles = MathOps::computeAngles( + directions, axis); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.86054803, angles[0], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.69914333, angles[1], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.27845478, angles[2], 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(2.61880376, angles[3], 1E-6); +} + +void MathOpsTest::testComputeBaseFromNormal() +{ + std::array normal = {1.0, 2.0, 3.0}; + std::array base = MathOps::computeBaseFromNormal(normal); + std::array baseRef = { + -0.53452248, 0.77454192, -0.33818712, + -0.80178373, -0.33818712, 0.49271932}; + for (size_t i = 0; i < 6; ++i) + { + std::ostringstream message; + message << "Mismatch at index " << i; + CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(message.str().c_str(), baseRef[i], base[i], 1E-6); + } +} diff --git a/src/ShapeRecogn/Test/MathOpsTest.hxx b/src/ShapeRecogn/Test/MathOpsTest.hxx new file mode 100644 index 000000000..2988a950b --- /dev/null +++ b/src/ShapeRecogn/Test/MathOpsTest.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 __MATHOPSTEST_HXX__ +#define __MATHOPSTEST_HXX__ + +#include +#include + +namespace MEDCoupling +{ + class MathOpsTest : public CppUnit::TestFixture + { + CPPUNIT_TEST_SUITE(MathOpsTest); + CPPUNIT_TEST(testLstsq); + CPPUNIT_TEST(testLstsq2); + CPPUNIT_TEST(testLstsqBig); + CPPUNIT_TEST(testComputeCov); + CPPUNIT_TEST(testComputePCAFirstAxis); + CPPUNIT_TEST(testComputeAngles); + CPPUNIT_TEST(testComputeBaseFromNormal); + CPPUNIT_TEST_SUITE_END(); + + public: + static void testLstsq(); + static void testLstsq2(); + static void testLstsqBig(); + static void testComputeCov(); + static void testComputePCAFirstAxis(); + static void testComputeAngles(); + static void testComputeBaseFromNormal(); + }; +}; + +#endif // __MATHOPSTEST_HXX__ diff --git a/src/ShapeRecogn/Test/PlaneTest.cxx b/src/ShapeRecogn/Test/PlaneTest.cxx new file mode 100644 index 000000000..6f47f161c --- /dev/null +++ b/src/ShapeRecogn/Test/PlaneTest.cxx @@ -0,0 +1,44 @@ +#include "PlaneTest.hxx" + +#include "ShapeRecognMeshBuilder.hxx" +#include "Areas.hxx" +#include "MathOps.hxx" +#include "TestInterpKernelUtils.hxx" // getResourceFile() + +using namespace MEDCoupling; + +void PlaneTest::setUp() +{ + std::string file = INTERP_TEST::getResourceFile("ShapeRecognPlane.med", 3); + srMesh = new ShapeRecognMeshBuilder(file); + srMesh->recognize(); + areas = srMesh->getAreas(); +} + +void PlaneTest::tearDown() +{ + if (srMesh != 0) + delete srMesh; + areas = 0; +} + +void PlaneTest::testArea() +{ + CPPUNIT_ASSERT_EQUAL(36, (int)areas->getNumberOfNodes(0)); + CPPUNIT_ASSERT_EQUAL(1, (int)areas->getNumberOfAreas()); + // Normal + std::array normal = areas->getNormal(0); + std::array normalRef = {0.781525, 0.310606, -0.541056}; + std::array affinePoint = areas->getAffinePoint(0); + double proportion0 = normal[0] / normalRef[0]; + double proportion1 = normal[1] / normalRef[1]; + double proportion2 = normal[2] / normalRef[2]; + double proportion3 = MathOps::dot(normal, affinePoint) / MathOps::dot(normalRef, affinePoint); + // Check proportions between the normal vectors of the two planes + CPPUNIT_ASSERT_DOUBLES_EQUAL(proportion0, proportion1, 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(proportion1, proportion2, 1E-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(proportion2, proportion3, 1E-6); + // Check the angle + double angle = MathOps::computeAngle(normal, normalRef); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, angle, 1E-6); +} \ No newline at end of file diff --git a/src/ShapeRecogn/Test/PlaneTest.hxx b/src/ShapeRecogn/Test/PlaneTest.hxx new file mode 100644 index 000000000..af6eb4755 --- /dev/null +++ b/src/ShapeRecogn/Test/PlaneTest.hxx @@ -0,0 +1,49 @@ +// 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 __PLANETEST_HXX__ +#define __PLANETEST_HXX__ + +#include +#include + +namespace MEDCoupling +{ + class ShapeRecognMeshBuilder; + class Areas; + + class PlaneTest : public CppUnit::TestFixture + { + CPPUNIT_TEST_SUITE(PlaneTest); + CPPUNIT_TEST(testArea); + CPPUNIT_TEST_SUITE_END(); + + public: + void setUp() override; + void tearDown() override; + + void testArea(); + + private: + ShapeRecognMeshBuilder *srMesh = 0; + const Areas *areas; + }; +}; + +#endif // __PLANETEST_HXX__ diff --git a/src/ShapeRecogn/Test/SphereTest.cxx b/src/ShapeRecogn/Test/SphereTest.cxx new file mode 100644 index 000000000..4d7f64ecb --- /dev/null +++ b/src/ShapeRecogn/Test/SphereTest.cxx @@ -0,0 +1,36 @@ +#include "SphereTest.hxx" + +#include "ShapeRecognMeshBuilder.hxx" +#include "Areas.hxx" +#include "MathOps.hxx" +#include "TestInterpKernelUtils.hxx" // getResourceFile() + +using namespace MEDCoupling; + +void SphereTest::setUp() +{ + std::string file = INTERP_TEST::getResourceFile("ShapeRecognSphere.med", 3); + srMesh = new ShapeRecognMeshBuilder(file); + srMesh->recognize(); + areas = srMesh->getAreas(); +} + +void SphereTest::tearDown() +{ + if (srMesh != 0) + delete srMesh; + areas = 0; +} + +void SphereTest::testArea() +{ + CPPUNIT_ASSERT_EQUAL(1, (int)areas->getNumberOfAreas()); + // 8 double nodes so 147 - 6 nodes + CPPUNIT_ASSERT_EQUAL(141, (int)areas->getNumberOfNodes(0)); + CPPUNIT_ASSERT_EQUAL(PrimitiveType::Sphere, areas->getPrimitiveType(0)); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, areas->getRadius(0), 1E-2); + std::array centerRef = {5.3, -6.7, -9.02}; + std::array center = areas->getCenter(0); + for (size_t j = 0; j < 3; ++j) + CPPUNIT_ASSERT_DOUBLES_EQUAL(centerRef[j], center[j], 1E-2); +} diff --git a/src/ShapeRecogn/Test/SphereTest.hxx b/src/ShapeRecogn/Test/SphereTest.hxx new file mode 100644 index 000000000..6500def00 --- /dev/null +++ b/src/ShapeRecogn/Test/SphereTest.hxx @@ -0,0 +1,49 @@ +// 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 __SPHERETEST_HXX__ +#define __SPHERETEST_HXX__ + +#include +#include + +namespace MEDCoupling +{ + class ShapeRecognMeshBuilder; + class Areas; + + class SphereTest : public CppUnit::TestFixture + { + CPPUNIT_TEST_SUITE(SphereTest); + CPPUNIT_TEST(testArea); + CPPUNIT_TEST_SUITE_END(); + + public: + void setUp() override; + void tearDown() override; + + void testArea(); + + private: + ShapeRecognMeshBuilder *srMesh = 0; + const Areas *areas; + }; +}; + +#endif // __SPHERETEST_HXX__ diff --git a/src/ShapeRecogn/Test/TestShapeRecogn.cxx b/src/ShapeRecogn/Test/TestShapeRecogn.cxx new file mode 100644 index 000000000..b7cd2af14 --- /dev/null +++ b/src/ShapeRecogn/Test/TestShapeRecogn.cxx @@ -0,0 +1,34 @@ +// 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 "MathOpsTest.hxx" +#include "PlaneTest.hxx" +#include "CylinderTest.hxx" +#include "ConeTest.hxx" +#include "SphereTest.hxx" +#include "TorusTest.hxx" + +CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::MathOpsTest); +CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::PlaneTest); +CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::CylinderTest); +CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::ConeTest); +CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::SphereTest); +CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::TorusTest); + +#include "BasicMainTest.hxx" diff --git a/src/ShapeRecogn/Test/TorusTest.cxx b/src/ShapeRecogn/Test/TorusTest.cxx new file mode 100644 index 000000000..6e24d10ee --- /dev/null +++ b/src/ShapeRecogn/Test/TorusTest.cxx @@ -0,0 +1,38 @@ +#include "TorusTest.hxx" + +#include "ShapeRecognMeshBuilder.hxx" +#include "Areas.hxx" +#include "MathOps.hxx" +#include "TestInterpKernelUtils.hxx" // getResourceFile() + +using namespace MEDCoupling; + +void TorusTest::setUp() +{ + std::string file = INTERP_TEST::getResourceFile("ShapeRecognTorus.med", 3); + srMesh = new ShapeRecognMeshBuilder(file); + srMesh->recognize(); + areas = srMesh->getAreas(); +} + +void TorusTest::tearDown() +{ + if (srMesh != 0) + delete srMesh; + areas = 0; +} + +void TorusTest::testArea() +{ + CPPUNIT_ASSERT_EQUAL(275, (int)srMesh->getNodes()->getNbNodes()); + CPPUNIT_ASSERT_EQUAL(1, (int)areas->getNumberOfAreas()); + CPPUNIT_ASSERT_EQUAL(PrimitiveType::Torus, areas->getPrimitiveType(0)); + // Some nodes are unknown + CPPUNIT_ASSERT_EQUAL(272, (int)areas->getNumberOfNodes(0)); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.843297, areas->getMinorRadius(0), 1E-2); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.156428, areas->getRadius(0), 1E-1); + std::array centerRef = {7.687022, -3.726887, -9.02}; + std::array center = areas->getCenter(0); + for (size_t j = 0; j < 3; ++j) + CPPUNIT_ASSERT_DOUBLES_EQUAL(centerRef[j], center[j], 1E-2); +} diff --git a/src/ShapeRecogn/Test/TorusTest.hxx b/src/ShapeRecogn/Test/TorusTest.hxx new file mode 100644 index 000000000..07ae86d3a --- /dev/null +++ b/src/ShapeRecogn/Test/TorusTest.hxx @@ -0,0 +1,49 @@ +// 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 __TORUSTEST_HXX__ +#define __TORUSTEST_HXX__ + +#include +#include + +namespace MEDCoupling +{ + class ShapeRecognMeshBuilder; + class Areas; + + class TorusTest : public CppUnit::TestFixture + { + CPPUNIT_TEST_SUITE(TorusTest); + CPPUNIT_TEST(testArea); + CPPUNIT_TEST_SUITE_END(); + + public: + void setUp() override; + void tearDown() override; + + void testArea(); + + private: + ShapeRecognMeshBuilder *srMesh = 0; + const Areas *areas; + }; +}; + +#endif // __TORUSTEST_HXX__