diff --git a/ChangeLog.md b/ChangeLog.md index d5db097696..4d01f5694a 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,5 +1,8 @@ # DGtal 1.5beta +## Bug fixes +- *Geometry* + - Bug fix in ArithmeticalDSSComputerOnSurfels (Tristan Roussillon, [#1742](https://github.com/DGtal-team/DGtal/pull/1742)) # DGtal 1.4.1 diff --git a/src/DGtal/geometry/doc/moduleMaximalSegmentSliceEstimation.dox b/src/DGtal/geometry/doc/moduleMaximalSegmentSliceEstimation.dox index 861121ab93..d4a4b9ecbe 100644 --- a/src/DGtal/geometry/doc/moduleMaximalSegmentSliceEstimation.dox +++ b/src/DGtal/geometry/doc/moduleMaximalSegmentSliceEstimation.dox @@ -69,7 +69,7 @@ using Surfel = KSpace::SCell; using Container = std::vector; using Integer = int; short adjacency = 4; -using SegmentComputerOnSurfels = ArithmeticalDSSComputerOnSurfels; +using SegmentComputerOnSurfels = ArithmeticalDSSComputerOnSurfels; // Instantiation SegmentComputerOnSurfels computer(kspace, // A 3D Khalimsky space diff --git a/src/DGtal/geometry/surfaces/ArithmeticalDSSComputerOnSurfels.h b/src/DGtal/geometry/surfaces/ArithmeticalDSSComputerOnSurfels.h index 8cf64f820a..01b364d621 100644 --- a/src/DGtal/geometry/surfaces/ArithmeticalDSSComputerOnSurfels.h +++ b/src/DGtal/geometry/surfaces/ArithmeticalDSSComputerOnSurfels.h @@ -18,10 +18,10 @@ /** * @file ArithmeticalDSSComputerOnSurfels.h - * @author Jocelyn Meyron (\c - * jocelyn.meyron@liris.cnrs.fr ) Laboratoire d'InfoRmatique en - * Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, - * France + * @author Jocelyn Meyron (\c jocelyn.meyron@liris.cnrs.fr ) + * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France + * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr ) + * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France * * @date 2021/01/22 * @@ -46,7 +46,9 @@ #include #include "DGtal/base/Exceptions.h" #include "DGtal/base/Common.h" +#include "DGtal/kernel/SpaceND.h" #include "DGtal/kernel/PointVector.h" +#include "DGtal/kernel/BasicPointFunctors.h" #include "DGtal/kernel/CInteger.h" #include "DGtal/base/ReverseIterator.h" #include "DGtal/geometry/curves/ArithmeticalDSS.h" @@ -60,8 +62,10 @@ namespace DGtal // class ArithmeticalDSSComputerOnSurfels /** * \brief Aim: This class is a wrapper around ArithmeticalDSS that is devoted - * to the dynamic recognition of digital straight segments (DSS) along any - * sequence of 3D surfels. + * to the dynamic recognition of digital straight segments (DSS) along a + * sequence of surfels lying on a slice of the digital surface (i.e., the + * orthogonal direction of all surfels belong to a same plane, most pairs + * of consecutive surfels share a common linel). * * @tparam TKSpace type of Khalimsky space * @tparam TIterator type of iterator on 3d surfels, @@ -77,8 +81,7 @@ namespace DGtal * @see ArithmeticalDSS NaiveDSS8 StandardDSS4 */ template + typename TInteger = typename TKSpace::Space::Integer> class ArithmeticalDSSComputerOnSurfels { BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > )); @@ -97,11 +100,6 @@ namespace DGtal */ typedef typename KSpace::SCell SCell; - /** - * Type of unsigned cell - */ - typedef typename KSpace::Cell Cell; - /** * Type of iterator, at least readable and forward */ @@ -120,6 +118,11 @@ namespace DGtal */ typedef PointVector<2, TInteger> Point; + /** + * Type of 3d to 2d projector + */ + typedef functors::Projector > Projector; + /** * Type of coordinate */ @@ -135,7 +138,7 @@ namespace DGtal /** * Type of objects that represents DSSs */ - typedef ArithmeticalDSS DSS; + typedef ArithmeticalDSS DSS; //we expect that the iterator type returned DGtal points, used in the DSS representation BOOST_STATIC_ASSERT(( concepts::ConceptUtils::SameType< Point, typename DSS::Point >::value )); @@ -149,8 +152,32 @@ namespace DGtal */ typedef Point Vector; - typedef ArithmeticalDSSComputerOnSurfels Self; - typedef ArithmeticalDSSComputerOnSurfels,TInteger,adjacency> Reverse; + /** + * Alias of this class + */ + typedef ArithmeticalDSSComputerOnSurfels Self; + + /** + * Helpers used to extract relevant points from a pair of points + */ + struct DirectPairExtractor { + + virtual Point first(const std::pair& aPair) const { return aPair.first; } + virtual Point second(const std::pair& aPair) const { return aPair.second; } + + }; + struct IndirectPairExtractor : public DirectPairExtractor { + + Point first(const std::pair& aPair) const { return aPair.second; } + Point second(const std::pair& aPair) const { return aPair.first; } + + }; + typedef std::shared_ptr PairExtractor; + + /** + * Reversed version of this class (using reverse iterators) + */ + typedef ArithmeticalDSSComputerOnSurfels,TInteger> Reverse; // ----------------------- Standard services ------------------------------ public: @@ -164,11 +191,14 @@ namespace DGtal /** * Constructor. * @param aKSpace a Khalimsky space - * @param aDim1 the first direction to project - * @param aDim2 the second direction to project + * @param aDim1 a first direction that describes the projection plane + * @param aDim2 a second direction that describes the projection plane + * @param aFlagToReverse a boolean telling whether one has to reverse + * the orientation of the points associated to a surfel or not + * ('false' by default) */ - ArithmeticalDSSComputerOnSurfels(const KSpace& aKSpace, Dimension aDim1, Dimension aDim2); - + ArithmeticalDSSComputerOnSurfels(const KSpace& aKSpace, Dimension aDim1, Dimension aDim2, bool aFlagToReverse = false); + /** * Initialisation. * @param it an iterator on 3D surfels @@ -227,6 +257,9 @@ namespace DGtal * Tests whether the current DSS can be extended at the front. * * @return 'true' if yes, 'false' otherwise. + * + * @warning the caller must be sure that the iterator returned + * by 'end()' can be safely dereferenced. */ bool isExtendableFront(); @@ -240,7 +273,11 @@ namespace DGtal /** * Tests whether the current DSS can be extended at the front. * Computes the parameters of the extended DSS if yes. + * * @return 'true' if yes, 'false' otherwise. + * + * @warning the caller must be sure that the iterator returned + * by 'end()' can be safely dereferenced. */ bool extendFront(); @@ -265,6 +302,32 @@ namespace DGtal */ bool retractBack(); + /** + * Returns the ends of a unit segment corresponding + * to the projection of a given signed surfel. + * + * @param aSurfel any signed surfel. + * @return a pair of 2D points. + */ + std::pair getProjectedPointsFromSurfel(SCell const& aSurfel) const; + + /** + * Front end of the projection of a given surfel. + * + * @param aSurfel any signed surfel. + * @return the second 2D point. + * @see getProjectedPointsFromSurfel + */ + Point getNextProjectedPoint(SCell const& aSurfel) const; + + /** + * Back end of the projection of a given surfel. + * + * @param aSurfel any signed surfel. + * @return the second 2D point. + * @see getProjectedPointsFromSurfel + */ + Point getPreviousProjectedPoint(SCell const& aSurfel) const; // ------------------------- Accessors ------------------------------ /** @@ -329,92 +392,77 @@ namespace DGtal */ bool isValid() const; - /** - * @param aSCell a surfel. - * @return the pair of 2D points projected on the plane defined by \ref myProjection1, \ref myProjection2. - */ - std::pair projectSurfel(SCell const& aSCell) const; // ------------------------- Hidden services ------------------------------ private: - /** - * @param aSurfel1 the first unsigned surfel. - * @param aSurfel2 the second unsigned surfel. - * @param aLinel the common unsigned linel if it exists. - * @return 'true' if we found a common linel, 'false' otherwise. - */ - bool commonLinel (Cell const& aSurfel1, Cell const& aSurfel2, Cell& aLinel); /** - * @param aPoint a digital 3D point. - * @return the 2D orthogonal projection on the plane defined by \ref myProjection1, \ref myProjection2 + * Returns the ends of a unit segment corresponding + * to the projection of a given signed linel. + * + * @param aLinel any signed linel. + * @return a pair of 2D points. */ - Point projectInPlane (Point3 const& aPoint) const; - + std::pair getProjectedPointsFromLinel(SCell const& aLinel) const; + /** - * @param it an iterator on a surfel. - * @param aPoint the new 2D point of 'it' that is not common to the previous surfel (if the update is possible). - * @param aUpdatePrevious a boolean that indicates whether to update myPreviousSurfel or not. - * @param aIsFront a boolean indicating we want to update in the front or in the back direction. - * @return 'true' if the update is possible, 'false' otherwise. + * Returns the unique dimension in {0,1,2} \ {aDim1, aDim2}. + * + * @param aDim1 a dimension + * @param aDim2 a dimension */ - bool getOtherPointFromSurfel (ConstIterator const& it, Point& aPoint, bool aIsFront, bool aUpdatePrevious); - + Dimension dimNotIn(Dimension const& aDim1, Dimension const& aDim2) const; + // ------------------------- Protected Datas ------------------------------ protected: /** - * Khalimsky space + * (Pointer to) Khalimsky space */ const KSpace* myKSpace; /** - * The first projection dimension. + * A first direction that describes the projection plane */ - Dimension myDim1; - + Dimension mySliceAxis1; + /** - * The second projection dimension. + * A second direction that describes the projection plane */ - Dimension myDim2; - + Dimension mySliceAxis2; + /** - * The first projection direction. + * A direction along which the points are projected + * (and orthogonal to the projection plane) */ - Point3 myProjection1; + Dimension myProjectionAxis; /** - * The second projection direction. + * Functor that projects a 3D point to a 2D point along myProjectionAxis */ - Point3 myProjection2; + Projector my2DProjector; /** - * The direction that is orthogonal to the two projection directions. + * Smart pointer on an object used to extract relevant points + * from a pair of points */ - Point3 myProjectionNormal; - - /** - * We store the previous surfel in the 'front' direction to compute the common linel. - * Used in \ref getOtherPointFromSurfel. - */ - ConstIterator myPreviousSurfelFront; - - /** - * We store the previous surfel in the 'back' direction to compute the common linel. - * Used in \ref getOtherPointFromSurfel. - */ - ConstIterator myPreviousSurfelBack; - + PairExtractor myExtractor; + /** * DSS representation */ DSS myDSS; + /** * begin iterator */ ConstIterator myBegin; + /** * end iterator + * + * @warning the user must be sure that it can be safely dereferenced + * before calling 'isExtendableFront' and 'extendFront'. */ ConstIterator myEnd; @@ -438,9 +486,9 @@ namespace DGtal * @param object the object of class 'ArithmeticalDSSComputerOnSurfels' to write. * @return the output stream after the writing. */ -template +template std::ostream& -operator<< ( std::ostream & out, const ArithmeticalDSSComputerOnSurfels & object ) +operator<< ( std::ostream & out, const ArithmeticalDSSComputerOnSurfels & object ) { object.selfDisplay( out); return out; diff --git a/src/DGtal/geometry/surfaces/ArithmeticalDSSComputerOnSurfels.ih b/src/DGtal/geometry/surfaces/ArithmeticalDSSComputerOnSurfels.ih index d66c67b9ef..2cc96cb9c9 100644 --- a/src/DGtal/geometry/surfaces/ArithmeticalDSSComputerOnSurfels.ih +++ b/src/DGtal/geometry/surfaces/ArithmeticalDSSComputerOnSurfels.ih @@ -18,6 +18,8 @@ * @file ArithmeticalDSSComputerOnSurfels.ih * @author Jocelyn Meyron (\c jocelyn.meyron@liris.cnrs.fr ) * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France + * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr ) + * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France * * @date 2021/01/22 * @@ -47,103 +49,62 @@ // Implementation of inline methods // //----------------------------------------------------------------------------- -template +template inline -DGtal::ArithmeticalDSSComputerOnSurfels:: +DGtal::ArithmeticalDSSComputerOnSurfels:: ArithmeticalDSSComputerOnSurfels() - : myKSpace(nullptr), myDSS( Point(0,0) ), myBegin(), myEnd() -{ -} - -//----------------------------------------------------------------------------- -template -inline -DGtal::ArithmeticalDSSComputerOnSurfels:: -ArithmeticalDSSComputerOnSurfels(const KSpace& aKSpace, Dimension aDim1, Dimension aDim2) - : myKSpace(&aKSpace), myDim1(aDim1), myDim2(aDim2), myDSS( Point(0,0) ), myBegin(), myEnd() + : myKSpace(nullptr), mySliceAxis1(), mySliceAxis2(), myProjectionAxis(), + my2DProjector(), myExtractor(nullptr), myDSS( Point(0,0) ), myBegin(), myEnd() { - // Initialize projection vectors - myProjection1 = Point3::zero; - myProjection2 = Point3::zero; - myProjection1[myDim1] = 1; - myProjection2[myDim2] = 1; - - myProjectionNormal = myProjection1.crossProduct(myProjection2); } - //----------------------------------------------------------------------------- -template +template inline -void DGtal::ArithmeticalDSSComputerOnSurfels:: -init(const ConstIterator& it) +DGtal::ArithmeticalDSSComputerOnSurfels:: +ArithmeticalDSSComputerOnSurfels(const KSpace& aKSpace, Dimension aDim1, Dimension aDim2, bool aFlagToReverse) + : myKSpace(&aKSpace), mySliceAxis1(aDim1), mySliceAxis2(aDim2), myProjectionAxis(), + my2DProjector(), myExtractor(), myDSS( Point(0,0) ), myBegin(), myEnd() { - ASSERT(myKSpace != nullptr); + ASSERT(aDim1 != aDim2); + + myProjectionAxis = dimNotIn(aDim1, aDim2); - myBegin = it; - myEnd = it; - ++myEnd; - - SCell s = *it; - auto initialPoints = projectSurfel(s); - myPreviousSurfelFront = it; - myPreviousSurfelBack = it; - - // We ensure that the first two points are inserted in the correct order - Point firstPoint = initialPoints.first, secondPoint = initialPoints.second; - if (myKSpace->sIsValid(*myEnd)) { - auto nextPoints = projectSurfel(*myEnd); - if (nextPoints.first == initialPoints.first) { - secondPoint = initialPoints.first; - firstPoint = initialPoints.second; - } else if (nextPoints.first == initialPoints.second) { - secondPoint = initialPoints.second; - firstPoint = initialPoints.first; - } else if (nextPoints.second == initialPoints.first) { - secondPoint = initialPoints.first; - firstPoint = initialPoints.second; - } else if (nextPoints.second == initialPoints.second) { - secondPoint = initialPoints.second; - firstPoint = initialPoints.first; - } else { - ASSERT(false); - } - } + std::vector v = {aDim1, aDim2}; + my2DProjector.init(v.begin(),v.end()); - myDSS = DSS(firstPoint); - ASSERT(myDSS.isExtendableFront(secondPoint)); - myDSS.extendFront(secondPoint); + if (aFlagToReverse) { + myExtractor = PairExtractor(new IndirectPairExtractor()); + } else { + myExtractor = PairExtractor(new DirectPairExtractor()); + } } //----------------------------------------------------------------------------- -template +template inline -DGtal::ArithmeticalDSSComputerOnSurfels:: -ArithmeticalDSSComputerOnSurfels ( const ArithmeticalDSSComputerOnSurfels & other ) - : myKSpace(other.myKSpace), myDim1(other.myDim1), myDim2(other.myDim2), - myProjection1(other.myProjection1), myProjection2(other.myProjection2), myProjectionNormal(other.myProjectionNormal), - myPreviousSurfelFront(other.myPreviousSurfelFront), myPreviousSurfelBack(other.myPreviousSurfelBack), - myDSS(other.myDSS), myBegin(other.myBegin), myEnd(other.myEnd) +DGtal::ArithmeticalDSSComputerOnSurfels:: +ArithmeticalDSSComputerOnSurfels ( const ArithmeticalDSSComputerOnSurfels & other ) + : myKSpace(other.myKSpace), mySliceAxis1(other.mySliceAxis1), mySliceAxis2(other.mySliceAxis2), myProjectionAxis(other.myProjectionAxis), + my2DProjector(other.my2DProjector), myExtractor(other.myExtractor), myDSS(other.myDSS), myBegin(other.myBegin), myEnd(other.myEnd) { } //----------------------------------------------------------------------------- -template +template inline -typename DGtal::ArithmeticalDSSComputerOnSurfels& -DGtal::ArithmeticalDSSComputerOnSurfels:: -operator=( const ArithmeticalDSSComputerOnSurfels & other ) +typename DGtal::ArithmeticalDSSComputerOnSurfels& +DGtal::ArithmeticalDSSComputerOnSurfels:: +operator=( const ArithmeticalDSSComputerOnSurfels & other ) { if ( this != &other ) { myKSpace = other.myKSpace; - myDim1 = other.myDim1; - myDim2 = other.myDim2; - myProjection1 = other.myProjection1; - myProjection2 = other.myProjection2; - myProjectionNormal = other.myProjectionNormal; - myPreviousSurfelFront = other.myPreviousSurfelFront; - myPreviousSurfelBack = other.myPreviousSurfelBack; + mySliceAxis1 = other.mySliceAxis1; + mySliceAxis2 = other.mySliceAxis2; + myProjectionAxis = other.myProjectionAxis; + my2DProjector = other.my2DProjector; + myExtractor = other.myExtractor; myDSS = other.myDSS; myBegin = other.myBegin; myEnd = other.myEnd; @@ -152,31 +113,31 @@ operator=( const ArithmeticalDSSComputerOnSurfels +template inline -typename DGtal::ArithmeticalDSSComputerOnSurfels::Reverse -DGtal::ArithmeticalDSSComputerOnSurfels +typename DGtal::ArithmeticalDSSComputerOnSurfels::Reverse +DGtal::ArithmeticalDSSComputerOnSurfels ::getReverse() const { - return Reverse(*myKSpace, myDim1, myDim2); + return Reverse(*myKSpace, mySliceAxis1, mySliceAxis2, true); } //----------------------------------------------------------------------------- -template +template inline -typename DGtal::ArithmeticalDSSComputerOnSurfels::Self -DGtal::ArithmeticalDSSComputerOnSurfels +typename DGtal::ArithmeticalDSSComputerOnSurfels::Self +DGtal::ArithmeticalDSSComputerOnSurfels ::getSelf() const { - return Self(*myKSpace, myDim1, myDim2); + return Self(*myKSpace, mySliceAxis1, mySliceAxis2); } //----------------------------------------------------------------------------- -template +template inline bool -DGtal::ArithmeticalDSSComputerOnSurfels:: -operator==( const ArithmeticalDSSComputerOnSurfels& other ) const +DGtal::ArithmeticalDSSComputerOnSurfels:: +operator==( const ArithmeticalDSSComputerOnSurfels& other ) const { return ( (myBegin == other.myBegin) && (myEnd == other.myEnd) @@ -184,11 +145,11 @@ operator==( const ArithmeticalDSSComputerOnSurfels +template inline bool -DGtal::ArithmeticalDSSComputerOnSurfels:: -operator!=( const ArithmeticalDSSComputerOnSurfels & other ) const +DGtal::ArithmeticalDSSComputerOnSurfels:: +operator!=( const ArithmeticalDSSComputerOnSurfels & other ) const { return (!(*this == other)); } @@ -197,391 +158,337 @@ operator!=( const ArithmeticalDSSComputerOnSurfels +template inline -bool -DGtal::ArithmeticalDSSComputerOnSurfels::isExtendableFront() +bool +DGtal::ArithmeticalDSSComputerOnSurfels::isExtendableFront() { - Point p; - if (! getOtherPointFromSurfel(myEnd, p, true, false)) - return false; - - return myDSS.isExtendableFront( p ); + //the caller must be sure that 'myEnd' can be safely dereferenced + return myDSS.isExtendableFront( getNextProjectedPoint(*myEnd) ); } //-------------------------------------------------------------------- -template +template inline -bool -DGtal::ArithmeticalDSSComputerOnSurfels::isExtendableBack() +bool +DGtal::ArithmeticalDSSComputerOnSurfels::isExtendableBack() { - ConstIterator it = myBegin; - --it; - Point p; - if (! getOtherPointFromSurfel(it, p, false, false)) - return false; - - return myDSS.isExtendableBack( p ); + ConstIterator it = myBegin; + --it; + //the caller must be sure that 'it' can be safely dereferenced + return myDSS.isExtendableBack( getPreviousProjectedPoint(*it) ); } //----------------------------------------------------------------------------- -template +template inline -bool -DGtal::ArithmeticalDSSComputerOnSurfels::extendFront() +bool +DGtal::ArithmeticalDSSComputerOnSurfels::extendFront() { - Point p; - if (! getOtherPointFromSurfel(myEnd, p, true, true)) - return false; - - if (myDSS.extendFront(p)) + //the caller must be sure that 'myEnd' can be safely dereferenced + if (myDSS.extendFront(getNextProjectedPoint(*myEnd))) { ++myEnd; - return true; + return true; } - else - return false; + else + return false; } //-------------------------------------------------------------------- -template +template inline -bool -DGtal::ArithmeticalDSSComputerOnSurfels::extendBack() +bool +DGtal::ArithmeticalDSSComputerOnSurfels::extendBack() { - ConstIterator it = myBegin; - --it; - Point p; - if (! getOtherPointFromSurfel(it, p, false, true)) - return false; - - if (myDSS.extendBack(p)) + ConstIterator it = myBegin; + --it; + //the caller must be sure that 'it' can be safely dereferenced + if (myDSS.extendBack(getPreviousProjectedPoint(*it))) { myBegin = it; - return true; + return true; } - else - return false; + else + return false; } //-------------------------------------------------------------------- -template +template inline -bool -DGtal::ArithmeticalDSSComputerOnSurfels::retractFront() +bool +DGtal::ArithmeticalDSSComputerOnSurfels::retractFront() { if (myDSS.retractFront()) { - --myEnd; - return true; + --myEnd; + return true; } - else - return false; + else + return false; } //-------------------------------------------------------------------- -template +template inline -bool -DGtal::ArithmeticalDSSComputerOnSurfels::retractBack() +bool +DGtal::ArithmeticalDSSComputerOnSurfels::retractBack() { if (myDSS.retractBack()) { ++myBegin; - return true; + return true; } - else - return false; + else + return false; } /////////////////////////////////////////////////////////////////////////////// // Accessors // /////////////////////////////////////////////////////////////////////////////// //------------------------------------------------------------------------- -template +template inline -const typename DGtal::ArithmeticalDSSComputerOnSurfels::Primitive& -DGtal::ArithmeticalDSSComputerOnSurfels::primitive() const +const typename DGtal::ArithmeticalDSSComputerOnSurfels::Primitive& +DGtal::ArithmeticalDSSComputerOnSurfels::primitive() const { return myDSS; } //------------------------------------------------------------------------- -template +template inline TInteger -DGtal::ArithmeticalDSSComputerOnSurfels::a() const +DGtal::ArithmeticalDSSComputerOnSurfels::a() const { return myDSS.a(); } //------------------------------------------------------------------------- -template +template inline TInteger -DGtal::ArithmeticalDSSComputerOnSurfels::b() const +DGtal::ArithmeticalDSSComputerOnSurfels::b() const { return myDSS.b(); } //------------------------------------------------------------------------- -template +template inline TInteger -DGtal::ArithmeticalDSSComputerOnSurfels::mu() const +DGtal::ArithmeticalDSSComputerOnSurfels::mu() const { return myDSS.mu(); } //------------------------------------------------------------------------- -template +template inline TInteger -DGtal::ArithmeticalDSSComputerOnSurfels::omega() const +DGtal::ArithmeticalDSSComputerOnSurfels::omega() const { return myDSS.omega(); } //------------------------------------------------------------------------- -template +template inline -typename DGtal::ArithmeticalDSSComputerOnSurfels::Point -DGtal::ArithmeticalDSSComputerOnSurfels::Uf() const +typename DGtal::ArithmeticalDSSComputerOnSurfels::Point +DGtal::ArithmeticalDSSComputerOnSurfels::Uf() const { return myDSS.Uf(); } //------------------------------------------------------------------------- -template +template inline -typename DGtal::ArithmeticalDSSComputerOnSurfels::Point -DGtal::ArithmeticalDSSComputerOnSurfels::Ul() const +typename DGtal::ArithmeticalDSSComputerOnSurfels::Point +DGtal::ArithmeticalDSSComputerOnSurfels::Ul() const { return myDSS.Ul(); } //------------------------------------------------------------------------- -template +template inline -typename DGtal::ArithmeticalDSSComputerOnSurfels::Point -DGtal::ArithmeticalDSSComputerOnSurfels::Lf() const +typename DGtal::ArithmeticalDSSComputerOnSurfels::Point +DGtal::ArithmeticalDSSComputerOnSurfels::Lf() const { return myDSS.Lf(); } //------------------------------------------------------------------------- -template +template inline -typename DGtal::ArithmeticalDSSComputerOnSurfels::Point -DGtal::ArithmeticalDSSComputerOnSurfels::Ll() const +typename DGtal::ArithmeticalDSSComputerOnSurfels::Point +DGtal::ArithmeticalDSSComputerOnSurfels::Ll() const { return myDSS.Ll(); } //------------------------------------------------------------------------- -template +template inline -typename DGtal::ArithmeticalDSSComputerOnSurfels::Point -DGtal::ArithmeticalDSSComputerOnSurfels::back() const +typename DGtal::ArithmeticalDSSComputerOnSurfels::Point +DGtal::ArithmeticalDSSComputerOnSurfels::back() const { return myDSS.back(); } //------------------------------------------------------------------------- -template +template inline -typename DGtal::ArithmeticalDSSComputerOnSurfels::Point -DGtal::ArithmeticalDSSComputerOnSurfels::front() const +typename DGtal::ArithmeticalDSSComputerOnSurfels::Point +DGtal::ArithmeticalDSSComputerOnSurfels::front() const { return myDSS.front(); } //------------------------------------------------------------------------- -template +template inline TIterator -DGtal::ArithmeticalDSSComputerOnSurfels::begin() const +DGtal::ArithmeticalDSSComputerOnSurfels::begin() const { return myBegin; } //------------------------------------------------------------------------- -template +template inline TIterator -DGtal::ArithmeticalDSSComputerOnSurfels::end() const +DGtal::ArithmeticalDSSComputerOnSurfels::end() const { return myEnd; } //----------------------------------------------------------------- -template +template inline bool -DGtal::ArithmeticalDSSComputerOnSurfels::isValid() const +DGtal::ArithmeticalDSSComputerOnSurfels::isValid() const { return ( (myDSS.isValid())&&(isNotEmpty(myBegin,myEnd)) ); } //----------------------------------------------------------------- -template +template inline void -DGtal::ArithmeticalDSSComputerOnSurfels::selfDisplay ( std::ostream & out) const +DGtal::ArithmeticalDSSComputerOnSurfels::selfDisplay ( std::ostream & out) const { out << "[ArithmeticalDSSComputerOnSurfels] " << myDSS; } -//----------------------------------------------------------------- -template -inline -bool -DGtal::ArithmeticalDSSComputerOnSurfels::commonLinel (Cell const& aSurfel1, - Cell const& aSurfel2, - Cell& aLinel) -{ - ASSERT(myKSpace != nullptr); - - typename KSpace::DirIterator q1_1 = myKSpace->uDirs(aSurfel1), q1_2 = q1_1; - ++q1_2; - typename KSpace::DirIterator q2_1 = myKSpace->uDirs(aSurfel2), q2_2 = q2_1; - ++q2_2; - - std::set linels1 = { - myKSpace->uIncident(aSurfel1, *q1_1, true), - myKSpace->uIncident(aSurfel1, *q1_1, false), - myKSpace->uIncident(aSurfel1, *q1_2, true), - myKSpace->uIncident(aSurfel1, *q1_2, false), - }; - - std::set linels2 = { - myKSpace->uIncident(aSurfel2, *q2_1, true), - myKSpace->uIncident(aSurfel2, *q2_1, false), - myKSpace->uIncident(aSurfel2, *q2_2, true), - myKSpace->uIncident(aSurfel2, *q2_2, false), - }; - - std::vector inter; - std::set_intersection(linels1.begin(), linels1.end(), - linels2.begin(), linels2.end(), - std::back_inserter(inter)); - - if (inter.size() == 1) - { - // The two surfels intersect on one linel - aLinel = inter[0]; - return true; - } - - // The surfels don't intersect or are the same - return false; -} +/////////////////////////////////////////////////////////////////////////////// +// Initialization // +/////////////////////////////////////////////////////////////////////////////// -//----------------------------------------------------------------- -template +//----------------------------------------------------------------------------- +template inline -bool -DGtal::ArithmeticalDSSComputerOnSurfels::getOtherPointFromSurfel(ConstIterator const& it, - Point& aPoint, - bool aIsFront, - bool aUpdatePrevious) +void DGtal::ArithmeticalDSSComputerOnSurfels:: +init(const ConstIterator& it) { - ASSERT(myKSpace != nullptr); - - SCell surfel = *it; - Point p1, p2; - std::tie(p1, p2) = projectSurfel(surfel); + ASSERT(myKSpace != nullptr); - ConstIterator& previousSurfel = aIsFront ? myPreviousSurfelFront : myPreviousSurfelBack; + myBegin = it; + myEnd = it; + ++myEnd; - // Find the common unsigned linel between surfel and previousSurfel - Cell linel; - if (! commonLinel(myKSpace->unsigns(surfel), myKSpace->unsigns(*previousSurfel), linel)) - { - return false; - } + SCell surfel = *it; + auto pair = getProjectedPointsFromSurfel(surfel); + auto p = myExtractor->first(pair); + myDSS = DSS( p ); + auto q = myExtractor->second(pair); + ASSERT(myDSS.isExtendableFront( q )); + myDSS.extendFront( q ); +} - // For the next point, choose the point that is not part of the common linel - typename KSpace::DirIterator q_linel = myKSpace->uDirs(linel); - Point linel1 = projectInPlane(myKSpace->uCoords(myKSpace->uIncident(linel, *q_linel, true))), - linel2 = projectInPlane(myKSpace->uCoords(myKSpace->uIncident(linel, *q_linel, false))); +/////////////////////////////////////////////////////////////////////////////// +// Projection // +/////////////////////////////////////////////////////////////////////////////// - if (p1 == linel1) - { - aPoint = p2; - } - else if (p1 == linel2) - { - aPoint = p2; - } - else if (p2 == linel1) - { - aPoint = p1; - } - else if (p2 == linel2) - { - aPoint = p1; - } - else - { - ASSERT(false); - } +//----------------------------------------------------------------- +template +inline +std::pair::Point, + typename DGtal::ArithmeticalDSSComputerOnSurfels::Point> +DGtal::ArithmeticalDSSComputerOnSurfels::getProjectedPointsFromSurfel(SCell const& aSurfel) const +{ - if (aUpdatePrevious) - { - previousSurfel = it; - } + SCell linel1; + //this convention has been chosen so that + //linels always stand at the same side + if (myKSpace->sSign(aSurfel) == myKSpace->POS) { + linel1 = myKSpace->sDirectIncident(aSurfel, myProjectionAxis); + } else { + linel1 = myKSpace->sIndirectIncident(aSurfel, myProjectionAxis); + } - return true; + return getProjectedPointsFromLinel(linel1); } //----------------------------------------------------------------- -template +template inline -typename DGtal::ArithmeticalDSSComputerOnSurfels::Point -DGtal::ArithmeticalDSSComputerOnSurfels::projectInPlane(Point3 const& aPoint) const +typename DGtal::ArithmeticalDSSComputerOnSurfels::Point +DGtal::ArithmeticalDSSComputerOnSurfels::getNextProjectedPoint(SCell const& aSurfel) const { - static const Point3 aOrigin = Point3::zero; - - // Orthogonal projection on the plane with a given unit normal - Point3 h = (aPoint - aOrigin) - myProjectionNormal * (aPoint - aOrigin).dot(myProjectionNormal); - - // We simply project the point on the plane defined by - // the two directions 'myProjection1' and 'myProjection2' passing through the origin point 'aOrigin' - return Point(h.dot(myProjection1), h.dot(myProjection2)); + return myExtractor->second(getProjectedPointsFromSurfel(aSurfel)); } //----------------------------------------------------------------- -template +template inline -std::pair::Point, - typename DGtal::ArithmeticalDSSComputerOnSurfels::Point> -DGtal::ArithmeticalDSSComputerOnSurfels::projectSurfel(SCell const& aSCell) const +typename DGtal::ArithmeticalDSSComputerOnSurfels::Point +DGtal::ArithmeticalDSSComputerOnSurfels::getPreviousProjectedPoint(SCell const& aSurfel) const { - ASSERT(myKSpace != nullptr); + return myExtractor->first(getProjectedPointsFromSurfel(aSurfel)); +} - typename KSpace::DirIterator q1 = myKSpace->sDirs(aSCell), q2 = q1; - ++q2; +/////////////////////////////////////////////////////////////////////////////// +// Private methods // +/////////////////////////////////////////////////////////////////////////////// - // We pick 2 linels of the surfel - SCell linel1 = myKSpace->sIncident(aSCell, *q1, true), - linel2 = myKSpace->sIncident(aSCell, *q1, false); - // 4 points of the surfel - Point3 p1_1 = myKSpace->sCoords(myKSpace->sIncident(linel1, *q2, false)), - p1_2 = myKSpace->sCoords(myKSpace->sIncident(linel1, *q2, true)), - p2_1 = myKSpace->sCoords(myKSpace->sIncident(linel2, *q2, false)), - p2_2 = myKSpace->sCoords(myKSpace->sIncident(linel2, *q2, true)); +//----------------------------------------------------------------- +template +inline +std::pair::Point, + typename DGtal::ArithmeticalDSSComputerOnSurfels::Point> +DGtal::ArithmeticalDSSComputerOnSurfels::getProjectedPointsFromLinel(SCell const& aLinel) const +{ + auto dim = *myKSpace->sDirs(aLinel); + + auto pointel1 = myKSpace->sIndirectIncident(aLinel, dim); + auto p1 = my2DProjector(myKSpace->sCoords(pointel1)); - std::set points; - points.insert(projectInPlane(p1_1)); - points.insert(projectInPlane(p1_2)); - points.insert(projectInPlane(p2_1)); - points.insert(projectInPlane(p2_2)); + auto pointel2 = myKSpace->sDirectIncident(aLinel, dim); + auto p2 = my2DProjector(myKSpace->sCoords(pointel2)); - ASSERT(points.size() == 2); + return {p1, p2}; +} - Point p1 = *points.begin(), p2 = *(++points.begin()); +//----------------------------------------------------------------- +template +inline +DGtal::Dimension +DGtal::ArithmeticalDSSComputerOnSurfels::dimNotIn(Dimension const& aDim1, Dimension const& aDim2) const +{ + ASSERT( KSpace::dimension == 3 ); + ASSERT( aDim1 != aDim2 ); + + bool marks[3] = {false, false, false}; + marks[aDim1] = true; + marks[aDim2] = true; + + int i = 0; + while (marks[i] == true) { + i++; + } - return { p1, p2 }; + ASSERT( (marks[i] == false) && (i < 3) ); + return i; } diff --git a/src/DGtal/geometry/surfaces/estimation/MaximalSegmentSliceEstimation.h b/src/DGtal/geometry/surfaces/estimation/MaximalSegmentSliceEstimation.h index 8bae427c16..1aec8c7c6c 100644 --- a/src/DGtal/geometry/surfaces/estimation/MaximalSegmentSliceEstimation.h +++ b/src/DGtal/geometry/surfaces/estimation/MaximalSegmentSliceEstimation.h @@ -94,7 +94,7 @@ namespace DGtal using RealPoint2 = PointVector<2, double>; using Container = std::deque; using Iterator = typename Container::const_iterator; - using DSSComputer = ArithmeticalDSSComputerOnSurfels, Integer, 4>; + using DSSComputer = ArithmeticalDSSComputerOnSurfels, Integer>; // ----------------------- Standard services ------------------------------ public: diff --git a/tests/geometry/surfaces/testArithmeticalDSSComputerOnSurfels.cpp b/tests/geometry/surfaces/testArithmeticalDSSComputerOnSurfels.cpp index a679b9c807..6c038139e9 100644 --- a/tests/geometry/surfaces/testArithmeticalDSSComputerOnSurfels.cpp +++ b/tests/geometry/surfaces/testArithmeticalDSSComputerOnSurfels.cpp @@ -49,7 +49,7 @@ using KSpace = Z3i::KSpace; using SH3 = Shortcuts; using Surfel = KSpace::SCell; -using SegmentComputerOnSurfels = ArithmeticalDSSComputerOnSurfels::const_iterator, int, 4>; +using SegmentComputerOnSurfels = ArithmeticalDSSComputerOnSurfels::const_iterator, int>; using SegmentationSurfels = SaturatedSegmentation; using SegmentComputer = ArithmeticalDSSComputer::const_iterator, int, 4>; @@ -79,7 +79,8 @@ std::pair getSlice (std::string const& shape = "ellipsoid", doubl Surfel surfel = Surfaces::findABel(kspace, *binary_image, 10000); KSpace::DirIterator q1 = kspace.sDirs(surfel); - Dimension dim1 = *q1, dim2 = kspace.sOrthDir(surfel); + KSpace::DirIterator q2 = kspace.sOrthDirs(surfel); + Dimension dim1 = *q1, dim2 = *q2; auto tracker = surface->container().newTracker(surfel); SurfaceSlice surfaceSlice(tracker, dim1); delete tracker; @@ -95,42 +96,15 @@ std::vector extractPoints (SegmentComputerOnSurfels const& sc, Slice { std::vector points; - auto initialPoints = sc.projectSurfel(slice.start); + auto initialPoints = sc.getProjectedPointsFromSurfel(slice.start); points.push_back(initialPoints.first); points.push_back(initialPoints.second); - int currentIdx = 0; - bool firstIt = true; for (auto sit = slice.contour.begin() + 1; sit != slice.contour.end(); ++sit) { Surfel s = *sit; - auto projectedPoints = sc.projectSurfel(s); - - if (firstIt) { - if (projectedPoints.first == points[currentIdx]) { - points.push_back(projectedPoints.second); - } else if (projectedPoints.first == points[currentIdx + 1]) { - points.push_back(projectedPoints.second); - } else if (projectedPoints.second == points[currentIdx]) { - points.push_back(projectedPoints.first); - } else if (projectedPoints.second == points[currentIdx + 1]) { - points.push_back(projectedPoints.first); - } else { - assert(false); - } - - firstIt = false; - } else { - if (projectedPoints.first == points[currentIdx]) { - points.push_back(projectedPoints.second); - } else if (projectedPoints.second == points[currentIdx]) { - points.push_back(projectedPoints.first); - } else { - assert(false); - } - } - - currentIdx = (int)points.size() - 1; + auto pt = sc.getNextProjectedPoint(s); + points.push_back(pt); } return points; @@ -160,10 +134,10 @@ TEST_CASE("Testing ArithmeticalDSSComputerOnSurfels") auto segIt = segmentation.begin(); auto segSurfelIt = segmentationSurfels.begin(); while (segIt != segmentation.end() && segSurfelIt != segmentationSurfels.end()) { - allEqual = allEqual && (segIt->primitive() == segSurfelIt->primitive()); - - ++segIt; - ++segSurfelIt; + + allEqual = allEqual && (segIt->primitive() == segSurfelIt->primitive()); + ++segIt; + ++segSurfelIt; } REQUIRE(allEqual);