DGtal 1.3.0
Loading...
Searching...
No Matches
ArithmeticalDSSComputerOnSurfels.ih
1/**
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.
6 *
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.
11 *
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/>.
14 *
15 **/
16
17/**
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
21 *
22 * @date 2021/01/22
23 *
24 * Implementation of inline methods defined in ArithmeticalDSSComputerOnSurfels.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29///////////////////////////////////////////////////////////////////////////////
30// IMPLEMENTATION of inline methods.
31///////////////////////////////////////////////////////////////////////////////
32
33//////////////////////////////////////////////////////////////////////////////
34#include <cstdlib>
35#include <algorithm>
36#include <boost/version.hpp>
37#if BOOST_VERSION < 105800
38#include <boost/math/common_factor_rt.hpp>
39#else
40#include <boost/integer/common_factor_rt.hpp>
41#endif
42
43//////////////////////////////////////////////////////////////////////////////
44
45
46///////////////////////////////////////////////////////////////////////////////
47// Implementation of inline methods //
48
49//-----------------------------------------------------------------------------
50template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
51inline
52DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
53ArithmeticalDSSComputerOnSurfels()
54 : myKSpace(nullptr), myDSS( Point(0,0) ), myBegin(), myEnd()
55{
56}
57
58//-----------------------------------------------------------------------------
59template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
60inline
61DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
62ArithmeticalDSSComputerOnSurfels(const KSpace& aKSpace, Dimension aDim1, Dimension aDim2)
63 : myKSpace(&aKSpace), myDim1(aDim1), myDim2(aDim2), myDSS( Point(0,0) ), myBegin(), myEnd()
64{
65 // Initialize projection vectors
66 myProjection1 = Point3::zero;
67 myProjection2 = Point3::zero;
68 myProjection1[myDim1] = 1;
69 myProjection2[myDim2] = 1;
70
71 myProjectionNormal = myProjection1.crossProduct(myProjection2);
72}
73
74
75//-----------------------------------------------------------------------------
76template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
77inline
78void DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
79init(const ConstIterator& it)
80{
81 ASSERT(myKSpace != nullptr);
82
83 myBegin = it;
84 myEnd = it;
85 ++myEnd;
86
87 SCell s = *it;
88 auto initialPoints = projectSurfel(s);
89 myPreviousSurfelFront = it;
90 myPreviousSurfelBack = it;
91
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;
108 } else {
109 ASSERT(false);
110 }
111 }
112
113 myDSS = DSS(firstPoint);
114 ASSERT(myDSS.isExtendableFront(secondPoint));
115 myDSS.extendFront(secondPoint);
116}
117
118//-----------------------------------------------------------------------------
119template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
120inline
121DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
122ArithmeticalDSSComputerOnSurfels ( 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)
127{
128}
129
130//-----------------------------------------------------------------------------
131template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
132inline
133typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>&
134DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
135operator=( const ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency> & other )
136{
137 if ( this != &other )
138 {
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;
147 myDSS = other.myDSS;
148 myBegin = other.myBegin;
149 myEnd = other.myEnd;
150 }
151 return *this;
152}
153
154//-----------------------------------------------------------------------------
155template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
156inline
157typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Reverse
158DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>
159::getReverse() const
160{
161 return Reverse(*myKSpace, myDim1, myDim2);
162}
163
164//-----------------------------------------------------------------------------
165template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
166inline
167typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Self
168DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>
169::getSelf() const
170{
171 return Self(*myKSpace, myDim1, myDim2);
172}
173
174//-----------------------------------------------------------------------------
175template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
176inline
177bool
178DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
179operator==( const ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>& other ) const
180{
181 return ( (myBegin == other.myBegin)
182 && (myEnd == other.myEnd)
183 && (myDSS == other.myDSS) );
184}
185
186//-----------------------------------------------------------------------------
187template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
188inline
189bool
190DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::
191operator!=( const ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency> & other ) const
192{
193 return (!(*this == other));
194}
195
196///////////////////////////////////////////////////////////////////////////////
197// Update methods //
198///////////////////////////////////////////////////////////////////////////////
199//--------------------------------------------------------------------
200template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
201inline
202bool
203DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::isExtendableFront()
204{
205 Point p;
206 if (! getOtherPointFromSurfel(myEnd, p, true, false))
207 return false;
208
209 return myDSS.isExtendableFront( p );
210}
211
212//--------------------------------------------------------------------
213template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
214inline
215bool
216DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::isExtendableBack()
217{
218 ConstIterator it = myBegin;
219 --it;
220 Point p;
221 if (! getOtherPointFromSurfel(it, p, false, false))
222 return false;
223
224 return myDSS.isExtendableBack( p );
225}
226
227//-----------------------------------------------------------------------------
228template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
229inline
230bool
231DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::extendFront()
232{
233 Point p;
234 if (! getOtherPointFromSurfel(myEnd, p, true, true))
235 return false;
236
237 if (myDSS.extendFront(p))
238 {
239 ++myEnd;
240 return true;
241 }
242 else
243 return false;
244}
245
246//--------------------------------------------------------------------
247template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
248inline
249bool
250DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::extendBack()
251{
252 ConstIterator it = myBegin;
253 --it;
254 Point p;
255 if (! getOtherPointFromSurfel(it, p, false, true))
256 return false;
257
258 if (myDSS.extendBack(p))
259 {
260 myBegin = it;
261 return true;
262 }
263 else
264 return false;
265}
266
267//--------------------------------------------------------------------
268template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
269inline
270bool
271DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::retractFront()
272{
273 if (myDSS.retractFront())
274 {
275 --myEnd;
276 return true;
277 }
278 else
279 return false;
280}
281
282//--------------------------------------------------------------------
283template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
284inline
285bool
286DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::retractBack()
287{
288 if (myDSS.retractBack())
289 {
290 ++myBegin;
291 return true;
292 }
293 else
294 return false;
295}
296
297///////////////////////////////////////////////////////////////////////////////
298// Accessors //
299///////////////////////////////////////////////////////////////////////////////
300//-------------------------------------------------------------------------
301template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
302inline
303const typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Primitive&
304DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::primitive() const
305{
306 return myDSS;
307}
308
309//-------------------------------------------------------------------------
310template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
311inline
312TInteger
313DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::a() const
314{
315 return myDSS.a();
316}
317
318//-------------------------------------------------------------------------
319template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
320inline
321TInteger
322DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::b() const
323{
324 return myDSS.b();
325}
326
327//-------------------------------------------------------------------------
328template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
329inline
330TInteger
331DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::mu() const
332{
333 return myDSS.mu();
334}
335
336//-------------------------------------------------------------------------
337template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
338inline
339TInteger
340DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::omega() const
341{
342 return myDSS.omega();
343}
344
345//-------------------------------------------------------------------------
346template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
347inline
348typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
349DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Uf() const
350{
351 return myDSS.Uf();
352}
353
354//-------------------------------------------------------------------------
355template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
356inline
357typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
358DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Ul() const
359{
360 return myDSS.Ul();
361}
362
363//-------------------------------------------------------------------------
364template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
365inline
366typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
367DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Lf() const
368{
369 return myDSS.Lf();
370}
371
372//-------------------------------------------------------------------------
373template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
374inline
375typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
376DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Ll() const
377{
378 return myDSS.Ll();
379}
380
381//-------------------------------------------------------------------------
382template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
383inline
384typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
385DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::back() const
386{
387 return myDSS.back();
388}
389
390//-------------------------------------------------------------------------
391template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
392inline
393typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
394DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::front() const
395{
396 return myDSS.front();
397}
398
399//-------------------------------------------------------------------------
400template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
401inline
402TIterator
403DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::begin() const
404{
405 return myBegin;
406}
407
408//-------------------------------------------------------------------------
409template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
410inline
411TIterator
412DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::end() const
413{
414 return myEnd;
415}
416
417//-----------------------------------------------------------------
418template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
419inline
420bool
421DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::isValid() const
422{
423 return ( (myDSS.isValid())&&(isNotEmpty(myBegin,myEnd)) );
424}
425
426//-----------------------------------------------------------------
427template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
428inline
429void
430DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::selfDisplay ( std::ostream & out) const
431{
432 out << "[ArithmeticalDSSComputerOnSurfels] " << myDSS;
433}
434
435//-----------------------------------------------------------------
436template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
437inline
438bool
439DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::commonLinel (Cell const& aSurfel1,
440 Cell const& aSurfel2,
441 Cell& aLinel)
442{
443 ASSERT(myKSpace != nullptr);
444
445 typename KSpace::DirIterator q1_1 = myKSpace->uDirs(aSurfel1), q1_2 = q1_1;
446 ++q1_2;
447 typename KSpace::DirIterator q2_1 = myKSpace->uDirs(aSurfel2), q2_2 = q2_1;
448 ++q2_2;
449
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),
455 };
456
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),
462 };
463
464 std::vector<Cell> inter;
465 std::set_intersection(linels1.begin(), linels1.end(),
466 linels2.begin(), linels2.end(),
467 std::back_inserter(inter));
468
469 if (inter.size() == 1)
470 {
471 // The two surfels intersect on one linel
472 aLinel = inter[0];
473 return true;
474 }
475
476 // The surfels don't intersect or are the same
477 return false;
478}
479
480//-----------------------------------------------------------------
481template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
482inline
483bool
484DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::getOtherPointFromSurfel(ConstIterator const& it,
485 Point& aPoint,
486 bool aIsFront,
487 bool aUpdatePrevious)
488{
489 ASSERT(myKSpace != nullptr);
490
491 SCell surfel = *it;
492 Point p1, p2;
493 std::tie(p1, p2) = projectSurfel(surfel);
494
495 ConstIterator& previousSurfel = aIsFront ? myPreviousSurfelFront : myPreviousSurfelBack;
496
497 // Find the common unsigned linel between surfel and previousSurfel
498 Cell linel;
499 if (! commonLinel(myKSpace->unsigns(surfel), myKSpace->unsigns(*previousSurfel), linel))
500 {
501 return false;
502 }
503
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)));
508
509 if (p1 == linel1)
510 {
511 aPoint = p2;
512 }
513 else if (p1 == linel2)
514 {
515 aPoint = p2;
516 }
517 else if (p2 == linel1)
518 {
519 aPoint = p1;
520 }
521 else if (p2 == linel2)
522 {
523 aPoint = p1;
524 }
525 else
526 {
527 ASSERT(false);
528 }
529
530 if (aUpdatePrevious)
531 {
532 previousSurfel = it;
533 }
534
535 return true;
536}
537
538//-----------------------------------------------------------------
539template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
540inline
541typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point
542DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::projectInPlane(Point3 const& aPoint) const
543{
544 static const Point3 aOrigin = Point3::zero;
545
546 // Orthogonal projection on the plane with a given unit normal
547 Point3 h = (aPoint - aOrigin) - myProjectionNormal * (aPoint - aOrigin).dot(myProjectionNormal);
548
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));
552}
553
554//-----------------------------------------------------------------
555template <typename TKSpace, typename TIterator, typename TInteger, unsigned short adjacency>
556inline
557std::pair<typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point,
558 typename DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::Point>
559DGtal::ArithmeticalDSSComputerOnSurfels<TKSpace,TIterator,TInteger,adjacency>::projectSurfel(SCell const& aSCell) const
560{
561 ASSERT(myKSpace != nullptr);
562
563 typename KSpace::DirIterator q1 = myKSpace->sDirs(aSCell), q2 = q1;
564 ++q2;
565
566 // We pick 2 linels of the surfel
567 SCell linel1 = myKSpace->sIncident(aSCell, *q1, true),
568 linel2 = myKSpace->sIncident(aSCell, *q1, false);
569
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));
575
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));
581
582 ASSERT(points.size() == 2);
583
584 Point p1 = *points.begin(), p2 = *(++points.begin());
585
586 return { p1, p2 };
587}