levelset: miscellaneous fixes

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<double,3> &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<char>::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);
                     }
                 }
 

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<double, 3> &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<double, 3>
 
     // Evaluate unsigned distance
     double unsignedDistance = norm2(pointProjectionVector);
+    if (!signedDistance) {
+        return unsignedDistance;
+    }
 
     // Signed distance
     //
@@ -1762,7 +1767,7 @@ short LevelSetSegmentationObject::_evalSign(const std::array<double,3> &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<double,3> &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<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);
+    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<double,3> LevelSetSegmentationObject::_evalGradient(const std::array<
     }
 
     // Evaluate levelset gradient
-    std::array<double,3> gradient = m_surfaceInfo->evalDistanceVector(point, supportItr);
-    if (signedLevelSet) {
-        gradient /= distance;
-    } else {
-        gradient /= std::abs(distance);
-    }
+    std::array<double,3> gradient = m_surfaceInfo->evalDistanceVector(point, supportItr) / distance;
 
     return gradient;
 }

src/levelset/levelSetSegmentationObject.hpp

@@ -70,7 +70,7 @@ class LevelSetSegmentationSurfaceInfo {
 
     double getFeatureAngle() const;
 
-    double evalDistance(const std::array<double, 3> &point, const SegmentConstIterator &segmentItr) const;
+    double evalDistance(const std::array<double, 3> &point, const SegmentConstIterator &segmentItr, bool signedDistance) const;
     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;