2 * This program is free software: you can redistribute it and/or modify
3 * it under the terms of the GNU Lesser General Public License as
4 * published by the Free Software Foundation, either version 3 of the
5 * License, or (at your option) any later version.
7 * This program is distributed in the hope that it will be useful,
8 * but WITHOUT ANY WARRANTY; without even the implied warranty of
9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10 * GNU General Public License for more details.
12 * You should have received a copy of the GNU General Public License
13 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 * @file ArithmeticalDSSComputerOnSurfels.ih
19 * @author Jocelyn Meyron (\c jocelyn.meyron@liris.cnrs.fr )
20 * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
24 * Implementation of inline methods defined in ArithmeticalDSSComputerOnSurfels.h
26 * This file is part of the DGtal library.
29 ///////////////////////////////////////////////////////////////////////////////
30 // IMPLEMENTATION of inline methods.
31 ///////////////////////////////////////////////////////////////////////////////
33 //////////////////////////////////////////////////////////////////////////////
36 #include <boost/version.hpp>
37 #if BOOST_VERSION < 105800
38 #include <boost/math/common_factor_rt.hpp>
40 #include <boost/integer/common_factor_rt.hpp>
43 //////////////////////////////////////////////////////////////////////////////
46 ///////////////////////////////////////////////////////////////////////////////
47 // Implementation of inline methods //
49 //-----------------------------------------------------------------------------
50 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
52 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
53 ArithmeticalDSSComputerOnSurfels()
54 : myKSpace(nullptr), myDSS( Point(0,0) ), myBegin(), myEnd()
58 //-----------------------------------------------------------------------------
59 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
61 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
62 ArithmeticalDSSComputerOnSurfels(const KSpace& aKSpace, Dimension aDim1, Dimension aDim2)
63 : myKSpace(&aKSpace), myDim1(aDim1), myDim2(aDim2), myDSS( Point(0,0) ), myBegin(), myEnd()
65 // Initialize projection vectors
66 myProjection1 = Point3::zero;
67 myProjection2 = Point3::zero;
68 myProjection1[myDim1] = 1;
69 myProjection2[myDim2] = 1;
71 myProjectionNormal = myProjection1.crossProduct(myProjection2);
75 //-----------------------------------------------------------------------------
76 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
78 void DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
79 init(const ConstIterator& it)
81 ASSERT(myKSpace != nullptr);
88 auto initialPoints = projectSurfel(s);
89 myPreviousSurfelFront = it;
90 myPreviousSurfelBack = it;
92 // We ensure that the first two points are inserted in the correct order
93 Point firstPoint = initialPoints.first, secondPoint = initialPoints.second;
94 if (myKSpace->sIsValid(*myEnd)) {
95 auto nextPoints = projectSurfel(*myEnd);
96 if (nextPoints.first == initialPoints.first) {
97 secondPoint = initialPoints.first;
98 firstPoint = initialPoints.second;
99 } else if (nextPoints.first == initialPoints.second) {
100 secondPoint = initialPoints.second;
101 firstPoint = initialPoints.first;
102 } else if (nextPoints.second == initialPoints.first) {
103 secondPoint = initialPoints.first;
104 firstPoint = initialPoints.second;
105 } else if (nextPoints.second == initialPoints.second) {
106 secondPoint = initialPoints.second;
107 firstPoint = initialPoints.first;
113 myDSS = DSS(firstPoint);
114 ASSERT(myDSS.isExtendableFront(secondPoint));
115 myDSS.extendFront(secondPoint);
118 //-----------------------------------------------------------------------------
119 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
121 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
122 ArithmeticalDSSComputerOnSurfels ( const ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency> & other )
123 : myKSpace(other.myKSpace), myDim1(other.myDim1), myDim2(other.myDim2),
124 myProjection1(other.myProjection1), myProjection2(other.myProjection2), myProjectionNormal(other.myProjectionNormal),
125 myPreviousSurfelFront(other.myPreviousSurfelFront), myPreviousSurfelBack(other.myPreviousSurfelBack),
126 myDSS(other.myDSS), myBegin(other.myBegin), myEnd(other.myEnd)
130 //-----------------------------------------------------------------------------
131 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
133 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>&
134 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
135 operator=( const ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency> & other )
137 if ( this != &other )
139 myKSpace = other.myKSpace;
140 myDim1 = other.myDim1;
141 myDim2 = other.myDim2;
142 myProjection1 = other.myProjection1;
143 myProjection2 = other.myProjection2;
144 myProjectionNormal = other.myProjectionNormal;
145 myPreviousSurfelFront = other.myPreviousSurfelFront;
146 myPreviousSurfelBack = other.myPreviousSurfelBack;
148 myBegin = other.myBegin;
154 //-----------------------------------------------------------------------------
155 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
157 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Reverse
158 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>
161 return Reverse(*myKSpace, myDim1, myDim2);
164 //-----------------------------------------------------------------------------
165 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
167 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Self
168 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>
171 return Self(*myKSpace, myDim1, myDim2);
174 //-----------------------------------------------------------------------------
175 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
178 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
179 operator==( const ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>& other ) const
181 return ( (myBegin == other.myBegin)
182 && (myEnd == other.myEnd)
183 && (myDSS == other.myDSS) );
186 //-----------------------------------------------------------------------------
187 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
190 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
191 operator!=( const ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency> & other ) const
193 return (!(*this == other));
196 ///////////////////////////////////////////////////////////////////////////////
198 ///////////////////////////////////////////////////////////////////////////////
199 //--------------------------------------------------------------------
200 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
203 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::isExtendableFront()
206 if (! getOtherPointFromSurfel(myEnd, p, true, false))
209 return myDSS.isExtendableFront( p );
212 //--------------------------------------------------------------------
213 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
216 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::isExtendableBack()
218 ConstIterator it = myBegin;
221 if (! getOtherPointFromSurfel(it, p, false, false))
224 return myDSS.isExtendableBack( p );
227 //-----------------------------------------------------------------------------
228 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
231 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::extendFront()
234 if (! getOtherPointFromSurfel(myEnd, p, true, true))
237 if (myDSS.extendFront(p))
246 //--------------------------------------------------------------------
247 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
250 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::extendBack()
252 ConstIterator it = myBegin;
255 if (! getOtherPointFromSurfel(it, p, false, true))
258 if (myDSS.extendBack(p))
267 //--------------------------------------------------------------------
268 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
271 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::retractFront()
273 if (myDSS.retractFront())
282 //--------------------------------------------------------------------
283 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
286 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::retractBack()
288 if (myDSS.retractBack())
297 ///////////////////////////////////////////////////////////////////////////////
299 ///////////////////////////////////////////////////////////////////////////////
300 //-------------------------------------------------------------------------
301 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
303 const typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Primitive&
304 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::primitive() const
309 //-------------------------------------------------------------------------
310 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
313 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::a() const
318 //-------------------------------------------------------------------------
319 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
322 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::b() const
327 //-------------------------------------------------------------------------
328 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
331 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::mu() const
336 //-------------------------------------------------------------------------
337 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
340 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::omega() const
342 return myDSS.omega();
345 //-------------------------------------------------------------------------
346 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
348 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
349 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Uf() const
354 //-------------------------------------------------------------------------
355 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
357 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
358 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Ul() const
363 //-------------------------------------------------------------------------
364 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
366 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
367 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Lf() const
372 //-------------------------------------------------------------------------
373 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
375 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
376 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Ll() const
381 //-------------------------------------------------------------------------
382 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
384 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
385 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::back() const
390 //-------------------------------------------------------------------------
391 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
393 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
394 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::front() const
396 return myDSS.front();
399 //-------------------------------------------------------------------------
400 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
403 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::begin() const
408 //-------------------------------------------------------------------------
409 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
412 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::end() const
417 //-----------------------------------------------------------------
418 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
421 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::isValid() const
423 return ( (myDSS.isValid())&&(isNotEmpty(myBegin,myEnd)) );
426 //-----------------------------------------------------------------
427 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
430 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::selfDisplay ( std::ostream & out) const
432 out << "[ArithmeticalDSSComputerOnSurfels] " << myDSS;
435 //-----------------------------------------------------------------
436 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
439 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::commonLinel (Cell const& aSurfel1,
440 Cell const& aSurfel2,
443 ASSERT(myKSpace != nullptr);
445 typename KSpace::DirIterator q1_1 = myKSpace->uDirs(aSurfel1), q1_2 = q1_1;
447 typename KSpace::DirIterator q2_1 = myKSpace->uDirs(aSurfel2), q2_2 = q2_1;
450 std::set<Cell> linels1 = {
451 myKSpace->uIncident(aSurfel1, *q1_1, true),
452 myKSpace->uIncident(aSurfel1, *q1_1, false),
453 myKSpace->uIncident(aSurfel1, *q1_2, true),
454 myKSpace->uIncident(aSurfel1, *q1_2, false),
457 std::set<Cell> linels2 = {
458 myKSpace->uIncident(aSurfel2, *q2_1, true),
459 myKSpace->uIncident(aSurfel2, *q2_1, false),
460 myKSpace->uIncident(aSurfel2, *q2_2, true),
461 myKSpace->uIncident(aSurfel2, *q2_2, false),
464 std::vector<Cell> inter;
465 std::set_intersection(linels1.begin(), linels1.end(),
466 linels2.begin(), linels2.end(),
467 std::back_inserter(inter));
469 if (inter.size() == 1)
471 // The two surfels intersect on one linel
476 // The surfels don't intersect or are the same
480 //-----------------------------------------------------------------
481 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
484 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::getOtherPointFromSurfel(ConstIterator const& it,
487 bool aUpdatePrevious)
489 ASSERT(myKSpace != nullptr);
493 std::tie(p1, p2) = projectSurfel(surfel);
495 ConstIterator& previousSurfel = aIsFront ? myPreviousSurfelFront : myPreviousSurfelBack;
497 // Find the common unsigned linel between surfel and previousSurfel
499 if (! commonLinel(myKSpace->unsigns(surfel), myKSpace->unsigns(*previousSurfel), linel))
504 // For the next point, choose the point that is not part of the common linel
505 typename KSpace::DirIterator q_linel = myKSpace->uDirs(linel);
506 Point linel1 = projectInPlane(myKSpace->uCoords(myKSpace->uIncident(linel, *q_linel, true))),
507 linel2 = projectInPlane(myKSpace->uCoords(myKSpace->uIncident(linel, *q_linel, false)));
513 else if (p1 == linel2)
517 else if (p2 == linel1)
521 else if (p2 == linel2)
538 //-----------------------------------------------------------------
539 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
541 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
542 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::projectInPlane(Point3 const& aPoint) const
544 static const Point3 aOrigin = Point3::zero;
546 // Orthogonal projection on the plane with a given unit normal
547 Point3 h = (aPoint - aOrigin) - myProjectionNormal * (aPoint - aOrigin).dot(myProjectionNormal);
549 // We simply project the point on the plane defined by
550 // the two directions 'myProjection1' and 'myProjection2' passing through the origin point 'aOrigin'
551 return Point(h.dot(myProjection1), h.dot(myProjection2));
554 //-----------------------------------------------------------------
555 template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
557 std::pair<typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point,
558 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point>
559 DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::projectSurfel(SCell const& aSCell) const
561 ASSERT(myKSpace != nullptr);
563 typename KSpace::DirIterator q1 = myKSpace->sDirs(aSCell), q2 = q1;
566 // We pick 2 linels of the surfel
567 SCell linel1 = myKSpace->sIncident(aSCell, *q1, true),
568 linel2 = myKSpace->sIncident(aSCell, *q1, false);
570 // 4 points of the surfel
571 Point3 p1_1 = myKSpace->sCoords(myKSpace->sIncident(linel1, *q2, false)),
572 p1_2 = myKSpace->sCoords(myKSpace->sIncident(linel1, *q2, true)),
573 p2_1 = myKSpace->sCoords(myKSpace->sIncident(linel2, *q2, false)),
574 p2_2 = myKSpace->sCoords(myKSpace->sIncident(linel2, *q2, true));
576 std::set<Point> points;
577 points.insert(projectInPlane(p1_1));
578 points.insert(projectInPlane(p1_2));
579 points.insert(projectInPlane(p2_1));
580 points.insert(projectInPlane(p2_2));
582 ASSERT(points.size() == 2);
584 Point p1 = *points.begin(), p2 = *(++points.begin());