Skip to content

Commit

Permalink
levelset: micro-optimization in the evaluatiion of gradient
Browse files Browse the repository at this point in the history
  • Loading branch information
ksamouchos authored and Konstantinos committed Apr 12, 2024
1 parent f6bb3ac commit b3817e5
Show file tree
Hide file tree
Showing 2 changed files with 32 additions and 3 deletions.
34 changes: 31 additions & 3 deletions src/levelset/levelSetSegmentationObject.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -299,6 +299,19 @@ std::array<double, 3> 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<double,3> 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
Expand Down Expand Up @@ -2872,19 +2885,34 @@ std::array<double,3> 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);

int nSegmentVertices = supportItr->getVertexCount();
BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
std::array<double, 3> projectionPoint = m_surfaceInfo->evalProjection(point, supportItr, lambda);
std::array<double, 3> gradient = point - projectionPoint;

double unsignedDistance = norm2(gradient);

// Early return if the point lies on the surface
if (evalValueSign(distance) == 0) {
double distanceTolerance = m_surfaceInfo->getSurface().getTol();
if (utils::DoubleFloatingEqual()(unsignedDistance, 0., distanceTolerance, distanceTolerance)) {
if (signedLevelSet) {
return m_surfaceInfo->evalNormal(point, supportItr);
} else {
return {{0., 0., 0.}};
}
}
gradient /= unsignedDistance;

// Evaluate levelset gradient
std::array<double,3> gradient = m_surfaceInfo->evalDistanceVector(point, supportItr) / distance;
if (signedLevelSet) {
std::array<double, 3> pseudoNormal = m_surfaceInfo->evalPseudoNormal(supportItr, lambda);
double normalGradient = dotProduct(gradient, pseudoNormal);
if (utils::DoubleFloatingEqual()(normalGradient, 0., distanceTolerance, distanceTolerance)) {
throw std::runtime_error("Unable to evaluate gradient: the point lies on the normal plane!");
}
gradient *= sign(normalGradient);
}

return gradient;
}
Expand Down
1 change: 1 addition & 0 deletions src/levelset/levelSetSegmentationObject.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@ class LevelSetSegmentationSurfaceInfo {
std::array<double, 3> evalDistanceVector(const std::array<double, 3> &point, const SegmentConstIterator &segmentItr) const;

std::array<double, 3> evalNormal(const std::array<double, 3> &point, const SegmentConstIterator &segmentItr) const;
std::array<double,3> evalPseudoNormal(const SurfUnstructured::CellConstIterator &segmentIterator, const double *lambda) const;

void evalBarycentricCoordinates(const std::array<double, 3> &point, const SegmentConstIterator &segmentItr, double *lambda) const;
void evalProjection(const std::array<double, 3> &point, const SegmentConstIterator &segmentItr, std::array<double, 3> *projectionPoint, std::array<double, 3> *projectionNormal) const;
Expand Down

0 comments on commit b3817e5

Please sign in to comment.