ENH: Extract curvature and internal and external closeness fields.

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

@@ -0,0 +1,203 @@
+#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) 2011-2011 OpenCFD Ltd.
+     \\/     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 typename HalfedgeDS::Traits     Traits;
+    typedef typename 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) 2011-2011 OpenCFD Ltd.
+     \\/     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/Make/files

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

applications/utilities/surface/surfaceFeatureExtract/Make/options

@@ -1,9 +1,25 @@
+EXE_FROUNDING_MATH = -frounding-math
+USE_F2C = -DCGAL_USE_F2C
+
+include $(GENERAL_RULES)/CGAL
+
 EXE_INC = \
+    ${EXE_FROUNDING_MATH} \
+    ${USE_F2C} \
+    ${CGAL_INC} \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/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
+    -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) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,6 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-
 #include "triangle.H"
 #include "triSurface.H"
 #include "argList.H"
@@ -39,6 +38,15 @@ Description
 #include "treeBoundBox.H"
 #include "meshTools.H"
 #include "OFstream.H"
+#include "triSurfaceMesh.H"
+#include "vtkSurfaceWriter.H"
+#include "triSurfaceFields.H"
+
+#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>
 
 using namespace Foam;
 
@@ -89,6 +97,154 @@ void deleteBox
 }
 
 
+scalarField curvature(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::endl;
+
+        // std::cout<< "k1 " << monge_form.principal_curvatures(0) << std::endl;
+        // std::cout<< "k2 " << monge_form.principal_curvatures(1) << std::endl;
+
+        // std::vector<Point_3>::iterator itbp = in_points.begin();
+        // std::vector<Point_3>::iterator itep = in_points.end();
+
+        // std::cout << "in_points list : " << std::endl;
+
+        // for (; itbp != itep; itbp++)
+        // {
+        //     std::cout << *itbp << std::endl;
+        // }
+
+        // std::cout << "--- vertex " <<  vertI
+        //     << " : " << v->point() << std::endl
+        //     << "number of points used : " << in_points.size()
+        //     << std::endl;
+
+        k[vertI++] = Foam::sqrt
+        (
+            sqr(monge_form.principal_curvatures(0))
+          + sqr(monge_form.principal_curvatures(1))
+        );
+    }
+
+    return k;
+}
+
+
 // Main program:
 
 int main(int argc, char *argv[])
@@ -137,6 +293,22 @@ int main(int argc, char *argv[])
         "((x0 y0 z0)(x1 y1 z1))",
         "remove edges within specified bounding box"
     );
+    argList::addBoolOption
+    (
+        "writeObj",
+        "write featureEdgeMesh obj files"
+    );
+    argList::addOption
+    (
+        "closeness",
+        "scalar",
+        "span to look for surface closeness"
+    );
+    argList::addBoolOption
+    (
+        "writeVTK",
+        "write surface property VTK files"
+    );
 
 #   include "setRootCase.H"
 #   include "createTime.H"
@@ -163,7 +335,6 @@ int main(int argc, char *argv[])
     Info<< endl;
 
 
-
     // Either construct features from surface&featureangle or read set.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -319,9 +490,338 @@ int main(int argc, char *argv[])
 
     feMesh.write();
 
-    feMesh.writeObj(sFeatFileName);
+    if (args.optionFound("writeObj"))
+    {
+        feMesh.writeObj(runTime.constant()/"featureEdgeMesh"/sFeatFileName);
+    };
 
