diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index e631083aa0..8bec6c4ae8 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -316,6 +316,8 @@ int main(int argc, char *argv[])
     Info<< "Feature line extraction is only valid on closed manifold surfaces."
         << endl;
 
+    bool writeVTK = args.optionFound("writeVTK");
+    bool writeObj = args.optionFound("writeObj");
 
     const fileName surfFileName = args[1];
     const fileName outFileName  = args[2];
@@ -324,6 +326,8 @@ int main(int argc, char *argv[])
         << "Output feature set : " << outFileName << nl
         << endl;
 
+    fileName sFeatFileName = surfFileName.lessExt().name();
+
 
     // Read
     // ~~~~
@@ -334,6 +338,13 @@ int main(int argc, char *argv[])
     surf.writeStats(Info);
     Info<< endl;
 
+    faceList faces(surf.size());
+
+    forAll(surf, fI)
+    {
+        faces[fI] = surf[fI].triFaceFace();
+    }
+
 
     // Either construct features from surface&featureangle or read set.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -467,13 +478,42 @@ 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."
         << endl;
 
-    fileName sFeatFileName = surfFileName.lessExt().name();
-
     featureEdgeMesh feMesh
     (
         IOobject
@@ -490,7 +530,7 @@ int main(int argc, char *argv[])
 
     feMesh.write();
 
-    if (args.optionFound("writeObj"))
+    if (writeObj)
     {
         feMesh.writeObj(runTime.constant()/"featureEdgeMesh"/sFeatFileName);
     };
@@ -517,13 +557,6 @@ int main(int argc, char *argv[])
 
     const vectorField& normals = searchSurf.faceNormals();
 
-    faceList faces(surf.size());
-
-    forAll(surf, fI)
-    {
-        faces[fI] = surf[fI].triFaceFace();
-    }
-
     scalar span = searchSurf.bounds().mag();
 
     args.optionReadIfPresent("closeness", span);
@@ -794,7 +827,7 @@ int main(int argc, char *argv[])
 
     kField.write();
 
-    if (args.optionFound("writeVTK"))
+    if (writeVTK)
     {
         vtkSurfaceWriter<scalar>().write
         (
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
index 050b183737..9a2eb5e770 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -546,7 +546,7 @@ void Foam::surfaceFeatures::findFeatures(const scalar includedAngle)
 
 // Remove small strings of edges. First assigns string number to
 // every edge and then works out the length of them.
-void Foam::surfaceFeatures::trimFeatures
+Foam::labelList Foam::surfaceFeatures::trimFeatures
 (
     const scalar minLen,
     const label minElems
@@ -679,6 +679,8 @@ void Foam::surfaceFeatures::trimFeatures
 
     // Convert back to edge labels
     setFromStatus(edgeStat);
+
+    return featLines;
 }
 
 
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
index a086752175..6f668ff76b 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -257,7 +257,7 @@ public:
 
             //- Delete small sets of edges. Edges are stringed up and any
             //  string of length < minLen (or nElems < minElems) is deleted.
-            void trimFeatures(const scalar minLen, const label minElems);
+            labelList trimFeatures(const scalar minLen, const label minElems);
 
             //- From member feature edges to status per edge.
             List<edgeStatus> toStatus() const;