ENH: surfaceFeatureExtract now builds with CGAL support for curvature

extendedFeatureEdgeMesh: Add function to find all feature points within a sphere
treeDataPoint: Add support for point overlap test

applications/utilities/surface/surfaceFeatureExtract/Allwmake

@@ -0,0 +1,33 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+set -x
+
+if [ ! -e "Make/files" ] || [ ! -e "Make/options" ]
+then
+    mkdir -p Make
+
+    if [ -n "$CGAL_ARCH_PATH" ]
+    then
+        cp -r MakeWithCGAL/* Make
+
+        echo
+        echo Compiling surfaceFeatureExtract with CGAL support for curvature
+        echo
+
+        wmake
+    else
+        cp -r MakeWithoutCGAL/* Make
+
+        echo
+        echo Compiling surfaceFeatureExtract without CGAL support for curvature
+        echo
+
+        wmake
+    fi
+else
+    echo surfaceFeatureExtract already has a Make folder
+fi
+
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H

@@ -0,0 +1,228 @@
+#ifndef CGAL_PSURF_RINGS_H_
+#define CGAL_PSURF_RINGS_H_
+
+// This file adapted from
+//     CGAL-3.7/examples/Jet_fitting_3//PolyhedralSurf_rings.h
+// Licensed under CGAL-3.7/LICENSE.FREE_USE
+
+// Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
+// Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
+// Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
+// Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
+// Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
+// (Israel).  All rights reserved.
+
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+
+#include <cassert>
+
+using namespace std;
+
+template
+<
+    class TPoly,
+    class VertexPropertyMap
+>
+
+class T_PolyhedralSurf_rings
+{
+
+protected:
+
+    // Polyhedron
+    typedef typename TPoly::Vertex Vertex;
+    typedef typename TPoly::Halfedge Halfedge;
+    typedef typename TPoly::Facet Facet;
+    typedef typename TPoly::Halfedge_around_vertex_circulator
+        Halfedge_around_vertex_circulator;
+    typedef typename TPoly::Vertex_iterator Vertex_iterator;
+
+    // vertex indices are initialised to -1
+    static void reset_ring_indices
+    (
+        std::vector< Vertex* >& vces,
+        VertexPropertyMap& vpm
+    );
+
+    // i >= 1; from a start vertex on the current i-1 ring, push non-visited
+    // neighbors of start in the nextRing and set indices to i. Also add these
+    // vertices in all.
+    static void push_neighbours_of
+    (
+        Vertex* start,
+        int ith,
+        std::vector< Vertex* >& nextRing,
+        std::vector< Vertex* >& all,
+        VertexPropertyMap& vpm
+    );
+
+    // i >= 1, from a currentRing i-1, collect all neighbors, set indices
+    // to i and store them in nextRing and all.
+    static void collect_ith_ring
+    (
+        int ith,
+        std::vector< Vertex* >& currentRing,
+        std::vector< Vertex* >& nextRing,
+        std::vector< Vertex* >& all,
+        VertexPropertyMap& vpm
+    );
+
+
+public:
+
+    // collect i>=1 rings : all neighbours up to the ith ring,
+    static void collect_i_rings
+    (
+        Vertex* v,
+        int ring_i,
+        std::vector< Vertex* >& all,
+        VertexPropertyMap& vpm
+    );
+
+    //collect enough rings (at least 1), to get at least min_nb of neighbors
+    static void collect_enough_rings
+    (
+        Vertex* v,
+        unsigned int min_nb,
+        std::vector< Vertex* >& all,
+        VertexPropertyMap& vpm
+    );
+};
+
+
+////IMPLEMENTATION////////////////////////////////////////////////////
+
+template < class TPoly , class VertexPropertyMap>
+void T_PolyhedralSurf_rings <TPoly, VertexPropertyMap>::
+push_neighbours_of(Vertex* start, int ith,
+                   std::vector< Vertex* >& nextRing,
+                   std::vector< Vertex* >& all,
+                   VertexPropertyMap& vpm)
+{
+    Vertex *v;
+    Halfedge_around_vertex_circulator
+        hedgeb = start->vertex_begin(), hedgee = hedgeb;
+
+    CGAL_For_all(hedgeb, hedgee)
+    {
+        v = &*(hedgeb->opposite()->vertex());
+        if (get(vpm, v) != -1) continue; //if visited: next
+
+        put(vpm, v, ith);
+        nextRing.push_back(v);
+        all.push_back(v);
+    }
+}
+
+
+template <class TPoly, class VertexPropertyMap>
+void T_PolyhedralSurf_rings <TPoly, VertexPropertyMap>::
+collect_ith_ring(int ith, std::vector< Vertex* >& currentRing,
+                 std::vector< Vertex* >& nextRing,
+                 std::vector< Vertex* >& all,
+                 VertexPropertyMap& vpm)
+{
+    typename std::vector< Vertex* >::iterator
+        itb = currentRing.begin(), ite = currentRing.end();
+
+    CGAL_For_all(itb, ite)
+    {
+        push_neighbours_of(*itb, ith, nextRing, all, vpm);
+    }
+}
+
+
+template <class TPoly, class VertexPropertyMap>
+void T_PolyhedralSurf_rings <TPoly, VertexPropertyMap>::
+reset_ring_indices(std::vector< Vertex* >& vces,
+                   VertexPropertyMap& vpm)
+{
+    typename std::vector< Vertex* >::iterator
+        itb = vces.begin(), ite = vces.end();
+
+    CGAL_For_all(itb, ite)
+    {
+        put(vpm, *itb, -1);
+    }
+}
+
+
+template <class TPoly, class VertexPropertyMap>
+void T_PolyhedralSurf_rings <TPoly, VertexPropertyMap>::
+collect_i_rings(Vertex* v,
+                int ring_i,
+                std::vector< Vertex* >& all,
+                VertexPropertyMap& vpm)
+{
+    std::vector<Vertex*> current_ring, next_ring;
+    std::vector<Vertex*> *p_current_ring, *p_next_ring;
+    assert(ring_i >= 1);
+
+    //initialize
+    p_current_ring = &current_ring;
+    p_next_ring = &next_ring;
+    put(vpm, v, 0);
+    current_ring.push_back(v);
+    all.push_back(v);
+
+    for (int i=1; i<=ring_i; i++)
+    {
+      collect_ith_ring(i, *p_current_ring, *p_next_ring, all, vpm);
+      //next round must be launched from p_nextRing...
+      p_current_ring->clear();
+      std::swap(p_current_ring, p_next_ring);
+    }
+
+    //clean up
+    reset_ring_indices(all, vpm);
+}
+
+
+template <class TPoly, class VertexPropertyMap>
+void T_PolyhedralSurf_rings <TPoly, VertexPropertyMap>::
+collect_enough_rings(Vertex* v,
+                     unsigned int min_nb,
+                     std::vector< Vertex* >& all,
+                     VertexPropertyMap& vpm)
+{
+    std::vector<Vertex*> current_ring, next_ring;
+    std::vector<Vertex*> *p_current_ring, *p_next_ring;
+
+    //initialize
+    p_current_ring = &current_ring;
+    p_next_ring = &next_ring;
+    put(vpm, v, 0);
+    current_ring.push_back(v);
+    all.push_back(v);
+
+    int i = 1;
+
+    while ( (all.size() < min_nb) &&  (p_current_ring->size() != 0) )
+    {
+        collect_ith_ring(i, *p_current_ring, *p_next_ring, all, vpm);
+        //next round must be launched from p_nextRing...
+        p_current_ring->clear();
+        std::swap(p_current_ring, p_next_ring);
+        i++;
+    }
+
+    //clean up
+    reset_ring_indices(all, vpm);
+}
+
+
+#endif

applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "buildCGALPolyhedron.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::buildCGALPolyhedron::buildCGALPolyhedron
+(
+    const Foam::triSurface& surf
+)
+:
+    CGAL::Modifier_base<HalfedgeDS>(),
+    surf_(surf)
+{}
+
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::buildCGALPolyhedron::~buildCGALPolyhedron()
+{}
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+void Foam::buildCGALPolyhedron::operator()
+(
+    HalfedgeDS& hds
+)
+{
+    typedef HalfedgeDS::Traits     Traits;
+    typedef Traits::Point_3        Point;
+
+    // Postcondition: `hds' is a valid polyhedral surface.
+    CGAL::Polyhedron_incremental_builder_3<HalfedgeDS> B(hds, false);
+
+    B.begin_surface
+    (
+        surf_.points().size(),      // n points
+        surf_.size(),               // n facets
+        2*surf_.edges().size()      // n halfedges
+    );
+
+    forAll(surf_.points(), pI)
+    {
+        const Foam::point& p = surf_.points()[pI];
+
+        B.add_vertex(Point(p.x(), p.y(), p.z()));
+    }
+
+    forAll(surf_, fI)
+    {
+        B.begin_facet();
+
+        B.add_vertex_to_facet(surf_[fI][0]);
+        B.add_vertex_to_facet(surf_[fI][1]);
+        B.add_vertex_to_facet(surf_[fI][2]);
+
+        B.end_facet();
+    }
+
+    B.end_surface();
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::buildCGALPolyhedron
+
+Description
+    Convert a triSurface into a CGAL Polyhedron
+
+SourceFiles
+    buildCGALPolyhedron.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef buildCGALPolyhedron_H
+#define buildCGALPolyhedron_H
+
+#include "triSurface.H"
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Polyhedron_incremental_builder_3.h>
+#include <CGAL/Polyhedron_3.h>
+
+typedef CGAL::Simple_cartesian<double>               Kernel;
+typedef CGAL::Polyhedron_3<Kernel>                   Polyhedron;
+typedef Polyhedron::HalfedgeDS                       HalfedgeDS;
+typedef Polyhedron::Vertex                           Vertex;
+typedef Polyhedron::Vertex_iterator                  Vertex_iterator;
+typedef Kernel::Point_3                              Point_3;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class buildCGALPolyhedron Declaration
+\*---------------------------------------------------------------------------*/
+
+class buildCGALPolyhedron
+:
+    public CGAL::Modifier_base<HalfedgeDS>
+{
+    // Private data
+
+        //- Reference to triSurface to convert
+        const Foam::triSurface& surf_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        buildCGALPolyhedron(const buildCGALPolyhedron&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const buildCGALPolyhedron&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct with reference to triSurface
+        buildCGALPolyhedron(const triSurface& surf);
+
+
+    //- Destructor
+    ~buildCGALPolyhedron();
+
+
+    // Member Operators
+
+        //- operator() of this `modifier' called by delegate function of
+        //  Polyhedron
+        void operator()(HalfedgeDS& hds);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files

@@ -0,0 +1,4 @@
+surfaceFeatureExtract.C
+CGALPolyhedron/buildCGALPolyhedron.C
+
+EXE = $(FOAM_APPBIN)/surfaceFeatureExtract

applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options

@@ -0,0 +1,29 @@
+EXE_FROUNDING_MATH = -frounding-math
+EXE_NDEBUG = -DNDEBUG
+USE_F2C = -DCGAL_USE_F2C
+include $(GENERAL_RULES)/CGAL
+
+EXE_INC = \
+    -DENABLE_CURVATURE \
+    ${EXE_FROUNDING_MATH} \
+    ${EXE_NDEBUG} \
+    ${USE_F2C} \
+    ${CGAL_INC} \
+    -ICGALPolyhedron \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    $(CGAL_LIBS) \
+    -L$(CGAL_ARCH_PATH)/lib \
+    -llapack \
+    -lblas \
+    -lCGAL \
+    -lmeshTools \
+    -ledgeMesh \
+    -ltriSurface \
+    -lsampling

applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files

@@ -0,0 +1,3 @@
+surfaceFeatureExtract.C
+
+EXE = $(FOAM_APPBIN)/surfaceFeatureExtract

applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options

@@ -0,0 +1,13 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -ledgeMesh \
+    -ltriSurface \
+    -lsampling

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,10 +46,263 @@ Description
 #include "unitConversion.H"
 #include "plane.H"
 
+#ifdef ENABLE_CURVATURE
+#include "buildCGALPolyhedron.H"
+#include "CGALPolyhedronRings.H"
+#include <CGAL/Monge_via_jet_fitting.h>
+#include <CGAL/Lapack/Linear_algebra_lapack.h>
+#include <CGAL/property_map.h>
+#endif
+
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef ENABLE_CURVATURE
+scalarField calcCurvature(const triSurface& surf)
+{
+    scalarField k(surf.points().size(), 0);
+
+    Polyhedron P;
+
+    buildCGALPolyhedron convert(surf);
+    P.delegate(convert);
+
+    // Info<< "Created CGAL Polyhedron with " << label(P.size_of_vertices())
+    //     << " vertices and " << label(P.size_of_facets())
+    //     << " facets. " << endl;
+
+    // The rest of this function adapted from
+    //     CGAL-3.7/examples/Jet_fitting_3/Mesh_estimation.cpp
+    // Licensed under CGAL-3.7/LICENSE.FREE_USE
+
+    // Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
+    // Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
+    // Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
+    // Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
+    // Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
+    // (Israel).  All rights reserved.
+
+    // Permission is hereby granted, free of charge, to any person obtaining a
+    // copy of this software and associated documentation files (the
+    // "Software"), to deal in the Software without restriction, including
+    // without limitation the rights to use, copy, modify, merge, publish,
+    // distribute, sublicense, and/or sell copies of the Software, and to permit
+    // persons to whom the Software is furnished to do so, subject to the
+    // following conditions:
+
+    // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+    // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+    // IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+    // CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT
+    // OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
+    // THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+     //Vertex property map, with std::map
+    typedef std::map<Vertex*, int> Vertex2int_map_type;
+    typedef boost::associative_property_map< Vertex2int_map_type >
+        Vertex_PM_type;
+    typedef T_PolyhedralSurf_rings<Polyhedron, Vertex_PM_type > Poly_rings;
+
+    typedef CGAL::Monge_via_jet_fitting<Kernel>         Monge_via_jet_fitting;
+    typedef Monge_via_jet_fitting::Monge_form           Monge_form;
+
+    std::vector<Point_3> in_points;  //container for data points
+
+    // default parameter values and global variables
+    unsigned int d_fitting = 2;
+    unsigned int d_monge = 2;
+    unsigned int min_nb_points = (d_fitting + 1)*(d_fitting + 2)/2;
+
+    //initialize the tag of all vertices to -1
+    Vertex_iterator vitb = P.vertices_begin();
+    Vertex_iterator vite = P.vertices_end();
+
+    Vertex2int_map_type vertex2props;
+    Vertex_PM_type vpm(vertex2props);
+
+    CGAL_For_all(vitb, vite)
+    {
+        put(vpm, &(*vitb), -1);
+    }
+
+    vite = P.vertices_end();
+
+    label vertI = 0;
+
+    for (vitb = P.vertices_begin(); vitb != vite; vitb++)
+    {
+        //initialize
+        Vertex* v = &(*vitb);
+
+        //gather points around the vertex using rings
+        // From: gather_fitting_points(v, in_points, vpm);
+        {
+            std::vector<Vertex*> gathered;
+            in_points.clear();
+
+            Poly_rings::collect_enough_rings(v, min_nb_points, gathered, vpm);
+
+            //store the gathered points
+            std::vector<Vertex*>::iterator itb = gathered.begin();
+            std::vector<Vertex*>::iterator ite = gathered.end();
+
+            CGAL_For_all(itb, ite)
+            {
+                in_points.push_back((*itb)->point());
+            }
+        }
+
+        //skip if the nb of points is to small
+        if ( in_points.size() < min_nb_points )
+        {
+            std::cerr
+                << "not enough pts for fitting this vertex"
+                << in_points.size()
+                << std::endl;
+
+            continue;
+        }
+
+        // perform the fitting
+        Monge_via_jet_fitting monge_fit;
+
+        Monge_form monge_form = monge_fit
+        (
+            in_points.begin(),
+            in_points.end(),
+            d_fitting,
+            d_monge
+        );
+
+//        std::cout<< monge_form;;
+//        std::cout<< "condition number : "
+//                 << monge_fit.condition_number() << nl << std::endl;
+
+        // Use the maximum curvature to give smaller cell sizes later.
+        k[vertI++]
+            = max
+              (
+                   mag(monge_form.principal_curvatures(0)),
+                   mag(monge_form.principal_curvatures(1))
+              );
+    }
+
+    return k;
+}
+#endif
+
+
+bool edgesConnected(const edge& e1, const edge& e2)
+{
+    if
+    (
+        e1.start() == e2.start() || e1.start() == e2.end()
+     || e1.end() == e2.start() || e1.end() == e2.end()
+    )
+    {
+        return true;
+    }
+
+    return false;
+}
+
+
+scalar calcProximityOfFeaturePoints
+(
+    const List<pointIndexHit>& hitList,
+    const scalar defaultCellSize
+)
+{
+    scalar minDist = defaultCellSize;
+
+    for
+    (
+        label hI1 = 0;
+        hI1 < hitList.size() - 1;
+        ++hI1
+    )
+    {
+        const pointIndexHit& pHit1 = hitList[hI1];
+
+        if (pHit1.hit())
+        {
+            for
+            (
+                label hI2 = hI1 + 1;
+                hI2 < hitList.size();
+                ++hI2
+            )
+            {
+                const pointIndexHit& pHit2 = hitList[hI2];
+
+                if (pHit2.hit())
+                {
+                    scalar curDist = mag(pHit1.hitPoint() - pHit2.hitPoint());
+
+                    minDist = min(curDist, minDist);
+                }
+            }
+        }
+    }
+
+    return minDist;
+}
+
+
+scalar calcProximityOfFeatureEdges
+(
+    const extendedFeatureEdgeMesh& efem,
+    const List<pointIndexHit>& hitList,
+    const scalar defaultCellSize
+)
+{
+    scalar minDist = defaultCellSize;
+
+    for
+    (
+        label hI1 = 0;
+        hI1 < hitList.size() - 1;
+        ++hI1
+    )
+    {
+        const pointIndexHit& pHit1 = hitList[hI1];
+
+        if (pHit1.hit())
+        {
+            const edge& e1 = efem.edges()[pHit1.index()];
+
+            for
+            (
+                label hI2 = hI1 + 1;
+                hI2 < hitList.size();
+                ++hI2
+            )
+            {
+                const pointIndexHit& pHit2 = hitList[hI2];
+
+                if (pHit2.hit())
+                {
+                    const edge& e2 = efem.edges()[pHit2.index()];
+
+                    // Don't refine if the edges are connected to each other
+                    if (!edgesConnected(e1, e2))
+                    {
+                        scalar curDist
+                            = mag(pHit1.hitPoint() - pHit2.hitPoint());
+
+                        minDist = min(curDist, minDist);
+                    }
+                }
+            }
+        }
+    }
+
+    return minDist;
+}
+
+
 void dumpBox(const treeBoundBox& bb, const fileName& fName)
 {
     OFstream str(fName);
@@ -299,10 +552,9 @@ int main(int argc, char *argv[])
         "writeVTK",
         "write extendedFeatureEdgeMesh vtk files"
     );
-    argList::addOption
+    argList::addBoolOption
     (
         "closeness",
-        "scalar",
         "span to look for surface closeness"
     );
     argList::addOption
@@ -331,7 +583,7 @@ int main(int argc, char *argv[])
 #   ifdef ENABLE_CURVATURE
     argList::addBoolOption
     (
-        "calcCurvature",
+        "curvature",
         "calculate curvature and closeness fields"
     );
 #   endif
@@ -527,6 +779,7 @@ int main(int argc, char *argv[])
             << " will be included as feature edges."<< endl;
     }
 
+
     surfaceFeatures newSet(surf);
     newSet.setFromStatus(edgeStat);
 
@@ -564,7 +817,6 @@ int main(int argc, char *argv[])
 
     feMesh.write();
 
-
     // Write a featureEdgeMesh for backwards compatibility
     {
         featureEdgeMesh bfeMesh
@@ -595,7 +847,7 @@ int main(int argc, char *argv[])
         (
             sFeatFileName + ".closeness",
             runTime.constant(),
-            "extendedFeatureEdgeMesh",
+            "triSurface",
             runTime,
             IOobject::NO_READ,
             IOobject::NO_WRITE
@@ -603,13 +855,6 @@ int main(int argc, char *argv[])
         surf
     );
 
-    if (!curvature)
-    {
-        Info<< "End\n" << endl;
-
-        return 0;
-    }
-
     // Find close features
 
     // // Dummy trim operation to mark features
@@ -677,8 +922,10 @@ int main(int argc, char *argv[])
     //     )
     // );
 
-    Info<< "Examine curvature, feature proximity and internal and "
+    if (args.optionFound("closeness"))
-        << "external closeness." << endl;
+    {
+        Info<< nl << "Extracting internal and external closeness of surface."
+            << endl;
 
         // Internal and external closeness
 
@@ -688,7 +935,7 @@ int main(int argc, char *argv[])
 
         scalar span = searchSurf.bounds().mag();
 
-    args.optionReadIfPresent("closeness", span);
+        //args.optionReadIfPresent("closeness", span);
 
         scalar externalAngleTolerance = 10;
         scalar externalToleranceCosAngle = Foam::cos
@@ -842,7 +1089,7 @@ int main(int argc, char *argv[])
             (
                 sFeatFileName + ".internalCloseness",
                 runTime.constant(),
-            "extendedFeatureEdgeMesh",
+                "triSurface",
                 runTime,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
@@ -860,7 +1107,7 @@ int main(int argc, char *argv[])
             (
                 sFeatFileName + ".externalCloseness",
                 runTime.constant(),
-            "extendedFeatureEdgeMesh",
+                "triSurface",
                 runTime,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
@@ -872,39 +1119,6 @@ int main(int argc, char *argv[])
 
         externalClosenessField.write();
 
-
-#ifdef ENABLE_CURVATURE
-    scalarField k = calcCurvature(surf);
-
-    // Modify the curvature values on feature edges and points to be zero.
-
-    forAll(newSet.featureEdges(), fEI)
-    {
-        const edge& e = surf.edges()[newSet.featureEdges()[fEI]];
-
-        k[surf.meshPoints()[e.start()]] = 0.0;
-        k[surf.meshPoints()[e.end()]] = 0.0;
-    }
-
-    triSurfacePointScalarField kField
-    (
-        IOobject
-        (
-            sFeatFileName + ".curvature",
-            runTime.constant(),
-            "extendedFeatureEdgeMesh",
-            runTime,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        surf,
-        dimLength,
-        k
-    );
-
-    kField.write();
-#endif
-
         if (writeVTK)
         {
             vtkSurfaceWriter().write
@@ -930,8 +1144,47 @@ int main(int argc, char *argv[])
                 false,                              // isNodeValues
                 true                                // verbose
             );
+        }
+    }
+
 
 #ifdef ENABLE_CURVATURE
+    if (args.optionFound("curvature"))
+    {
+        Info<< nl << "Extracting curvature of surface at the points." << endl;
+
+        scalarField k = calcCurvature(surf);
+
+        // Modify the curvature values on feature edges and points to be zero.
+
+    //    forAll(newSet.featureEdges(), fEI)
+    //    {
+    //        const edge& e = surf.edges()[newSet.featureEdges()[fEI]];
+    //
+    //        k[surf.meshPoints()[e.start()]] = 0.0;
+    //        k[surf.meshPoints()[e.end()]] = 0.0;
+    //    }
+
+        triSurfacePointScalarField kField
+        (
+            IOobject
+            (
+                sFeatFileName + ".curvature",
+                runTime.constant(),
+                "triSurface",
+                runTime,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            surf,
+            dimLength,
+            k
+        );
+
+        kField.write();
+
+        if (writeVTK)
+        {
             vtkSurfaceWriter().write
             (
                 runTime.constant()/"triSurface",    // outputDir
@@ -943,7 +1196,82 @@ int main(int argc, char *argv[])
                 true,                               // isNodeValues
                 true                                // verbose
             );
+        }
+    }
 #endif
+
+
+    if (args.optionFound("featureProximity"))
+    {
+        Info<< nl << "Extracting proximity of close feature points and edges "
+            << "to the surface" << endl;
+
+        const scalar searchDistance =
+            args.optionRead<scalar>("featureProximity");
+
+        const scalar radiusSqr = sqr(searchDistance);
+
+        scalarField featureProximity(surf.size(), searchDistance);
+
+        forAll(surf, fI)
+        {
+            const triPointRef& tri = surf[fI].tri(surf.points());
+            const point& triCentre = tri.circumCentre();
+
+            List<pointIndexHit> hitList;
+
+            feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList);
+
+            featureProximity[fI] =
+                calcProximityOfFeatureEdges
+                (
+                    feMesh,
+                    hitList,
+                    featureProximity[fI]
+                );
+
+            feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList);
+
+            featureProximity[fI] =
+                calcProximityOfFeaturePoints
+                (
+                    hitList,
+                    featureProximity[fI]
+                );
+        }
+
+        triSurfaceScalarField featureProximityField
+        (
+            IOobject
+            (
+                sFeatFileName + ".featureProximity",
+                runTime.constant(),
+                "triSurface",
+                runTime,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            surf,
+            dimLength,
+            featureProximity
+        );
+
+        featureProximityField.write();
+
+        if (writeVTK)
+        {
+            vtkSurfaceWriter().write
+            (
+                runTime.constant()/"triSurface",    // outputDir
+                sFeatFileName,                      // surfaceName
+                surf.points(),
+                faces,
+                "featureProximity",                // fieldName
+                featureProximity,
+                false,                              // isNodeValues
+                true                                // verbose
+            );
+        }
     }
 
     Info<< "End\n" << endl;

src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C

@@ -819,6 +819,39 @@ void Foam::extendedFeatureEdgeMesh::nearestFeatureEdgeByType
 }
 
 
+void Foam::extendedFeatureEdgeMesh::allNearestFeaturePoints
+(
+    const point& sample,
+    scalar searchRadiusSqr,
+    List<pointIndexHit>& info
+) const
+{
+    DynamicList<pointIndexHit> dynPointHit;
+
+    // Pick up all the feature points that intersect the search sphere
+    labelList elems = pointTree().findSphere
+    (
+        sample,
+        searchRadiusSqr
+    );
+
+    forAll(elems, elemI)
+    {
+        label index = elems[elemI];
+        label ptI = pointTree().shapes().pointLabels()[index];
+        const point& pt = points()[ptI];
+
+        pointIndexHit nearHit;
+
+        nearHit = pointIndexHit(true, pt, index);
+
+        dynPointHit.append(nearHit);
+    }
+
+    info.transfer(dynPointHit);
+}
+
+
 void Foam::extendedFeatureEdgeMesh::allNearestFeatureEdges
 (
     const point& sample,

src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H

@@ -283,6 +283,14 @@ public:
                 List<pointIndexHit>& info
             ) const;
 
+            //- Find all the feature points within searchDistSqr of sample
+            void allNearestFeaturePoints
+            (
+                const point& sample,
+                scalar searchRadiusSqr,
+                List<pointIndexHit>& info
+            ) const;
+
             //- Find all the feature edges within searchDistSqr of sample
             void allNearestFeatureEdges
             (
@@ -291,6 +299,7 @@ public:
                 List<pointIndexHit>& info
             ) const;
 
+
         // Access
 
             //- Return the index of the start of the convex feature points

src/meshTools/indexedOctree/treeDataPoint.C

@@ -93,6 +93,25 @@ bool Foam::treeDataPoint::overlaps
 }
 
 
+// Check if any point on shape is inside sphere.
+bool Foam::treeDataPoint::overlaps
+(
+    const label index,
+    const point& centre,
+    const scalar radiusSqr
+) const
+{
+    label pointI = (useSubset_ ? pointLabels_[index] : index);
+
+    if (magSqr(points_[pointI] - centre) <= radiusSqr)
+    {
+        return true;
+    }
+
+    return false;
+}
+
+
 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
 // nearestPoint.
 void Foam::treeDataPoint::findNearest

src/meshTools/indexedOctree/treeDataPoint.H

@@ -128,6 +128,14 @@ public:
                 const treeBoundBox& sampleBb
             ) const;
 
+            //- Does shape at index overlap the sphere
+            bool overlaps
+            (
+                const label index,
+                const point& centre,
+                const scalar radiusSqr
+            ) const;
+
             //- Calculates nearest (to sample) point in shape.
             //  Returns actual point and distance (squared)
             void findNearest