DGtal  1.4.beta
FP.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 FP.ih
19  * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr )
20  * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21  *
22  * @date 2011/01/26
23  *
24  * @brief Implementation of inline methods defined in FP.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 //////////////////////////////////////////////////////////////////////////////
33 
34 ///////////////////////////////////////////////////////////////////////////////
35 // IMPLEMENTATION of inline methods.
36 ///////////////////////////////////////////////////////////////////////////////
37 
38 ///////////////////////////////////////////////////////////////////////////////
39 // ----------------------- Standard services ------------------------------
40 
41 template <typename TIterator, typename TInteger, int connectivity>
42 inline
43 DGtal::FP<TIterator,TInteger,connectivity>::
44 FP(const TIterator& itb,
45  const TIterator& ite )
46 {
47  algorithm(itb, ite);
48 }
49 
50 template <typename TIterator, typename TInteger, int connectivity>
51 inline
52 void
53 DGtal::FP<TIterator,TInteger,connectivity>::
54 algorithm(const TIterator& itb,
55  const TIterator& ite )
56 {
57  if (isNotEmpty(itb,ite))
58  {
59  typedef typename IteratorCirculatorTraits<TIterator>::Type Type;
60  algorithm( itb, ite, Type() );
61  }
62 }
63 
64 template <typename TIterator, typename TInteger, int connectivity>
65 inline
66 void
67 DGtal::FP<TIterator,TInteger,connectivity>::
68 algorithm(const TIterator& itb,
69  const TIterator& ite, IteratorType )
70 {
71  ASSERT(isNotEmpty(itb,ite));
72 
73  /////////////////////////////////////// open
74  //list of successive upper (U) and lower (L)
75  // leaning points.
76  std::list<Point> vTmpU, vTmpL;
77  vTmpU.push_back(*itb);
78  vTmpL.push_back(*itb);
79 
80  //longest DSS
81  DSSComputer longestDSS;
82  longestDSS.init(itb);
83 
84  while ( (longestDSS.end() != ite)&&(longestDSS.extendFront()) )
85  {
86  //store the first upper leaning point
87  if (longestDSS.Uf() != vTmpU.back()) {
88  vTmpU.push_back(longestDSS.Uf());
89  }
90  //store the first lower leaning point
91  if (longestDSS.Lf() != vTmpL.back()) {
92  vTmpL.push_back(longestDSS.Lf());
93  }
94  }
95 
96  //std::cerr << "First MS" << std::endl << longestDSS << std::endl;
97 
98  typedef detail::DSSDecorator<DSSComputer> DSSDecorator;
99  DSSDecorator* adapter = 0; //adapter for the current DSS
100 
101  if (longestDSS.end() != ite)
102  {
103  //if the curve is not straight
104  //it is either locally convex or concave at longestDSS.end()
105  adapter = initConvexityConcavity<DSSDecorator>( longestDSS );
106 
107  ASSERT(adapter);
108  if ( adapter->isInConvexPart() ) myPolygon = vTmpU;
109  else myPolygon = vTmpL;
110  ASSERT(myPolygon.size() > 0);
111 
112  bool go = true;
113  while ( go )
114  { //main loop
115  if ( !removingStep( adapter ) )
116  { //disconnected digital curve
117  throw InputException();
118  }
119  go = addingStep( adapter, ite );
120  if ( go )
121  {
122  delete (adapter);
123  adapter = initConvexityConcavity<DSSDecorator>( longestDSS );
124  }
125  }
126 
127  }
128  else
129  {
130  //the curve is assumed to be convex
131  //if it is straight
132  myPolygon = vTmpU;
133  adapter = new detail::DSSDecorator4ConvexPart<DSSComputer>(longestDSS);
134  ASSERT(adapter);
135  }
136 
137  //store the last leaning point
138  //if the first and last leaning points
139  //of the MS are not confounded
140  if (adapter->firstLeaningPoint() != adapter->lastLeaningPoint()) {
141  myPolygon.push_back(adapter->lastLeaningPoint());
142  }
143 
144  //last removing step
145  while ( longestDSS.retractBack() ) {
146  //store the last leaning point
147  if (adapter->lastLeaningPoint() != myPolygon.back()) {
148  myPolygon.push_back(adapter->lastLeaningPoint());
149  }
150  }
151 
152  ASSERT(adapter);//had to be allocated
153  delete (adapter);
154 }
155 
156 template <typename TIterator, typename TInteger, int connectivity>
157 inline
158 void
159 DGtal::FP<TIterator,TInteger,connectivity>::
160 algorithm(const TIterator& itb,
161  const TIterator& ite, CirculatorType )
162 {
163  ASSERT(isNotEmpty(itb,ite)); boost::ignore_unused_variable_warning(ite);
164 
165  /////////////////////////// closed
166  //first maximal DSS
167  DSSComputer currentDSS;
168  currentDSS.init( itb );
169  //forward extension
170  while ( currentDSS.extendFront() ) {}
171  //backward extension
172  while ( currentDSS.extendBack() ) {}
173 
174  //local convexity
175  typedef detail::DSSDecorator<DSSComputer> DSSDecorator;
176  DSSDecorator* adapter = 0; //adapter for the current DSS
177  adapter = initConvexityConcavity<DSSDecorator>( currentDSS );
178 
179  if (adapter->firstLeaningPoint() == adapter->lastLeaningPoint()) {
180  myPolygon.push_back( adapter->lastLeaningPoint() );
181  }//stored in removingStep function otherwise
182 
183  //first MS
184  DSSComputer firstMS = currentDSS;
185 
186  //since there are more than 2 MS (closed case)
187  //but no more than 2 MS sharing a leaning point
188  //at the same position
189  if ( !removingStep( adapter ) )
190  { //disconnected digital curve
191  throw InputException();
192  }
193  addingStep( adapter );
194  bool go = true;
195  while ( go )
196  { //main loop
197  delete (adapter);
198  adapter = initConvexityConcavity<DSSDecorator>( currentDSS );
199 
200  if ( !removingStep( adapter ) )
201  { //disconnected digital curve
202  throw InputException();
203  }
204  addingStep( adapter );
205  if (currentDSS == firstMS)
206  { //stop and remove the last point to avoid a repetition
207  if (adapter->firstLeaningPoint() == myPolygon.front()) {
208  ASSERT( myPolygon.size() > 0 );
209  myPolygon.pop_back();
210  }
211  go = false;
212  }
213  }
214 
215  delete (adapter);
216 }
217 
218 template <typename TIterator, typename TInteger, int connectivity>
219 template<typename Adapter>
220 inline
221 Adapter*
222 DGtal::FP<TIterator,TInteger,connectivity>
223 ::initConvexityConcavity( typename Adapter::DSS &aDSS )
224 {
225 
226  ASSERT( aDSS.isExtendableFront() == false );
227 
228  if ( aDSS.remainder( aDSS.end() ) < (aDSS.mu()) )
229  { //concave part
230  return new detail
231  ::DSSDecorator4ConcavePart<typename Adapter::DSS>(aDSS);
232  }
233  else
234  { //convex part
235  return new detail
236  ::DSSDecorator4ConvexPart<typename Adapter::DSS>(aDSS);
237  }
238 }
239 
240 template <typename TIterator, typename TInteger, int connectivity>
241 template<typename Adapter>
242 inline
243 bool
244 DGtal::FP<TIterator,TInteger,connectivity>
245 ::removingStep( Adapter* adapter )
246 {
247 
248  ASSERT( adapter->isExtendableFront() == false );
249 
250  //store the last leaning point
251  //if the first and last leaning points
252  //of the MS are not confounded
253  if (adapter->firstLeaningPoint() != adapter->lastLeaningPoint()) {
254  myPolygon.push_back(adapter->lastLeaningPoint());
255  }
256 
257  //removing step
258  while ( !adapter->isExtendableFront() ) {
259  //remove a point from the back
260  if ( adapter->retractBack() ) {
261  //store the last leaning point
262  ASSERT( myPolygon.size() > 0 );
263  if (adapter->lastLeaningPoint() != myPolygon.back()) {
264  myPolygon.push_back(adapter->lastLeaningPoint());
265  }
266  } else {
267  return false;
268  }
269  }
270 
271  return true;
272 }
273 
274 template <typename TIterator, typename TInteger, int connectivity>
275 template<typename Adapter>
276 inline
277 bool
278 DGtal::FP<TIterator,TInteger,connectivity>
279 ::addingStep( Adapter* adapter,
280  const typename Adapter::DSS::ConstIterator& itEnd )
281 {
282 
283  ASSERT( adapter->isExtendableFront() == true );
284 
285  //remove the last leaning point
286  //if the first and last leaning points
287  //of the current DSS are not confounded
288  if (adapter->firstLeaningPoint() != adapter->lastLeaningPoint()) {
289  ASSERT( myPolygon.size() > 0 );
290  myPolygon.pop_back();
291  }
292 
293  //adding step
294  while ( (adapter->end() != itEnd)&&(adapter->extendFront()) ) {
295  //store the first leaning point
296  ASSERT( myPolygon.size() > 0 );
297  if (adapter->firstLeaningPoint() != myPolygon.back()) {
298  myPolygon.push_back(adapter->firstLeaningPoint());
299  }
300  }
301 
302  return (adapter->end() != itEnd);
303 }
304 
305 template <typename TIterator, typename TInteger, int connectivity>
306 template<typename Adapter>
307 inline
308 void
309 DGtal::FP<TIterator,TInteger,connectivity>
310 ::addingStep( Adapter* adapter )
311 {
312 
313  ASSERT( adapter->isExtendableFront() == true );
314 
315  //remove the last leaning point
316  //if the first and last leaning points
317  //of the current DSS are not confounded
318  if (adapter->firstLeaningPoint() != adapter->lastLeaningPoint()) {
319  ASSERT( myPolygon.size() > 0 );
320  myPolygon.pop_back();
321  }
322 
323  //adding step
324  while ( adapter->extendFront() ) {
325  //store the first leaning point
326  ASSERT( myPolygon.size() > 0 );
327  if (adapter->firstLeaningPoint() != myPolygon.back()) {
328  myPolygon.push_back(adapter->firstLeaningPoint());
329  }
330  }
331 
332 }
333 
334 template <typename TIterator, typename TInteger, int connectivity>
335 inline
336 DGtal::FP<TIterator,TInteger,connectivity>::~FP()
337 {
338 }
339 
340 ///////////////////////////////////////////////////////////////////////////////
341 // Interface - public :
342 
343 template <typename TIterator, typename TInteger, int connectivity>
344 inline
345 const typename DGtal::FP<TIterator,TInteger,connectivity>::Polygon&
346 DGtal::FP<TIterator,TInteger,connectivity>::polygon() const
347 {
348  return myPolygon;
349 }
350 
351 template <typename TIterator, typename TInteger, int connectivity>
352 inline
353 bool
354 DGtal::FP<TIterator,TInteger,connectivity>::isClosed() const
355 {
356  //since the input digital curve has to be connected
357  return IsCirculator<TIterator>::value;
358 }
359 
360 template <typename TIterator, typename TInteger, int connectivity>
361 inline
362 bool
363 DGtal::FP<TIterator,TInteger,connectivity>
364 ::isValid(const Point& a,const Point& b, const Point& c) const {
365 
366  Vector e1 = b - a; //previous edge
367  Vector e2 = c - b; //next edge
368 
369  if ( (quadrant(e1,1))&&(quadrant(e2,1)) ) {
370  return true;
371  } else if ( (quadrant(e1,2))&&(quadrant(e2,2)) ) {
372  return true;
373  } else if ( (quadrant(e1,3))&&(quadrant(e2,3)) ) {
374  return true;
375  } else if ( (quadrant(e1,4))&&(quadrant(e2,4)) ) {
376  return true;
377  } else {
378  return false;
379  }
380 
381 }
382 
383 template <typename TIterator, typename TInteger, int connectivity>
384 inline
385 bool
386 DGtal::FP<TIterator,TInteger,connectivity>::isValid() const
387 {
388  typedef typename Polygon::const_iterator I;
389 
390  //not two consecutive confounded points
391  bool flag1 = true;
392  I it = myPolygon.begin();
393  if (it != myPolygon.end())
394  {
395  I itPrev = it; ++it;
396  if (it != myPolygon.end())
397  {
398  while ( (it != myPolygon.end())&&(flag1) )
399  {
400  flag1 = (*itPrev != *it);
401  itPrev = it;
402  ++it;
403  }
404  }//if only one point ok
405  }//if no point ok
406 
407  //valid quadrant for all consecutive pairs of edges
408  bool flag2 = true;
409  if (myPolygon.size() >= 3) {
410 
411  I i, j, k;
412  i = myPolygon.begin();
413  j = i; ++j;
414  k = j; ++k;
415 
416  if ( isClosed() ) { ///////// closed
417 
418  //first point
419  flag2 = isValid(myPolygon.back(),*i,*j);
420  //middle points
421  while ( (k != myPolygon.end())&&(flag2) ) {
422  flag2 = isValid(*i,*j,*k);
423  ++i; ++j; ++k;
424  }
425  //last point
426  if (flag2)
427  flag2 = isValid(*i,*j,myPolygon.front());
428 
429  } else { ////////////////////// open
430 
431  //middle points
432  while ( (k != myPolygon.end())&&(flag2) ) {
433  flag2 = isValid(*i,*j,*k);
434  ++i; ++j; ++k;
435  }
436 
437  }
438  }//special case of less than 3 points ok
439 
440  return flag1 && flag2;
441 }
442 
443 template <typename TIterator, typename TInteger, int connectivity>
444 inline
445 typename DGtal::FP<TIterator,TInteger,connectivity>::Polygon::size_type
446 DGtal::FP<TIterator,TInteger,connectivity>::size() const
447 {
448  return myPolygon.size();
449 }
450 
451 
452 template <typename TIterator, typename TInteger, int connectivity>
453 template <typename OutputIterator>
454 inline
455 OutputIterator
456 DGtal::FP<TIterator,TInteger,connectivity>::copyFP(OutputIterator result) const {
457  typename Polygon::const_iterator i = myPolygon.begin();
458  while ( i != myPolygon.end() ) {
459  *result++ = *i++;
460  }
461  return result;
462 }
463 
464 template <typename TIterator, typename TInteger, int connectivity>
465 template <typename OutputIterator>
466 inline
467 OutputIterator
468 DGtal::FP<TIterator,TInteger,connectivity>::copyMLP(OutputIterator result) const {
469 
470  ASSERT( isValid() );
471 
472  unsigned int n = (unsigned int) myPolygon.size();
473  if (n < 3) { //special case < 3 points
474 
475  typename Polygon::const_iterator i = myPolygon.begin();
476  for ( ; i!= myPolygon.end() ; ++i) {
477  *result++ = RealPoint( *i );
478  }
479 
480  } else { //standard case
481 
482  typename Polygon::const_iterator i, j, k;
483  i = myPolygon.begin();
484  j = i; ++j;
485  k = j; ++k;
486 
487  if ( isClosed() ) { ///////// closed
488 
489  //first point
490  *result++ = getRealPoint(myPolygon.back(),*i,*j);
491  //middle points
492  while ( k != myPolygon.end() ) {
493  *result++ = getRealPoint(*i,*j,*k);
494  ++i; ++j; ++k;
495  }
496  //last point
497  *result++ = getRealPoint(*i,*j,myPolygon.front());
498 
499  } else { ////////////////////// open
500 
501  //first point
502  *result++ = RealPoint( myPolygon.front() );
503  //middle points
504  while ( k != myPolygon.end() ) {
505  *result++ = getRealPoint(*i,*j,*k);
506  ++i; ++j; ++k;
507  }
508  //last point
509  *result++ = RealPoint( myPolygon.back() );
510 
511  }
512 
513  }
514  return result;
515 }
516 
517 template <typename TIterator, typename TInteger, int connectivity>
518 inline
519 DGtal::PointVector<2,double>
520 DGtal::FP<TIterator,TInteger,connectivity>
521 ::getRealPoint (const Point& a,const Point& b, const Point& c) const {
522 
523  ASSERT( isValid() );
524 
525  RealVector shift;
526 
527  Vector e1 = b - a; //previous edge
528  Vector e2 = c - b; //next edge
529 
530  if ( (e1[0]*e2[1]-e1[1]*e2[0]) <= 0 ) {
531 
532  //convex turn
533  if ( (quadrant(e1,1))&&(quadrant(e2,1)) ) {
534  shift = RealVector(0.5,-0.5);
535  } else if ( (quadrant(e1,2))&&(quadrant(e2,2)) ) {
536  shift = RealVector(-0.5,-0.5);
537  } else if ( (quadrant(e1,3))&&(quadrant(e2,3)) ) {
538  shift = RealVector(-0.5,0.5);
539  } else if ( (quadrant(e1,4))&&(quadrant(e2,4)) ) {
540  shift = RealVector(0.5,0.5);
541  } else {
542  ASSERT(false && "DGtal::FP<TIterator,TInteger,connectivity>::getRealPoint: not valid polygon" );
543  }
544 
545  } else {
546 
547  //concave turn
548  if ( (quadrant(e1,1))&&(quadrant(e2,1)) ) {
549  shift = RealVector(-0.5,0.5);
550  } else if ( (quadrant(e1,2))&&(quadrant(e2,2)) ) {
551  shift = RealVector(0.5,0.5);
552  } else if ( (quadrant(e1,3))&&(quadrant(e2,3)) ) {
553  shift = RealVector(0.5,-0.5);
554  } else if ( (quadrant(e1,4))&&(quadrant(e2,4)) ) {
555  shift = RealVector(-0.5,-0.5);
556  } else {
557  ASSERT(false && "DGtal::FP<TIterator,TInteger,connectivity>::getRealPoint: not valid polygon" );
558  }
559 
560  }
561 
562  return ( RealPoint(b) + shift );
563 }
564 
565 template <typename TIterator, typename TInteger, int connectivity>
566 inline
567 bool
568 DGtal::FP<TIterator,TInteger,connectivity>
569 ::quadrant (const Vector& v, const int& q) const {
570 
571  if (q == 1) {
572  return ( (v[0]>=0)&&(v[1]>=0) );
573  } else if (q == 2) {
574  return ( (v[0]>=0)&&(v[1]<=0) );
575  } else if (q == 3) {
576  return ( (v[0]<=0)&&(v[1]<=0) );
577  } else if (q == 4) {
578  return ( (v[0]<=0)&&(v[1]>=0) );
579  } else {
580  ASSERT(false &&
581  "DGtal::FP<TIterator,TInteger,connectivity>::quadrant: quadrant number should be 0,1,2 or 3" );
582  return false;
583  }
584 }
585 
586 ///////////////////////////////////////////////////////////////////////////////
587 // Display :
588 
589 template <typename TIterator, typename TInteger, int connectivity>
590 inline
591 std::string
592 DGtal::FP<TIterator,TInteger,connectivity>::className() const
593 {
594  return "FP";
595 }
596 
597 
598 template <typename TIterator, typename TInteger, int connectivity>
599 inline
600 void
601 DGtal::FP<TIterator,TInteger,connectivity>::selfDisplay ( std::ostream & out ) const
602 {
603  out << "[FP]" << std::endl;
604  typename Polygon::const_iterator i = myPolygon.begin();
605  for ( ;i != myPolygon.end();++i)
606  {
607  out << "\t " << (*i) << std::endl;
608  }
609  out << "[end FP]" << std::endl;
610 }
611 
612 
613 
614 
615 
616 ///////////////////////////////////////////////////////////////////////////////
617 // Implementation of inline functions //
618 
619 template <typename TIterator, typename TInteger, int connectivity>
620 inline
621 std::ostream&
622 DGtal::operator<< ( std::ostream & out,
623  const FP<TIterator,TInteger,connectivity> & object )
624 {
625  object.selfDisplay( out );
626  return out;
627 }
628 
629 // //
630 ///////////////////////////////////////////////////////////////////////////////
631 
632