diff --git a/src/levelset/levelSetObject.cpp b/src/levelset/levelSetObject.cpp index a7e46d52c8..f0111b5273 100644 --- a/src/levelset/levelSetObject.cpp +++ b/src/levelset/levelSetObject.cpp @@ -903,6 +903,14 @@ bool LevelSetObject::isCellInNarrowBand(long id)const * * The value of the levelset is evaluated and compared with the specified narrow band size. * + * If a point is inside a cell that belongs to the narrow band because it is a neighbour of a + * cell with a different levelset sign, its function may identify the point as outside the + * narrow band. That's because it will only compare the levelset of the point with the narrow + * band size. The extreme case is when the narrow band size is set to zero and being a neighbour + * of a cell with a different levelset sign is the only criterion to identify cells inside the + * narrow band. In this situation, this function will identify as inside the narrow band only + * the points that lie on the levelset-zero iso surface. + * * \param point are the coordinates of the point * \result Return true if the cell is in the narrow band, false otherwise. */ @@ -914,7 +922,7 @@ bool LevelSetObject::isInNarrowBand(const std::array &point)const } // Early return if no narrow band was set - if (m_narrowBandSize != LEVELSET_NARROW_BAND_UNLIMITED) { + if (m_narrowBandSize == LEVELSET_NARROW_BAND_UNLIMITED) { return true; } @@ -1111,9 +1119,11 @@ void LevelSetObject::fillCellPropagatedSignCache() // Set the sign of the cell // // If the sign of the cell has already been set, we check if it matches - // the seed sign and then we skip the cell. - // - // If the cell is a seed, we remove it from the seed list. + // the seed sign and then we skip the cell, otherwise we set the sign of + // the cell and we process its neighoburs. Once the sign of a cell is + // set, we can remove it form the seed list (there is no need to start + // again the porpagation from that cell, beause it will not reach any + // new portions of the domain). if (cellId != seedId) { CellCacheCollection::ValueCache::Entry cellSignEntry = propagatedSignCache->findEntry(cellId); if (cellSignEntry.isValid()) { @@ -1125,10 +1135,9 @@ void LevelSetObject::fillCellPropagatedSignCache() } else { continue; } - - seedIds.erase(cellId); } else { propagatedSignCache->insertEntry(cellId, seedSign); + seedIds.erase(cellId); } } diff --git a/src/levelset/levelSetSegmentationObject.cpp b/src/levelset/levelSetSegmentationObject.cpp index d9977a85fd..8ddd25899c 100644 --- a/src/levelset/levelSetSegmentationObject.cpp +++ b/src/levelset/levelSetSegmentationObject.cpp @@ -183,14 +183,16 @@ double LevelSetSegmentationSurfaceInfo::getFeatureAngle() const { } /*! - * Evaluate the signed distance function at the specified point. + * 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 - * @return The signed distance function at the specified point. + * @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) const + const SegmentConstIterator &segmentItr, + bool signedDistance) const { // Project the point on the surface and evaluate the point-projection vector int nSegmentVertices = segmentItr->getVertexCount(); @@ -200,6 +202,9 @@ double LevelSetSegmentationSurfaceInfo::evalDistance(const std::array // Evaluate unsigned distance double unsignedDistance = norm2(pointProjectionVector); + if (!signedDistance) { + return unsignedDistance; + } // Signed distance // @@ -1762,7 +1767,7 @@ short LevelSetSegmentationObject::_evalSign(const std::array &point, l LevelSetSegmentationSurfaceInfo::SegmentConstIterator supportItr = getSurface().getCellConstIterator(support); - return evalValueSign(m_surfaceInfo->evalDistance(point, supportItr)); + return evalValueSign(m_surfaceInfo->evalDistance(point, supportItr, true)); } /*! @@ -1789,7 +1794,7 @@ double LevelSetSegmentationObject::_evalValue(const std::array &point, // Evaluate the distance of the point from the surface LevelSetSegmentationSurfaceInfo::SegmentConstIterator supportItr = getSurface().getCellConstIterator(support); - double distance = m_surfaceInfo->evalDistance(point, supportItr); + double distance = m_surfaceInfo->evalDistance(point, supportItr, signedLevelSet); // Early return if the point lies on the surface if (evalValueSign(distance) == 0) { @@ -1831,7 +1836,7 @@ 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); + double distance = m_surfaceInfo->evalDistance(point, supportItr, signedLevelSet); // Early return if the point lies on the surface if (evalValueSign(distance) == 0) { @@ -1843,12 +1848,7 @@ std::array LevelSetSegmentationObject::_evalGradient(const std::array< } // Evaluate levelset gradient - std::array gradient = m_surfaceInfo->evalDistanceVector(point, supportItr); - if (signedLevelSet) { - gradient /= distance; - } else { - gradient /= std::abs(distance); - } + std::array gradient = m_surfaceInfo->evalDistanceVector(point, supportItr) / distance; return gradient; } diff --git a/src/levelset/levelSetSegmentationObject.hpp b/src/levelset/levelSetSegmentationObject.hpp index 464e0281dc..2368083a79 100644 --- a/src/levelset/levelSetSegmentationObject.hpp +++ b/src/levelset/levelSetSegmentationObject.hpp @@ -70,7 +70,7 @@ class LevelSetSegmentationSurfaceInfo { double getFeatureAngle() const; - double evalDistance(const std::array &point, const SegmentConstIterator &segmentItr) 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;