From dec13b2b4f71fd377102e17f0ca2da7c0a5461e1 Mon Sep 17 00:00:00 2001 From: Tristan Roussillon Date: Tue, 17 Sep 2024 13:16:30 +0200 Subject: [PATCH 1/9] adding plane-probing L-algorithm --- src/DGtal/doc/global.bib | 18 + .../helpers/PlaneProbingEstimatorHelper.h | 208 ++++++++ .../estimation/PlaneProbingLNeighborhood.h | 369 ++++++++++++++ .../estimation/PlaneProbingLNeighborhood.ih | 465 ++++++++++++++++++ .../estimation/PlaneProbingNeighborhood.h | 4 +- .../PlaneProbingTetrahedronEstimator.h | 5 + .../PlaneProbingTetrahedronEstimator.ih | 13 + .../testPlaneProbingTetrahedronEstimator.cpp | 34 ++ 8 files changed, 1114 insertions(+), 2 deletions(-) create mode 100644 src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h create mode 100644 src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih 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/helpers/PlaneProbingEstimatorHelper.h b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h index fc361e8645..6b258893f5 100644 --- a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h +++ b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h @@ -201,6 +201,214 @@ namespace DGtal */ template < typename Integer > std::ostream& operator<< (std::ostream& aOs, PointOnProbingRay const& aRay); + + ///////////////////////////////////////////////////////////////////////////// + // template class GridPoint + /** + * Description of template class 'GridPoint'

