diff --git a/ChangeLog.md b/ChangeLog.md index 78547b9e52..50a37f2cb8 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,12 +1,16 @@ # DGtal 1.5beta +## New features + +- *Geometry* + - Implementation of the plane-probing L-algorithm (Tristan Roussillon, [#1744](https://github.com/DGtal-team/DGtal/pull/1744)) + ## Bug fixes - *General* - Fix cmake CGAL 6.0 Breaking change. (David Coeurjolly, [#1745](https://github.com/DGtal-team/DGtal/pull/1745)) - Adding a new `DGTAL_REMOVE_UNINSTALL` cmake option to disable the `uninstall` target. (David Coeurjolly, [#1746](https://github.com/DGtal-team/DGtal/pull/1746) - - *Geometry* - Bug fix in ArithmeticalDSSComputerOnSurfels (Tristan Roussillon, [#1742](https://github.com/DGtal-team/DGtal/pull/1742)) diff --git a/examples/geometry/surfaces/examplePlaneProbingTetrahedronEstimator.cpp b/examples/geometry/surfaces/examplePlaneProbingTetrahedronEstimator.cpp index 0e178180db..55faa252c4 100644 --- a/examples/geometry/surfaces/examplePlaneProbingTetrahedronEstimator.cpp +++ b/examples/geometry/surfaces/examplePlaneProbingTetrahedronEstimator.cpp @@ -51,7 +51,7 @@ int main(void) //! [PlaneProbingTetrahedronEstimatorConstruction] // The general form is ProbingEstimator where // - Predicate is a model of concepts::PointPredicate, see DigitalPlanePredicate or DigitalSurfacePredicate for instance, - // - mode specifies the candidate set, it is one of { ProbingMode::H, ProbingMode::R, ProbingMode::R1 }. + // - mode specifies the candidate set, it is one of { ProbingMode::H, ProbingMode::R, ProbingMode::R1, ProbingMode::L }. using DigitalPlane = DigitalPlanePredicate; using Estimator = PlaneProbingTetrahedronEstimator; diff --git a/src/DGtal/doc/global.bib b/src/DGtal/doc/global.bib index 63e3de10ed..d8f4ae94f2 100644 --- a/src/DGtal/doc/global.bib +++ b/src/DGtal/doc/global.bib @@ -1329,6 +1329,24 @@ @article{LMRJMIV2020 HAL_VERSION = {v1}, } +@InProceedings{Lu2022, +author="Lu, Jui-Ting +and Roussillon, Tristan +and Coeurjolly, David", +editor="Baudrier, {\'E}tienne +and Naegel, Beno{\^i}t +and Kr{\"a}henb{\"u}hl, Adrien +and Tajine, Mohamed", +title="A New Lattice-Based Plane-Probing Algorithm", +booktitle="Discrete Geometry and Mathematical Morphology", +year="2022", +publisher="Springer International Publishing", +address="Cham", +pages="366--381", +abstract="Plane-probing algorithms have become fundamental tools to locally capture arithmetical and geometrical properties of digital surfaces (boundaries of a connected set of voxels), and especially normal vector information. On a digital plane, the overall idea is to consider a local pattern, a triangle, that is expanded starting from a point of interest using simple probes of the digital plane with a predicate ``Is a point x in the digital plane?''. Challenges in plane-probing methods are to design an algorithm that terminates on a triangle with several geometrical properties: its normal vector should match with the expected one for digital plane (correctness), the triangle should be as compact as possible (acute or right angles only), and probes should be as close as possible to the source point (locality property). In addition, we also wish to minimize the number of iterations or probes during the computations. Existing methods provide correct outputs but only experimental evidence for these properties. In this paper, we present a new plane-probing algorithm that is theoretically correct on digital planes, and with better experimental compactness and locality than existing solutions. Additional properties of this new approach also suggest that theoretical proofs of the aforementioned geometrical properties could be achieved.", +isbn="978-3-031-19897-7" +} + @INPROCEEDINGS{Lachaud03c, AUTHOR = {J.-O. Lachaud and A. Vialard}, TITLE = {Geometric measures on arbitrary dimensional digital surfaces}, diff --git a/src/DGtal/geometry/doc/modulePlaneProbing.dox b/src/DGtal/geometry/doc/modulePlaneProbing.dox index 31f16a52e7..9693a35867 100644 --- a/src/DGtal/geometry/doc/modulePlaneProbing.dox +++ b/src/DGtal/geometry/doc/modulePlaneProbing.dox @@ -38,7 +38,7 @@ testPlaneProbingParallelepipedEstimator.cpp. \section sectPlaneProbing1 Introduction to plane-probing algorithms -A plane-probing algorithm (see @cite LPRJMIV2017, @cite RLDGCI2019 and @cite LMRJMIV2020) +A plane-probing algorithm (see @cite LPRJMIV2017, @cite RLDGCI2019, @cite LMRJMIV2020 and @cite Lu2022) computes the normal vector of a set of digital points from a starting point and a predicate \b InPlane: "Is a point x in the set of digital points?". This predicate is used to probe the set as locally as possible diff --git a/src/DGtal/geometry/doc/packageGeometry.dox b/src/DGtal/geometry/doc/packageGeometry.dox index 3cdb630dde..95dc39cfba 100644 --- a/src/DGtal/geometry/doc/packageGeometry.dox +++ b/src/DGtal/geometry/doc/packageGeometry.dox @@ -71,7 +71,7 @@ of arbitrary dimension, by the means of separable and incremental distance trans - \subpage moduleIntegralInvariant (Jérémy Levallois, David Coeurjolly, Jacques-Olivier Lachaud) - \subpage LocalEstimatorsFromSurfel (David Coeurjolly) - \subpage moduleVCM (Louis Cuel, Jacques-Olivier Lachaud, Quentin Mérigot, Boris Thibert) - - \subpage modulePlaneProbing (Jacques-Olivier Lachaud, Jocelyn Meyron, Tristan Roussillon) + - \subpage modulePlaneProbing (Jacques-Olivier Lachaud, Jui-Ting Lu, Jocelyn Meyron, Tristan Roussillon) - \subpage moduleMaximalSegmentSliceEstimation (Jocelyn Meyron, Tristan Roussillon) - Mesh geometric estimators diff --git a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h index fc361e8645..0b2b236885 100644 --- a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h +++ b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h @@ -23,6 +23,11 @@ * * @date 2020/12/04 * + * @author Tristan Roussilllon (\c tristan.roussillon@liris.cnrs.fr ) + * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France + * + * @date 2024/09/16 + * * Helper functions for plane-probing algorithms. * * This file is part of the DGtal library. @@ -94,7 +99,7 @@ namespace DGtal // template class PointOnProbingRay /** * Description of template class 'PointOnProbingRay'

