DGtal  1.5.beta
LatticePolytope2D.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 LatticePolytope2D.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21  * @author Emilie Charrier
22  *
23  * @date 2012/04/19
24  *
25  * Implementation of inline methods defined in LatticePolytope2D.h
26  *
27  * This file is part of the DGtal library.
28  */
29 
30 
31 //////////////////////////////////////////////////////////////////////////////
32 #include <cstdlib>
33 #include <iterator>
34 #include "DGtal/kernel/NumberTraits.h"
35 #include "DGtal/kernel/sets/CDigitalSet.h"
36 //////////////////////////////////////////////////////////////////////////////
37 
38 ///////////////////////////////////////////////////////////////////////////////
39 // IMPLEMENTATION of inline methods.
40 ///////////////////////////////////////////////////////////////////////////////
41 
42 ///////////////////////////////////////////////////////////////////////////////
43 // ----------------------- Standard services ------------------------------
44 
45 //-----------------------------------------------------------------------------
46 template <typename TSpace, typename TSequence>
47 inline
48 DGtal::LatticePolytope2D<TSpace,TSequence>::~LatticePolytope2D()
49 { // Nothing to do.
50 }
51 //-----------------------------------------------------------------------------
52 template <typename TSpace, typename TSequence>
53 inline
54 DGtal::LatticePolytope2D<TSpace,TSequence>::LatticePolytope2D()
55 { // Nothing to do.
56 }
57 //-----------------------------------------------------------------------------
58 template <typename TSpace, typename TSequence>
59 inline
60 DGtal::LatticePolytope2D<TSpace,TSequence>::LatticePolytope2D
61 ( const Self & other )
62  : myVertices( other.myVertices )
63 { // Nothing to do.
64 }
65 //-----------------------------------------------------------------------------
66 template <typename TSpace, typename TSequence>
67 inline
68 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Self &
69 DGtal::LatticePolytope2D<TSpace,TSequence>::operator=
70 ( const Self & other )
71 { // Nothing to do.
72  if ( this != &other )
73  myVertices = other.myVertices;
74  return *this;
75 }
76 //-----------------------------------------------------------------------------
77 template <typename TSpace, typename TSequence>
78 inline
79 typename DGtal::LatticePolytope2D<TSpace,TSequence>::ConstIterator
80 DGtal::LatticePolytope2D<TSpace,TSequence>::
81 begin() const
82 {
83  return myVertices.begin();
84 }
85 //-----------------------------------------------------------------------------
86 template <typename TSpace, typename TSequence>
87 inline
88 typename DGtal::LatticePolytope2D<TSpace,TSequence>::ConstIterator
89 DGtal::LatticePolytope2D<TSpace,TSequence>::
90 end() const
91 {
92  return myVertices.end();
93 }
94 //-----------------------------------------------------------------------------
95 template <typename TSpace, typename TSequence>
96 inline
97 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Iterator
98 DGtal::LatticePolytope2D<TSpace,TSequence>::
99 begin()
100 {
101  return myVertices.begin();
102 }
103 //-----------------------------------------------------------------------------
104 template <typename TSpace, typename TSequence>
105 inline
106 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Iterator
107 DGtal::LatticePolytope2D<TSpace,TSequence>::
108 end()
109 {
110  return myVertices.end();
111 }
112 //-----------------------------------------------------------------------------
113 template <typename TSpace, typename TSequence>
114 inline
115 bool
116 DGtal::LatticePolytope2D<TSpace,TSequence>::
117 empty() const
118 {
119  return myVertices.empty();
120 }
121 //-----------------------------------------------------------------------------
122 template <typename TSpace, typename TSequence>
123 inline
124 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Size
125 DGtal::LatticePolytope2D<TSpace,TSequence>::
126 size() const
127 {
128  return myVertices.size();
129 }
130 //-----------------------------------------------------------------------------
131 template <typename TSpace, typename TSequence>
132 inline
133 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Size
134 DGtal::LatticePolytope2D<TSpace,TSequence>::
135 max_size() const
136 {
137  return myVertices.max_size();
138 }
139 //-----------------------------------------------------------------------------
140 template <typename TSpace, typename TSequence>
141 inline
142 void
143 DGtal::LatticePolytope2D<TSpace,TSequence>::
144 clear()
145 {
146  myVertices.clear();
147 }
148 //-----------------------------------------------------------------------------
149 template <typename TSpace, typename TSequence>
150 inline
151 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Iterator
152 DGtal::LatticePolytope2D<TSpace,TSequence>::
153 erase( Iterator it )
154 {
155  return myVertices.erase( it );
156 }
157 //-----------------------------------------------------------------------------
158 template <typename TSpace, typename TSequence>
159 inline
160 void
161 DGtal::LatticePolytope2D<TSpace,TSequence>::
162 swap( LatticePolytope2D & other )
163 {
164  myVertices.swap( other.myVertices );
165 }
166 
167 //-----------------------------------------------------------------------------
168 template <typename TSpace, typename TSequence>
169 inline
170 typename DGtal::LatticePolytope2D<TSpace,TSequence>::MyIntegerComputer &
171 DGtal::LatticePolytope2D<TSpace,TSequence>::
172 ic() const
173 {
174  return _ic;
175 }
176 
177 //-----------------------------------------------------------------------------
178 template <typename TSpace, typename TSequence>
179 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Domain
180 DGtal::LatticePolytope2D<TSpace,TSequence>::
181 boundingBoxDomain() const
182 {
183  ConstIterator it = begin();
184  ConstIterator it_end = end();
185  Point infimum = *it;
186  Point supremum = *it;
187  for ( ++it; it != it_end; ++it )
188  {
189  infimum = infimum.inf( *it );
190  supremum = supremum.sup( *it );
191  }
192  return Domain( infimum, supremum );
193 }
194 //-----------------------------------------------------------------------------
195 template <typename TSpace, typename TSequence>
196 void
197 DGtal::LatticePolytope2D<TSpace,TSequence>::
198 purge()
199 {
200  Iterator it = begin();
201  Iterator it_end = end();
202  Iterator it_prev;
203  while ( it != it_end )
204  {
205  _A = *it;
206  it_prev = it;
207  ++it;
208  while ( ( it != it_end ) && ( _A == *it ) )
209  it = erase( it );
210  }
211  // Checks case where first vertex is also last vertex.
212  if ( ( size() > 1 ) && ( * begin() == *it_prev ) )
213  erase( begin() );
214 }
215 //-----------------------------------------------------------------------------
216 template <typename TSpace, typename TSequence>
217 inline
218 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Iterator
219 DGtal::LatticePolytope2D<TSpace,TSequence>::
220 insertBefore( const Iterator & pos, const Point & K )
221 {
222  return myVertices.insert( pos, K );
223 }
224 //-----------------------------------------------------------------------------
225 template <typename TSpace, typename TSequence>
226 inline
227 void
228 DGtal::LatticePolytope2D<TSpace,TSequence>::
229 pushBack( const Point & K )
230 {
231  myVertices.push_back( K );
232 }
233 //-----------------------------------------------------------------------------
234 template <typename TSpace, typename TSequence>
235 inline
236 void
237 DGtal::LatticePolytope2D<TSpace,TSequence>::
238 pushFront( const Point & K )
239 {
240  myVertices.push_front( K );
241 }
242 //-----------------------------------------------------------------------------
243 template <typename TSpace, typename TSequence>
244 inline
245 void
246 DGtal::LatticePolytope2D<TSpace,TSequence>::
247 push_back( const Point & K )
248 {
249  myVertices.push_back( K );
250 }
251 //-----------------------------------------------------------------------------
252 template <typename TSpace, typename TSequence>
253 inline
254 void
255 DGtal::LatticePolytope2D<TSpace,TSequence>::
256 push_front( const Point & K )
257 {
258  myVertices.push_front( K );
259 }
260 //-----------------------------------------------------------------------------
261 template <typename TSpace, typename TSequence>
262 inline
263 const typename DGtal::LatticePolytope2D<TSpace,TSequence>::Integer &
264 DGtal::LatticePolytope2D<TSpace,TSequence>::
265 twiceArea() const
266 {
267  _a = NumberTraits<Integer>::ZERO;
268  ConstIterator it = begin();
269  ConstIterator it_end = end();
270  ConstIterator it_next = it; ++it_next;
271  while ( it_next != it_end )
272  {
273  _a += _ic.crossProduct( *it, *it_next );
274  it = it_next; ++it_next;
275  }
276  _a += _ic.crossProduct( *it, *(begin()) );
277  return _a;
278 }
279 //-----------------------------------------------------------------------------
280 template <typename TSpace, typename TSequence>
281 inline
282 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Point3I
283 DGtal::LatticePolytope2D<TSpace,TSequence>::
284 centroid() const
285 {
286  _a = twiceArea();
287  return centroid( _a );
288 }
289 //-----------------------------------------------------------------------------
290 template <typename TSpace, typename TSequence>
291 inline
292 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Point3I
293 DGtal::LatticePolytope2D<TSpace,TSequence>::
294 centroid( const Integer & twice_area ) const
295 {
296  _a = _b = NumberTraits<Integer>::ZERO;
297  _A = Point( NumberTraits<Integer>::ZERO, NumberTraits<Integer>::ZERO );
298  ConstIterator it_begin = begin();
299  ConstIterator it = it_begin;
300  ConstIterator it_end = end();
301  if( twice_area > NumberTraits<Integer>::ZERO )
302  {
303  _den = 3 * twice_area;
304  ConstIterator it_next = it; ++it_next;
305  while ( it_next != it_end )
306  {
307  _ic.getCrossProduct( _c, *it, *it_next );
308  _B = *it + *it_next;
309  _B *= _c;
310  _A += _B;
311  it = it_next; ++it_next;
312  }
313  _ic.getCrossProduct( _c, *it, *it_begin );
314  _B = *it + *it_begin;
315  _B *= _c;
316  _A += _B;
317  }
318  else
319  {
320  _den = NumberTraits<Integer>::ZERO;
321  for ( ; it != it_end; ++it )
322  {
323  _A += *it;
324  ++_den;
325  }
326  }
327  return Point3I( _A[ 0 ], _A[ 1 ], _den );
328 }
329 //-----------------------------------------------------------------------------
330 template <typename TSpace, typename TSequence>
331 inline
332 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Integer
333 DGtal::LatticePolytope2D<TSpace,TSequence>::
334 numberBoundaryPoints() const
335 {
336  _a = NumberTraits<Integer>::ZERO;
337  ConstIterator it = begin();
338  ConstIterator it_end = end();
339  ConstIterator it_next = it; ++it_next;
340  while ( it_next != it_end )
341  {
342  _N = *it_next - *it;
343  ic().getGcd( _g, _N[ 0 ], _N[ 1 ] );
344  _a += _g;
345  it = it_next; ++it_next;
346  }
347  _N = *it - *(begin());
348  ic().getGcd( _g, _N[ 0 ], _N[ 1 ] );
349  _a += _g;
350  return _a;
351 }
352 //-----------------------------------------------------------------------------
353 template <typename TSpace, typename TSequence>
354 inline
355 typename DGtal::LatticePolytope2D<TSpace,TSequence>::Integer
356 DGtal::LatticePolytope2D<TSpace,TSequence>::
357 numberInteriorPoints() const
358 {
359  _c1 = numberBoundaryPoints();
360  _c3 = twiceArea();
361  _c = _c3 - _c1 + 2;
362  ASSERT( ic().isZero( _c % 2 ) );
363  return _c / 2;
364 }
365 //-----------------------------------------------------------------------------
366 template <typename TSpace, typename TSequence>
367 inline
368 typename DGtal::LatticePolytope2D<TSpace,TSequence>::SizeCouple
369 DGtal::LatticePolytope2D<TSpace,TSequence>::
370 findCut( Iterator & it_next_is_outside, Iterator & it_next_is_inside,
371  const HalfSpace & hs )
372 {
373  Size nbWithin = (Size) 0; // JOL: problem on macos. unsigned long has no NumberTraits. NumberTraits<Size>::ZERO;
374  Size nb = (Size) 0; // JOL: problem on macos. unsigned long has no NumberTraits. NumberTraits<Size>::ZERO;
375  Iterator it = begin();
376  Iterator it_prev = it;
377  Iterator it_end = end();
378  it_next_is_outside = it_next_is_inside = it_end;
379  if ( it == it_end )
380  return std::make_pair( (Size) 0, (Size) 0 ); // JOL: problem on macos. unsigned long has no NumberTraits. NumberTraits<Size>::ZERO;
381  bool visibility_begin_vtx; // visibility of begin vertex.
382  bool visibility_prev_vtx; // visibility of previous vertex when visiting the list.
383  bool visibility_vtx; // visibility of current vertex when visiting the list.
384  ++nb;
385  if ( ( visibility_begin_vtx = hs( *it++ ) ) ) ++nbWithin; // Assignment
386  visibility_prev_vtx = visibility_begin_vtx;
387  visibility_vtx = visibility_begin_vtx;
388  for ( ; it != it_end; ++it, ++nb )
389  {
390  if ( ( visibility_vtx = hs( *it ) ) ) ++nbWithin; // Assignment
391  if ( visibility_vtx != visibility_prev_vtx )
392  {
393  if ( visibility_prev_vtx ) it_next_is_outside = it_prev;
394  else it_next_is_inside = it_prev;
395  }
396  visibility_prev_vtx = visibility_vtx;
397  it_prev = it;
398  }
399  if ( visibility_vtx != visibility_begin_vtx )
400  {
401  if ( visibility_vtx ) it_next_is_outside = it_prev;
402  else it_next_is_inside = it_prev;
403  }
404  ASSERT( ( nbWithin == 0 ) || ( nbWithin == size() )
405  || ( ( it_next_is_outside != it_end ) && ( it_next_is_inside != it_end )
406  && ( it_next_is_inside != it_next_is_outside ) ) );
407  return std::make_pair( nbWithin, nb );
408 }
409 
410 //-----------------------------------------------------------------------------
411 template <typename TSpace, typename TSequence>
412 inline
413 bool
414 DGtal::LatticePolytope2D<TSpace,TSequence>::
415 cut( const HalfSpace & hs )
416 {
417  Iterator it_next_is_outside;
418  Iterator it_next_is_inside;
419  SizeCouple nbs = findCut( it_next_is_outside, it_next_is_inside, hs );
420 
421  // Take care of easy cases.
422  if ( nbs.first == nbs.second ) { return false; }
423  if ( nbs.first == (Size) 0 ) { clear(); return true; } // JOL see findCut
424 
425  // Otherwise, determines A1B1 and A2B2.
426  twiceArea(); // result in _a;
427  HalfSpace hs1 = halfSpace( it_next_is_outside );
428  HalfSpace hs3 = halfSpace( it_next_is_inside );
429  //hs3.negate();
430  _A1 = *it_next_is_outside;
431  ++it_next_is_outside;
432  if ( it_next_is_outside == end() ) it_next_is_outside = begin();
433  _B1 = *it_next_is_outside;
434  _B2 = *it_next_is_inside;
435  ++it_next_is_inside;
436  if ( it_next_is_inside == end() ) it_next_is_inside = begin();
437  _A2 = *it_next_is_inside;
438 
439  // Erases outside vertices.
440  while ( it_next_is_outside != it_next_is_inside )
441  {
442  it_next_is_outside = erase( it_next_is_outside );
443  if ( it_next_is_outside == end() ) it_next_is_outside = begin();
444  }
445  // Both iterators point on the right place.
446  if ( _a > NumberTraits<Integer>::ZERO )
447  { //convex not reduced to a straight line segment
448  std::insert_iterator<ClockwiseVertexSequence> itOut
449  = std::inserter<ClockwiseVertexSequence>( myVertices, it_next_is_outside );
450  computeConvexHullBorder( itOut, _A1, _A2, hs1, hs, hs3 );
451  }
452  else //convex reduced to a straight line segment
453  {
454  //compute the new extremity of the straight line segment
455  _v = _B1 - _A1;
456  _ic.reduce( _v );
457  _a = ( hs.c - hs.N.dot( _A1 ) ) / ( hs.N.dot( _v ) );
458  _A1 += _v * _a;
459  insertBefore( it_next_is_outside, _A1 );
460  }
461  purge(); // O(n)
462  return true;
463 }
464 //-----------------------------------------------------------------------------
465 template <typename TSpace, typename TSequence>
466 inline
467 typename DGtal::LatticePolytope2D<TSpace,TSequence>::HalfSpace
468 DGtal::LatticePolytope2D<TSpace,TSequence>::
469 halfSpace( ConstIterator it ) const
470 {
471  Point A( *it ); ++it;
472  if ( it == end() ) it = begin();
473  Point B( *it ); ++it;
474  if ( it == end() ) it = begin();
475  _N[ 0 ] = A[ 1 ] - B[ 1 ];
476  _N[ 1 ] = B[ 0 ] - A[ 0 ];
477  _ic.getDotProduct( _c, _N, A );
478  _ic.getDotProduct( _c1, _N, *it );
479  if ( _c1 > _c )
480  {
481  _N.negate();
482  _c = -_c;
483  }
484  //simplification of the constraint
485  _ic.getGcd( _g, _N[ 0 ], _N[ 1 ] );
486  _N /= _g;
487  return HalfSpace( _N, _ic.floorDiv( _c, _g) );
488 }
489 //-----------------------------------------------------------------------------
490 template <typename TSpace, typename TSequence>
491 inline
492 typename DGtal::LatticePolytope2D<TSpace,TSequence>::HalfSpace
493 DGtal::LatticePolytope2D<TSpace,TSequence>::
494 halfSpace( const Point & A, const Point & B, const Point & inP ) const
495 {
496  _N[ 0 ] = A[ 1 ] - B[ 1 ];
497  _N[ 1 ] = B[ 0 ] - A[ 0 ];
498  _ic.getDotProduct( _c, _N, A );
499  _ic.getDotProduct( _c1, _N, inP );
500  if ( _c1 > _c )
501  {
502  _N.negate();
503  _c = -_c;
504  }
505  //simplification of the constraint
506  _ic.getGcd( _g, _N[ 0 ], _N[ 1 ] );
507  _N /= _g;
508  return HalfSpace( _N, _ic.floorDiv( _c, _g) );
509 }
510 
511  /**
512  Computes the set \a aSet all the digital points that belongs to this polygon.
513 
514  @param aSet (returns) the set that contains as output all the
515  digital points of this polygon.
516 
517  @todo this method is for now not efficient.
518  */
519 //-----------------------------------------------------------------------------
520 template <typename TSpace, typename TSequence>
521 template <typename DigitalSet>
522 void
523 DGtal::LatticePolytope2D<TSpace,TSequence>::
524 getIncludedDigitalPoints( DigitalSet & aSet ) const
525 {
526  BOOST_CONCEPT_ASSERT(( concepts::CDigitalSet< DigitalSet > ));
527  typedef typename DigitalSet::Domain DigitalSetDomain;
528  BOOST_STATIC_ASSERT(( concepts::ConceptUtils::SameType< DigitalSetDomain, Domain >::value ));
529  typedef typename DigitalSetDomain::ConstIterator DomainConstIterator;
530  typedef typename DigitalSet::Iterator DigitalSetIterator;
531  aSet.clear();
532  // case 1 : empty polygon.
533  if ( empty() ) return;
534  // case 2 : one vertex.
535  size_t s = size();
536  if ( s == 1 )
537  {
538  aSet.insert( *begin() );
539  return;
540  }
541  // case 3 : 2 vertices
542  if ( s == 2 )
543  {
544  ConstIterator it_vtx = begin();
545  ConstIterator it_vtx2 = it_vtx; ++it_vtx2;
546  HalfSpace hs = halfSpace( it_vtx );
547  Vector orthN( -hs.N[ 1 ], hs.N[ 0 ] );
548  HalfSpace hs_orth( orthN, _ic.dotProduct( orthN, *it_vtx ) );
549  if ( ! hs_orth( *it_vtx2 ) ) hs_orth.negate();
550  HalfSpace hs2 = halfSpace( it_vtx2 );
551  Vector orthN2( -hs2.N[ 1 ], hs2.N[ 0 ] );
552  HalfSpace hs2_orth( orthN2, _ic.dotProduct( orthN2, *it_vtx2 ) );
553  if ( ! hs2_orth( *it_vtx ) ) hs2_orth.negate();
554 
555  for ( DomainConstIterator it = aSet.domain().begin(),
556  it_end = aSet.domain().end(); it != it_end; ++it )
557  if ( hs( *it ) ) aSet.insert( *it );
558  for ( DigitalSetIterator it_set = aSet.begin(), it_set_end = aSet.end();
559  it_set != it_set_end; )
560  {
561  DigitalSetIterator it_cur = it_set; ++it_set;
562  if ( ! ( hs_orth( *it_cur )
563  && hs2( *it_cur )
564  && hs2_orth( *it_cur ) ) )
565  aSet.erase( it_cur );
566  }
567  return;
568  }
569  // case 4 : >= 3 vertices
570  ConstIterator it_vtx = begin();
571  ConstIterator it_vtx_end = end();
572  HalfSpace hs = halfSpace( it_vtx );
573  for ( DomainConstIterator it = aSet.domain().begin(),
574  it_end = aSet.domain().end(); it != it_end; ++it )
575  if ( hs( *it ) ) aSet.insert( *it );
576  ++it_vtx;
577  for ( ; it_vtx != it_vtx_end; ++it_vtx )
578  {
579  hs = halfSpace( it_vtx );
580  for ( DigitalSetIterator it_set = aSet.begin(), it_set_end = aSet.end();
581  it_set != it_set_end; )
582  {
583  DigitalSetIterator it_cur = it_set; ++it_set;
584  if ( ! hs( *it_cur ) ) aSet.erase( it_cur );
585  }
586  }
587 }
588 //-----------------------------------------------------------------------------
589 template <typename TSpace, typename TSequence>
590 bool
591 DGtal::LatticePolytope2D<TSpace,TSequence>::
592 getFirstPointsOfHull
593 ( Vector & v,
594  Point & inPt, // must belong to hs1.
595  Point & outPt,
596  const HalfSpace & hs1,
597  const HalfSpace & hs2 ) const
598 {
599  ASSERT( hs1.isOnBoundary( inPt ) );
600  bool exactIntersection;
601 
602  // initialize vector directionVector (not definitive)
603  _DV = hs1.tangent();
604 
605  //compute the intersection of ray (inPt, _DV) with constraint C2
606  _ic.getCoefficientIntersection( _fl, _ce,
607  inPt, _DV, hs2.N, hs2.c );
608 
609  // uses the intersection to compute the first vertex of the upper
610  // convex hull (inside convex hull), i.e. the grid point closest to
611  // C2 and satisfying C2
612  _ic.getDotProduct( _a, inPt, hs2.N );
613  // if ( ( _a > hs2.c ) && ( _fl == NumberTraits<Integer>::ZERO ) )
614  // {
615  // inPt += _DV * _ce;
616  // _DV.neg();
617  // }
618  // else if ( ( ( _a <= hs2.c ) && ( _fl >= NumberTraits<Integer>::ZERO ) )
619  // || ( ( _a > hs2.c2 ) && ( _fl <= NumberTraits<Integer>::ZERO ) ) )
620  if ( ( ( _a <= hs2.c ) && ( _fl >= NumberTraits<Integer>::ZERO ) )
621  || ( ( _a > hs2.c ) && ( _fl < NumberTraits<Integer>::ZERO ) ) )
622  {
623  inPt += _DV * _fl;
624  }
625  else
626  {
627  inPt += _DV * _ce;
628  _DV.negate();
629  }
630 
631  // compute the first vertex of the lower convex hull
632  if ( _fl == _ce ) //integer intersection
633  {
634  outPt = inPt;
635  exactIntersection = true;
636  }
637  else
638  {
639  outPt = inPt + _DV;
640  //initialization of v: valid Bezout vector of u
641  _ic.getValidBezout( v, inPt, _DV, hs1.N, hs1.c, hs2.N, hs2.c, true );
642  exactIntersection = false;
643  }
644 
645 #ifdef DEBUG_LatticePolytope2D
646  std::cerr << "[CIP::getFirstPointsOfHull] in=" << inPt
647  << " out=" << outPt << " v=" << v << std::endl;
648 #endif //DEBUG_LatticePolytope2D
649  ASSERT( hs1.isOnBoundary( inPt ) );
650  ASSERT( hs1.isOnBoundary( outPt ) );
651  ASSERT( hs2( inPt ) );
652  ASSERT( ( exactIntersection && hs2( outPt ) ) || ( ! exactIntersection && ! hs2( outPt ) ) );
653  return exactIntersection;
654 }
655 
656 /**
657  Computes the border of the upper and of the lower convex hull
658  from the starting points inPts[0] (up) and outPts[0]
659  down, along the constraint N2.p <= c2 while the vertices
660  satisfy the constraint N3.p <= c3. The vertices of the two
661  borders are stored at the end of inPts and outPts.
662 
663  @param inPts (in, out) as input, contains the first point, as
664  output the sequence of points satisfying \a hs2 and \a hs3.
665 
666  @param outPts (in, out) as input, contains the first point, as
667  output the sequence of points not satisfying \a hs2 and satisfying
668  \a hs3.
669 
670  @param BV the Bezout vector of the vector between inPts[ 0 ] and outPts[ 0 ].
671 
672  @param hs2 the half-space that is approached by the two sequences of points.
673  @param hs3 the limiting half-space which defines the bounds of the approximation.
674  */
675 //-----------------------------------------------------------------------------
676 template <typename TSpace, typename TSequence>
677 void
678 DGtal::LatticePolytope2D<TSpace,TSequence>::
679 getAllPointsOfHull( std::vector<Point> & inPts,
680  std::vector<Point> & outPts,
681  const Vector & BV,
682  const HalfSpace & hs2,
683  const HalfSpace & hs3 ) const
684 {
685  ASSERT( hs2.N != hs3.N );
686 
687  // _A and _B represents the last computed vertex above and below the constraint resp.
688  // _u, _v pair of Bezout vectors.
689  _v = BV;
690  // initialize A and B
691  _A = inPts[0];
692  _B = outPts[0];
693 
694  // while A and B do not lie on the supporting line of hs2
695  // and satisfies hs3 and while v is not parallel to hs2
696  while ( ( ! hs2.isOnBoundary( _A ) ) // stops if A reaches hs2.
697  && ( hs3( _A ) ) // stops if A does not satisfy hs3.
698  && ( ! hs2.isOnBoundary( _B ) ) // stops if B reaches hs2.
699  && ( hs3( _B ) ) // stops if B does not satisfy hs3.
700  && ( _ic.dotProduct( _v, hs2.N ) != NumberTraits<Integer>::ZERO ) )
701  {
702 #ifdef DEBUG_LatticePolytope2D
703  std::cerr << "[CIP::getAllPointsOfHull] A=" << _A
704  << " B=" << _B << " v= " << _v << std::endl;
705 #endif // DEBUG_LatticePolytope2D
706  if ( _ic.dotProduct( _v, hs2.N ) < NumberTraits<Integer>::ZERO )
707  { //------ second configuration, we find a new B --------------------
708  _ic.getCoefficientIntersection( _fl, _ce,
709  _B, _v, hs2.N, hs2.c );
710  _u = _B; _u += _v *_fl;
711  if ( hs3( _u ) )
712  { // Point is still within bounds.
713  _B = _u;
714  outPts.push_back( _B );
715 #ifdef DEBUG_LatticePolytope2D
716  std::cerr << "[CIP::getAllPointsOfHull] add B=" << _B << std::endl;
717 #endif // ifdef DEBUG_LatticePolytope2D
718  }
719  else
720  { // Point is too far away. Check intersection with hs3 instead.
721  _ic.getCoefficientIntersection( _fl, _ce,
722  _B, _v, hs3.N, hs3.c );
723  _B += _v * _fl;
724  outPts.push_back( _B );
725 #ifdef DEBUG_LatticePolytope2D
726  std::cerr << "[CIP::getAllPointsOfHull] add B=" << _B << std::endl;
727 #endif //ifdef DEBUG_LatticePolytope2D
728  ASSERT( hs3( _B ) );
729  break; // necessarily the last point.
730  }
731  }
732  else
733  { //----- first configuration, we find a new A -----------------------
734  _ic.getCoefficientIntersection( _fl, _ce,
735  _A, _v, hs2.N, hs2.c );
736  _u = _A; _u += _v *_fl;
737  if ( hs3( _u ) )
738  { // Point is still within bounds.
739  _A = _u;
740  inPts.push_back( _A );
741 #ifdef DEBUG_LatticePolytope2D
742  std::cerr << "[CIP::getAllPointsOfHull] add A=" << _A << std::endl;
743 #endif //ifdef DEBUG_LatticePolytope2D
744  }
745  else
746  { // Point is too far away. Check intersection with hs3 instead.
747  _ic.getCoefficientIntersection( _fl, _ce,
748  _A, _v, hs3.N, hs3.c );
749  _A += _v * _fl;
750  inPts.push_back( _A );
751 #ifdef DEBUG_LatticePolytope2D
752  std::cerr << "[CIP::getAllPointsOfHull] add A=" << _A << std::endl;
753 #endif //ifdef DEBUG_LatticePolytope2D
754  ASSERT( hs3( _A ) );
755  break; // necessarily the last point.
756  }
757  }
758  // update u and v
759  _u = _B;
760  _u -= _A;
761  // From _ic.getValidBezout( _A, _u, _v, N2, c2, N2, c2, 0);
762  // JOL: to check.
763  _ic.getValidBezout( _v, _A, _u, hs2.N, hs2.c, hs2.N, hs2.c, false );
764  }
765  // when the loop finishes, we have to complete the computation
766  // of each convex hull
767  // if ( ! hs3( _A ) ) // A does not satisfy C3, we remove it.
768  // inPts.pop_back();
769  // else if ( ! hs3( _B ) ) // B does not satisfy C3, we remove it
770  // outPts.pop_back();
771  // else
772  if ( hs2.isOnBoundary( _A ) ) // A lies on C2 so A is also a vertex of outPts
773  outPts.push_back( _A );
774  else if ( hs2.isOnBoundary( _B ) ) // B lies on C2 so B is also a vertex of inPts
775  inPts.push_back( _B );
776 }
777 
778 /**
779  * compute the convex hull of grid points satisfying the
780  * constraints N1.P<=c1, N2.P<=c2 and N3.P>=c3.
781  *
782  * N2.P<=c2 corresponds to the cut two parts of computation: from
783  * constraint 1 to constraint 3 and from constraint 3 to
784  * constraint 1.
785  *
786  * The computed vertices are inserted at position [pos] in some list.
787  *
788  * @param pointRefC1 and pointRefC3 corresponds to grid point lying on
789  * the supporting lines of C1 and of C3 resp.
790  *
791  * @param itOut corresponds to an iterator in the list of vertices
792  * of the convex, to add the next new vertices
793  *
794  * NB: the method also computes grid point satisfying N1.P<=c1 and
795  * N3.P>=c3 but not satisfying N2.P<=c2. They are stored in
796  * "resultdown" of size "nbverticesdown". the algorithm uses
797  * these points that's why they appear in the code.
798  */
799 //-----------------------------------------------------------------------------
800 template <typename TSpace, typename TSequence>
801 template <typename OutputIterator>
802 OutputIterator
803 DGtal::LatticePolytope2D<TSpace,TSequence>::
804 computeConvexHullBorder
805 ( OutputIterator itOut,
806  const Point & pointRefC1,
807  const Point & pointRefC3,
808  const HalfSpace & hs1,
809  const HalfSpace & hs2,
810  const HalfSpace & hs3 ) const
811 {
812  // _u, _v: vectors u and v to determine the next vertex
813  bool exactIntersection;
814  // initializes A, B, u, v and the two first vertices of resultup and
815  // resultdown.
816  _inPts.resize( 1 ); //to store half convex hull border
817  _outPts.resize( 1 ); //to store half convex hull border
818  _inPts[ 0 ] = pointRefC1;
819 
820  // exactIntersection is equal to one when the intersection of the
821  // supporting lines of C1 and C2 corresponds to an integer point
822  exactIntersection = getFirstPointsOfHull( _v, _inPts[ 0 ], _outPts[ 0 ],
823  hs1, hs2 );
824  if ( ! exactIntersection ) // not integer intersection
825  {
826  //computation of the first part of the border
827  getAllPointsOfHull( _inPts, _outPts, _v, hs2, hs3 );
828  }
829  // fill in convexup
830  for( Size i = 0; i < _inPts.size(); ++i )
831  *itOut++ = _inPts[ i ];
832 
833  // second part
834  // initializes A, B, u, v and the two first vertices of resultup and
835  // resultdown.
836  _inPts.resize( 1 );
837  _outPts.resize( 1 );
838  _inPts[0] = pointRefC3;
839  // exactIntersection is equal to one when the intersection of the
840  // supporting lines of C3 and C2 corresponds to an integer point
841  exactIntersection = getFirstPointsOfHull( _v, _inPts[ 0 ], _outPts[ 0 ],
842  hs3, hs2 );
843  if ( ! exactIntersection ) // not integer intersection
844  {
845  //computation of the second part of the border
846  getAllPointsOfHull( _inPts, _outPts, _v, hs2, hs1 ); // check not hs1.
847  }
848  //fill in convexup
849  for( Size i = _inPts.size(); i != 0; --i )
850  {
851  *itOut++ = _inPts[ i - 1 ];
852  }
853  return itOut;
854 }
855 
856 
857 
858 ///////////////////////////////////////////////////////////////////////////////
859 // Interface - public :
860 
861 /**
862  * Writes/Displays the object on an output stream.
863  * @param out the output stream where the object is written.
864  */
865 template <typename TSpace, typename TSequence>
866 inline
867 void
868 DGtal::LatticePolytope2D<TSpace,TSequence>::selfDisplay ( std::ostream & out ) const
869 {
870  out << "[LatticePolytope2D #Vertices=" << size() << "]";
871 }
872 
873 /**
874  * Checks the validity/consistency of the object.
875  * @return 'true' if the object is valid, 'false' otherwise.
876  */
877 template <typename TSpace, typename TSequence>
878 inline
879 bool
880 DGtal::LatticePolytope2D<TSpace,TSequence>::isValid() const
881 {
882  return true;
883 }
884 //-----------------------------------------------------------------------------
885 template <typename TSpace, typename TSequence>
886 inline
887 std::string
888 DGtal::LatticePolytope2D<TSpace,TSequence>::className
889 () const
890 {
891  return "LatticePolytope2D";
892 }
893 
894 
895 
896 ///////////////////////////////////////////////////////////////////////////////
897 // Implementation of inline functions //
898 
899 template <typename TSpace, typename TSequence>
900 inline
901 std::ostream&
902 DGtal::operator<< ( std::ostream & out,
903  const LatticePolytope2D<TSpace,TSequence> & object )
904 {
905  object.selfDisplay( out );
906  return out;
907 }
908 
909 // //
910 ///////////////////////////////////////////////////////////////////////////////