LpMetric in R^d

ChangeLog.md

@@ -2,6 +2,7 @@
 # DGtal 1.0
 
 ## New Features / Critical Changes
+
 - *Base package*
   - Adding FunctorHolder to transform any callable object (e.g. function,
     functor, lambda function,...) into a valid DGtal functor.
@@ -37,6 +38,10 @@
   - Remove the internal object from VoxelComplex, improving performance
     (Pablo Hernandez, [#1369](https://github.com/DGtal-team/DGtal/pull/1369))
 
+- *Geometry*
+  - New LpMetric class (model of CMetricSpace) for distance computations in R^n.
+    (David Coeurjolly,  [#1388](https://github.com/DGtal-team/DGtal/pull/1388))
+
 ## Bug Fixes
 - *Configuration/General*
   - Continuous integration AppVeyor fix

examples/geometry/surfaces/exampleEstimatorFromSurfelFunctors.cpp

@@ -44,6 +44,7 @@
 #include "DGtal/shapes/Shapes.h"
 
 #include "DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h"
+#include "DGtal/geometry/volumes/distance/LpMetric.h"
 #ifdef WITH_CGAL
 #include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingGaussianCurvatureEstimator.h"
 #include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingMeanCurvatureEstimator.h"
@@ -109,16 +110,16 @@ int main(  )
   //constant convolution functor
   typedef functors::ConstValue<double> ConstFunctor;
 
-  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, Z3i::L2Metric, FunctorGaussian, ConstFunctor> ReporterK;
-  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, Z3i::L2Metric, FunctorMean, ConstFunctor> ReporterH;
-  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, Z3i::L2Metric, FunctorNormal, ConstFunctor> ReporterNormal;
-  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, Z3i::L2Metric, FunctorNormalLeast, ConstFunctor> ReporterNormalLeast;
+  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, LpMetric<Z3i::Space>, FunctorGaussian, ConstFunctor> ReporterK;
+  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, LpMetric<Z3i::Space>, FunctorMean, ConstFunctor> ReporterH;
+  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, LpMetric<Z3i::Space>, FunctorNormal, ConstFunctor> ReporterNormal;
+  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, LpMetric<Z3i::Space>, FunctorNormalLeast, ConstFunctor> ReporterNormalLeast;
 #endif
 
   ///For Elmentary convolution, we specify a Gaussian convolution
   ///kernel from the BasicFunctors.h file
   typedef functors::ElementaryConvolutionNormalVectorEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorNormalElem;
-  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, Z3i::L2Metric,
+  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, LpMetric<Z3i::Space>,
                                                  FunctorNormalElem, DGtal::functors::GaussianKernel> ReporterNormalElem;
   //! [SurfelFunctorsType]
 
@@ -140,10 +141,11 @@ int main(  )
   ReporterNormalLeast reporterL;
 #endif
 
+  LpMetric<Z3i::Space> l2;
   FunctorNormalElem estimatorNormalElem(embedder,1.0);
   ///sigma = 2.0 for the gaussian smoothing
   DGtal::functors::GaussianKernel gaussian(2.0);
-  ReporterNormalElem reporterElem(surface, l2Metric,
+  ReporterNormalElem reporterElem(surface, l2,
                                   estimatorNormalElem, gaussian);
   //! [SurfelFunctorsInstances]
 
@@ -159,10 +161,10 @@ int main(  )
   reporterN.init(1, surface.begin(), surface.end());
   reporterL.init(1, surface.begin(), surface.end());
 
-  reporterK.setParams(l2Metric, estimatorK, constFunctor, 5.0);
-  reporterH.setParams(l2Metric, estimatorH, constFunctor, 5.0);
-  reporterN.setParams(l2Metric, estimatorN, constFunctor, 5.0);
-  reporterL.setParams(l2Metric, estimatorL, constFunctor, 5.0);
+  reporterK.setParams(l2, estimatorK, constFunctor, 5.0);
+  reporterH.setParams(l2, estimatorH, constFunctor, 5.0);
+  reporterN.setParams(l2, estimatorN, constFunctor, 5.0);
+  reporterL.setParams(l2, estimatorL, constFunctor, 5.0);
 
   FunctorGaussian::Quantity valK = reporterK.eval( surface.begin());
   FunctorMean::Quantity valH = reporterH.eval( surface.begin());
@@ -171,7 +173,7 @@ int main(  )
 #endif
 
   reporterElem.attach(surface);
-  reporterElem.setParams(l2Metric,
+  reporterElem.setParams(l2,
                          estimatorNormalElem, gaussian, 5.0);
   reporterElem.init(1.0, surface.begin(), surface.end());
   FunctorNormalElem::Quantity valElem = reporterElem.eval( surface.begin());

src/DGtal/geometry/doc/moduleMetrics.dox

@@ -62,12 +62,14 @@ In the following, we restrict our attention to metric spaces even if some of the
 In DGtal, the concept of CMetricSpace is associated with metric space
 definition. A first refinement of this concept is the concept of
 CDigitalMetricSpace whose models implement distance functions on
-integer numbers (i.e. @f$(\mathbb{Z}^n, \mathbb{Z}, d)@f$)). For
+integer numbers (i.e. @f$(\mathbb{Z}^n, \mathbb{Z}, d)@f$). For
 example,
 
+- LpMetric are models of CMetricSpace on @f$(\mathbb{R}^n, \mathbb{R}, d)@f$ spaces (i.e. all computations are performed
+on C++ "double").
 - @f$ l_p@f$ metrics (ExactPredicateLpSeparableMetric or
-  InexactPredicateLpSeparableMetric) are models of CMetricSpace for
-  @f$p>1@f$.
+  InexactPredicateLpSeparableMetric) are models of CMetricSpace on @f$(\mathbb{Z}^n, \mathbb{Z}, d)@f$ for
+  @f$p>1@f$ with specific separability methods.
 - Chamfer norms with axis symmetric unit ball (see experimental::ChamferNorm2D in dimension 2).
 - DigitalMetricAdapter which adapts any metric (model of CMetricSpace)
   to digital space using the ceil function. This induces a digital

src/DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h

@@ -93,7 +93,7 @@ namespace DGtal
    * canonical embedding of surfel elements (cf CanonicSCellEmbedder).
    *
    *  @tparam TDigitalSurfaceContainer any model of digital surface container concept (CDigitalSurfaceContainer)
-   *  @tparam TMetric any model of CMetricSpace to be used in the neighborhood construction.
+   *  @tparam TMetric any model of CMetricSpace to be used in the neighborhood construction (e.g. LpMetric)
    *  @tparam TFunctorOnSurfel an estimator on surfel set (model of CLocalEstimatorFromSurfelFunctor)
    *  @tparam TConvolutionFunctor type of  functor on double
    *  [0,1]->[0,1] to implement the response of a symmetric convolution kernel.

src/DGtal/geometry/volumes/distance/CMetricSpace.h

@@ -114,7 +114,7 @@ squares of Point::Coordiante values.
 
 ### Models
 
-ExactPredicateLpSeparableMetric, InexactPredicateLpSeparableMetric.
+LpMetric, ExactPredicateLpSeparableMetric, InexactPredicateLpSeparableMetric.
 
 ### Notes
 

src/DGtal/geometry/volumes/distance/LpMetric.h

@@ -0,0 +1,236 @@
+/**
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as
+ *  published by the Free Software Foundation, either version 3 of the
+ *  License, or  (at your option) any later version.
+ *
+ *  This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ **/
+
+#pragma once
+
+/**
+ * @file LpMetric.h
+ * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
+ * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
+ *
+ * @date 2019/27/01
+ *
+ *
+ * This file is part of the DGtal library.
+ */
+
+#if defined(LpMetric_RECURSES)
+#error Recursive header files inclusion detected in LpMetric.h
+#else // defined(LpMetric_RECURSES)
+/** Prevents recursive inclusion of headers. */
+#define LpMetric_RECURSES
+
+#if !defined LpMetric_h
+/** Prevents repeated inclusion of headers. */
+#define LpMetric_h
+
+//////////////////////////////////////////////////////////////////////////////
+// Inclusions
+#include <iostream>
+#include <cmath>
+#include "DGtal/base/Common.h"
+#include "DGtal/kernel/CSpace.h"
+#include "DGtal/kernel/CEuclideanRing.h"
+//////////////////////////////////////////////////////////////////////////////
+
+namespace DGtal
+{
+
+/////////////////////////////////////////////////////////////////////////////
+// template class LpMetric
+/**
+ * Description of template class 'LpMetric' <p>
+ * \brief Aim: implements  l_p metrics.
+ *
+ * Given a parameter p, the class implements classical l_p
+ * metric as a model of CMetricSpace. Hence, given two points
+ * @f$ x=(x_0...x_{n-1})@f$, @f$ y=(y_0...y_{n-1})@f$, we define a metric as:
+ *
+ * @f$ distance(x,y)= \left(
+ * \sum_{i=0}^{n-1} |x_i-y_i|^p\right)^{1/p}@f$
+ *
+ * This class performs all computations on C++ double converting the digital
+ * points to Space::RealPoint
+ *
+ * The exponent @a p is specifed at the constructor.
+ *
+ * @tparam TSpace the model of CSpace on which the metric is
+ * defined.
+ */
+  template <typename TSpace>
+  class LpMetric
+  {
+    // ----------------------- Standard services ------------------------------
+  public:
+
+    ///Copy the space type
+    typedef TSpace Space;
+    BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
+
+    ///Type for points (RealPoint for this class)
+    typedef typename Space::RealPoint Point;
+    ///Type for distance values
+    typedef double Value;
+    ///Type for raw distance values
+    typedef double RawValue;
+
+    /**
+     * Default Constructor (p=2.0)
+     *
+     */
+    LpMetric(): myExponent(2.0) {}
+
+
+      /**
+     * Constructor.
+     *
+     * @param anExponent the exponent (p) of the lp metric.
+     */
+    LpMetric( const double anExponent): myExponent(anExponent)
+    {}
+
+
+    /**
+     * Destructor.
+     */
+    ~LpMetric()
+    {}
+
+    /**
+     * Copy constructor.
+     * @param other the object to clone.
+     */
+    LpMetric ( const LpMetric & other )
+    {
+      myExponent = other.myExponent;
+    }
+
+    /**
+     * Assignment.
+     * @param other the object to copy.
+     * @return a reference on 'this'.
+     */
+    LpMetric & operator= ( const LpMetric & other )
+    {
+      myExponent = other.myExponent;
+      return *this;
+    }
+
+    // ----------------------- Interface --------------------------------------
+  public:
+
+    // ----------------------- CMetric --------------------------------------
+    /**
+     * Compute the distance between @a aP and @a aQ.
+     *
+     * @param aP a first point.
+     * @param aQ a second point.
+     *
+     * @return the distance between aP and aQ.
+     */
+    Value operator()(const Point & aP, const Point &aQ) const
+    {
+      return std::pow( rawDistance(aP,aQ), 1.0/myExponent);
+    }
+
+    /**
+     * Compute the raw distance between @a aP and @a aQ.
+     *
+     * @param aP a first point.
+     * @param aQ a second point.
+     *
+     * @return the distance between aP and aQ.
+     */
+    RawValue rawDistance(const Point & aP, const Point &aQ) const
+    {
+      RawValue tmp=0.0;
+      for(auto i = 0; i < aP.size(); ++i)
+        tmp += std::pow(NumberTraits<typename Point::Coordinate>::castToDouble(std::abs(aP[i] - aQ[i])),
+                        myExponent);
+      return tmp;
+    }
+
+    /**
+     * Given an origin and two points, this method decides which one
+     * is closest to the origin. This method should be faster than
+     * comparing distance values.
+     *
+     * @param origin the origin
+     * @param first  the first point
+     * @param second the second point
+     *
+     * @return a Closest enum: FIRST, SECOND or BOTH.
+     */
+    Closest closest(const Point &origin,
+                    const Point &first,
+                    const Point &second) const
+    {
+      auto dfirst = rawDistance(origin,first);
+      auto dsecond = rawDistance(origin,second);
+      if (dfirst < dsecond)
+        return ClosestFIRST;
+      else
+        if (dfirst>dsecond)
+          return ClosestSECOND;
+
+      return ClosestBOTH;
+    }
+
+    /**
+     * Writes/Displays the object on an output stream.
+     * @param out the output stream where the object is written.
+     */
+    void selfDisplay ( std::ostream & out ) const
+    {
+      out << "[LpMetric] Lp Metric exponent=" << myExponent ;
+    }
+
+    /**
+     * Checks the validity/consistency of the object.
+     * @return 'true' if the object is valid, 'false' otherwise.
+     */
+    bool isValid() const
+    {
+      return true;
+    }
+
+    // ------------------------- Private Datas --------------------------------
+  private:
+
+    ///Exponent value
+    Value myExponent;
+
+  }; // end of class LpMetric
+
+  /**
+   * Overloads 'operator<<' for displaying objects of class 'LpMetric'.
+   * @param out the output stream where the object is written.
+   * @param object the object of class 'LpMetric' to write.
+   * @return the output stream after the writing.
+   */
+  template <typename T>
+  std::ostream&
+  operator<< ( std::ostream & out, const LpMetric<T> & object );
+
+} // namespace DGtal
+
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#endif // !defined LpMetric_h
+
+#undef LpMetric_RECURSES
+#endif // else defined(LpMetric_RECURSES)