diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 756128d12e..af2f6d645b 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -535,37 +535,6 @@ int main(int argc, char *argv[])
         << "        internal edges : " << newSet.nInternalEdges() << nl
         << endl;
 
-    // Dummy trim operation to mark features
-    labelList featureEdgeIndexing = newSet.trimFeatures(-GREAT, 0);
-
-    scalarField surfacePtFeatureIndex(surf.points().size(), -1);
-
-    forAll(newSet.featureEdges(), eI)
-    {
-        const edge& e = surf.edges()[newSet.featureEdges()[eI]];
-
-        surfacePtFeatureIndex[surf.meshPoints()[e.start()]] =
-            featureEdgeIndexing[newSet.featureEdges()[eI]];
-
-        surfacePtFeatureIndex[surf.meshPoints()[e.end()]] =
-            featureEdgeIndexing[newSet.featureEdges()[eI]];
-    }
-
-    if (writeVTK)
-    {
-        vtkSurfaceWriter<scalar>().write
-        (
-            runTime.constant()/"triSurface",    // outputDir
-            sFeatFileName,                      // surfaceName
-            surf.points(),
-            faces,
-            "surfacePtFeatureIndex",            // fieldName
-            surfacePtFeatureIndex,
-            true,                               // isNodeValues
-            true                                // verbose
-        );
-    }
-
     // Extracting and writing a featureEdgeMesh
 
     Pout<< nl << "Writing featureEdgeMesh to constant/featureEdgeMesh."
@@ -608,30 +577,61 @@ int main(int argc, char *argv[])
 
     // Find close features
 
-    Random rndGen(343267);
+    // // Dummy trim operation to mark features
+    // labelList featureEdgeIndexing = newSet.trimFeatures(-GREAT, 0);
 
-    treeBoundBox surfBB
-    (
-        treeBoundBox(searchSurf.bounds()).extend(rndGen, 1e-4)
-    );
+    // scalarField surfacePtFeatureIndex(surf.points().size(), -1);
 
-    surfBB.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-    surfBB.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    // forAll(newSet.featureEdges(), eI)
+    // {
+    //     const edge& e = surf.edges()[newSet.featureEdges()[eI]];
 
-    indexedOctree<treeDataEdge> ftEdTree
-    (
-        treeDataEdge
-        (
-            false,
-            surf.edges(),
-            surf.localPoints(),
-            newSet.featureEdges()
-        ),
-        surfBB,
-        8,      // maxLevel
-        10,     // leafsize
-        3.0     // duplicity
-    );
+    //     surfacePtFeatureIndex[surf.meshPoints()[e.start()]] =
+    //     featureEdgeIndexing[newSet.featureEdges()[eI]];
+
+    //     surfacePtFeatureIndex[surf.meshPoints()[e.end()]] =
+    //     featureEdgeIndexing[newSet.featureEdges()[eI]];
+    // }
+
+    // if (writeVTK)
+    // {
+    //     vtkSurfaceWriter<scalar>().write
+    //     (
+    //         runTime.constant()/"triSurface",    // outputDir
+    //         sFeatFileName,                      // surfaceName
+    //         surf.points(),
+    //         faces,
+    //         "surfacePtFeatureIndex",            // fieldName
+    //         surfacePtFeatureIndex,
+    //         true,                               // isNodeValues
+    //         true                                // verbose
+    //     );
+    // }
+
+    // Random rndGen(343267);
+
+    // treeBoundBox surfBB
+    // (
+    //     treeBoundBox(searchSurf.bounds()).extend(rndGen, 1e-4)
+    // );
+
+    // surfBB.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    // surfBB.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+    // indexedOctree<treeDataEdge> ftEdTree
+    // (
+    //     treeDataEdge
+    //     (
+    //         false,
+    //         surf.edges(),
+    //         surf.localPoints(),
+    //         newSet.featureEdges()
+    //     ),
+    //     surfBB,
+    //     8,      // maxLevel
+    //     10,     // leafsize
+    //     3.0     // duplicity
+    // );
 
     // labelList nearPoints = ftEdTree.findBox
     // (
@@ -905,6 +905,8 @@ int main(int argc, char *argv[])
         );
     }
 
+    Info<< nl << "End" << endl;
+
     return 0;
 }
 
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
index 24fcd501f7..011b9c19a3 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
@@ -263,6 +263,16 @@ Foam::surfaceFeatures::labelScalar Foam::surfaceFeatures::walkSegment
 
     label vertI = startPointI;
 
+    scalar visitedLength = 0.0;
+
+    label nVisited = 0;
+
+    if (findIndex(featurePoints_, startPointI) >= 0)
+    {
+        // Do not walk across feature points
+
+        return labelScalar(nVisited, visitedLength);
+    }
 
     //
     // Now we have:
@@ -284,11 +294,6 @@ Foam::surfaceFeatures::labelScalar Foam::surfaceFeatures::walkSegment
         unsetVal = currentFeatI;
     }
 
-
-    scalar visitedLength = 0.0;
-
-    label nVisited = 0;
-
     do
     {
         // Step to next feature edge with value unsetVal