From 7c98677a6b2c577669ee87886aba3c145ee6bea6 Mon Sep 17 00:00:00 2001 From: Konstantinos Samouchos Date: Fri, 1 Mar 2024 18:08:11 +0100 Subject: [PATCH] levelset: micro-optimization in the evaluation of gradient --- src/levelset/levelSetSegmentationObject.cpp | 44 ++++++++++++++++++--- src/levelset/levelSetSegmentationObject.hpp | 2 + 2 files changed, 40 insertions(+), 6 deletions(-) diff --git a/src/levelset/levelSetSegmentationObject.cpp b/src/levelset/levelSetSegmentationObject.cpp index 70f18a3001..661550ef97 100644 --- a/src/levelset/levelSetSegmentationObject.cpp +++ b/src/levelset/levelSetSegmentationObject.cpp @@ -225,20 +225,22 @@ bitpit::LevelSetSurfaceSmoothing LevelSetSegmentationSurfaceInfo::getSurfaceSmoo * @param[in] point are the coordinates of point * @param[in] segmentItr is an iterator pointing to the closest segment * @param[in] signedDistance controls if the signed or unsigned distance will be evaluated + * @param[out] the distance function at the specified point * @return The distance function at the specified point. */ double LevelSetSegmentationSurfaceInfo::evalDistance(const std::array &point, const SegmentConstIterator &segmentItr, - bool signedDistance) const + bool signedDistance, + std::array *distanceVector) const { // Project the point on the surface and evaluate the point-projection vector int nSegmentVertices = segmentItr->getVertexCount(); BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES); std::array pointProjection = evalProjection(point, segmentItr, lambda); - std::array pointProjectionVector = point - pointProjection; + (*distanceVector) = point - pointProjection; // Evaluate unsigned distance - double unsignedDistance = norm2(pointProjectionVector); + double unsignedDistance = norm2(*distanceVector); if (!signedDistance) { return unsignedDistance; } @@ -249,7 +251,7 @@ double LevelSetSegmentationSurfaceInfo::evalDistance(const std::array // plane. This case is not supported, because it would require to evaluate the sign taking // into account the the curvature of the surface. std::array pseudoNormal = computePseudoNormal(segmentItr, lambda); - double pointProjectionNormalComponent = dotProduct(pointProjectionVector, pseudoNormal); + double pointProjectionNormalComponent = dotProduct(*distanceVector, pseudoNormal); double distanceTolerance = m_surface->getTol(); if (utils::DoubleFloatingEqual()(pointProjectionNormalComponent, 0., distanceTolerance, distanceTolerance)) { @@ -262,6 +264,22 @@ double LevelSetSegmentationSurfaceInfo::evalDistance(const std::array return sign(pointProjectionNormalComponent) * unsignedDistance; } +/*! + * Evaluate the distance function at the specified point. + * + * @param[in] point are the coordinates of point + * @param[in] segmentItr is an iterator pointing to the closest segment + * @param[in] signedDistance controls if the signed or unsigned distance will be evaluated + * @return The distance function at the specified point. + */ +double LevelSetSegmentationSurfaceInfo::evalDistance(const std::array &point, + const SegmentConstIterator &segmentItr, + bool signedDistance) const +{ + std::array distanceVector; + return evalDistance(point, segmentItr, signedDistance, &distanceVector); +} + /*! * Evaluate the distance vector function at the specified point. * @@ -299,6 +317,19 @@ std::array LevelSetSegmentationSurfaceInfo::evalNormal(const std::arr return projectionNormal; } +/*! + * Compute the pseudo-normal at specified point of the given segment. + * + * @param[in] segmentItr is an iterator pointing to the closest segment + * @param[in] lambda are the barycentric coordinates of the point + * @return the pseudo-normal at specified point of the given segment + */ +std::array LevelSetSegmentationSurfaceInfo::evalPseudoNormal(const SegmentConstIterator &segmentItr, + const double *lambda ) const +{ + return computePseudoNormal(segmentItr, lambda); +} + /*! * Evaluate the projection of the given point on the surface created based on * the points representing the specified segment. The surface passes from these @@ -2880,7 +2911,8 @@ std::array LevelSetSegmentationObject::_evalGradient(const std::array< // Evaluate the distance of the point from the surface LevelSetSegmentationSurfaceInfo::SegmentConstIterator supportItr = getSurface().getCellConstIterator(support); - double distance = m_surfaceInfo->evalDistance(point, supportItr, signedLevelSet); + std::array distanceVector; + double distance = m_surfaceInfo->evalDistance(point, supportItr, signedLevelSet, &distanceVector); // Early return if the point lies on the surface if (evalValueSign(distance) == 0) { @@ -2892,7 +2924,7 @@ std::array LevelSetSegmentationObject::_evalGradient(const std::array< } // Evaluate levelset gradient - std::array gradient = m_surfaceInfo->evalDistanceVector(point, supportItr) / distance; + std::array gradient = distanceVector / distance; return gradient; } diff --git a/src/levelset/levelSetSegmentationObject.hpp b/src/levelset/levelSetSegmentationObject.hpp index 91e83f102c..38cb0b24a6 100644 --- a/src/levelset/levelSetSegmentationObject.hpp +++ b/src/levelset/levelSetSegmentationObject.hpp @@ -73,10 +73,12 @@ class LevelSetSegmentationSurfaceInfo { double getFeatureAngle() const; LevelSetSurfaceSmoothing getSurfaceSmoothing() const; + double evalDistance(const std::array &point, const SegmentConstIterator &segmentItr, bool signedDistance, std::array *distanceVector) const; double evalDistance(const std::array &point, const SegmentConstIterator &segmentItr, bool signedDistance) const; std::array evalDistanceVector(const std::array &point, const SegmentConstIterator &segmentItr) const; std::array evalNormal(const std::array &point, const SegmentConstIterator &segmentItr) const; + std::array evalPseudoNormal(const SurfUnstructured::CellConstIterator &segmentIterator, const double *lambda) const; void evalBarycentricCoordinates(const std::array &point, const SegmentConstIterator &segmentItr, double *lambda) const; void evalProjection(const std::array &point, const SegmentConstIterator &segmentItr, std::array *projectionPoint, std::array *projectionNormal) const;