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 //-----------------------------------------------------------------------------
45 template <typename TFraction>
47 DGtal::StandardDSLQ0<TFraction>::
51 //-----------------------------------------------------------------------------
52 template <typename TFraction>
54 DGtal::StandardDSLQ0<TFraction>::
59 //-----------------------------------------------------------------------------
60 template <typename TFraction>
62 DGtal::StandardDSLQ0<TFraction>::
63 StandardDSLQ0 ( const Self & other )
64 : myPattern( other.myPattern ), myMu( other.myMu )
67 //-----------------------------------------------------------------------------
68 template <typename TFraction>
70 DGtal::StandardDSLQ0<TFraction> &
71 DGtal::StandardDSLQ0<TFraction>::
72 operator= ( const Self & other )
74 myPattern = other.myPattern;
78 //-----------------------------------------------------------------------------
79 template <typename TFraction>
81 DGtal::StandardDSLQ0<TFraction>::
82 StandardDSLQ0( Fraction aSlope, IntegerParamType aMu )
83 : myPattern( aSlope ), myMu( aMu )
86 //-----------------------------------------------------------------------------
87 template <typename TFraction>
89 DGtal::StandardDSLQ0<TFraction>::
90 StandardDSLQ0( IntegerParamType a1, IntegerParamType b1, IntegerParamType mu1 )
91 : myPattern( a1, b1 ), myMu( mu1 )
94 //-----------------------------------------------------------------------------
95 template <typename TFraction>
98 DGtal::StandardDSLQ0<TFraction>::
99 operator()( const Point & p ) const
101 if ( slope().null() ) return false;
103 return ( mu() <= _r ) && ( _r < ( mu() + pattern().length() ) );
105 //-----------------------------------------------------------------------------
106 template <typename TFraction>
108 typename DGtal::StandardDSLQ0<TFraction>::Integer
109 DGtal::StandardDSLQ0<TFraction>::
110 r( const Point & p ) const
112 ASSERT( ! slope().null() );
113 return a() * p[ 0 ] - b() * p[ 1 ];
115 //-----------------------------------------------------------------------------
116 template <typename TFraction>
118 typename DGtal::StandardDSLQ0<TFraction>::Fraction
119 DGtal::StandardDSLQ0<TFraction>::
122 return pattern().slope();
124 //-----------------------------------------------------------------------------
125 template <typename TFraction>
127 const DGtal::Pattern<TFraction> &
128 DGtal::StandardDSLQ0<TFraction>::
133 //-----------------------------------------------------------------------------
134 template <typename TFraction>
136 const typename DGtal::StandardDSLQ0<TFraction>::Integer &
137 DGtal::StandardDSLQ0<TFraction>::
142 //-----------------------------------------------------------------------------
143 template <typename TFraction>
145 typename DGtal::StandardDSLQ0<TFraction>::Integer
146 DGtal::StandardDSLQ0<TFraction>::
149 return myMu + pattern().length() - NumberTraits<Integer>::ONE;
151 //-----------------------------------------------------------------------------
152 template <typename TFraction>
154 typename DGtal::StandardDSLQ0<TFraction>::Integer
155 DGtal::StandardDSLQ0<TFraction>::
160 //-----------------------------------------------------------------------------
161 template <typename TFraction>
163 typename DGtal::StandardDSLQ0<TFraction>::Integer
164 DGtal::StandardDSLQ0<TFraction>::
169 //-----------------------------------------------------------------------------
170 template <typename TFraction>
172 typename DGtal::StandardDSLQ0<TFraction>::Vector2I
173 DGtal::StandardDSLQ0<TFraction>::
176 return pattern().v();
178 //-----------------------------------------------------------------------------
179 template <typename TFraction>
181 typename DGtal::StandardDSLQ0<TFraction>::Point
182 DGtal::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 //-----------------------------------------------------------------------------
192 template <typename TFraction>
194 typename DGtal::StandardDSLQ0<TFraction>::Point
195 DGtal::StandardDSLQ0<TFraction>::
198 return Point( U() + pattern().L( NumberTraits<Quotient>::ZERO ) );
200 //-----------------------------------------------------------------------------
201 template <typename TFraction>
203 typename DGtal::StandardDSLQ0<TFraction>::Point
204 DGtal::StandardDSLQ0<TFraction>::
205 lowestY( IntegerParamType x ) const
207 Integer q = a() * x - mup();
208 return Point( x, ic.ceilDiv( q, b() ) );
210 //-----------------------------------------------------------------------------
211 template <typename TFraction>
213 typename DGtal::StandardDSLQ0<TFraction>::Point
214 DGtal::StandardDSLQ0<TFraction>::
215 uppermostY( IntegerParamType x ) const
217 Integer q = a() * x - mu();
218 return Point( x, ic.floorDiv( q, b() ) );
220 //-----------------------------------------------------------------------------
221 template <typename TFraction>
223 typename DGtal::StandardDSLQ0<TFraction>::Point
224 DGtal::StandardDSLQ0<TFraction>::
225 lowestX( IntegerParamType y ) const
227 Integer q = mu() + b() * y;
228 return Point( ic.ceilDiv( q, a() ), y );
230 //-----------------------------------------------------------------------------
231 template <typename TFraction>
233 typename DGtal::StandardDSLQ0<TFraction>::Point
234 DGtal::StandardDSLQ0<TFraction>::
235 uppermostX( IntegerParamType y ) const
237 Integer q = mup() + b() * y;
238 return Point( ic.floorDiv( q, a() ), y );
240 //-----------------------------------------------------------------------------
241 template <typename TFraction>
244 DGtal::StandardDSLQ0<TFraction>::
245 before( const Point & p1, const Point & p2 ) const
247 return ( p1[ 0 ] < p2[ 0 ] )
248 || ( ( p1[ 0 ] == p2[ 0 ] ) && ( p1[ 1 ] < p2[ 1 ] ) );
250 //-----------------------------------------------------------------------------
251 template <typename TFraction>
254 DGtal::StandardDSLQ0<TFraction>::
255 beforeOrEqual( const Point & p1, const Point & p2 ) const
257 return ( p1[ 0 ] < p2[ 0 ] )
258 || ( ( p1[ 0 ] == p2[ 0 ] ) && ( p1[ 1 ] <= p2[ 1 ] ) );
260 //-----------------------------------------------------------------------------
261 template <typename TFraction>
262 typename DGtal::StandardDSLQ0<TFraction>::ConstIterator
263 DGtal::StandardDSLQ0<TFraction>::
264 begin( Point p ) const
266 return ConstIterator( *this, p );
268 //-----------------------------------------------------------------------------
269 template <typename TFraction>
270 typename DGtal::StandardDSLQ0<TFraction>::ConstIterator
271 DGtal::StandardDSLQ0<TFraction>::
274 ConstIterator it( *this, p );
278 //-----------------------------------------------------------------------------
279 template <typename TFraction>
280 typename DGtal::StandardDSLQ0<TFraction>::Self
281 DGtal::StandardDSLQ0<TFraction>::
282 reversedSmartDSS( 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 //-----------------------------------------------------------------------------
295 template <typename TFraction>
296 typename DGtal::StandardDSLQ0<TFraction>::Self
297 DGtal::StandardDSLQ0<TFraction>::
298 reversedSmartDSS( 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 //-----------------------------------------------------------------------------
383 template <typename TFraction>
384 typename DGtal::StandardDSLQ0<TFraction>::Self
385 DGtal::StandardDSLQ0<TFraction>::
386 DSSWithinTwoPatterns( 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 //-----------------------------------------------------------------------------
500 template <typename TFraction>
501 typename DGtal::StandardDSLQ0<TFraction>::Self
502 DGtal::StandardDSLQ0<TFraction>::
503 smartDSS( 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 //-----------------------------------------------------------------------------
591 template <typename TFraction>
593 typename DGtal::StandardDSLQ0<TFraction>::Fraction
594 DGtal::StandardDSLQ0<TFraction>::
595 deepest( Fraction f1, Fraction f2, Fraction f3 )
597 return deepest( f1, deepest( f2, f3 ) );
599 //-----------------------------------------------------------------------------
600 template <typename TFraction>
602 typename DGtal::StandardDSLQ0<TFraction>::Fraction
603 DGtal::StandardDSLQ0<TFraction>::
604 deepest( 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.
618 template <typename TFraction>
621 DGtal::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.
631 template <typename TFraction>
634 DGtal::StandardDSLQ0<TFraction>::isValid() const
641 ///////////////////////////////////////////////////////////////////////////////
642 // Implementation of inline functions //
644 template <typename TFraction>
647 DGtal::operator<< ( std::ostream & out,
648 const StandardDSLQ0<TFraction> & object )
650 object.selfDisplay( out );
655 ///////////////////////////////////////////////////////////////////////////////