- * A ray consists of a permutation \f$ \sigma \f$ and an integer index \f$ \lambda \f$ (position on the ray). + * \brief A ray consists of a permutation \f$ \sigma \f$ and an integer index \f$ \lambda \f$ (position on the ray). * For a triplet of vectors \f$ (m_k)_{0 \leq k \leq 2} \f$ and a point \f$ q \f$, a point on the ray is defined as: * \f$ q - m_{\sigma(0)} + m_{\sigma(1)} + \lambda m_{\sigma(2)} \f$. \f$ q - m_{\sigma(0)} + m_{\sigma(1)} \f$ is called the \e base point. * @@ -103,14 +108,14 @@ namespace DGtal * * @tparam Integer the integer type, model of concepts::CInteger. */ - template < typename Integer = int > + template < typename Integer = int, typename Index = std::size_t > class PointOnProbingRay { BOOST_CONCEPT_ASSERT(( concepts::CInteger ) ); // ----------------------- Public types ------------------------------ public: - using Permutation = std::array; + using Permutation = std::array; public: /** @@ -122,9 +127,9 @@ namespace DGtal * Constructs a ray with a permutation and an index. * * @param aSigma a permutation. - * @param aIndex an index. + * @param aInt an integer that determines a point along the ray. */ - PointOnProbingRay (Permutation const& aSigma, Integer const& aIndex = Integer(0)); + PointOnProbingRay (Permutation const& aSigma, Integer const& aInt = Integer(0)); /** * @return the base point of the ray (with index 0). @@ -140,13 +145,23 @@ namespace DGtal * @param aIndex an index between 0 and 2. * @return the i-th element of the permutation that defines the ray. */ - int sigma (int aIndex) const; + Index sigma (Index const& aIndex) const; /** - * @return index of the current point on the ray. + * @return integer that locates the current point on the ray. */ - Integer const& index () const; + Integer const& position () const; + /** + * Returns the geometric realization of this grid point. + * + * @param aM an array of three points. + * @tparam Point a type for points. + * @return the computed point. + */ + template < typename Point > + Point relativePoint (std::array const& aM) const; + /** * Equality test between two rays: the internal permutations and * indices must be the same. @@ -189,7 +204,7 @@ namespace DGtal private: Permutation mySigma; /**< The permutation. */ - Integer myIndex; /**< The index. */ + Integer myPosition; /**< The index. */ }; // end of class PointOnProbingRay /** @@ -199,8 +214,316 @@ namespace DGtal * @param aRay the probing ray to display. * @return the output stream after the writing. */ - template < typename Integer > - std::ostream& operator<< (std::ostream& aOs, PointOnProbingRay const& aRay); + template < typename Integer, typename Index > + std::ostream& operator<< (std::ostream& aOs, PointOnProbingRay const& aRay); + + ///////////////////////////////////////////////////////////////////////////// + // template class GridPoint + /** + * Description of template class 'GridPoint'

+ * \brief + * A grid point consists of a couple of nonnegative coordinates \f$ (x,y) \f$ and an integer index \f$ k \f$ + * that determines a point used as origin. + * For a triplet of vectors \f$ (m_k)_{0 \leq k \leq 2} \f$ and a point \f$ q \f$, a grid point is defined as: + * \f$ q - m_{k} + x m_{(k+1)\bmod 3} + y m_{(k+2)\bmod 3} \f$. \f$ q - m_{k} \f$, called base point, is used + * as origin. + * + * This class is used to represent candidate points for the plane-probing L-algorithm. In practice, + * the point \f$ q \f$ is the fixed point and the three vectors \f$ (m_k)_{0 \leq k \leq 2} \f$ are + * the vectors defining the current probing frame. + * + * @tparam Integer the integer type, model of concepts::CInteger. + */ + template < typename Integer = int, typename Index = std::size_t > + class GridPoint + { + BOOST_CONCEPT_ASSERT(( concepts::CInteger ) ); + + public: + + /** + * Default constructor. + */ + GridPoint () = default; + + /** + * Constructs a grid point from a couple of coordinates and + * the index of the point used as origin. + * + * @param aDir a pair of nonnegative integers. + * @param aK an index in {0,1,2}. + */ + GridPoint (std::pair const& aDir, Index const& aK ) : myDir(aDir), myK(aK) {} + + /** + * Constructs a grid point from a couple of coordinates and + * the index of the point used as origin. + * + * @param aX first coordinate. + * @param aY second coordinate. + * @param aK an index in {0,1,2}. + */ + GridPoint (Integer const& aX, Integer const& aY, Index const& aK ) : myDir(std::make_pair(aX,aY)), myK(aK) {} + + /** + * Returns the couple of coordinates, i.e., + * the direction going from the origin to + * the grid point. + * + * @return the direction. + */ + std::pair direction() const { + return myDir; + } + + /** + * Returns the index of the point used as origin. + * + * @return the index. + */ + Index k() const { + return myK; + } + + /** + * Returns the vector going from the base point to the + * geometric realization of this object. + * + * @param aM an array of three vectors. + * @tparam Vector a type for vectors. + * @return the computed vector. + */ + template < typename Vector > + Vector directionVector (std::array const& aM) const { + return aM[(myK+1)%3]*myDir.first + aM[(myK+2)%3]*myDir.second; + } + + /** + * Returns the geometric realization of this grid point. + * + * @param aM an array of three points. + * @tparam Point a type for points. + * @return the computed point. + */ + template < typename Point > + Point relativePoint (std::array const& aM) const { + return -aM[myK] + directionVector(aM); + } + + /** + * Tells whether this grid point is equal to another or not. + * Two grid points are equal iff their members are equal. + * + * @param other another grid point. + * @return 'true' if equal, 'false' otherwise. + */ + bool operator== (GridPoint const& other) const { + return (myDir == other.myDir) && (myK == other.myK); + } + + /** + * Tells whether this grid point is different from another or not. + * + * @param other another grid point. + * @return 'true' if different, i.e., not equal, 'false' otherwise. + */ + bool operator!= (GridPoint const& other) const { + return !(*this == other); + } + + /** + * Returns a grid point given a couple of coordinates. + * + * @param aDir a couple of coordinates. + * @return the resulting grid point. + */ + GridPoint getOnSameGrid(const std::pair& aDir) const { + return GridPoint(aDir,myK); + } + + /** + * Adds a grid point to this one (as if they were vectors). + * + * @param other another grid point. + * @return the resulting point. + */ + GridPoint operator+(const GridPoint & other) const { + ASSERT(myK == other.myK); + std::pair d = std::make_pair(myDir.first+other.myDir.first, + myDir.second+other.myDir.second); + return getOnSameGrid(d); + } + + /** + * Scales this grid point by a scalar (as if it was vector). + * + * @param aValue a scalar value + * @return the resulting point. + */ + GridPoint operator*(Integer aValue) const { + std::pair d = std::make_pair(myDir.first*aValue, + myDir.second*aValue); + return getOnSameGrid(d); + } + + /** + * Checks whether the representation of the grid point is valid, + * i.e., the coordinates are nonnegative coordinates, not both + * equal to zero, and the index is in {0,1,2}. + * + * @return 'true' if valid, 'false' otherwise. + */ + bool isValid() const { + if ( (myDir.first != 0) || (myDir.second != 0) ) { //not both null + return ( (myDir.first >= 0) && (myDir.second >= 0) + && (myK >= 0) && (myK <= 2) + ); + } else { + return false; + } + } + + private: + + std::pair myDir; /**< Couple of coordinates giving a direction */ + Index myK; /**< Index of a point used as origin */ + + }; //end of class GridPoint + + /** + * Display a grid point on the standard output. + * + * @param aOs the output stream where the object is written. + * @param aGridPoint the grid point to display. + * @return the output stream after the writing. + */ + template < typename Integer, typename Index > + std::ostream& operator<< (std::ostream& aOs, GridPoint const& aGridPoint) { + aOs << "GridPoint[k=" << aGridPoint.k() + << ", a=" << aGridPoint.direction().first + << ", b=" << aGridPoint.direction().second + << "]"; + return aOs; + } + + ///////////////////////////////////////////////////////////////////////////// + // template class GridPointOnProbingRay + /** + * Description of template class 'GridPointOnProbingRay'

+ * \brief Aim: Represents a grid point along a discrete ray defined on a grid. + * + * More precisely, a ray consists of a starting point (represented as an instance of 'GridPoint') + * and a direction (respresented as a couple of coordinates for the basis of the underlying grid). + * A grid point along that ray is determined by an index, 0 being the starting point of the ray. + * The class provides several methods to compare and move points along the ray. + * + * @tparam Integer the integer type, model of concepts::CInteger. + */ + template < typename Integer = int, typename Index = std::size_t > + class GridPointOnProbingRay + { + BOOST_CONCEPT_ASSERT(( concepts::CInteger ) ); + + public: + /** + * Default constructor. + */ + GridPointOnProbingRay () = default; + + /** + * Constructor. + * + * @param aGridPoint starting point of the ray + * @param aDirection direction of the ray + * @param aIdx index of the grid point along the ray + */ + GridPointOnProbingRay (const GridPoint& aGridPoint, + const std::pair& aDirection, + const Integer& aIdx = 0) + : myOrigin(aGridPoint), myDirection(aDirection), myIdx(aIdx) {} + + /** + * Equality test. The two objects are equal iff + * the underlying rays are the same and + * the indices are the same. + * + * @param other another instance of GridPointOnProbingRay + * @return 'true' if equal, 'false' otherwise + */ + bool operator== (GridPointOnProbingRay const& other) const { + return ( (myOrigin == other.myOrigin) && + (myDirection == other.myDirection) && + (myIdx == other.myIdx) ); + } + + /** + * Difference test. + * + * @param other another instance of GridPointOnProbingRay + * @return 'true' if different, i.e. not equal, 'false' otherwise + */ + bool operator!= (GridPointOnProbingRay const& other) const { + return !(*this == other); + } + + /** + * Returns a grid point lying after this one along the ray. + * The distance is given as an input parameter. + * + * @param aInc an increment. + * @return a new grid point on the same ray, with index, + * the current index incremented by 'aInc'. + */ + GridPointOnProbingRay next(const Integer& aInc) const { + return GridPointOnProbingRay(myOrigin, myDirection, myIdx+aInc); + } + + /** + * Returns a grid point lying before this one along the ray. + * The distance is given as an input parameter. + * + * @param aDec a decrement. + * @return a new grid point on the same ray, with index, + * the current index decremented by 'aDec'. + */ + GridPointOnProbingRay previous(const Integer& aDec) const { + return GridPointOnProbingRay(myOrigin, myDirection, myIdx-aDec); + } + + /** + * @return index of the current grid point on the ray. + */ + Integer index() const { + return myIdx; + } + + /** + * @return the current grid point as an instance of GridPoint. + */ + GridPoint gridPoint() const { + return myOrigin + myOrigin.getOnSameGrid(myDirection)*myIdx; + } + + /** + * Returns the geometric realization of this grid point. + * + * @param aM an array of three points. + * @tparam Point a type for points. + * @return the computed point. + */ + template < typename Point > + Point relativePoint (std::array const& aM) const { + return gridPoint().relativePoint(aM); + } + + + private: + GridPoint myOrigin; /**< starting point of the ray */ + std::pair myDirection; /**< direction of the ray */ + Integer myIdx; /**< index of the point along the ray */ + + }; //end of class GridPointOnProbingRay + } // namespace detail } // namespace DGtal diff --git a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.ih b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.ih index 4a91f5c9ef..619055842b 100644 --- a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.ih +++ b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.ih @@ -119,107 +119,115 @@ DGtal::detail::isBasisReduced (Point const& aU, Point const& aV) // ----------------------- Standard services ------------------------------ // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > inline -DGtal::detail::PointOnProbingRay:: -PointOnProbingRay (Permutation const& aSigma, Integer const& aIndex) - : mySigma(aSigma), myIndex(aIndex) +DGtal::detail::PointOnProbingRay:: +PointOnProbingRay (Permutation const& aSigma, Integer const& aPosition) + : mySigma(aSigma), myPosition(aPosition) { } // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > inline -DGtal::detail::PointOnProbingRay -DGtal::detail::PointOnProbingRay::getBase () const +DGtal::detail::PointOnProbingRay +DGtal::detail::PointOnProbingRay::getBase () const { return PointOnProbingRay(mySigma, 0); } // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > inline -typename DGtal::detail::PointOnProbingRay::Permutation const& -DGtal::detail::PointOnProbingRay::sigma () const +typename DGtal::detail::PointOnProbingRay::Permutation const& +DGtal::detail::PointOnProbingRay::sigma () const { return mySigma; } // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > inline -int -DGtal::detail::PointOnProbingRay::sigma (int aIndex) const +Index +DGtal::detail::PointOnProbingRay::sigma (Index const& aIndex) const { assert(aIndex >= 0 && aIndex <= 2); return mySigma[aIndex]; } // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > inline Integer const& -DGtal::detail::PointOnProbingRay::index () const +DGtal::detail::PointOnProbingRay::position () const { - return myIndex; + return myPosition; } // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > +template < typename Point > +inline +Point +DGtal::detail::PointOnProbingRay::relativePoint (std::array const& aM) const { + return -aM[mySigma[0]] + aM[mySigma[1]] + aM[mySigma[2]] * myPosition; +} + +// ------------------------------------------------------------------------ +template < typename Integer, typename Index > inline bool -DGtal::detail::PointOnProbingRay::operator== (PointOnProbingRay const& aRay) const +DGtal::detail::PointOnProbingRay::operator== (PointOnProbingRay const& aRay) const { - return (mySigma == aRay.mySigma) && (myIndex == aRay.index()); + return (mySigma == aRay.mySigma) && (myPosition == aRay.position()); } - // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > inline bool -DGtal::detail::PointOnProbingRay::operator!= (PointOnProbingRay const& aRay) const +DGtal::detail::PointOnProbingRay::operator!= (PointOnProbingRay const& aRay) const { return !(*this == aRay); } // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > inline bool -DGtal::detail::PointOnProbingRay::operator<= (PointOnProbingRay const& aRay) const +DGtal::detail::PointOnProbingRay::operator<= (PointOnProbingRay const& aRay) const { - return (mySigma == aRay.mySigma) && (myIndex <= aRay.index()); + return (mySigma == aRay.mySigma) && (myPosition <= aRay.position()); } // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > inline -DGtal::detail::PointOnProbingRay -DGtal::detail::PointOnProbingRay::next (Integer const& aInc) const +DGtal::detail::PointOnProbingRay +DGtal::detail::PointOnProbingRay::next (Integer const& aInc) const { - return PointOnProbingRay(mySigma, myIndex + aInc); + return PointOnProbingRay(mySigma, myPosition + aInc); } // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > inline -DGtal::detail::PointOnProbingRay -DGtal::detail::PointOnProbingRay::previous (Integer const& aDec) const +DGtal::detail::PointOnProbingRay +DGtal::detail::PointOnProbingRay::previous (Integer const& aDec) const { - return PointOnProbingRay(mySigma, myIndex - aDec); + return PointOnProbingRay(mySigma, myPosition - aDec); } // ------------------------------------------------------------------------ -template < typename Integer > +template < typename Integer, typename Index > inline std::ostream& -DGtal::detail::operator<< (std::ostream& aOs, PointOnProbingRay const& aRay) +DGtal::detail::operator<< (std::ostream& aOs, PointOnProbingRay const& aRay) { aOs << "sigma=(" << aRay.sigma(0) << ", " << aRay.sigma(1) << ", " << - aRay.sigma(2) << "); i=" << aRay.index(); + aRay.sigma(2) << "); i=" << aRay.position(); return aOs; } diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h new file mode 100644 index 0000000000..ba0947efa7 --- /dev/null +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h @@ -0,0 +1,301 @@ +/** + * 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 . + * + **/ + +#pragma once + +/** + * @file + * @author Tristan Roussilllon (\c tristan.roussillon@liris.cnrs.fr ) + * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France + * + * @date 2024/09/16 + * + * Header file for module PlaneProbingLNeighborhood.cpp + * + * This file is part of the DGtal library. + */ + +#if defined(PlaneProbingLNeighborhood_RECURSES) +#error Recursive header files inclusion detected in PlaneProbingLNeighborhood.h +#else // defined(PlaneProbingLNeighborhood_RECURSES) +/** Prevents recursive inclusion of headers. */ +#define PlaneProbingLNeighborhood_RECURSES + +#if !defined PlaneProbingLNeighborhood_h +/** Prevents repeated inclusion of headers. */ +#define PlaneProbingLNeighborhood_h + +////////////////////////////////////////////////////////////////////////////// +// Inclusions +#include +#include +#include +#include "DGtal/base/Common.h" +#include "DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h" +#include "DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.h" +#include "DGtal/kernel/CPointPredicate.h" +////////////////////////////////////////////////////////////////////////////// + +namespace DGtal +{ + + ///////////////////////////////////////////////////////////////////////////// + // template class PlaneProbingLNeighborhood + /** + * Description of template class 'PlaneProbingLNeighborhood'

+ * \brief Aim: Represents a way to probe the L-neighborhood, + * see \cite Lu2022 for details. + * + * \tparam TPredicate the probing predicate, a model of concepts::CPointPredicate. + */ + template + class PlaneProbingLNeighborhood: public DGtal::PlaneProbingRNeighborhood + { + BOOST_CONCEPT_ASSERT((DGtal::concepts::CPointPredicate)); + + // ----------------------- Public types ------------------------------ + public: + using Predicate = TPredicate; + using Point = typename Predicate::Point; + using Vector = Point; + using Integer = typename Point::Coordinate; + using Triangle = std::array; + + using HexagonState = typename PlaneProbingNeighborhood::HexagonState; + using UpdateOperation = typename PlaneProbingNeighborhood::UpdateOperation; + + using Index = typename PlaneProbingNeighborhood::Index; + using PointOnProbingRay = typename PlaneProbingNeighborhood::PointOnProbingRay; + using GridPoint = typename detail::GridPoint; + using GridPointOnProbingRay = typename detail::GridPointOnProbingRay; + + // ----------------------- Internal type ------------------------------- + private: + /** + * Description of data structure 'ClosestGridPoint'

+ * \brief Aim: Used to store the closest grid point associated + * to a vertex of the triangle and two extra boolean values + * about the local configuration at that vertex. + * + * More precisely, given a triplet of vectors \f$ (m_k)_{0 \leq k \leq 2} \f$ + * and a point \f$ q \f$, let us denote \f$ v \f$ the vertex equal to + * \f$ q - m_k \f$. The first boolean value is 'true' iff + * the predicate returns 'true' on \f$ v - m_{k+1} \f$ and, + * similarly, the second boolean value is 'true' iff + * the predicate returns 'true' on \f$ v - m_{k+2} \f$ + * (the indices are taken modulo 3). + */ + struct ClosestGridPoint + { + /** + * Default constructor. + */ + ClosestGridPoint () = default; + + /** + * Constructor + * + * @param aGridPoint the closest grid point (possibly invalid) + * @param aFirst the first boolean value + * @param aSecond the second boolean value + */ + ClosestGridPoint (const GridPoint& aGridPoint, + const bool& aFirst, const bool& aSecond ) + : myGridPoint(aGridPoint), myPair(std::make_pair(aFirst,aSecond)) {} + + GridPoint myGridPoint; /**< a grid point, which can be invalid + if no grid point belong to the underlying surface */ + + std::pair myPair; /**< pair of boolean values that encode + the local configuration */ + }; + + // ----------------------- Standard services ------------------------------ + public: + /** + * Default constructor. + */ + PlaneProbingLNeighborhood() = delete; + + /** + * Constructor. + * + * @param aPredicate a probing predicate. + * @param aQ the fixed point 'q'. + * @param aM a frame composed of the three vectors. + */ + PlaneProbingLNeighborhood(Predicate const& aPredicate, Point const& aQ, Triangle const& aM); + + /** + * Destructor. + */ + ~PlaneProbingLNeighborhood(); + + /** + * Copy constructor. + * @param other the object to clone. + */ + PlaneProbingLNeighborhood ( const PlaneProbingLNeighborhood & other ) = delete; + + /** + * Move constructor. + * @param other the object to move. + */ + PlaneProbingLNeighborhood ( PlaneProbingLNeighborhood && other ) = delete; + + /** + * Copy assignment operator. + * @param other the object to copy. + * @return a reference on 'this'. + */ + PlaneProbingLNeighborhood & operator= ( const PlaneProbingLNeighborhood & other ) = delete; + + /** + * Move assignment operator. + * @param other the object to move. + * @return a reference on 'this'. + */ + PlaneProbingLNeighborhood & operator= ( PlaneProbingLNeighborhood && other ) = delete; + + // ----------------------- Plane-Probing services ------------------------------ + public: + + /** + * Computes the current state of the neighborhood. + * This is the function that is overloaded for the different probing modes. + * + * @return the hexagon state, see HexagonState. + */ + HexagonState hexagonState () override; + + /** + * Computes the closest candidate point, used for updating a frame in a plane probing based estimator. + * + * @return the update operation to apply. + */ + UpdateOperation closestCandidate () override; + + // ------------------------- Protected Datas ------------------------------ + protected: + + std::vector myGrids; /**< closest point and additional useful data stored at each vertex. */ + + // ------------------------- Helpers to find a closest point -------------- + protected: + + /** + * Computes the closest candidate point in a given grid identified by the + * index of the associated vertex. + * + * @param aIdx + * @return an instance of ClosestGridPoint. + */ + ClosestGridPoint closestInGrid (const Index& aIdx) const; + + /** + * Computes the candidate grid points lying in a cone given by two grid points. + * + * @param y1 a first grid point + * @param y2 a second grid point + * @param out an output iterator on grid points + */ + void candidatesInGrid (const GridPoint& y1, const GridPoint& y2, + std::back_insert_iterator > out) const; + + /** + * Finds a closest point on a given ray using a linear search. + * + * @param aRay a ray. + * @param aBound a bound that limits the search range. + * @return a closest point on the ray. + */ + GridPointOnProbingRay closestOnBoundedRayLinearWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const; + + /** + * Finds a closest point on a given ray using a binary search. + * + * @param aRay a ray. + * @param aBound a bound that limits the search range. + * @return a closest point on the ray. + */ + GridPointOnProbingRay closestOnBoundedRayLogWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const; + + /** + * Constructs an update operation from the closest candidate point. + * + * @param aClosest the closest candidate point. + * @return the update operation. + */ + UpdateOperation getOperationFromGridPoint (GridPoint const& aClosest) const; + + /** + * Update a grid after a triangle update. This procedure is called at the + * beginning of every call to hexagonState, which must prepare the computations. + * + * @param aIdx + */ + void updateGrid (const Index& aIdx); + + /** + * Returns the vector from the base to a grid point. + * + * @param aP a point on a grid. + * @return the vector. + */ + Point direction (GridPoint const& aP) const; + + // ----------------------- Interface -------------------------------------- + public: + + /** + * Writes/Displays the object on an output stream. + * @param out the output stream where the object is written. + */ + void selfDisplay ( std::ostream & out ) const; + + /** + * Checks the validity/consistency of the object. + * @return 'true' if the object is valid, 'false' otherwise. + */ + bool isValid() const; + + }; // end of class PlaneProbingLNeighborhood + + + /** + * Overloads 'operator<<' for displaying objects of class 'PlaneProbingLNeighborhood'. + * @param out the output stream where the object is written. + * @param object the object of class 'PlaneProbingLNeighborhood' to write. + * @return the output stream after the writing. + */ + template + std::ostream& + operator<< ( std::ostream & out, const PlaneProbingLNeighborhood & object ); + +} // namespace DGtal + + + /////////////////////////////////////////////////////////////////////////////// + // Includes inline functions. +#include "DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih" + +// // +/////////////////////////////////////////////////////////////////////////////// + +#endif // !defined PlaneProbingLNeighborhood_h + +#undef PlaneProbingLNeighborhood_RECURSES +#endif // else defined(PlaneProbingLNeighborhood_RECURSES) diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih new file mode 100644 index 0000000000..9857233b0d --- /dev/null +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih @@ -0,0 +1,386 @@ +/** + * 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 . + * + **/ + +/** + * @file + * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr ) + * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France + * + * @date 2024/09/16 + * + * Implementation of inline methods defined in PlaneProbingLNeighborhood.h + * + * This file is part of the DGtal library. + */ + + +////////////////////////////////////////////////////////////////////////////// +#include +////////////////////////////////////////////////////////////////////////////// + + +/////////////////////////////////////////////////////////////////////////////// +// IMPLEMENTATION of inline methods. +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// ----------------------- Standard services ------------------------------ + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +DGtal::PlaneProbingLNeighborhood:: +PlaneProbingLNeighborhood(Predicate const& aPredicate, Point const& aQ, Triangle const& aM) + : DGtal::PlaneProbingRNeighborhood(aPredicate, aQ, aM) +{ + for (int k = 0; k < 3; k++) { + myGrids.push_back(closestInGrid(k)); + } +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +DGtal::PlaneProbingLNeighborhood::~PlaneProbingLNeighborhood() +{} + +/////////////////////////////////////////////////////////////////////////////// +// ----------------------- Plane Probing services ------------------------------ + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::HexagonState +DGtal::PlaneProbingLNeighborhood::hexagonState () +{ + for (int k = 0; k < 3; k++) { + updateGrid(k); + } + + std::array state({ false, false, false, false, false, false }); + for (int k = 0; k < 3; k++) { + state[2*k] = myGrids[k].myPair.first; + state[2*k+1] = myGrids[k].myPair.second; + } + + return this->classify(state); +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::UpdateOperation +DGtal::PlaneProbingLNeighborhood:: +closestCandidate () +{ + + std::vector validGridPoints; + for (int k = 0; k < 3; k++) { + GridPoint gp = myGrids[k].myGridPoint; + if (gp.isValid()) + validGridPoints.push_back(gp); + } + + if (validGridPoints.size() == 1) { + + return getOperationFromGridPoint(validGridPoints[0]); + + } else { + // One should call hexagonState before closestCandidate, and check the return value + // to ensure that there is at least one point in the plane in the H-neighbhorhood + ASSERT(validGridPoints.size() == 2); + + if (this->isSmallest(this->relativePoint(validGridPoints[0]), + this->relativePoint(validGridPoints[1]))) { + return getOperationFromGridPoint(validGridPoints[1]); + } else { + return getOperationFromGridPoint(validGridPoints[0]); + } + } + +} + +/////////////////////////////////////////////////////////////////////////////// +// ----------------------- Helpers to find a closest point ---------------- + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::ClosestGridPoint +DGtal::PlaneProbingLNeighborhood:: +closestInGrid (const Index& aIdx) const +{ + std::vector candidates; + + Index k = aIdx; + GridPoint y1(1,0,k); //y1 = vk + m1 + GridPoint y2(0,1,k); //y2 = vk + m2 + + if (this->myPredicate(this->absolutePoint(y1))) { + if (this->myPredicate(this->absolutePoint(y2))) { + //at least 2 candidates + if (this->isSmallest(this->relativePoint(y1), this->relativePoint(y2))) { + candidates.push_back(y2); + } else { + candidates.push_back(y1); + } + + ASSERT(candidates.size() == 1); + candidatesInGrid(y1, y2, std::back_inserter(candidates)); + GridPoint closest = this->closestPointInList(candidates); + return ClosestGridPoint( closest, true, true ); + + } else { + //1 candidate + return ClosestGridPoint( y1, true, false ); + } + } else { + if (this->myPredicate(this->absolutePoint(y2))) { + //1 candidate + return ClosestGridPoint( y2, false, true ); + } else { + //no candidate + return ClosestGridPoint( GridPoint(0, 0, k), false, false ); + } + } +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +void +DGtal::PlaneProbingLNeighborhood:: +candidatesInGrid (const GridPoint& y1, const GridPoint& y2, + std::back_insert_iterator > out) const +{ + using DGtal::detail::squaredNorm; + + ASSERT( (this->myPredicate(this->absolutePoint(y1)) && + (this->myPredicate(this->absolutePoint(y2)))) ); + + Integer a = direction(y1).dot(direction(y2)); + if (a < 0) { + + GridPoint y = y1 + y2; + Integer a1 = direction(y1).dot(direction(y)); + Integer a2 = direction(y2).dot(direction(y)); + if ( (a1 < 0)||(a2 < 0) ) { + + //if a2 < 0 + GridPoint u = y1; + GridPoint w = y2; + if (a1 < 0) { + std::swap(u, w); + } + ASSERT(squaredNorm(direction(u)) > squaredNorm(direction(w))); + + Integer gamma = (-a)/squaredNorm(direction(w)); + GridPoint w2 = u + w*gamma; + GridPoint w2Next = u + w*(gamma+1); + + if (direction(w2).dot(direction(w2Next)) < 0) { + + if (this->myPredicate(this->absolutePoint(w2))) { + candidatesInGrid(w, w2, out); + } else { + //we look for a closest point on the ray u +cw, c= 0 ); + ASSERT( direction(w).dot(direction(w2Next)) >= 0 ); + + //whatever w2Next is in plane or not, + //we look for a closest point on the ray u +cw, c<=gamma+1 + GridPointOnProbingRay gr = closestOnBoundedRayLogWithPredicate(GridPointOnProbingRay(u,w.direction()),(gamma+2)); + ASSERT( gr == closestOnBoundedRayLinearWithPredicate(GridPointOnProbingRay(u,w.direction()),(gamma+2)) ); + *out++ = gr.gridPoint(); + } + + } else { //a1 and a2 are both acute (or right), + //only y has yet to be considered (in addition to y1, y2) + //(because the angle between y1 and y2 is obtuse) + if (this->myPredicate(this->absolutePoint(y))) { + *out++ = y; + } + } + } //if a >= 0, we have an acute or right angle + //and no other point has to be considered. + //(if right angle, then case of cospherical points) +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::GridPointOnProbingRay +DGtal::PlaneProbingLNeighborhood::closestOnBoundedRayLogWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const +{ + GridPointOnProbingRay Xk = aRay, Xl = aRay.next(aRay.index()+aBound); + ASSERT(this->myPredicate(this->absolutePoint(Xk))); + + // Binary search + Integer d = Xl.index() - Xk.index(); + while (d > 4) { + ASSERT(this->myPredicate(this->absolutePoint(Xk))); + + GridPointOnProbingRay Xalpha = Xk.next(Integer(d / 4)), + Xbeta = Xk.next(Integer(d / 2)), + Xgamma = Xk.next(Integer(3*d/4)); + + ASSERT(Xk.index() < Xalpha.index() && Xalpha.index() < Xbeta.index() && + Xbeta.index() < Xgamma.index() && Xgamma.index() < Xl.index()); + + if (this->myPredicate(this->absolutePoint(Xbeta)) && + this->isSmallest(this->relativePoint(Xbeta), this->relativePoint(Xgamma))) { + Xk = Xbeta; + } else if (! this->myPredicate(this->absolutePoint(Xalpha)) || + this->isSmallest(this->relativePoint(Xbeta), this->relativePoint(Xalpha))) { + Xl = Xbeta; + } else { + Xk = Xalpha; + Xl = Xgamma; + } + + d = Xl.index() - Xk.index(); + } + + return closestOnBoundedRayLinearWithPredicate(Xk, d); +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::GridPointOnProbingRay +DGtal::PlaneProbingLNeighborhood::closestOnBoundedRayLinearWithPredicate (GridPointOnProbingRay const& aRay, const Integer& aBound) const +{ + ASSERT(this->myPredicate(this->absolutePoint(aRay))); + + Integer counter = 0; + GridPointOnProbingRay previousX = aRay, currentX = previousX.next(1); + while (this->myPredicate(this->absolutePoint(currentX)) && + this->isSmallest(this->relativePoint(previousX), this->relativePoint(currentX)) && + counter < aBound) { + previousX = currentX; + currentX = previousX.next(1); + } + + ASSERT(this->myPredicate(this->absolutePoint(previousX))); + + return previousX; +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::UpdateOperation +DGtal::PlaneProbingLNeighborhood:: +getOperationFromGridPoint (GridPoint const& aP) const +{ + auto k = aP.k(); + auto d = aP.direction(); + return { { k, (k+1)%3, (k+2)%3 }, //permutation + { 1, -d.first, -d.second } //coefficients + }; +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +void +DGtal::PlaneProbingLNeighborhood:: +updateGrid (const Index& aIdx) +{ + //Let r1 and r2 be the two rays bound to the vertex of index 'aIdx' + PointOnProbingRay r1 = PointOnProbingRay({aIdx, (aIdx+1)%3, (aIdx+2)%3}); + PointOnProbingRay r2 = PointOnProbingRay({aIdx, (aIdx+2)%3, (aIdx+1)%3}); + //Check: do they are in the allowed rays? + if (this->isNeighbor(r1)) { + if (this->isNeighbor(r2)) { + //if both, + myGrids[aIdx] = closestInGrid(aIdx); + } else { + //if only r1, + myGrids[aIdx] = ClosestGridPoint( GridPoint(1, 0, aIdx), true, false ); + } + } else { + if (this->isNeighbor(r2)) { + //if only r2, + myGrids[aIdx] = ClosestGridPoint( GridPoint(0, 1, aIdx), false, true ); + } else { + //if none of them, no candidate + myGrids[aIdx] = ClosestGridPoint( GridPoint(0, 0, aIdx), false, false ); + } + } +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::Point +DGtal::PlaneProbingLNeighborhood:: +direction (GridPoint const& aP) const +{ + return aP.directionVector(this->myM); +} + +/////////////////////////////////////////////////////////////////////////////// +// Interface - public : + +/** + * Writes/Displays the object on an output stream. + * @param out the output stream where the object is written. + */ +template +inline +void +DGtal::PlaneProbingLNeighborhood::selfDisplay ( std::ostream & out ) const +{ + out << "[PlaneProbingLNeighborhood]"; +} + +/** + * Checks the validity/consistency of the object. + * @return 'true' if the object is valid, 'false' otherwise. + */ +template +inline bool DGtal::PlaneProbingLNeighborhood::isValid() const +{ + return true; +} + + + +/////////////////////////////////////////////////////////////////////////////// +// Implementation of inline functions // + +template +inline +std::ostream& +DGtal::operator<< ( std::ostream & out, + const DGtal::PlaneProbingLNeighborhood & object ) +{ + object.selfDisplay( out ); + return out; +} + +// // +/////////////////////////////////////////////////////////////////////////////// + + diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h index c4d367439f..890f7f68a3 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h @@ -73,7 +73,8 @@ namespace DGtal using Vector = Point; using Integer = typename Point::Coordinate; using Triangle = std::array; - using PointOnProbingRay = detail::PointOnProbingRay; + using Index = std::size_t; + using PointOnProbingRay = detail::PointOnProbingRay; /** * Represents a configuration of the H-neighborhood. @@ -95,7 +96,7 @@ namespace DGtal */ struct UpdateOperation { - std::array sigma; + std::array sigma; std::array coeffs; }; @@ -166,7 +167,7 @@ namespace DGtal * * @return the update operation to apply. */ - UpdateOperation closestCandidate (); + virtual UpdateOperation closestCandidate (); /** * Constructs an update operation from the closest candidate point. @@ -174,7 +175,7 @@ namespace DGtal * @param aClosest the closest candidate point. * @return the update operation. */ - virtual UpdateOperation getOperation (PointOnProbingRay const& aClosest) const; + UpdateOperation getOperation (PointOnProbingRay const& aClosest) const; /** * Classify the state of the H-neighborhood encoded as an array of 6 booleans. @@ -218,9 +219,11 @@ namespace DGtal * Computes the closest point, among a list of candidates, using a Delaunay-based criterion. * * @param aPoints the list of points. + * @tparam TPointAdapter a type for the input list of points * @return the closest point. */ - PointOnProbingRay closestPointInList (std::vector const& aPoints) const; + template + TPointAdapter closestPointInList (std::vector const& aPoints) const; /** * Tests whether a ray should be probed or not, according to the current @@ -242,20 +245,26 @@ namespace DGtal bool isSmallest (Point const& aX, Point const& aY) const; /** - * Returns the vector from the point q to the current point on the ray. + * Returns the vector from the point q to the current point, + * represented by the input object. * - * @param aRay a point on a ray. - * @return the vector from the fixed point 'q' to the current point on the ray. + * @param aPoint a representation of the current point. + * @tparam TPointAdapter a type for the representation of the input point. + * @return the vector from the fixed point 'q' to the current point. */ - Point relativePoint (PointOnProbingRay const& aRay) const; + template + Point relativePoint (TPointAdapter const& aPoint) const; /** - * Returns the current point on the ray. + * Returns the geometric realization of the input representation + * of the current point on the ray. * - * @param aRay a point on a ray. - * @return the current point on the ray. + * @param aPoint a representation of the current point. + * @tparam TPointAdapter a type for the representation of the input point. + * @return the current point. */ - Point absolutePoint (PointOnProbingRay const& aRay) const; + template + Point absolutePoint (TPointAdapter const& aPoint) const; // ------------------------- Internals ------------------------------------ private: diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih index aae16af8ee..57669ffce8 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih @@ -70,23 +70,22 @@ DGtal::PlaneProbingNeighborhood::~PlaneProbingNeighborhood() // ------------------------------------------------------------------------ template < typename TPredicate > +template < typename TPointAdapter > inline -typename DGtal::PlaneProbingNeighborhood::PointOnProbingRay -DGtal::PlaneProbingNeighborhood::closestPointInList (std::vector const& aPoints) const +TPointAdapter +DGtal::PlaneProbingNeighborhood::closestPointInList (std::vector const& aPoints) const { - const auto N = aPoints.size(); - if (N == 1) { - return aPoints[0]; + auto b = aPoints.begin(); + auto e = aPoints.end(); + ASSERT(b != e); + auto minPoint = *b; + + for (auto it = ++b; it != e; ++it) { + if (isSmallest(relativePoint(minPoint), relativePoint(*it))) { + minPoint = *it; } - - PointOnProbingRay minPoint = aPoints[N-1]; - for (Dimension k = 0; k < (Dimension)((int)N-1); ++k) { - if (isSmallest(relativePoint(minPoint), relativePoint(aPoints[k]))) { - minPoint = aPoints[k]; - } - } - - return minPoint; + } + return minPoint; } // ------------------------------------------------------------------------ @@ -177,7 +176,7 @@ getOperation (PointOnProbingRay const& aClosest) const { return { aClosest.sigma(), - { 1, -1, -aClosest.index() }, + { 1, -1, -aClosest.position() }, }; } @@ -224,22 +223,24 @@ isSmallest (Point const& aX, Point const& aY) const // ------------------------------------------------------------------------ template < typename TPredicate > +template < typename TPointAdapter > inline typename DGtal::PlaneProbingNeighborhood::Point DGtal::PlaneProbingNeighborhood:: -relativePoint (PointOnProbingRay const& aRay) const +relativePoint (TPointAdapter const& aPoint) const { - return -myM[aRay.sigma(0)] + myM[aRay.sigma(1)] + myM[aRay.sigma(2)] * aRay.index(); + return aPoint.relativePoint(myM); } // ------------------------------------------------------------------------ template < typename TPredicate > +template < typename TPointAdapter > inline typename DGtal::PlaneProbingNeighborhood::Point DGtal::PlaneProbingNeighborhood:: -absolutePoint (PointOnProbingRay const& aRay) const +absolutePoint (TPointAdapter const& aPoint) const { - return myQ + relativePoint(aRay); + return myQ + relativePoint(aPoint); } diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.h index 1e08496cf9..a19896a693 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.h @@ -73,6 +73,8 @@ namespace DGtal using Integer = typename Point::Coordinate; using HexagonState = typename PlaneProbingRNeighborhood::HexagonState; + using Index = typename PlaneProbingRNeighborhood::Index; + // ----------------------- Standard services ------------------------------ public: /** @@ -160,7 +162,7 @@ namespace DGtal * @param index an integer between 0 and 2. * @return a pair of a point on a ray and a ray. */ - std::pair candidateRay (int index) const; + std::pair candidateRay (Index const& index) const; /** * @param aPoint a point on a ray. diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.ih index 86ce29d713..c605fbe66f 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.ih @@ -124,7 +124,9 @@ DGtal::PlaneProbingR1Neighborhood::hexagonState () case 35: { std::pair candidate = candidateRay(0); - PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 2, 1, 0 }, 0) }); + std::vector shortList = { candidate.second, PointOnProbingRay({ 2, 1, 0 }, 0) }; + PointOnProbingRay closest = this->closestPointInList(shortList); + //PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 2, 1, 0 }, 0) }); this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest))); } break; @@ -132,7 +134,8 @@ DGtal::PlaneProbingR1Neighborhood::hexagonState () case 7: { std::pair candidate = candidateRay(0); - PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 1, 2, 0 }, 0) }); + std::vector shortList = { candidate.second, PointOnProbingRay({ 1, 2, 0 }, 0) }; + PointOnProbingRay closest = this->closestPointInList(shortList); this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest))); } break; @@ -140,7 +143,8 @@ DGtal::PlaneProbingR1Neighborhood::hexagonState () case 14: { std::pair candidate = candidateRay(1); - PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 0, 2, 1 }, 0) }); + std::vector shortList = { candidate.second, PointOnProbingRay({ 0, 2, 1 }, 0) }; + PointOnProbingRay closest = this->closestPointInList(shortList); this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest))); } break; @@ -148,7 +152,8 @@ DGtal::PlaneProbingR1Neighborhood::hexagonState () case 28: { std::pair candidate = candidateRay(1); - PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 2, 0, 1 }, 0) }); + std::vector shortList = { candidate.second, PointOnProbingRay({ 2, 0, 1 }, 0) }; + PointOnProbingRay closest = this->closestPointInList(shortList); this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest))); } break; @@ -156,7 +161,8 @@ DGtal::PlaneProbingR1Neighborhood::hexagonState () case 56: { std::pair candidate = candidateRay(2); - PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 1, 0, 2 }, 0) }); + std::vector shortList = { candidate.second, PointOnProbingRay({ 1, 0, 2 }, 0) }; + PointOnProbingRay closest = this->closestPointInList(shortList); this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest))); } break; @@ -164,7 +170,8 @@ DGtal::PlaneProbingR1Neighborhood::hexagonState () case 49: { std::pair candidate = candidateRay(2); - PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 0, 1, 2 }, 0) }); + std::vector shortList = { candidate.second, PointOnProbingRay({ 0, 1, 2 }, 0) }; + PointOnProbingRay closest = this->closestPointInList(shortList); this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest))); } break; @@ -208,9 +215,9 @@ template < typename TPredicate > inline std::pair::PointOnProbingRay, typename DGtal::PlaneProbingR1Neighborhood::PointOnProbingRay> -DGtal::PlaneProbingR1Neighborhood::candidateRay (int index) const +DGtal::PlaneProbingR1Neighborhood::candidateRay (Index const& index) const { - int indexM1 = (index - 1 + 3) % 3, indexM2 = (index - 2 + 3) % 3; + Index indexM1 = (index + 2) % 3, indexM2 = (index + 1) % 3; PointOnProbingRay r1({ index, indexM1, indexM2 }, 0), r2({ index, indexM2, indexM1 }, 0); @@ -316,7 +323,7 @@ DGtal::PlaneProbingR1Neighborhood::isValidIntersectSphereRay (PointO } else { int n = 15; PointOnProbingRay currentX = startingPoint; - while (res && currentX.index() < n) { + while (res && currentX.position() < n) { res = !this->isSmallest(refPoint, this->relativePoint(currentX)); currentX = currentX.next(1); } diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.h index bcb724f4fd..a81b87e653 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.h @@ -152,7 +152,8 @@ namespace DGtal * @param aRay a ray. * @return the closest point on the ray. */ - PointOnProbingRay closestPointOnRayLogWithPredicate (PointOnProbingRay const& aRay) const; + template + TPointAdapter closestPointOnRayLogWithPredicate (TPointAdapter const& aRay) const; /** * Finds the closest point on a given ray using a linear search. @@ -160,7 +161,8 @@ namespace DGtal * @param aRay a ray. * @return the closest point on the ray. */ - PointOnProbingRay closestPointOnRayLinearWithPredicate (PointOnProbingRay const& aRay) const; + template + TPointAdapter closestPointOnRayLinearWithPredicate (TPointAdapter const& aRay) const; // ------------------------- Internals ------------------------------------ private: diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.ih index ff45c5eba5..27c8230eef 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.ih @@ -90,33 +90,34 @@ DGtal::PlaneProbingRNeighborhood::hexagonState () // ------------------------------------------------------------------------ template < typename TPredicate > +template < typename TPointAdapter > inline -typename DGtal::PlaneProbingRNeighborhood::PointOnProbingRay -DGtal::PlaneProbingRNeighborhood::closestPointOnRayLogWithPredicate (PointOnProbingRay const& aRay) const +TPointAdapter +DGtal::PlaneProbingRNeighborhood::closestPointOnRayLogWithPredicate (TPointAdapter const& aRay) const { assert(this->myPredicate(this->absolutePoint(aRay))); // Exponential march - PointOnProbingRay Xk = aRay, Xl = aRay.next(1); + TPointAdapter Xk = aRay, Xl = aRay.next(1); while (this->myPredicate(this->absolutePoint(Xk)) && this->isSmallest(this->relativePoint(Xk), this->relativePoint(Xl))) { - Integer d = Xl.index() - Xk.index(); + Integer d = Xl.position() - Xk.position(); Xk = Xl; Xl = Xl.next(2 * d); } - Xk = Xk.previous(Integer((Xl.index() - Xk.index()) / 2)); + Xk = Xk.previous(Integer((Xl.position() - Xk.position()) / 2)); // Binary search - Integer d = Xl.index() - Xk.index(); + Integer d = Xl.position() - Xk.position(); while (d > 4) { assert(this->myPredicate(this->absolutePoint(Xk))); - PointOnProbingRay Xalpha = Xk.next(Integer(d / 4)), + TPointAdapter Xalpha = Xk.next(Integer(d / 4)), Xbeta = Xk.next(Integer(d / 2)), Xgamma = Xk.next(Integer(3*d/4)); - assert(Xk.index() < Xalpha.index() && Xalpha.index() < Xbeta.index() && - Xbeta.index() < Xgamma.index() && Xgamma.index() < Xl.index()); + assert(Xk.position() < Xalpha.position() && Xalpha.position() < Xbeta.position() && + Xbeta.position() < Xgamma.position() && Xgamma.position() < Xl.position()); if (this->myPredicate(this->absolutePoint(Xbeta)) && this->isSmallest(this->relativePoint(Xbeta), this->relativePoint(Xgamma))) { @@ -129,7 +130,7 @@ DGtal::PlaneProbingRNeighborhood::closestPointOnRayLogWithPredicate Xl = Xgamma; } - d = Xl.index() - Xk.index(); + d = Xl.position() - Xk.position(); } return closestPointOnRayLinearWithPredicate(Xk); @@ -137,13 +138,14 @@ DGtal::PlaneProbingRNeighborhood::closestPointOnRayLogWithPredicate // ------------------------------------------------------------------------ template < typename TPredicate > +template < typename TPointAdapter > inline -typename DGtal::PlaneProbingRNeighborhood::PointOnProbingRay -DGtal::PlaneProbingRNeighborhood::closestPointOnRayLinearWithPredicate (PointOnProbingRay const& aRay) const +TPointAdapter +DGtal::PlaneProbingRNeighborhood::closestPointOnRayLinearWithPredicate (TPointAdapter const& aRay) const { assert(this->myPredicate(this->absolutePoint(aRay))); - PointOnProbingRay previousX = aRay, currentX = previousX.next(1); + TPointAdapter previousX = aRay, currentX = previousX.next(1); while (this->myPredicate(this->absolutePoint(currentX)) && this->isSmallest(this->relativePoint(previousX), this->relativePoint(currentX))) { previousX = currentX; diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingTetrahedronEstimator.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingTetrahedronEstimator.h index 72f2fd88b4..c652f7af2c 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingTetrahedronEstimator.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingTetrahedronEstimator.h @@ -60,6 +60,7 @@ namespace DGtal H, /**< The H-neighborhood composed of 6 points on an hexagon, see \ref PlaneProbingHNeighborhood.*/ R, /**< The R-neighborhood composed of 6 rays, see \ref PlaneProbingRNeighborhood. */ R1, /**< The R-neighborhood but with an optimization to reduce the number of calls to the predicate, see \ref PlaneProbingR1Neighborhood. */ + L, /**< The L-neighborhood composed of three lattices, see \ref PlaneProbingLNeighborhood. */ }; /** @@ -85,6 +86,10 @@ namespace DGtal case ProbingMode::R1: aOs << "R1"; break; + + case ProbingMode::L: + aOs << "L"; + break; } return aOs; diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingTetrahedronEstimator.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingTetrahedronEstimator.ih index f2642646a8..77d50e7d16 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingTetrahedronEstimator.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingTetrahedronEstimator.ih @@ -36,6 +36,7 @@ #include "DGtal/geometry/surfaces/estimation/PlaneProbingHNeighborhood.h" #include "DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.h" #include "DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.h" +#include "DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h" ////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// @@ -93,6 +94,18 @@ namespace DGtal } }; + template < typename TPredicate > + struct PlaneProbingNeighborhoodSelector + { + static DGtal::PlaneProbingNeighborhood* + select(TPredicate const& aPredicate, + typename DGtal::PlaneProbingNeighborhood::Point const& aQ, + typename DGtal::PlaneProbingNeighborhood::Triangle const& aM) + { + return new DGtal::PlaneProbingLNeighborhood(aPredicate, aQ, aM); + } + }; + } // namespace detail } // namespace DGtal diff --git a/tests/geometry/surfaces/testPlaneProbingParallelepipedEstimator.cpp b/tests/geometry/surfaces/testPlaneProbingParallelepipedEstimator.cpp index a0f2d3d753..d75f5c2862 100644 --- a/tests/geometry/surfaces/testPlaneProbingParallelepipedEstimator.cpp +++ b/tests/geometry/surfaces/testPlaneProbingParallelepipedEstimator.cpp @@ -140,6 +140,32 @@ TEST_CASE( "Testing PlaneProbingParallelepipedEstimator" ) REQUIRE(nbNormals == nbOk); } + SECTION("L-algorithm should return the correct normal and a reduced basis") + { + int nbNormals = 0; + int nbOk = 0; + + for (const auto& n: NORMALS) { + for (int height = 0; height < min(int(n.normInfinity()), MAX_HEIGHT); ++height) { + ++nbNormals; + + TestPlaneProbingParallelepipedEstimator::compute + (n, height, + [&] (TestPlaneProbingParallelepipedEstimator::Estimator& estimator) { + auto estimated = estimator.compute(); + bool isReduced = estimator.isReduced(); + + if (estimated == n && isReduced) + { + nbOk++; + } + }); + } + } + + REQUIRE(nbNormals == nbOk); + } + #ifdef WITH_GMP SECTION("H-algorithm should return the correct normal with BigInteger") { diff --git a/tests/geometry/surfaces/testPlaneProbingTetrahedronEstimator.cpp b/tests/geometry/surfaces/testPlaneProbingTetrahedronEstimator.cpp index 2dc3293128..cf0ba99414 100644 --- a/tests/geometry/surfaces/testPlaneProbingTetrahedronEstimator.cpp +++ b/tests/geometry/surfaces/testPlaneProbingTetrahedronEstimator.cpp @@ -60,7 +60,7 @@ static const Z3i::Vector NORMALS[100] = { }; static const Z3i::Vector NORMALS_BIG[2] = { - {1, 59438, 82499}, {2071, 8513, 6444}, + {1, 59438, 82499}, {2071, 8513, 6444} }; template < typename Integer, ProbingMode mode > @@ -140,14 +140,48 @@ TEST_CASE("Testing PlaneProbingTetrahedronEstimator") REQUIRE(nbOk == 100); } + + SECTION("R and L algorithms should return the correct, same normal and the final basis should be reduced") + { + using Point = PointVector<3, int>; + + int nbOk = 0; + Point estimatedR, estimatedL; + bool isReducedR = false, isReducedL = false; + + for (const auto& n: NORMALS) { + TestPlaneProbingTetrahedronEstimator::compute + (n, + [&] (TestPlaneProbingTetrahedronEstimator::Estimator& estimator) { + estimatedR = estimator.compute(); + isReducedR = estimator.isReduced(); + }); + + TestPlaneProbingTetrahedronEstimator::compute + (n, + [&] (TestPlaneProbingTetrahedronEstimator::Estimator& estimator) { + estimatedL = estimator.compute(); + isReducedL = estimator.isReduced(); + }); + + if (estimatedR == n && estimatedL == estimatedR && + isReducedR && isReducedL) + { + nbOk++; + } + } + + REQUIRE(nbOk == 100); + } + #ifdef WITH_GMP - SECTION("H and R algorithm should return the correct normal, R-algorithm a reduced basis with BigInteger") + SECTION("H , R and L algorithms should return the correct normal, R and L algorithms a reduced basis with BigInteger") { using Point = PointVector<3, BigInteger>; int nbOk = 0; - Point estimatedH, estimatedR; - bool isReducedH = false, isReducedR = false; + Point estimatedH, estimatedR, estimatedL; + bool isReducedH = false, isReducedR = false, isReducedL = false; for (const auto& n: NORMALS_BIG) { TestPlaneProbingTetrahedronEstimator::compute @@ -160,11 +194,19 @@ TEST_CASE("Testing PlaneProbingTetrahedronEstimator") TestPlaneProbingTetrahedronEstimator::compute (n, [&] (TestPlaneProbingTetrahedronEstimator::Estimator& estimator) { - estimatedR = estimator.compute(); - isReducedR = estimator.isReduced(); + estimatedR = estimator.compute(); + isReducedR = estimator.isReduced(); }); - - if (estimatedH == n && estimatedR == estimatedH && !isReducedH && isReducedR) + + TestPlaneProbingTetrahedronEstimator::compute + (n, + [&] (TestPlaneProbingTetrahedronEstimator::Estimator& estimator) { + estimatedL = estimator.compute(); + isReducedL = estimator.isReduced(); + }); + + if (estimatedH == n && estimatedR == estimatedH && estimatedL == estimatedH + && !isReducedH && isReducedR && isReducedL) { nbOk++; }