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 StandardDSLQ0.ih
19 * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20 * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
24 * Implementation of inline methods defined in StandardDSLQ0.h
26 * This file is part of the DGtal library.
30//////////////////////////////////////////////////////////////////////////////
32//////////////////////////////////////////////////////////////////////////////
34///////////////////////////////////////////////////////////////////////////////
35// IMPLEMENTATION of inline methods.
36///////////////////////////////////////////////////////////////////////////////
38///////////////////////////////////////////////////////////////////////////////
39// ----------------------- Standard services ------------------------------
44//-----------------------------------------------------------------------------
45template <typename TFraction>
47DGtal::StandardDSLQ0<TFraction>::
51//-----------------------------------------------------------------------------
52template <typename TFraction>
54DGtal::StandardDSLQ0<TFraction>::
59//-----------------------------------------------------------------------------
60template <typename TFraction>
62DGtal::StandardDSLQ0<TFraction>::
63StandardDSLQ0 ( const Self & other )
64 : myPattern( other.myPattern ), myMu( other.myMu )
67//-----------------------------------------------------------------------------
68template <typename TFraction>
70DGtal::StandardDSLQ0<TFraction> &
71DGtal::StandardDSLQ0<TFraction>::
72operator= ( const Self & other )
74 myPattern = other.myPattern;
78//-----------------------------------------------------------------------------
79template <typename TFraction>
81DGtal::StandardDSLQ0<TFraction>::
82StandardDSLQ0( Fraction aSlope, IntegerParamType aMu )
83 : myPattern( aSlope ), myMu( aMu )
86//-----------------------------------------------------------------------------
87template <typename TFraction>
89DGtal::StandardDSLQ0<TFraction>::
90StandardDSLQ0( IntegerParamType a1, IntegerParamType b1, IntegerParamType mu1 )
91 : myPattern( a1, b1 ), myMu( mu1 )
94//-----------------------------------------------------------------------------
95template <typename TFraction>
98DGtal::StandardDSLQ0<TFraction>::
99operator()( const Point & p ) const
101 if ( slope().null() ) return false;
103 return ( mu() <= _r ) && ( _r < ( mu() + pattern().length() ) );
105//-----------------------------------------------------------------------------
106template <typename TFraction>
108typename DGtal::StandardDSLQ0<TFraction>::Integer
109DGtal::StandardDSLQ0<TFraction>::
110r( const Point & p ) const
112 ASSERT( ! slope().null() );
113 return a() * p[ 0 ] - b() * p[ 1 ];
115//-----------------------------------------------------------------------------
116template <typename TFraction>
118typename DGtal::StandardDSLQ0<TFraction>::Fraction
119DGtal::StandardDSLQ0<TFraction>::
122 return pattern().slope();
124//-----------------------------------------------------------------------------
125template <typename TFraction>
127const DGtal::Pattern<TFraction> &
128DGtal::StandardDSLQ0<TFraction>::
133//-----------------------------------------------------------------------------
134template <typename TFraction>
136const typename DGtal::StandardDSLQ0<TFraction>::Integer &
137DGtal::StandardDSLQ0<TFraction>::
142//-----------------------------------------------------------------------------
143template <typename TFraction>
145typename DGtal::StandardDSLQ0<TFraction>::Integer
146DGtal::StandardDSLQ0<TFraction>::
149 return myMu + pattern().length() - NumberTraits<Integer>::ONE;
151//-----------------------------------------------------------------------------
152template <typename TFraction>
154typename DGtal::StandardDSLQ0<TFraction>::Integer
155DGtal::StandardDSLQ0<TFraction>::
160//-----------------------------------------------------------------------------
161template <typename TFraction>
163typename DGtal::StandardDSLQ0<TFraction>::Integer
164DGtal::StandardDSLQ0<TFraction>::
169//-----------------------------------------------------------------------------
170template <typename TFraction>
172typename DGtal::StandardDSLQ0<TFraction>::Vector2I
173DGtal::StandardDSLQ0<TFraction>::
176 return pattern().v();
178//-----------------------------------------------------------------------------
179template <typename TFraction>
181typename DGtal::StandardDSLQ0<TFraction>::Point
182DGtal::StandardDSLQ0<TFraction>::
185 Vector2I u = pattern().bezout();
186 Integer c = ( mu() * u[ 0 ] ) / b();
187 return Point( mu() > NumberTraits<Integer>::ZERO
188 ? v() * ( c + 1 ) - u * mu()
189 : u * mu() - v() * c );
191//-----------------------------------------------------------------------------
192template <typename TFraction>
194typename DGtal::StandardDSLQ0<TFraction>::Point
195DGtal::StandardDSLQ0<TFraction>::
198 return Point( U() + pattern().L( NumberTraits<Quotient>::ZERO ) );
200//-----------------------------------------------------------------------------
201template <typename TFraction>
203typename DGtal::StandardDSLQ0<TFraction>::Point
204DGtal::StandardDSLQ0<TFraction>::
205lowestY( IntegerParamType x ) const
207 Integer q = a() * x - mup();
208 return Point( x, ic.ceilDiv( q, b() ) );
210//-----------------------------------------------------------------------------
211template <typename TFraction>
213typename DGtal::StandardDSLQ0<TFraction>::Point
214DGtal::StandardDSLQ0<TFraction>::
215uppermostY( IntegerParamType x ) const
217 Integer q = a() * x - mu();
218 return Point( x, ic.floorDiv( q, b() ) );
220//-----------------------------------------------------------------------------
221template <typename TFraction>
223typename DGtal::StandardDSLQ0<TFraction>::Point
224DGtal::StandardDSLQ0<TFraction>::
225lowestX( IntegerParamType y ) const
227 Integer q = mu() + b() * y;
228 return Point( ic.ceilDiv( q, a() ), y );
230//-----------------------------------------------------------------------------
231template <typename TFraction>
233typename DGtal::StandardDSLQ0<TFraction>::Point
234DGtal::StandardDSLQ0<TFraction>::
235uppermostX( IntegerParamType y ) const
237 Integer q = mup() + b() * y;
238 return Point( ic.floorDiv( q, a() ), y );
240//-----------------------------------------------------------------------------
241template <typename TFraction>
244DGtal::StandardDSLQ0<TFraction>::
245before( const Point & p1, const Point & p2 ) const
247 return ( p1[ 0 ] < p2[ 0 ] )
248 || ( ( p1[ 0 ] == p2[ 0 ] ) && ( p1[ 1 ] < p2[ 1 ] ) );
250//-----------------------------------------------------------------------------
251template <typename TFraction>
254DGtal::StandardDSLQ0<TFraction>::
255beforeOrEqual( const Point & p1, const Point & p2 ) const
257 return ( p1[ 0 ] < p2[ 0 ] )
258 || ( ( p1[ 0 ] == p2[ 0 ] ) && ( p1[ 1 ] <= p2[ 1 ] ) );
260//-----------------------------------------------------------------------------
261template <typename TFraction>
262typename DGtal::StandardDSLQ0<TFraction>::ConstIterator
263DGtal::StandardDSLQ0<TFraction>::
264begin( Point p ) const
266 return ConstIterator( *this, p );
268//-----------------------------------------------------------------------------
269template <typename TFraction>
270typename DGtal::StandardDSLQ0<TFraction>::ConstIterator
271DGtal::StandardDSLQ0<TFraction>::
274 ConstIterator it( *this, p );
278//-----------------------------------------------------------------------------
279template <typename TFraction>
280typename DGtal::StandardDSLQ0<TFraction>::Self
281DGtal::StandardDSLQ0<TFraction>::
282reversedSmartDSS( const Point & A, const Point & B ) const
285 Integer cA = ic.floorDiv( A[ 0 ] - _U[ 0 ], v()[ 0 ] );
286 Point U1 = _U + v() * cA;
287 Integer cB = ic.ceilDiv( B[ 0 ] - _U[ 0 ], v()[ 0 ] );
288 Point U2 = _U + v() * cB;
289 if ( before( A, U1 ) ) U1 -= v();
290 if ( before( U2, B ) ) U2 += v();
291 return reversedSmartDSS( U1, U2, A, B );
294//-----------------------------------------------------------------------------
295template <typename TFraction>
296typename DGtal::StandardDSLQ0<TFraction>::Self
297DGtal::StandardDSLQ0<TFraction>::
298reversedSmartDSS( Point U1, Point U2,
299 const Point & A, const Point & B ) const
302 std::cerr << "[reversedSmartDSS] " << (*this)
303 << " " << pattern().rE() << std::endl
304 << " +- U1=" << U1 << " A=" << A
305 << " B=" << B << " U2=" << U2 << std::endl
306 << " v()=" << pattern().v()
307 << " u()=" << pattern().bezout()
308 << " r(U())=" << r(U())
309 << " mu=" << mu() << " r(U1)=" << r(U1)
310 << " r(U2)=" << r(U2)
314 << " DSS(A)=" << this->operator()( A )
315 << " DSS(B)=" << this->operator()( B ) << std::endl;
317 ASSERT( ! slope().null() );
318 ASSERT( r( U1 ) == mu() && r( U2 ) == mu()
319 && this->operator()( A ) && this->operator()( B ) );
320 ASSERT( beforeOrEqual( U1, A ) );
321 ASSERT( before( A, B ) );
322 ASSERT( beforeOrEqual( B, U2 ) );
323 if ( A[ 0 ] == B[ 0 ] ) return Self( 1, 0, A[ 0 ] );
324 if ( A[ 1 ] == B[ 1 ] ) return Self( 0, 1, -A[ 1 ] );
325 Integer dU = U2[ 0 ] - U1[ 0 ];
328 std::cerr << " +- dU=" << dU << std::endl;
330 if ( ( dU >= (3*b()) )
331 || ( ( dU == (2*b()) ) && ( A == U1 || B == U2 ) )
332 || ( A == U1 && B == U2 ) )
335 std::cerr << "[reversedSmartDSS] 3 patterns || ..." << std::endl;
339 // [A,B] is included in two patterns of this DSL.
343 std::cerr << "[reversedSmartDSS] 2 patterns" << std::endl;
345 return DSSWithinTwoPatterns( U1, U2, A, B );
347 // [A,B] is included in one patterns of this DSL.
348 Pattern<Fraction> subpattern;
351 Integer posA = ( A - U1 ).norm1();
352 Integer posB = ( B - U1 ).norm1();
354 std::cerr << "[reversedSmartDSS] 1 pattern" << std::endl;
356 bool m = pattern().getSmallestCoveringSubpattern
357 ( subpattern, nb, startPos, posA, posB );
359 std::cerr << " - smallest:" << subpattern.rE() << " at " << startPos << endl;
362 { // smallest covering subpattern is the pattern itself.
363 bool n = pattern().getGreatestIncludedSubpattern
364 ( subpattern, nb, startPos, posA, posB );
366 std::cerr << " - greatest:" << subpattern.rE() << " at " << startPos << endl;
368 ASSERT( n ); boost::ignore_unused_variable_warning(n);
369 Point NU( U1 + startPos );
370 Integer nmu = subpattern.slope().p() * NU[ 0 ]
371 - subpattern.slope().q() * NU[ 1 ];
372 return Self( subpattern.slope(), nmu );
374 // last case, recursive call.
375 Point NU1( U1 + startPos );
376 Point NU2( NU1 + subpattern.v()*nb );
377 Integer nmu = subpattern.slope().p() * NU1[ 0 ]
378 - subpattern.slope().q() * NU1[ 1 ];
379 StandardDSLQ0 ndsl( subpattern.slope(), nmu );
380 return ndsl.reversedSmartDSS( NU1, NU2, A, B );
382//-----------------------------------------------------------------------------
383template <typename TFraction>
384typename DGtal::StandardDSLQ0<TFraction>::Self
385DGtal::StandardDSLQ0<TFraction>::
386DSSWithinTwoPatterns( Point U1, Point U2,
387 const Point & A, const Point & B ) const
390 Point Um = U1 + pattern().v();
391 ASSERT( Um + pattern().v() == U2 );
392 ASSERT( before( A, B ) );
393 ASSERT( before( A, Um ) );
394 ASSERT( before( Um, B ) );
395 ASSERT( r( U1 ) == mu() && r( U2 ) == mu() );
396 bool readyLU = false;
397 bool readyRU = false;
399 Point L1 = U1 + pattern().L( 0 );
400 Point L2 = L1 + pattern().v();
401 Pattern<Fraction> subpattern;
402 Pattern<Fraction> patDeepest;
403 Pattern<Fraction> patLU = pattern();
404 Pattern<Fraction> patRU = pattern();
405 Pattern<Fraction> patL = pattern();
409 std::cerr << "[DSSWithinTwoPatterns] " << (*this)
410 << " " << pattern().rE() << std::endl
411 << " +- U1=" << U1 << " A=" << A
412 << " B=" << B << " U2=" << U2
413 << " L1=" << L1 << std::endl;
415 while( true ) //for ( Quotient i = NumberTraits<Quotient>::ZERO; i <= pattern().slope().k(); ++i )
419 bool mLU = patLU.getSmallestCoveringSubpattern
420 ( subpattern, nb, startPos,
421 ( A - U1 ).norm1(), ( Um - U1 ).norm1(), false );
422 if ( ! mLU || nb > 1 )
424 bool n = patLU.getGreatestIncludedSubpattern
425 ( subpattern, nb, startPos,
426 ( A - U1 ).norm1(), ( Um - U1 ).norm1(), false );
427 ASSERT( n ); boost::ignore_unused_variable_warning(n);
435 bool mRU = patRU.getSmallestCoveringSubpattern
436 ( subpattern, nb, startPos,
437 0, ( B - Um ).norm1(), false );
438 if ( ! mRU || nb > 1 )
440 bool n = patRU.getGreatestIncludedSubpattern
441 ( subpattern, nb, startPos,
442 0, ( B - Um ).norm1(), false );
443 ASSERT( n ); boost::ignore_unused_variable_warning(n);
447 ASSERT( ! patRU.slope().null() );
448 U2 = Um + patRU.v() * nb;
452 posA = L1[ 0 ] <= A[ 0 ] ? ( A - L1 ).norm1() : 0;
453 posB = L2[ 0 ] > B[ 0 ] ? ( B - L1 ).norm1() : patL.length();
454 bool mL = patL.getSmallestCoveringSubpattern
455 ( subpattern, nb, startPos, posA, posB, true );
458 bool n = patL.getGreatestIncludedSubpattern
459 ( subpattern, nb, startPos, posA, posB, true );
460 ASSERT( n ); boost::ignore_unused_variable_warning(n);
465 { // One has to keep the pertinent pattern
466 // It is the reversed pattern around Um.
468 L2 = Um + patL.L( 0 );
471 // patL = subpattern;
473 // L2 = L1 + patL.v() * nb;
476 std::cerr << " (*) " << (readyLU ? '*' : '-')
477 << "LU=" << patLU.rE() << " at " << U1 << std::endl;
478 std::cerr << " (*) " << (readyRU ? '*' : '-')
479 << "RU=" << patRU.rE() << " til " << U2 << std::endl;
480 std::cerr << " (*) " << (readyL ? '*' : '-')
481 << "L =" << patL.rE() << " at " << L1 << std::endl;
483 if ( readyLU || readyRU || readyL )
485 patDeepest = deepest( patLU.slope(), patRU.slope(), patL.slope() );
487 std::cerr << " => deepest is " << patDeepest.rE() << std::endl;
490 if ( ( readyLU && patDeepest.slope().q() == patLU.slope().q() )
491 || ( readyRU && patDeepest.slope().q() == patRU.slope().q() )
492 || ( readyL && patDeepest.slope().q() == patL.slope().q() ) )
495 Integer nmu = patDeepest.slope().p() * Um[ 0 ]
496 - patDeepest.slope().q() * Um[ 1 ];
497 return StandardDSLQ0( patDeepest.slope(), nmu );
499//-----------------------------------------------------------------------------
500template <typename TFraction>
501typename DGtal::StandardDSLQ0<TFraction>::Self
502DGtal::StandardDSLQ0<TFraction>::
503smartDSS( const Point & A, const Point & B ) const
506 std::cerr << "[smartDSS] " << (*this)
507 << " " << pattern().rE() << std::endl
508 << " A=" << A << " B=" << B << std::endl;
510 ASSERT( ! slope().null() );
511 ASSERT( this->operator()( A ) && this->operator()( B ) );
512 ASSERT( before( A, B ) );
513 Fraction f10( 1, 0 );
514 Pattern<Fraction> p( 0, 1 );
520 Point2I _Up = _U + Point2I(0,1);
521 Point2I _Lp = _L + Point2I(1,-1);
522 UnsignedInteger AB1 = (B-A).norm1();
523 while ( ( (_Up - A).norm1() <= AB1 )
524 && this->operator()( _Up ) )
527 std::cerr << "Vertical" << std::endl;
529 p = Pattern<Fraction>( f10 );
540 while ( p.slope() != this->slope() )
543 std::cerr << "[smartDSS] v=" << p.v()
544 << " bez=" << p.bezout()
545 << " U=(" << _U[0] << "," << _U[1] << ")"
546 << " L=(" << _L[0] << "," << _L[1] << ")"
547 << " Up=(" << _Up[0] << "," << _Up[1] << ")"
548 << " Lp=(" << _Lp[0] << "," << _Lp[1] << ")"
551 ASSERT( p.v()[1]*p.bezout()[0] - p.v()[0]*p.bezout()[1] == -1 );
552 if ( ( (_Up - A).norm1() > AB1 ) && ( (_Lp - A).norm1() > AB1 ) ) break;
553 else if ( _Up[ 1 ] <= B[ 1 ] && this->operator()( _Up ) )
555 Fraction np = p.slope().right();
556 for ( Quotient i = 1; i < delta; ++i )
558 _L = _Lp + p.bezout() - p.v();
559 if ( ! lul ) _L -= p.v();
560 p = Pattern<Fraction>( np );
561 ASSERT( p.v()[1]*p.bezout()[0] - p.v()[0]*p.bezout()[1] == -1 );
562 _Up = _U + p.v() + p.bezout();
563 _Lp = _L + p.v() + p.v() - p.bezout();
564 delta = 1; ulu = true; lul = false;
566 else if ( _Lp[ 0 ] <= B[ 0 ] && this->operator()( _Lp ) )
568 Fraction np = p.slope().left();
569 for ( Quotient i = 1; i < delta; ++i )
571 _U = p.slope() == f10 ? _Up - Point2I( 0,1 ) : _Up - p.bezout();
572 if ( ! ulu ) _U -= p.v();
573 p = Pattern<Fraction>( np );
574 ASSERT( p.v()[1]*p.bezout()[0] - p.v()[0]*p.bezout()[1] == -1 );
575 _Up = _U + p.v() + p.bezout();
576 _Lp = _L + p.v() + p.v() - p.bezout();
577 delta = 1; ulu = false; lul = true;
586 Integer nmu = p.slope().p() * _U[ 0 ] - p.slope().q() * _U[ 1 ];
587 return StandardDSLQ0( p.slope(), nmu );
590//-----------------------------------------------------------------------------
591template <typename TFraction>
593typename DGtal::StandardDSLQ0<TFraction>::Fraction
594DGtal::StandardDSLQ0<TFraction>::
595deepest( Fraction f1, Fraction f2, Fraction f3 )
597 return deepest( f1, deepest( f2, f3 ) );
599//-----------------------------------------------------------------------------
600template <typename TFraction>
602typename DGtal::StandardDSLQ0<TFraction>::Fraction
603DGtal::StandardDSLQ0<TFraction>::
604deepest( Fraction f1, Fraction f2 )
606 return ( ( f1.k() > f2.k() )
607 || ( ( f1.k() == f2.k() ) && ( f1.u() >= f2.u() ) ) )
611///////////////////////////////////////////////////////////////////////////////
612// Interface - public :
615 * Writes/Displays the object on an output stream.
616 * @param out the output stream where the object is written.
618template <typename TFraction>
621DGtal::StandardDSLQ0<TFraction>::selfDisplay ( std::ostream & out ) const
623 out << "[StandardDSLQ0"
624 << " a=" << a() << ", b=" << b() << ", mu=" << mu() << "]";
628 * Checks the validity/consistency of the object.
629 * @return 'true' if the object is valid, 'false' otherwise.
631template <typename TFraction>
634DGtal::StandardDSLQ0<TFraction>::isValid() const
641///////////////////////////////////////////////////////////////////////////////
642// Implementation of inline functions //
644template <typename TFraction>
647DGtal::operator<< ( std::ostream & out,
648 const StandardDSLQ0<TFraction> & object )
650 object.selfDisplay( out );
655///////////////////////////////////////////////////////////////////////////////