-    Info<< "End\n" << endl;
+    // Examine curvature, feature proximity and internal and external closeness.
+
+    // Internal and external closeness
+
+    triSurfaceMesh searchSurf
+    (
+        IOobject
+        (
+            sFeatFileName + ".closeness",
+            runTime.constant(),
+            "featureEdgeMesh",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf
+    );
+
+    // Prepare start and end points for intersection tests
+
+    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);
+
+    Info<< "span " << span << endl;
+
+    pointField start = searchSurf.faceCentres() - span*normals;
+    pointField end = searchSurf.faceCentres() + span*normals;
+    const pointField& faceCentres = searchSurf.faceCentres();
+
+    // {
+    //     label testI = 192;
+
+    //     List<List<pointIndexHit> > testAllHitInfo;
+
+    //     searchSurf.findLineAll
+    //     (
+    //         pointField(1, start[testI]),
+    //         pointField(1, end[testI]),
+    //         testAllHitInfo
+    //     );
+
+    //     Info<< testAllHitInfo << endl;
+
+    //     // return 0;
+    // }
+
+    List<List<pointIndexHit> > allHitInfo;
+
+    // Find all intersections (in order)
+    searchSurf.findLineAll(start, end, allHitInfo);
+
+    scalarField internalCloseness(start.size(), GREAT);
+    scalarField externalCloseness(start.size(), GREAT);
+
+    forAll(allHitInfo, fI)
+    {
+        const List<pointIndexHit>& hitInfo = allHitInfo[fI];
+
+        if (hitInfo.size() < 1)
+        {
+            Info<< nl << "# fI " << fI
+                << nl << "# start " << start[fI]
+                << nl << "# f centre " << faceCentres[fI]
+                << nl << "# end " << end[fI]
+                << endl;
+
+            meshTools::writeOBJ(Info, start[fI]);
+            meshTools::writeOBJ(Info, faceCentres[fI]);
+            meshTools::writeOBJ(Info, end[fI]);
+
+            Info<< "l 1 2 3" << endl;
+
+            meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
+            meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
+            meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
+
+            Info<< "f 4 5 6" << endl;
+
+            FatalErrorIn(args.executable())
+                << "findLineAll did not hit its own face."
+                << exit(FatalError);
+        }
+        else if (hitInfo.size() == 1)
+        {
+            if (!hitInfo[0].hit())
+            {
+                FatalErrorIn(args.executable())
+                    << "findLineAll did not hit any face."
+                    << exit(FatalError);
+            }
+            else if (hitInfo[0].index() != fI)
+            {
+                Info<< nl << "# findLineAll did not hit its own face."
+                    << nl << "# fI " << fI
+                    << nl << "# start " << start[fI]
+                    << nl << "# f centre " << faceCentres[fI]
+                    << nl << "# end " << end[fI]
+                    << nl << "# hitInfo " << hitInfo
+                    << endl;
+
+                meshTools::writeOBJ(Info, start[fI]);
+                meshTools::writeOBJ(Info, faceCentres[fI]);
+                meshTools::writeOBJ(Info, end[fI]);
+
+                Info<< "l 1 2 3" << endl;
+
+                meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
+                meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
+                meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
+
+                Info<< "f 4 5 6" << endl;
+
+                forAll(hitInfo, hI)
+                {
+                    label hFI = hitInfo[hI].index();
+
+                    meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
+                    meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
+                    meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
+
+                    Info<< "f "
+                        << 3*hI + 7 << " "
+                        << 3*hI + 8 << " "
+                        << 3*hI + 9
+                        << endl;
+                }
+                // FatalErrorIn(args.executable())
+                //     << "findLineAll did not hit its own face."
+                //     << exit(FatalError);
+            }
+        }
+        else
+        {
+            label ownHitI = -1;
+
+            forAll(hitInfo, hI)
+            {
+                // Find the hit on the triangle that launched the ray
+
+                if (hitInfo[hI].index() == fI)
+                {
+                    ownHitI = hI;
+
+                    break;
+                }
+            }
+
+            if (ownHitI < 0)
+            {
+                Info<< nl << "# findLineAll did not hit its own face."
+                    << nl << "# fI " << fI
+                    << nl << "# start " << start[fI]
+                    << nl << "# f centre " << faceCentres[fI]
+                    << nl << "# end " << end[fI]
+                    << nl << "# hitInfo " << hitInfo
+                    << endl;
+
+                meshTools::writeOBJ(Info, start[fI]);
+                meshTools::writeOBJ(Info, faceCentres[fI]);
+                meshTools::writeOBJ(Info, end[fI]);
+
+                Info<< "l 1 2 3" << endl;
+
+                meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
+                meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
+                meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
+
+                Info<< "f 4 5 6" << endl;
+
+                forAll(hitInfo, hI)
+                {
+                    label hFI = hitInfo[hI].index();
+
+                    meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
+                    meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
+                    meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
+
+                    Info<< "f "
+                        << 3*hI + 7 << " "
+                        << 3*hI + 8 << " "
+                        << 3*hI + 9
+                        << endl;
+                }
+
+                // FatalErrorIn(args.executable())
+                //     << "findLineAll did not hit its own face."
+                //     << exit(FatalError);
+            }
+            else if (ownHitI == 0)
+            {
+                // There are no internal hits, the first hit is the closest
+                // external hit
+
+                externalCloseness[fI] = mag
+                (
+                    faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
+                );
+            }
+            else if (ownHitI == hitInfo.size() - 1)
+            {
+                // There are no external hits, the last but one hit is the
+                // closest internal hit
+
+                internalCloseness[fI] = mag
+                (
+                    faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
+                );
+            }
+            else
+            {
+                externalCloseness[fI] = mag
+                (
+                    faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
+                );
+
+                internalCloseness[fI] = mag
+                (
+                    faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
+                );
+            }
+        }
+    }
+
+    triSurfaceScalarField internalClosenessField
+    (
+        IOobject
+        (
+            sFeatFileName + ".internalCloseness",
+            runTime.constant(),
+            "featureEdgeMesh",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf,
+        dimLength,
+        internalCloseness
+    );
+
+    internalClosenessField.write();
+
+    triSurfaceScalarField externalClosenessField
+    (
+        IOobject
+        (
+            sFeatFileName + ".externalCloseness",
+            runTime.constant(),
+            "featureEdgeMesh",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf,
+        dimLength,
+        externalCloseness
+    );
+
+    externalClosenessField.write();
+
+    scalarField k = curvature(surf);
+
+    triSurfacePointScalarField kField
+    (
+        IOobject
+        (
+            sFeatFileName + ".curvature",
+            runTime.constant(),
+            "featureEdgeMesh",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf,
+        dimLength,
+        k
+    );
+
+    kField.write();
+
+    if (args.optionFound("writeObj"))
+    {
+        vtkSurfaceWriter<scalar>().write
+        (
+            runTime.constant()/"triSurface",    // outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "internalCloseness",                // fieldName
+            internalCloseness,
+            false,                              // isNodeValues
+            true                                // verbose
+        );
+
+        vtkSurfaceWriter<scalar>().write
+        (
+            runTime.constant()/"triSurface",    // outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "externalCloseness",                // fieldName
+            externalCloseness,
+            false,                              // isNodeValues
+            true                                // verbose
+        );
+
+        vtkSurfaceWriter<scalar>().write
+        (
+            runTime.constant()/"triSurface",    // outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "curvature",                        // fieldName
+            k,
+            true,                               // isNodeValues
+            true                                // verbose
+        );
+    }
 
     return 0;
 }