+ * + * 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 > + 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, int 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 (int aX, int aY, int aK ) : myDir(std::make_pair(aX,aY)), myK(aK) {} + + // GridPoint (const GridPoint& other) : direction(other.direction), i(other.i) {} + + // GridPoint& operator=(const GridPoint& other) { + // if (this != &other) { + // direction = other.direction; + // i = other.i; + // } + // return *this; + // } + + // GridPoint getBase () const { + // return GridPoint(direction, 0); + // } + + /** + * 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. + */ + Integer 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 vectors. + * @tparam Point a type for vectors. + * @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); + } + + // bool operator<= (GridPoint const& other) const { + // return (myDir == other.myDir) && (i <= other.i); + // } + + /** + * 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 a scalar value + * @return the resulting point. + */ + GridPoint operator*(Integer value) const { + std::pair d = std::make_pair(myDir.first*value, + myDir.second*value); + 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 */ + Integer 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 > + std::ostream& operator<< (std::ostream& aOs, GridPoint const& aGridPoint) { + aOs << "GridPoint[k=" << aGridPoint.k() + << ", a=" << aGridPoint.direction().first + << ", b=" << aGridPoint.direction().second + << "]"; + return aOs; + } + } // namespace detail } // namespace DGtal diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h new file mode 100644 index 0000000000..c8738ce985 --- /dev/null +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h @@ -0,0 +1,369 @@ +/** + * 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/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::PlaneProbingNeighborhood + { + 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 PointOnProbingRay = typename PlaneProbingNeighborhood::PointOnProbingRay; + + using GridPoint = typename detail::GridPoint; + + // ----------------------- Inner types ------------------------------ + + // Data structure used to represent a ray on a grid. + struct GridRay + { + GridRay () = default; + + GridRay (const GridPoint& aGridPoint, + const std::pair& aDirection, + const Integer& aIdx = 0) + : myGridPoint(aGridPoint), myDirection(aDirection), myIdx(aIdx) {} + + GridRay (const GridRay& other) + : myGridPoint(other.myGridPoint), myDirection(other.myDirection), myIdx(other.myIdx) {} + + GridRay& operator=(const GridRay& other) { + if (this != &other) { + myGridPoint = other.myGridPoint; + myDirection = other.myDirection; + myIdx = other.myIdx; + } + return *this; + } + + bool operator== (GridRay const& other) const { + return ( (myGridPoint == other.myGridPoint) && + (myDirection == other.myDirection) && + (myIdx == other.myIdx) ); + } + + bool operator!= (GridRay const& other) const { + return !(*this == other); + } + + GridRay next(const int& offset) const { + return GridRay(myGridPoint, myDirection, myIdx+offset); + } + + Integer index() const { + return myIdx; + } + + GridPoint gridPoint() const { + return myGridPoint + myGridPoint.getOnSameGrid(myDirection)*myIdx; + } + + GridPoint myGridPoint; //grid point + std::pair myDirection; //direction + Integer myIdx; //index on the ray + }; + + // Data structure used to store the closest grid point associated with a vertex of the triangle. + struct ClosestGridPoint + { + ClosestGridPoint () = default; + + ClosestGridPoint (const GridPoint& aGridPoint, + const bool& first, const bool& second ) + : myGridPoint(aGridPoint), myPair(std::make_pair(first,second)) {} + + //a grid point, which can be not valid if no grid point belong to the plane + GridPoint myGridPoint; + //pair of bool values that encode the local configuration: + //at the position k, x if the predicate has returned x on v + mk, + //where v is the associated vertex + std::pair myPair; + }; + + // ----------------------- 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 int& 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; + + /** + * Computes the closest point, among a list of candidates, using a Delaunay-based criterion. + * + * @param aPoints the list of points. + * @return the closest point. + */ + GridPoint closestInList (std::vector const& aPoints) 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. + */ + GridRay closestOnRayLinearWithPredicate (GridRay 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. + */ + GridRay closestOnRayLogWithPredicate (GridRay 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 int& aIdx); + + // ----------------------- Helpers -------------------------------------- + protected: + + /** + * 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; + + /** + * Returns the vector from the point q to the current point on the grid. + * + * @param aP a point on a grid. + * @return the vector from the fixed point 'q' to the current point on the grid. + */ + Point relativePoint (GridPoint const& aP) const; + + /** + * Returns the current point on the grid. + * + * @param aP a point on a grid. + * @return the current point on the grid. + */ + Point absolutePoint (GridPoint const& aP) const; + + /** + * Returns the vector from the point q to the current point on the grid. + * + * @param aP a point on a ray. + * @return the vector from the fixed point 'q' to the current point on the grid. + */ + Point relativePoint (GridRay const& aP) const; + + /** + * Returns the current point on the grid. + * + * @param aP a point on a ray. + * @return the current point on the grid. + */ + Point absolutePoint (GridRay 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..26d72fff08 --- /dev/null +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih @@ -0,0 +1,465 @@ +/** + * 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::PlaneProbingNeighborhood(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 () +{ + // std::cout << "hexagon state" << std::endl; + 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(relativePoint(validGridPoints[0]), 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 int& aIdx) const +{ + std::vector candidates; + + int k = aIdx; + GridPoint y1(1,0,k); //y1 = vk + m1 + GridPoint y2(0,1,k); //y2 = vk + m2 + + if (this->myPredicate(absolutePoint(y1))) { + if (this->myPredicate(absolutePoint(y2))) { + //at least 2 candidates + if (this->isSmallest(relativePoint(y1), relativePoint(y2))) { + candidates.push_back(y2); + } else { + candidates.push_back(y1); + } + + ASSERT(candidates.size() == 1); + candidatesInGrid(y1, y2, std::back_inserter(candidates)); + GridPoint closest = closestInList(candidates); + return ClosestGridPoint( closest, true, true ); + + } else { + //1 candidate + return ClosestGridPoint( y1, true, false ); + } + } else { + if (this->myPredicate(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(absolutePoint(y1)) && + (this->myPredicate(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 = std::floor( (-((double) a) + /squaredNorm(direction(w))) ); + GridPoint w2 = u + w*gamma; + GridPoint w2Next = u + w*(gamma+1); + + // std::cout << "gamma=" << gamma << " "// << std::endl; + // << "w2" << direction(w2) << " " << squaredNorm(direction(w2)) << " " + // << "angle=" << (direction(w2).dot(direction(w2Next))) << " " + // << std::endl; + + if (direction(w2).dot(direction(w2Next)) < 0) { + //std::cout << "obtuse" << std::endl; + + if (this->myPredicate(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 + //std::cout << "search for a closest in plane" << std::endl; + GridRay gr = closestOnRayLogWithPredicate(GridRay(u,w.direction()),(gamma+2)); + ASSERT( gr == closestOnRayLinearWithPredicate(GridRay(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) + //std::cout << "acute" << std::endl; + if (this->myPredicate(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::GridPoint +DGtal::PlaneProbingLNeighborhood::closestInList (std::vector const& aPoints) const +{ + int N = aPoints.size(); + if (N == 1) { + return aPoints[0]; + } + + GridPoint minPoint = aPoints[N-1]; + for (int k = 0; k < N-1; ++k) { + if (this->isSmallest(relativePoint(minPoint), relativePoint(aPoints[k]))) { + minPoint = aPoints[k]; + } + } + + return minPoint; +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::GridRay +DGtal::PlaneProbingLNeighborhood::closestOnRayLogWithPredicate (GridRay const& aRay, Integer const& aBound) const +{ + GridRay Xk = aRay, Xl = aRay.next(aRay.index()+aBound); + ASSERT(this->myPredicate(absolutePoint(Xk))); + + // Binary search + Integer d = Xl.index() - Xk.index(); + while (d > 4) { + ASSERT(this->myPredicate(absolutePoint(Xk))); + + GridRay 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(absolutePoint(Xbeta)) && + this->isSmallest(relativePoint(Xbeta), relativePoint(Xgamma))) { + Xk = Xbeta; + } else if (! this->myPredicate(absolutePoint(Xalpha)) || + this->isSmallest(relativePoint(Xbeta), relativePoint(Xalpha))) { + Xl = Xbeta; + } else { + Xk = Xalpha; + Xl = Xgamma; + } + + d = Xl.index() - Xk.index(); + } + + return closestOnRayLinearWithPredicate(Xk, d); +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::GridRay +DGtal::PlaneProbingLNeighborhood::closestOnRayLinearWithPredicate (GridRay const& aRay, const Integer& aBound) const +{ + ASSERT(this->myPredicate(absolutePoint(aRay))); + + Integer counter = 0; + GridRay previousX = aRay, currentX = previousX.next(1); + while (this->myPredicate(absolutePoint(currentX)) && + this->isSmallest(relativePoint(previousX), relativePoint(currentX)) && + counter < aBound) { + previousX = currentX; + currentX = previousX.next(1); + } + + ASSERT(this->myPredicate(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 int& aIdx) +{ + // //TODO, we can be smarter if we do not recompute the closest point in some cases + // myGrids[aIdx] = closestInGrid(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 ); + } + } +} + +/////////////////////////////////////////////////////////////////////////////// +// ----------------------- Helpers ---------------------------------------- + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::Point +DGtal::PlaneProbingLNeighborhood:: +direction (GridPoint const& aP) const +{ + return aP.directionVector(this->myM); +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::Point +DGtal::PlaneProbingLNeighborhood:: +relativePoint (GridPoint const& aP) const +{ + return aP.relativePoint(this->myM); +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::Point +DGtal::PlaneProbingLNeighborhood:: +absolutePoint (GridPoint const& aP) const +{ + return this->myQ + aP.relativePoint(this->myM); +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::Point +DGtal::PlaneProbingLNeighborhood:: +relativePoint (GridRay const& aP) const +{ + return relativePoint(aP.gridPoint()); +} + +// ------------------------------------------------------------------------ +template < typename TPredicate > +inline +typename DGtal::PlaneProbingLNeighborhood::Point +DGtal::PlaneProbingLNeighborhood:: +absolutePoint (GridRay const& aP) const +{ + return absolutePoint(aP.gridPoint()); +} + + +/////////////////////////////////////////////////////////////////////////////// +// 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..997f3ff4b8 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h @@ -166,7 +166,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 +174,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. 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/testPlaneProbingTetrahedronEstimator.cpp b/tests/geometry/surfaces/testPlaneProbingTetrahedronEstimator.cpp index 2dc3293128..8712809659 100644 --- a/tests/geometry/surfaces/testPlaneProbingTetrahedronEstimator.cpp +++ b/tests/geometry/surfaces/testPlaneProbingTetrahedronEstimator.cpp @@ -140,6 +140,40 @@ 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") { From 1bb466b84ab45a0e6814a24d8e6f146fb76df809 Mon Sep 17 00:00:00 2001 From: Tristan Roussillon Date: Tue, 17 Sep 2024 18:35:55 +0200 Subject: [PATCH 2/9] moving a class from one file to another --- .../helpers/PlaneProbingEstimatorHelper.h | 144 ++++++++++++++---- .../estimation/PlaneProbingLNeighborhood.h | 103 +++++-------- .../estimation/PlaneProbingLNeighborhood.ih | 26 ++-- 3 files changed, 169 insertions(+), 104 deletions(-) diff --git a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h index 6b258893f5..6197d5aeda 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. * @@ -206,7 +211,7 @@ namespace DGtal // 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: @@ -238,7 +243,7 @@ namespace DGtal * @param aDir a pair of nonnegative integers. * @param aK an index in {0,1,2}. */ - GridPoint (std::pair const& aDir, int aK ) : myDir(aDir), myK(aK) {} + GridPoint (std::pair const& aDir, Integer aK ) : myDir(aDir), myK(aK) {} /** * Constructs a grid point from a couple of coordinates and @@ -248,22 +253,8 @@ namespace DGtal * @param aY second coordinate. * @param aK an index in {0,1,2}. */ - GridPoint (int aX, int aY, int aK ) : myDir(std::make_pair(aX,aY)), myK(aK) {} + GridPoint (Integer aX, Integer aY, Integer aK ) : myDir(std::make_pair(aX,aY)), myK(aK) {} - // GridPoint (const GridPoint& other) : direction(other.direction), i(other.i) {} - - // GridPoint& operator=(const GridPoint& other) { - // if (this != &other) { - // direction = other.direction; - // i = other.i; - // } - // return *this; - // } - - // GridPoint getBase () const { - // return GridPoint(direction, 0); - // } - /** * Returns the couple of coordinates, i.e., * the direction going from the origin to @@ -330,10 +321,6 @@ namespace DGtal return !(*this == other); } - // bool operator<= (GridPoint const& other) const { - // return (myDir == other.myDir) && (i <= other.i); - // } - /** * Returns a grid point given a couple of coordinates. * @@ -352,7 +339,7 @@ namespace DGtal */ GridPoint operator+(const GridPoint & other) const { ASSERT(myK == other.myK); - std::pair d = std::make_pair(myDir.first+other.myDir.first, + std::pair d = std::make_pair(myDir.first+other.myDir.first, myDir.second+other.myDir.second); return getOnSameGrid(d); } @@ -364,7 +351,7 @@ namespace DGtal * @return the resulting point. */ GridPoint operator*(Integer value) const { - std::pair d = std::make_pair(myDir.first*value, + std::pair d = std::make_pair(myDir.first*value, myDir.second*value); return getOnSameGrid(d); } @@ -388,8 +375,8 @@ namespace DGtal private: - std::pair myDir; /** Couple of coordinates giving a direction */ - Integer myK; /** Index of a point used as origin */ + std::pair myDir; /**< Couple of coordinates giving a direction */ + Integer myK; /**< Index of a point used as origin */ }; //end of class GridPoint @@ -409,6 +396,111 @@ namespace DGtal 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 > + 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; + } + + 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/surfaces/estimation/PlaneProbingLNeighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h index c8738ce985..ad12ae3c9e 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h @@ -76,76 +76,49 @@ namespace DGtal using HexagonState = typename PlaneProbingNeighborhood::HexagonState; using UpdateOperation = typename PlaneProbingNeighborhood::UpdateOperation; - using PointOnProbingRay = typename PlaneProbingNeighborhood::PointOnProbingRay; - - using GridPoint = typename detail::GridPoint; + using PointOnProbingRay = typename detail::PointOnProbingRay; + using GridPoint = typename detail::GridPoint; + using GridPointOnProbingRay = typename detail::GridPointOnProbingRay; - // ----------------------- Inner types ------------------------------ - - // Data structure used to represent a ray on a grid. - struct GridRay - { - GridRay () = default; - - GridRay (const GridPoint& aGridPoint, - const std::pair& aDirection, - const Integer& aIdx = 0) - : myGridPoint(aGridPoint), myDirection(aDirection), myIdx(aIdx) {} - - GridRay (const GridRay& other) - : myGridPoint(other.myGridPoint), myDirection(other.myDirection), myIdx(other.myIdx) {} - - GridRay& operator=(const GridRay& other) { - if (this != &other) { - myGridPoint = other.myGridPoint; - myDirection = other.myDirection; - myIdx = other.myIdx; - } - return *this; - } - - bool operator== (GridRay const& other) const { - return ( (myGridPoint == other.myGridPoint) && - (myDirection == other.myDirection) && - (myIdx == other.myIdx) ); - } - - bool operator!= (GridRay const& other) const { - return !(*this == other); - } - - GridRay next(const int& offset) const { - return GridRay(myGridPoint, myDirection, myIdx+offset); - } - - Integer index() const { - return myIdx; - } - - GridPoint gridPoint() const { - return myGridPoint + myGridPoint.getOnSameGrid(myDirection)*myIdx; - } - - GridPoint myGridPoint; //grid point - std::pair myDirection; //direction - Integer myIdx; //index on the ray - }; - - // Data structure used to store the closest grid point associated with a vertex of the triangle. + // ----------------------- 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& first, const bool& second ) : myGridPoint(aGridPoint), myPair(std::make_pair(first,second)) {} - //a grid point, which can be not valid if no grid point belong to the plane - GridPoint myGridPoint; - //pair of bool values that encode the local configuration: - //at the position k, x if the predicate has returned x on v + mk, - //where v is the associated vertex - std::pair myPair; + 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 ------------------------------ @@ -255,7 +228,7 @@ namespace DGtal * @param aBound a bound that limits the search range. * @return a closest point on the ray. */ - GridRay closestOnRayLinearWithPredicate (GridRay const& aRay, Integer const& aBound) const; + GridPointOnProbingRay closestOnRayLinearWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const; /** * Finds a closest point on a given ray using a binary search. @@ -264,7 +237,7 @@ namespace DGtal * @param aBound a bound that limits the search range. * @return a closest point on the ray. */ - GridRay closestOnRayLogWithPredicate (GridRay const& aRay, Integer const& aBound) const; + GridPointOnProbingRay closestOnRayLogWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const; /** * Constructs an update operation from the closest candidate point. @@ -315,7 +288,7 @@ namespace DGtal * @param aP a point on a ray. * @return the vector from the fixed point 'q' to the current point on the grid. */ - Point relativePoint (GridRay const& aP) const; + Point relativePoint (GridPointOnProbingRay const& aP) const; /** * Returns the current point on the grid. @@ -323,7 +296,7 @@ namespace DGtal * @param aP a point on a ray. * @return the current point on the grid. */ - Point absolutePoint (GridRay const& aP) const; + Point absolutePoint (GridPointOnProbingRay const& aP) const; // ----------------------- Interface -------------------------------------- public: diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih index 26d72fff08..2dfa24e6fe 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih @@ -205,8 +205,8 @@ candidatesInGrid (const GridPoint& y1, const GridPoint& y2, } else { //we look for a closest point on the ray u +cw, c::closestInList (std::vector inline -typename DGtal::PlaneProbingLNeighborhood::GridRay -DGtal::PlaneProbingLNeighborhood::closestOnRayLogWithPredicate (GridRay const& aRay, Integer const& aBound) const +typename DGtal::PlaneProbingLNeighborhood::GridPointOnProbingRay +DGtal::PlaneProbingLNeighborhood::closestOnRayLogWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const { - GridRay Xk = aRay, Xl = aRay.next(aRay.index()+aBound); + GridPointOnProbingRay Xk = aRay, Xl = aRay.next(aRay.index()+aBound); ASSERT(this->myPredicate(absolutePoint(Xk))); // Binary search @@ -272,7 +272,7 @@ DGtal::PlaneProbingLNeighborhood::closestOnRayLogWithPredicate (Grid while (d > 4) { ASSERT(this->myPredicate(absolutePoint(Xk))); - GridRay Xalpha = Xk.next(Integer(d / 4)), + GridPointOnProbingRay Xalpha = Xk.next(Integer(d / 4)), Xbeta = Xk.next(Integer(d / 2)), Xgamma = Xk.next(Integer(3*d/4)); @@ -299,13 +299,13 @@ DGtal::PlaneProbingLNeighborhood::closestOnRayLogWithPredicate (Grid // ------------------------------------------------------------------------ template < typename TPredicate > inline -typename DGtal::PlaneProbingLNeighborhood::GridRay -DGtal::PlaneProbingLNeighborhood::closestOnRayLinearWithPredicate (GridRay const& aRay, const Integer& aBound) const +typename DGtal::PlaneProbingLNeighborhood::GridPointOnProbingRay +DGtal::PlaneProbingLNeighborhood::closestOnRayLinearWithPredicate (GridPointOnProbingRay const& aRay, const Integer& aBound) const { ASSERT(this->myPredicate(absolutePoint(aRay))); Integer counter = 0; - GridRay previousX = aRay, currentX = previousX.next(1); + GridPointOnProbingRay previousX = aRay, currentX = previousX.next(1); while (this->myPredicate(absolutePoint(currentX)) && this->isSmallest(relativePoint(previousX), relativePoint(currentX)) && counter < aBound) { @@ -403,7 +403,7 @@ template < typename TPredicate > inline typename DGtal::PlaneProbingLNeighborhood::Point DGtal::PlaneProbingLNeighborhood:: -relativePoint (GridRay const& aP) const +relativePoint (GridPointOnProbingRay const& aP) const { return relativePoint(aP.gridPoint()); } @@ -413,7 +413,7 @@ template < typename TPredicate > inline typename DGtal::PlaneProbingLNeighborhood::Point DGtal::PlaneProbingLNeighborhood:: -absolutePoint (GridRay const& aP) const +absolutePoint (GridPointOnProbingRay const& aP) const { return absolutePoint(aP.gridPoint()); } From bb1b477c4205d4e27d0d210767c59d835e277d56 Mon Sep 17 00:00:00 2001 From: Tristan Roussillon Date: Wed, 18 Sep 2024 11:14:23 +0200 Subject: [PATCH 3/9] cleaning and refactoring with template methods --- .../helpers/PlaneProbingEstimatorHelper.h | 29 +++- .../estimation/PlaneProbingLNeighborhood.h | 86 ++++++----- .../estimation/PlaneProbingLNeighborhood.ih | 136 ++++-------------- .../estimation/PlaneProbingNeighborhood.h | 26 ++-- .../estimation/PlaneProbingNeighborhood.ih | 37 ++--- .../estimation/PlaneProbingR1Neighborhood.ih | 19 ++- .../estimation/PlaneProbingRNeighborhood.h | 6 +- .../estimation/PlaneProbingRNeighborhood.ih | 16 ++- 8 files changed, 158 insertions(+), 197 deletions(-) diff --git a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h index 6197d5aeda..3e11d893b2 100644 --- a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h +++ b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h @@ -152,6 +152,18 @@ namespace DGtal */ Integer const& index () 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 { + return -aM[mySigma[0]] + aM[mySigma[1]] + aM[mySigma[2]] * myIndex; + } + /** * Equality test between two rays: the internal permutations and * indices must be the same. @@ -291,8 +303,8 @@ namespace DGtal /** * Returns the geometric realization of this grid point. * - * @param aM an array of three vectors. - * @tparam Point a type for vectors. + * @param aM an array of three points. + * @tparam Point a type for points. * @return the computed point. */ template < typename Point > @@ -494,6 +506,19 @@ namespace DGtal 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 */ diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h index ad12ae3c9e..9287023718 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h @@ -45,6 +45,7 @@ #include #include "DGtal/base/Common.h" #include "DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h" +#include "DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.h" #include "DGtal/kernel/CPointPredicate.h" ////////////////////////////////////////////////////////////////////////////// @@ -61,7 +62,7 @@ namespace DGtal * \tparam TPredicate the probing predicate, a model of concepts::CPointPredicate. */ template - class PlaneProbingLNeighborhood: public DGtal::PlaneProbingNeighborhood + class PlaneProbingLNeighborhood: public DGtal::PlaneProbingRNeighborhood { BOOST_CONCEPT_ASSERT((DGtal::concepts::CPointPredicate)); @@ -212,14 +213,6 @@ namespace DGtal */ void candidatesInGrid (const GridPoint& y1, const GridPoint& y2, std::back_insert_iterator > out) const; - - /** - * Computes the closest point, among a list of candidates, using a Delaunay-based criterion. - * - * @param aPoints the list of points. - * @return the closest point. - */ - GridPoint closestInList (std::vector const& aPoints) const; /** * Finds a closest point on a given ray using a linear search. @@ -228,7 +221,7 @@ namespace DGtal * @param aBound a bound that limits the search range. * @return a closest point on the ray. */ - GridPointOnProbingRay closestOnRayLinearWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const; + GridPointOnProbingRay closestOnBoundedRayLinearWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const; /** * Finds a closest point on a given ray using a binary search. @@ -237,7 +230,7 @@ namespace DGtal * @param aBound a bound that limits the search range. * @return a closest point on the ray. */ - GridPointOnProbingRay closestOnRayLogWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const; + GridPointOnProbingRay closestOnBoundedRayLogWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const; /** * Constructs an update operation from the closest candidate point. @@ -254,9 +247,6 @@ namespace DGtal * @param aIdx */ void updateGrid (const int& aIdx); - - // ----------------------- Helpers -------------------------------------- - protected: /** * Returns the vector from the base to a grid point. @@ -265,38 +255,42 @@ namespace DGtal * @return the vector. */ Point direction (GridPoint const& aP) const; - - /** - * Returns the vector from the point q to the current point on the grid. - * - * @param aP a point on a grid. - * @return the vector from the fixed point 'q' to the current point on the grid. - */ - Point relativePoint (GridPoint const& aP) const; - - /** - * Returns the current point on the grid. - * - * @param aP a point on a grid. - * @return the current point on the grid. - */ - Point absolutePoint (GridPoint const& aP) const; - - /** - * Returns the vector from the point q to the current point on the grid. - * - * @param aP a point on a ray. - * @return the vector from the fixed point 'q' to the current point on the grid. - */ - Point relativePoint (GridPointOnProbingRay const& aP) const; - - /** - * Returns the current point on the grid. - * - * @param aP a point on a ray. - * @return the current point on the grid. - */ - Point absolutePoint (GridPointOnProbingRay const& aP) const; + + // // ----------------------- Helpers -------------------------------------- + // protected: + + + // /** + // * Returns the vector from the point q to the current point on the grid. + // * + // * @param aP a point on a grid. + // * @return the vector from the fixed point 'q' to the current point on the grid. + // */ + // Point relativePoint (GridPoint const& aP) const; + + // /** + // * Returns the current point on the grid. + // * + // * @param aP a point on a grid. + // * @return the current point on the grid. + // */ + // Point absolutePoint (GridPoint const& aP) const; + + // /** + // * Returns the vector from the point q to the current point on the grid. + // * + // * @param aP a point on a ray. + // * @return the vector from the fixed point 'q' to the current point on the grid. + // */ + // Point relativePoint (GridPointOnProbingRay const& aP) const; + + // /** + // * Returns the current point on the grid. + // * + // * @param aP a point on a ray. + // * @return the current point on the grid. + // */ + // Point absolutePoint (GridPointOnProbingRay const& aP) const; // ----------------------- Interface -------------------------------------- public: diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih index 2dfa24e6fe..720ac10aba 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih @@ -44,7 +44,7 @@ template < typename TPredicate > inline DGtal::PlaneProbingLNeighborhood:: PlaneProbingLNeighborhood(Predicate const& aPredicate, Point const& aQ, Triangle const& aM) - : DGtal::PlaneProbingNeighborhood(aPredicate, aQ, aM) + : DGtal::PlaneProbingRNeighborhood(aPredicate, aQ, aM) { for (int k = 0; k < 3; k++) { myGrids.push_back(closestInGrid(k)); @@ -66,7 +66,6 @@ inline typename DGtal::PlaneProbingLNeighborhood::HexagonState DGtal::PlaneProbingLNeighborhood::hexagonState () { - // std::cout << "hexagon state" << std::endl; for (int k = 0; k < 3; k++) { updateGrid(k); } @@ -104,7 +103,8 @@ closestCandidate () // to ensure that there is at least one point in the plane in the H-neighbhorhood ASSERT(validGridPoints.size() == 2); - if (this->isSmallest(relativePoint(validGridPoints[0]), relativePoint(validGridPoints[1]))) { + if (this->isSmallest(this->relativePoint(validGridPoints[0]), + this->relativePoint(validGridPoints[1]))) { return getOperationFromGridPoint(validGridPoints[1]); } else { return getOperationFromGridPoint(validGridPoints[0]); @@ -129,10 +129,10 @@ closestInGrid (const int& aIdx) const GridPoint y1(1,0,k); //y1 = vk + m1 GridPoint y2(0,1,k); //y2 = vk + m2 - if (this->myPredicate(absolutePoint(y1))) { - if (this->myPredicate(absolutePoint(y2))) { + if (this->myPredicate(this->absolutePoint(y1))) { + if (this->myPredicate(this->absolutePoint(y2))) { //at least 2 candidates - if (this->isSmallest(relativePoint(y1), relativePoint(y2))) { + if (this->isSmallest(this->relativePoint(y1), this->relativePoint(y2))) { candidates.push_back(y2); } else { candidates.push_back(y1); @@ -140,7 +140,7 @@ closestInGrid (const int& aIdx) const ASSERT(candidates.size() == 1); candidatesInGrid(y1, y2, std::back_inserter(candidates)); - GridPoint closest = closestInList(candidates); + GridPoint closest = this->closestPointInList(candidates); return ClosestGridPoint( closest, true, true ); } else { @@ -148,7 +148,7 @@ closestInGrid (const int& aIdx) const return ClosestGridPoint( y1, true, false ); } } else { - if (this->myPredicate(absolutePoint(y2))) { + if (this->myPredicate(this->absolutePoint(y2))) { //1 candidate return ClosestGridPoint( y2, false, true ); } else { @@ -168,8 +168,8 @@ candidatesInGrid (const GridPoint& y1, const GridPoint& y2, { using DGtal::detail::squaredNorm; - ASSERT( (this->myPredicate(absolutePoint(y1)) && - (this->myPredicate(absolutePoint(y2)))) ); + ASSERT( (this->myPredicate(this->absolutePoint(y1)) && + (this->myPredicate(this->absolutePoint(y2)))) ); Integer a = direction(y1).dot(direction(y2)); if (a < 0) { @@ -192,43 +192,33 @@ candidatesInGrid (const GridPoint& y1, const GridPoint& y2, GridPoint w2 = u + w*gamma; GridPoint w2Next = u + w*(gamma+1); - // std::cout << "gamma=" << gamma << " "// << std::endl; - // << "w2" << direction(w2) << " " << squaredNorm(direction(w2)) << " " - // << "angle=" << (direction(w2).dot(direction(w2Next))) << " " - // << std::endl; - if (direction(w2).dot(direction(w2Next)) < 0) { - //std::cout << "obtuse" << std::endl; - if (this->myPredicate(absolutePoint(w2))) { + 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 - //std::cout << "search for a closest in plane" << std::endl; - GridPointOnProbingRay gr = closestOnRayLogWithPredicate(GridPointOnProbingRay(u,w.direction()),(gamma+2)); - ASSERT( gr == closestOnRayLinearWithPredicate(GridPointOnProbingRay(u,w.direction()),(gamma+2)) ); + 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) - //std::cout << "acute" << std::endl; - if (this->myPredicate(absolutePoint(y))) { + if (this->myPredicate(this->absolutePoint(y))) { *out++ = y; } } @@ -237,40 +227,19 @@ candidatesInGrid (const GridPoint& y1, const GridPoint& y2, //(if right angle, then case of cospherical points) } -// ------------------------------------------------------------------------ -template < typename TPredicate > -inline -typename DGtal::PlaneProbingLNeighborhood::GridPoint -DGtal::PlaneProbingLNeighborhood::closestInList (std::vector const& aPoints) const -{ - int N = aPoints.size(); - if (N == 1) { - return aPoints[0]; - } - - GridPoint minPoint = aPoints[N-1]; - for (int k = 0; k < N-1; ++k) { - if (this->isSmallest(relativePoint(minPoint), relativePoint(aPoints[k]))) { - minPoint = aPoints[k]; - } - } - - return minPoint; -} - // ------------------------------------------------------------------------ template < typename TPredicate > inline typename DGtal::PlaneProbingLNeighborhood::GridPointOnProbingRay -DGtal::PlaneProbingLNeighborhood::closestOnRayLogWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const +DGtal::PlaneProbingLNeighborhood::closestOnBoundedRayLogWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const { GridPointOnProbingRay Xk = aRay, Xl = aRay.next(aRay.index()+aBound); - ASSERT(this->myPredicate(absolutePoint(Xk))); + ASSERT(this->myPredicate(this->absolutePoint(Xk))); // Binary search Integer d = Xl.index() - Xk.index(); while (d > 4) { - ASSERT(this->myPredicate(absolutePoint(Xk))); + ASSERT(this->myPredicate(this->absolutePoint(Xk))); GridPointOnProbingRay Xalpha = Xk.next(Integer(d / 4)), Xbeta = Xk.next(Integer(d / 2)), @@ -279,11 +248,11 @@ DGtal::PlaneProbingLNeighborhood::closestOnRayLogWithPredicate (Grid ASSERT(Xk.index() < Xalpha.index() && Xalpha.index() < Xbeta.index() && Xbeta.index() < Xgamma.index() && Xgamma.index() < Xl.index()); - if (this->myPredicate(absolutePoint(Xbeta)) && - this->isSmallest(relativePoint(Xbeta), relativePoint(Xgamma))) { + if (this->myPredicate(this->absolutePoint(Xbeta)) && + this->isSmallest(this->relativePoint(Xbeta), this->relativePoint(Xgamma))) { Xk = Xbeta; - } else if (! this->myPredicate(absolutePoint(Xalpha)) || - this->isSmallest(relativePoint(Xbeta), relativePoint(Xalpha))) { + } else if (! this->myPredicate(this->absolutePoint(Xalpha)) || + this->isSmallest(this->relativePoint(Xbeta), this->relativePoint(Xalpha))) { Xl = Xbeta; } else { Xk = Xalpha; @@ -293,27 +262,27 @@ DGtal::PlaneProbingLNeighborhood::closestOnRayLogWithPredicate (Grid d = Xl.index() - Xk.index(); } - return closestOnRayLinearWithPredicate(Xk, d); + return closestOnBoundedRayLinearWithPredicate(Xk, d); } // ------------------------------------------------------------------------ template < typename TPredicate > inline typename DGtal::PlaneProbingLNeighborhood::GridPointOnProbingRay -DGtal::PlaneProbingLNeighborhood::closestOnRayLinearWithPredicate (GridPointOnProbingRay const& aRay, const Integer& aBound) const +DGtal::PlaneProbingLNeighborhood::closestOnBoundedRayLinearWithPredicate (GridPointOnProbingRay const& aRay, const Integer& aBound) const { - ASSERT(this->myPredicate(absolutePoint(aRay))); + ASSERT(this->myPredicate(this->absolutePoint(aRay))); Integer counter = 0; GridPointOnProbingRay previousX = aRay, currentX = previousX.next(1); - while (this->myPredicate(absolutePoint(currentX)) && - this->isSmallest(relativePoint(previousX), relativePoint(currentX)) && + 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(absolutePoint(previousX))); + ASSERT(this->myPredicate(this->absolutePoint(previousX))); return previousX; } @@ -339,9 +308,6 @@ void DGtal::PlaneProbingLNeighborhood:: updateGrid (const int& aIdx) { - // //TODO, we can be smarter if we do not recompute the closest point in some cases - // myGrids[aIdx] = closestInGrid(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}); @@ -365,9 +331,6 @@ updateGrid (const int& aIdx) } } -/////////////////////////////////////////////////////////////////////////////// -// ----------------------- Helpers ---------------------------------------- - // ------------------------------------------------------------------------ template < typename TPredicate > inline @@ -378,47 +341,6 @@ direction (GridPoint const& aP) const return aP.directionVector(this->myM); } -// ------------------------------------------------------------------------ -template < typename TPredicate > -inline -typename DGtal::PlaneProbingLNeighborhood::Point -DGtal::PlaneProbingLNeighborhood:: -relativePoint (GridPoint const& aP) const -{ - return aP.relativePoint(this->myM); -} - -// ------------------------------------------------------------------------ -template < typename TPredicate > -inline -typename DGtal::PlaneProbingLNeighborhood::Point -DGtal::PlaneProbingLNeighborhood:: -absolutePoint (GridPoint const& aP) const -{ - return this->myQ + aP.relativePoint(this->myM); -} - -// ------------------------------------------------------------------------ -template < typename TPredicate > -inline -typename DGtal::PlaneProbingLNeighborhood::Point -DGtal::PlaneProbingLNeighborhood:: -relativePoint (GridPointOnProbingRay const& aP) const -{ - return relativePoint(aP.gridPoint()); -} - -// ------------------------------------------------------------------------ -template < typename TPredicate > -inline -typename DGtal::PlaneProbingLNeighborhood::Point -DGtal::PlaneProbingLNeighborhood:: -absolutePoint (GridPointOnProbingRay const& aP) const -{ - return absolutePoint(aP.gridPoint()); -} - - /////////////////////////////////////////////////////////////////////////////// // Interface - public : diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h index 997f3ff4b8..b21985d3f9 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h @@ -218,9 +218,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 +244,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..7e0050d450 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; b != e; ++b) { + 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; } // ------------------------------------------------------------------------ @@ -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.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.ih index 86ce29d713..cb0957b275 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; 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..3affa76230 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.ih @@ -90,14 +90,15 @@ 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(); @@ -111,7 +112,7 @@ DGtal::PlaneProbingRNeighborhood::closestPointOnRayLogWithPredicate 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)); @@ -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; From 9271b4e38c37d3df1c3c159c4cb11ab51f534a6b Mon Sep 17 00:00:00 2001 From: Tristan Roussillon Date: Wed, 18 Sep 2024 15:51:30 +0200 Subject: [PATCH 4/9] bug fix + test with BigInt --- .../helpers/PlaneProbingEstimatorHelper.h | 8 ++-- .../estimation/PlaneProbingLNeighborhood.h | 42 ++----------------- .../estimation/PlaneProbingLNeighborhood.ih | 3 +- .../estimation/PlaneProbingNeighborhood.ih | 2 +- .../testPlaneProbingTetrahedronEstimator.cpp | 24 +++++++---- 5 files changed, 25 insertions(+), 54 deletions(-) diff --git a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h index 3e11d893b2..aa855cfc68 100644 --- a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h +++ b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h @@ -255,7 +255,7 @@ namespace DGtal * @param aDir a pair of nonnegative integers. * @param aK an index in {0,1,2}. */ - GridPoint (std::pair const& aDir, Integer aK ) : myDir(aDir), myK(aK) {} + GridPoint (std::pair const& aDir, int aK ) : myDir(aDir), myK(aK) {} /** * Constructs a grid point from a couple of coordinates and @@ -265,7 +265,7 @@ namespace DGtal * @param aY second coordinate. * @param aK an index in {0,1,2}. */ - GridPoint (Integer aX, Integer aY, Integer aK ) : myDir(std::make_pair(aX,aY)), myK(aK) {} + GridPoint (Integer aX, Integer aY, int aK ) : myDir(std::make_pair(aX,aY)), myK(aK) {} /** * Returns the couple of coordinates, i.e., @@ -283,7 +283,7 @@ namespace DGtal * * @return the index. */ - Integer k() const { + int k() const { return myK; } @@ -388,7 +388,7 @@ namespace DGtal private: std::pair myDir; /**< Couple of coordinates giving a direction */ - Integer myK; /**< Index of a point used as origin */ + int myK; /**< Index of a point used as origin */ }; //end of class GridPoint diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h index 9287023718..3aaf2423c4 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h @@ -77,9 +77,9 @@ namespace DGtal using HexagonState = typename PlaneProbingNeighborhood::HexagonState; using UpdateOperation = typename PlaneProbingNeighborhood::UpdateOperation; - using PointOnProbingRay = typename detail::PointOnProbingRay; - using GridPoint = typename detail::GridPoint; - using GridPointOnProbingRay = typename detail::GridPointOnProbingRay; + using PointOnProbingRay = typename detail::PointOnProbingRay; + using GridPoint = typename detail::GridPoint; + using GridPointOnProbingRay = typename detail::GridPointOnProbingRay; // ----------------------- Internal type ------------------------------- private: @@ -255,42 +255,6 @@ namespace DGtal * @return the vector. */ Point direction (GridPoint const& aP) const; - - // // ----------------------- Helpers -------------------------------------- - // protected: - - - // /** - // * Returns the vector from the point q to the current point on the grid. - // * - // * @param aP a point on a grid. - // * @return the vector from the fixed point 'q' to the current point on the grid. - // */ - // Point relativePoint (GridPoint const& aP) const; - - // /** - // * Returns the current point on the grid. - // * - // * @param aP a point on a grid. - // * @return the current point on the grid. - // */ - // Point absolutePoint (GridPoint const& aP) const; - - // /** - // * Returns the vector from the point q to the current point on the grid. - // * - // * @param aP a point on a ray. - // * @return the vector from the fixed point 'q' to the current point on the grid. - // */ - // Point relativePoint (GridPointOnProbingRay const& aP) const; - - // /** - // * Returns the current point on the grid. - // * - // * @param aP a point on a ray. - // * @return the current point on the grid. - // */ - // Point absolutePoint (GridPointOnProbingRay const& aP) const; // ----------------------- Interface -------------------------------------- public: diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih index 720ac10aba..c3dc05324c 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih @@ -187,8 +187,7 @@ candidatesInGrid (const GridPoint& y1, const GridPoint& y2, } ASSERT(squaredNorm(direction(u)) > squaredNorm(direction(w))); - Integer gamma = std::floor( (-((double) a) - /squaredNorm(direction(w))) ); + Integer gamma = (-a)/squaredNorm(direction(w)); GridPoint w2 = u + w*gamma; GridPoint w2Next = u + w*(gamma+1); diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih index 7e0050d450..7684acfdb6 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih @@ -80,7 +80,7 @@ DGtal::PlaneProbingNeighborhood::closestPointInList (std::vector @@ -175,13 +175,13 @@ TEST_CASE("Testing PlaneProbingTetrahedronEstimator") } #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 @@ -194,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++; } From 19190814b3a6fde7428a22c6f9969919f68432fb Mon Sep 17 00:00:00 2001 From: Tristan Roussillon Date: Wed, 18 Sep 2024 16:20:25 +0200 Subject: [PATCH 5/9] checking tests/examples/doc for L-alg. --- ...xamplePlaneProbingTetrahedronEstimator.cpp | 2 +- src/DGtal/geometry/doc/modulePlaneProbing.dox | 2 +- src/DGtal/geometry/doc/packageGeometry.dox | 2 +- ...estPlaneProbingParallelepipedEstimator.cpp | 26 +++++++++++++++++++ 4 files changed, 29 insertions(+), 3 deletions(-) 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/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/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") { From df9e0cd81ddef378cbc4d8afb27c6b673627f669 Mon Sep 17 00:00:00 2001 From: Tristan Roussillon Date: Thu, 19 Sep 2024 17:33:16 +0200 Subject: [PATCH 6/9] changelog updated --- ChangeLog.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/ChangeLog.md b/ChangeLog.md index 4d01f5694a..fb39f511ff 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,6 +1,12 @@ # 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 + - *Geometry* - Bug fix in ArithmeticalDSSComputerOnSurfels (Tristan Roussillon, [#1742](https://github.com/DGtal-team/DGtal/pull/1742)) From 96fc91c633bffb539f327fdebc53a4ebae5d0f52 Mon Sep 17 00:00:00 2001 From: Tristan Roussillon Date: Fri, 20 Sep 2024 09:10:01 +0200 Subject: [PATCH 7/9] fixing doxygen warnings --- src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h | 8 ++++---- .../surfaces/estimation/PlaneProbingLNeighborhood.h | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h index aa855cfc68..db61f558c1 100644 --- a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h +++ b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h @@ -359,12 +359,12 @@ namespace DGtal /** * Scales this grid point by a scalar (as if it was vector). * - * @param a scalar value + * @param aValue a scalar value * @return the resulting point. */ - GridPoint operator*(Integer value) const { - std::pair d = std::make_pair(myDir.first*value, - myDir.second*value); + GridPoint operator*(Integer aValue) const { + std::pair d = std::make_pair(myDir.first*aValue, + myDir.second*aValue); return getOnSameGrid(d); } diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h index 3aaf2423c4..e248100230 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h @@ -112,8 +112,8 @@ namespace DGtal * @param aSecond the second boolean value */ ClosestGridPoint (const GridPoint& aGridPoint, - const bool& first, const bool& second ) - : myGridPoint(aGridPoint), myPair(std::make_pair(first,second)) {} + 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 */ From 753f611640566d28cb739f633142d2a4b8fdfba0 Mon Sep 17 00:00:00 2001 From: Tristan Roussillon Date: Thu, 3 Oct 2024 13:35:21 +0200 Subject: [PATCH 8/9] index type for the index value that identifies a triangle vertex --- .../helpers/PlaneProbingEstimatorHelper.h | 38 +++++---- .../helpers/PlaneProbingEstimatorHelper.ih | 80 ++++++++++--------- .../estimation/PlaneProbingLNeighborhood.h | 11 +-- .../estimation/PlaneProbingLNeighborhood.ih | 6 +- .../estimation/PlaneProbingNeighborhood.h | 5 +- .../estimation/PlaneProbingNeighborhood.ih | 2 +- .../estimation/PlaneProbingR1Neighborhood.h | 4 +- .../estimation/PlaneProbingR1Neighborhood.ih | 6 +- .../estimation/PlaneProbingRNeighborhood.ih | 12 +-- 9 files changed, 87 insertions(+), 77 deletions(-) diff --git a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h index db61f558c1..bd4bc36e41 100644 --- a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h +++ b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h @@ -108,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: /** @@ -127,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). @@ -145,12 +145,12 @@ 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. @@ -160,9 +160,7 @@ namespace DGtal * @return the computed point. */ template < typename Point > - Point relativePoint (std::array const& aM) const { - return -aM[mySigma[0]] + aM[mySigma[1]] + aM[mySigma[2]] * myIndex; - } + Point relativePoint (std::array const& aM) const; /** * Equality test between two rays: the internal permutations and @@ -206,7 +204,7 @@ namespace DGtal private: Permutation mySigma; /**< The permutation. */ - Integer myIndex; /**< The index. */ + Integer myPosition; /**< The index. */ }; // end of class PointOnProbingRay /** @@ -236,7 +234,7 @@ 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 GridPoint { BOOST_CONCEPT_ASSERT(( concepts::CInteger ) ); @@ -255,7 +253,7 @@ namespace DGtal * @param aDir a pair of nonnegative integers. * @param aK an index in {0,1,2}. */ - GridPoint (std::pair const& aDir, int aK ) : myDir(aDir), myK(aK) {} + GridPoint (std::pair const& aDir, Index const& aK ) : myDir(aDir), myK(aK) {} /** * Constructs a grid point from a couple of coordinates and @@ -265,7 +263,7 @@ namespace DGtal * @param aY second coordinate. * @param aK an index in {0,1,2}. */ - GridPoint (Integer aX, Integer aY, int aK ) : myDir(std::make_pair(aX,aY)), myK(aK) {} + 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., @@ -283,7 +281,7 @@ namespace DGtal * * @return the index. */ - int k() const { + Index k() const { return myK; } @@ -388,7 +386,7 @@ namespace DGtal private: std::pair myDir; /**< Couple of coordinates giving a direction */ - int myK; /**< Index of a point used as origin */ + Index myK; /**< Index of a point used as origin */ }; //end of class GridPoint @@ -421,7 +419,7 @@ 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 GridPointOnProbingRay { BOOST_CONCEPT_ASSERT(( concepts::CInteger ) ); @@ -439,7 +437,7 @@ namespace DGtal * @param aDirection direction of the ray * @param aIdx index of the grid point along the ray */ - GridPointOnProbingRay (const GridPoint& aGridPoint, + GridPointOnProbingRay (const GridPoint& aGridPoint, const std::pair& aDirection, const Integer& aIdx = 0) : myOrigin(aGridPoint), myDirection(aDirection), myIdx(aIdx) {} @@ -502,7 +500,7 @@ namespace DGtal /** * @return the current grid point as an instance of GridPoint. */ - GridPoint gridPoint() const { + GridPoint gridPoint() const { return myOrigin + myOrigin.getOnSameGrid(myDirection)*myIdx; } @@ -520,7 +518,7 @@ namespace DGtal private: - GridPoint myOrigin; /**< starting point of the ray */ + GridPoint myOrigin; /**< starting point of the ray */ std::pair myDirection; /**< direction of the ray */ Integer myIdx; /**< index of the point along the ray */ 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 index e248100230..ba0947efa7 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h @@ -77,9 +77,10 @@ namespace DGtal using HexagonState = typename PlaneProbingNeighborhood::HexagonState; using UpdateOperation = typename PlaneProbingNeighborhood::UpdateOperation; - using PointOnProbingRay = typename detail::PointOnProbingRay; - using GridPoint = typename detail::GridPoint; - using GridPointOnProbingRay = typename detail::GridPointOnProbingRay; + using Index = typename PlaneProbingNeighborhood::Index; + using PointOnProbingRay = typename PlaneProbingNeighborhood::PointOnProbingRay; + using GridPoint = typename detail::GridPoint; + using GridPointOnProbingRay = typename detail::GridPointOnProbingRay; // ----------------------- Internal type ------------------------------- private: @@ -202,7 +203,7 @@ namespace DGtal * @param aIdx * @return an instance of ClosestGridPoint. */ - ClosestGridPoint closestInGrid (const int& aIdx) const; + ClosestGridPoint closestInGrid (const Index& aIdx) const; /** * Computes the candidate grid points lying in a cone given by two grid points. @@ -246,7 +247,7 @@ namespace DGtal * * @param aIdx */ - void updateGrid (const int& aIdx); + void updateGrid (const Index& aIdx); /** * Returns the vector from the base to a grid point. diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih index c3dc05324c..9857233b0d 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.ih @@ -121,11 +121,11 @@ template < typename TPredicate > inline typename DGtal::PlaneProbingLNeighborhood::ClosestGridPoint DGtal::PlaneProbingLNeighborhood:: -closestInGrid (const int& aIdx) const +closestInGrid (const Index& aIdx) const { std::vector candidates; - int k = aIdx; + Index k = aIdx; GridPoint y1(1,0,k); //y1 = vk + m1 GridPoint y2(0,1,k); //y2 = vk + m2 @@ -305,7 +305,7 @@ template < typename TPredicate > inline void DGtal::PlaneProbingLNeighborhood:: -updateGrid (const int& aIdx) +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}); diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.h index b21985d3f9..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; }; diff --git a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih index 7684acfdb6..57669ffce8 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingNeighborhood.ih @@ -176,7 +176,7 @@ getOperation (PointOnProbingRay const& aClosest) const { return { aClosest.sigma(), - { 1, -1, -aClosest.index() }, + { 1, -1, -aClosest.position() }, }; } 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 cb0957b275..c605fbe66f 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.ih @@ -215,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); @@ -323,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.ih b/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.ih index 3affa76230..27c8230eef 100644 --- a/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.ih +++ b/src/DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.ih @@ -101,14 +101,14 @@ DGtal::PlaneProbingRNeighborhood::closestPointOnRayLogWithPredicate 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))); @@ -116,8 +116,8 @@ DGtal::PlaneProbingRNeighborhood::closestPointOnRayLogWithPredicate 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))) { @@ -130,7 +130,7 @@ DGtal::PlaneProbingRNeighborhood::closestPointOnRayLogWithPredicate Xl = Xgamma; } - d = Xl.index() - Xk.index(); + d = Xl.position() - Xk.position(); } return closestPointOnRayLinearWithPredicate(Xk); From e1c26b4cacdf9d83f28e18754d6ceaf855b7190b Mon Sep 17 00:00:00 2001 From: Tristan Roussillon Date: Thu, 3 Oct 2024 16:31:54 +0200 Subject: [PATCH 9/9] bug fix --- src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h index bd4bc36e41..0b2b236885 100644 --- a/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h +++ b/src/DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h @@ -214,8 +214,8 @@ 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 @@ -397,8 +397,8 @@ namespace DGtal * @param aGridPoint the grid point to display. * @return the output stream after the writing. */ - template < typename Integer > - std::ostream& operator<< (std::ostream& aOs, GridPoint const& aGridPoint) { + 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