DGtal  0.9.3
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  @return the number of inserted points.
518  @todo this method is for now not efficient.
519  */
520 //-----------------------------------------------------------------------------
521 template <typename TSpace, typename TSequence>
522 template <typename DigitalSet>
523 void
524 DGtal::LatticePolytope2D<TSpace,TSequence>::
525 getIncludedDigitalPoints( DigitalSet & aSet ) const
526 {
527  BOOST_CONCEPT_ASSERT(( concepts::CDigitalSet< DigitalSet > ));
528  typedef typename DigitalSet::Domain DigitalSetDomain;
529  BOOST_STATIC_ASSERT(( concepts::ConceptUtils::SameType< DigitalSetDomain, Domain >::value ));
530  typedef typename DigitalSetDomain::ConstIterator DomainConstIterator;
531  typedef typename DigitalSet::Iterator DigitalSetIterator;
532  aSet.clear();
533  // case 1 : empty polygon.
534  if ( empty() ) return;
535  // case 2 : one vertex.
536  size_t s = size();
537  if ( s == 1 )
538  {
539  aSet.insert( *begin() );
540  return;
541  }
542  // case 3 : 2 vertices
543  if ( s == 2 )
544  {
545  ConstIterator it_vtx = begin();
546  ConstIterator it_vtx2 = it_vtx; ++it_vtx2;
547  HalfSpace hs = halfSpace( it_vtx );
548  Vector orthN( -hs.N[ 1 ], hs.N[ 0 ] );
549  HalfSpace hs_orth( orthN, _ic.dotProduct( orthN, *it_vtx ) );
550  if ( ! hs_orth( *it_vtx2 ) ) hs_orth.negate();
551  HalfSpace hs2 = halfSpace( it_vtx2 );
552  Vector orthN2( -hs2.N[ 1 ], hs2.N[ 0 ] );
553  HalfSpace hs2_orth( orthN2, _ic.dotProduct( orthN2, *it_vtx2 ) );
554  if ( ! hs2_orth( *it_vtx ) ) hs2_orth.negate();
555 
556  for ( DomainConstIterator it = aSet.domain().begin(),
557  it_end = aSet.domain().end(); it != it_end; ++it )
558  if ( hs( *it ) ) aSet.insert( *it );
559  for ( DigitalSetIterator it_set = aSet.begin(), it_set_end = aSet.end();
560  it_set != it_set_end; )
561  {
562  DigitalSetIterator it_cur = it_set; ++it_set;
563  if ( ! ( hs_orth( *it_cur )
564  && hs2( *it_cur )
565  && hs2_orth( *it_cur ) ) )
566  aSet.erase( it_cur );
567  }
568  return;
569  }
570  // case 4 : >= 3 vertices
571  ConstIterator it_vtx = begin();
572  ConstIterator it_vtx_end = end();
573  HalfSpace hs = halfSpace( it_vtx );
574  for ( DomainConstIterator it = aSet.domain().begin(),
575  it_end = aSet.domain().end(); it != it_end; ++it )
576  if ( hs( *it ) ) aSet.insert( *it );
577  ++it_vtx;
578  for ( ; it_vtx != it_vtx_end; ++it_vtx )
579  {
580  hs = halfSpace( it_vtx );
581  for ( DigitalSetIterator it_set = aSet.begin(), it_set_end = aSet.end();
582  it_set != it_set_end; )
583  {
584  DigitalSetIterator it_cur = it_set; ++it_set;
585  if ( ! hs( *it_cur ) ) aSet.erase( it_cur );
586  }
587  }
588 }
589 //-----------------------------------------------------------------------------
590 template <typename TSpace, typename TSequence>
591 bool
592 DGtal::LatticePolytope2D<TSpace,TSequence>::
593 getFirstPointsOfHull
594 ( Vector & v,
595  Point & inPt, // must belong to hs1.
596  Point & outPt,
597  const HalfSpace & hs1,
598  const HalfSpace & hs2 ) const
599 {
600  ASSERT( hs1.isOnBoundary( inPt ) );
601  bool exactIntersection;
602 
603  // initialize vector directionVector (not definitive)
604  _DV = hs1.tangent();
605 
606  //compute the intersection of ray (inPt, _DV) with constraint C2
607  _ic.getCoefficientIntersection( _fl, _ce,
608  inPt, _DV, hs2.N, hs2.c );
609 
610  // uses the intersection to compute the first vertex of the upper
611  // convex hull (inside convex hull), i.e. the grid point closest to
612  // C2 and satisfying C2
613  _ic.getDotProduct( _a, inPt, hs2.N );
614  // if ( ( _a > hs2.c ) && ( _fl == NumberTraits<Integer>::ZERO ) )
615  // {
616  // inPt += _DV * _ce;
617  // _DV.neg();
618  // }
619  // else if ( ( ( _a <= hs2.c ) && ( _fl >= NumberTraits<Integer>::ZERO ) )
620  // || ( ( _a > hs2.c2 ) && ( _fl <= NumberTraits<Integer>::ZERO ) ) )
621  if ( ( ( _a <= hs2.c ) && ( _fl >= NumberTraits<Integer>::ZERO ) )
622  || ( ( _a > hs2.c ) && ( _fl < NumberTraits<Integer>::ZERO ) ) )
623  {
624  inPt += _DV * _fl;
625  }
626  else
627  {
628  inPt += _DV * _ce;
629  _DV.negate();
630  }
631 
632  // compute the first vertex of the lower convex hull
633  if ( _fl == _ce ) //integer intersection
634  {
635  outPt = inPt;
636  exactIntersection = true;
637  }
638  else
639  {
640  outPt = inPt + _DV;
641  //initialization of v: valid Bezout vector of u
642  _ic.getValidBezout( v, inPt, _DV, hs1.N, hs1.c, hs2.N, hs2.c, true );
643  exactIntersection = false;
644  }
645 
646 #ifdef DEBUG_LatticePolytope2D
647  std::cerr << "[CIP::getFirstPointsOfHull] in=" << inPt
648  << " out=" << outPt << " v=" << v << std::endl;
649 #endif //DEBUG_LatticePolytope2D
650  ASSERT( hs1.isOnBoundary( inPt ) );
651  ASSERT( hs1.isOnBoundary( outPt ) );
652  ASSERT( hs2( inPt ) );
653  ASSERT( ( exactIntersection && hs2( outPt ) ) || ( ! exactIntersection && ! hs2( outPt ) ) );
654  return exactIntersection;
655 }
656 
657 /**
658  Computes the border of the upper and of the lower convex hull
659  from the starting points inPts[0] (up) and outPts[0]
660  down, along the constraint N2.p <= c2 while the vertices
661  satisfy the constraint N3.p <= c3. The vertices of the two
662  borders are stored at the end of inPts and outPts.
663 
664  @param inPts (in, out) as input, contains the first point, as
665  output the sequence of points satisfying \a hs2 and \a hs3.
666 
667  @param outPts (in, out) as input, contains the first point, as
668  output the sequence of points not satisfying \a hs2 and satisfying
669  \a hs3.
670 
671  @param BV the Bezout vector of the vector between inPts[ 0 ] and outPts[ 0 ].
672 
673  @param hs2 the half-space that is approached by the two sequences of points.
674  @param hs3 the limiting half-space which defines the bounds of the approximation.
675  */
676 //-----------------------------------------------------------------------------
677 template <typename TSpace, typename TSequence>
678 void
679 DGtal::LatticePolytope2D<TSpace,TSequence>::
680 getAllPointsOfHull( std::vector<Point> & inPts,
681  std::vector<Point> & outPts,
682  const Vector & BV,
683  const HalfSpace & hs2,
684  const HalfSpace & hs3 ) const
685 {
686  ASSERT( hs2.N != hs3.N );
687 
688  // _A and _B represents the last computed vertex above and below the constraint resp.
689  // _u, _v pair of Bezout vectors.
690  _v = BV;
691  // initialize A and B
692  _A = inPts[0];
693  _B = outPts[0];
694 
695  // while A and B do not lie on the supporting line of hs2
696  // and satisfies hs3 and while v is not parallel to hs2
697  while ( ( ! hs2.isOnBoundary( _A ) ) // stops if A reaches hs2.
698  && ( hs3( _A ) ) // stops if A does not satisfy hs3.
699  && ( ! hs2.isOnBoundary( _B ) ) // stops if B reaches hs2.
700  && ( hs3( _B ) ) // stops if B does not satisfy hs3.
701  && ( _ic.dotProduct( _v, hs2.N ) != NumberTraits<Integer>::ZERO ) )
702  {
703 #ifdef DEBUG_LatticePolytope2D
704  std::cerr << "[CIP::getAllPointsOfHull] A=" << _A
705  << " B=" << _B << " v= " << _v << std::endl;
706 #endif // DEBUG_LatticePolytope2D
707  if ( _ic.dotProduct( _v, hs2.N ) < NumberTraits<Integer>::ZERO )
708  { //------ second configuration, we find a new B --------------------
709  _ic.getCoefficientIntersection( _fl, _ce,
710  _B, _v, hs2.N, hs2.c );
711  _u = _B; _u += _v *_fl;
712  if ( hs3( _u ) )
713  { // Point is still within bounds.
714  _B = _u;
715  outPts.push_back( _B );
716 #ifdef DEBUG_LatticePolytope2D
717  std::cerr << "[CIP::getAllPointsOfHull] add B=" << _B << std::endl;
718 #endif // ifdef DEBUG_LatticePolytope2D
719  }
720  else
721  { // Point is too far away. Check intersection with hs3 instead.
722  _ic.getCoefficientIntersection( _fl, _ce,
723  _B, _v, hs3.N, hs3.c );
724  _B += _v * _fl;
725  outPts.push_back( _B );
726 #ifdef DEBUG_LatticePolytope2D
727  std::cerr << "[CIP::getAllPointsOfHull] add B=" << _B << std::endl;
728 #endif //ifdef DEBUG_LatticePolytope2D
729  ASSERT( hs3( _B ) );
730  break; // necessarily the last point.
731  }
732  }
733  else
734  { //----- first configuration, we find a new A -----------------------
735  _ic.getCoefficientIntersection( _fl, _ce,
736  _A, _v, hs2.N, hs2.c );
737  _u = _A; _u += _v *_fl;
738  if ( hs3( _u ) )
739  { // Point is still within bounds.
740  _A = _u;
741  inPts.push_back( _A );
742 #ifdef DEBUG_LatticePolytope2D
743  std::cerr << "[CIP::getAllPointsOfHull] add A=" << _A << std::endl;
744 #endif //ifdef DEBUG_LatticePolytope2D
745  }
746  else
747  { // Point is too far away. Check intersection with hs3 instead.
748  _ic.getCoefficientIntersection( _fl, _ce,
749  _A, _v, hs3.N, hs3.c );
750  _A += _v * _fl;
751  inPts.push_back( _A );
752 #ifdef DEBUG_LatticePolytope2D
753  std::cerr << "[CIP::getAllPointsOfHull] add A=" << _A << std::endl;
754 #endif //ifdef DEBUG_LatticePolytope2D
755  ASSERT( hs3( _A ) );
756  break; // necessarily the last point.
757  }
758  }
759  // update u and v
760  _u = _B;
761  _u -= _A;
762  // From _ic.getValidBezout( _A, _u, _v, N2, c2, N2, c2, 0);
763  // JOL: to check.
764  _ic.getValidBezout( _v, _A, _u, hs2.N, hs2.c, hs2.N, hs2.c, false );
765  }
766  // when the loop finishes, we have to complete the computation
767  // of each convex hull
768  // if ( ! hs3( _A ) ) // A does not satisfy C3, we remove it.
769  // inPts.pop_back();
770  // else if ( ! hs3( _B ) ) // B does not satisfy C3, we remove it
771  // outPts.pop_back();
772  // else
773  if ( hs2.isOnBoundary( _A ) ) // A lies on C2 so A is also a vertex of outPts
774  outPts.push_back( _A );
775  else if ( hs2.isOnBoundary( _B ) ) // B lies on C2 so B is also a vertex of inPts
776  inPts.push_back( _B );
777 }
778 
779 /**
780  * compute the convex hull of grid points satisfying the
781  * constraints N1.P<=c1, N2.P<=c2 and N3.P>=c3.
782  *
783  * N2.P<=c2 corresponds to the cut two parts of computation: from
784  * constraint 1 to constraint 3 and from constraint 3 to
785  * constraint 1.
786  *
787  * The computed vertices are inserted at position [pos] in some list.
788  *
789  * @param pointRefC1 and pointRefC3 corresponds to grid point lying on
790  * the supporting lines of C1 and of C3 resp.
791  *
792  * @param pos corresponds to an iterator in the list of vertices
793  * of the convex, to add the next new vertices
794  *
795  * NB: the method also computes grid point satisfying N1.P<=c1 and
796  * N3.P>=c3 but not satisfying N2.P<=c2. They are stored in
797  * "resultdown" of size "nbverticesdown". the algorithm uses
798  * these points that's why they appear in the code.
799  */
800 //-----------------------------------------------------------------------------
801 template <typename TSpace, typename TSequence>
802 template <typename OutputIterator>
803 OutputIterator
804 DGtal::LatticePolytope2D<TSpace,TSequence>::
805 computeConvexHullBorder
806 ( OutputIterator itOut,
807  const Point & pointRefC1,
808  const Point & pointRefC3,
809  const HalfSpace & hs1,
810  const HalfSpace & hs2,
811  const HalfSpace & hs3 ) const
812 {
813  // _u, _v: vectors u and v to determine the next vertex
814  bool exactIntersection;
815  // initializes A, B, u, v and the two first vertices of resultup and
816  // resultdown.
817  _inPts.resize( 1 ); //to store half convex hull border
818  _outPts.resize( 1 ); //to store half convex hull border
819  _inPts[ 0 ] = pointRefC1;
820 
821  // exactIntersection is equal to one when the intersection of the
822  // supporting lines of C1 and C2 corresponds to an integer point
823  exactIntersection = getFirstPointsOfHull( _v, _inPts[ 0 ], _outPts[ 0 ],
824  hs1, hs2 );
825  if ( ! exactIntersection ) // not integer intersection
826  {
827  //computation of the first part of the border
828  getAllPointsOfHull( _inPts, _outPts, _v, hs2, hs3 );
829  }
830  // fill in convexup
831  for( Size i = 0; i < _inPts.size(); ++i )
832  *itOut++ = _inPts[ i ];
833 
834  // second part
835  // initializes A, B, u, v and the two first vertices of resultup and
836  // resultdown.
837  _inPts.resize( 1 );
838  _outPts.resize( 1 );
839  _inPts[0] = pointRefC3;
840  // exactIntersection is equal to one when the intersection of the
841  // supporting lines of C3 and C2 corresponds to an integer point
842  exactIntersection = getFirstPointsOfHull( _v, _inPts[ 0 ], _outPts[ 0 ],
843  hs3, hs2 );
844  if ( ! exactIntersection ) // not integer intersection
845  {
846  //computation of the second part of the border
847  getAllPointsOfHull( _inPts, _outPts, _v, hs2, hs1 ); // check not hs1.
848  }
849  //fill in convexup
850  for( Size i = _inPts.size(); i != 0; --i )
851  {
852  *itOut++ = _inPts[ i - 1 ];
853  }
854  return itOut;
855 }
856 
857 
858 
859 ///////////////////////////////////////////////////////////////////////////////
860 // Interface - public :
861 
862 /**
863  * Writes/Displays the object on an output stream.
864  * @param out the output stream where the object is written.
865  */
866 template <typename TSpace, typename TSequence>
867 inline
868 void
869 DGtal::LatticePolytope2D<TSpace,TSequence>::selfDisplay ( std::ostream & out ) const
870 {
871  out << "[LatticePolytope2D #Vertices=" << size() << "]";
872 }
873 
874 /**
875  * Checks the validity/consistency of the object.
876  * @return 'true' if the object is valid, 'false' otherwise.
877  */
878 template <typename TSpace, typename TSequence>
879 inline
880 bool
881 DGtal::LatticePolytope2D<TSpace,TSequence>::isValid() const
882 {
883  return true;
884 }
885 //-----------------------------------------------------------------------------
886 template <typename TSpace, typename TSequence>
887 inline
888 std::string
889 DGtal::LatticePolytope2D<TSpace,TSequence>::className
890 () const
891 {
892  return "LatticePolytope2D";
893 }
894 
895 
896 
897 ///////////////////////////////////////////////////////////////////////////////
898 // Implementation of inline functions //
899 
900 template <typename TSpace, typename TSequence>
901 inline
902 std::ostream&
903 DGtal::operator<< ( std::ostream & out,
904  const LatticePolytope2D<TSpace,TSequence> & object )
905 {
906  object.selfDisplay( out );
907  return out;
908 }
909 
910 // //
911 ///////////////////////////////////////////////////////////////////////////////
Generated on Sun Feb 12 2017 10:24:57 for DGtal by   doxygen 1.8.13