DGtal 1.3.0
Loading...
Searching...
No Matches
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//-----------------------------------------------------------------------------
46template <typename TSpace, typename TSequence>
47inline
48DGtal::LatticePolytope2D<TSpace,TSequence>::~LatticePolytope2D()
49{ // Nothing to do.
50}
51//-----------------------------------------------------------------------------
52template <typename TSpace, typename TSequence>
53inline
54DGtal::LatticePolytope2D<TSpace,TSequence>::LatticePolytope2D()
55{ // Nothing to do.
56}
57//-----------------------------------------------------------------------------
58template <typename TSpace, typename TSequence>
59inline
60DGtal::LatticePolytope2D<TSpace,TSequence>::LatticePolytope2D
61( const Self & other )
62 : myVertices( other.myVertices )
63{ // Nothing to do.
64}
65//-----------------------------------------------------------------------------
66template <typename TSpace, typename TSequence>
67inline
68typename DGtal::LatticePolytope2D<TSpace,TSequence>::Self &
69DGtal::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//-----------------------------------------------------------------------------
77template <typename TSpace, typename TSequence>
78inline
79typename DGtal::LatticePolytope2D<TSpace,TSequence>::ConstIterator
80DGtal::LatticePolytope2D<TSpace,TSequence>::
81begin() const
82{
83 return myVertices.begin();
84}
85//-----------------------------------------------------------------------------
86template <typename TSpace, typename TSequence>
87inline
88typename DGtal::LatticePolytope2D<TSpace,TSequence>::ConstIterator
89DGtal::LatticePolytope2D<TSpace,TSequence>::
90end() const
91{
92 return myVertices.end();
93}
94//-----------------------------------------------------------------------------
95template <typename TSpace, typename TSequence>
96inline
97typename DGtal::LatticePolytope2D<TSpace,TSequence>::Iterator
98DGtal::LatticePolytope2D<TSpace,TSequence>::
99begin()
100{
101 return myVertices.begin();
102}
103//-----------------------------------------------------------------------------
104template <typename TSpace, typename TSequence>
105inline
106typename DGtal::LatticePolytope2D<TSpace,TSequence>::Iterator
107DGtal::LatticePolytope2D<TSpace,TSequence>::
108end()
109{
110 return myVertices.end();
111}
112//-----------------------------------------------------------------------------
113template <typename TSpace, typename TSequence>
114inline
115bool
116DGtal::LatticePolytope2D<TSpace,TSequence>::
117empty() const
118{
119 return myVertices.empty();
120}
121//-----------------------------------------------------------------------------
122template <typename TSpace, typename TSequence>
123inline
124typename DGtal::LatticePolytope2D<TSpace,TSequence>::Size
125DGtal::LatticePolytope2D<TSpace,TSequence>::
126size() const
127{
128 return myVertices.size();
129}
130//-----------------------------------------------------------------------------
131template <typename TSpace, typename TSequence>
132inline
133typename DGtal::LatticePolytope2D<TSpace,TSequence>::Size
134DGtal::LatticePolytope2D<TSpace,TSequence>::
135max_size() const
136{
137 return myVertices.max_size();
138}
139//-----------------------------------------------------------------------------
140template <typename TSpace, typename TSequence>
141inline
142void
143DGtal::LatticePolytope2D<TSpace,TSequence>::
144clear()
145{
146 myVertices.clear();
147}
148//-----------------------------------------------------------------------------
149template <typename TSpace, typename TSequence>
150inline
151typename DGtal::LatticePolytope2D<TSpace,TSequence>::Iterator
152DGtal::LatticePolytope2D<TSpace,TSequence>::
153erase( Iterator it )
154{
155 return myVertices.erase( it );
156}
157//-----------------------------------------------------------------------------
158template <typename TSpace, typename TSequence>
159inline
160void
161DGtal::LatticePolytope2D<TSpace,TSequence>::
162swap( LatticePolytope2D & other )
163{
164 myVertices.swap( other.myVertices );
165}
166
167//-----------------------------------------------------------------------------
168template <typename TSpace, typename TSequence>
169inline
170typename DGtal::LatticePolytope2D<TSpace,TSequence>::MyIntegerComputer &
171DGtal::LatticePolytope2D<TSpace,TSequence>::
172ic() const
173{
174 return _ic;
175}
176
177//-----------------------------------------------------------------------------
178template <typename TSpace, typename TSequence>
179typename DGtal::LatticePolytope2D<TSpace,TSequence>::Domain
180DGtal::LatticePolytope2D<TSpace,TSequence>::
181boundingBoxDomain() 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//-----------------------------------------------------------------------------
195template <typename TSpace, typename TSequence>
196void
197DGtal::LatticePolytope2D<TSpace,TSequence>::
198purge()
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//-----------------------------------------------------------------------------
216template <typename TSpace, typename TSequence>
217inline
218typename DGtal::LatticePolytope2D<TSpace,TSequence>::Iterator
219DGtal::LatticePolytope2D<TSpace,TSequence>::
220insertBefore( const Iterator & pos, const Point & K )
221{
222 return myVertices.insert( pos, K );
223}
224//-----------------------------------------------------------------------------
225template <typename TSpace, typename TSequence>
226inline
227void
228DGtal::LatticePolytope2D<TSpace,TSequence>::
229pushBack( const Point & K )
230{
231 myVertices.push_back( K );
232}
233//-----------------------------------------------------------------------------
234template <typename TSpace, typename TSequence>
235inline
236void
237DGtal::LatticePolytope2D<TSpace,TSequence>::
238pushFront( const Point & K )
239{
240 myVertices.push_front( K );
241}
242//-----------------------------------------------------------------------------
243template <typename TSpace, typename TSequence>
244inline
245void
246DGtal::LatticePolytope2D<TSpace,TSequence>::
247push_back( const Point & K )
248{
249 myVertices.push_back( K );
250}
251//-----------------------------------------------------------------------------
252template <typename TSpace, typename TSequence>
253inline
254void
255DGtal::LatticePolytope2D<TSpace,TSequence>::
256push_front( const Point & K )
257{
258 myVertices.push_front( K );
259}
260//-----------------------------------------------------------------------------
261template <typename TSpace, typename TSequence>
262inline
263const typename DGtal::LatticePolytope2D<TSpace,TSequence>::Integer &
264DGtal::LatticePolytope2D<TSpace,TSequence>::
265twiceArea() 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//-----------------------------------------------------------------------------
280template <typename TSpace, typename TSequence>
281inline
282typename DGtal::LatticePolytope2D<TSpace,TSequence>::Point3I
283DGtal::LatticePolytope2D<TSpace,TSequence>::
284centroid() const
285{
286 _a = twiceArea();
287 return centroid( _a );
288}
289//-----------------------------------------------------------------------------
290template <typename TSpace, typename TSequence>
291inline
292typename DGtal::LatticePolytope2D<TSpace,TSequence>::Point3I
293DGtal::LatticePolytope2D<TSpace,TSequence>::
294centroid( 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//-----------------------------------------------------------------------------
330template <typename TSpace, typename TSequence>
331inline
332typename DGtal::LatticePolytope2D<TSpace,TSequence>::Integer
333DGtal::LatticePolytope2D<TSpace,TSequence>::
334numberBoundaryPoints() 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//-----------------------------------------------------------------------------
353template <typename TSpace, typename TSequence>
354inline
355typename DGtal::LatticePolytope2D<TSpace,TSequence>::Integer
356DGtal::LatticePolytope2D<TSpace,TSequence>::
357numberInteriorPoints() 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//-----------------------------------------------------------------------------
366template <typename TSpace, typename TSequence>
367inline
368typename DGtal::LatticePolytope2D<TSpace,TSequence>::SizeCouple
369DGtal::LatticePolytope2D<TSpace,TSequence>::
370findCut( 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//-----------------------------------------------------------------------------
411template <typename TSpace, typename TSequence>
412inline
413bool
414DGtal::LatticePolytope2D<TSpace,TSequence>::
415cut( 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//-----------------------------------------------------------------------------
465template <typename TSpace, typename TSequence>
466inline
467typename DGtal::LatticePolytope2D<TSpace,TSequence>::HalfSpace
468DGtal::LatticePolytope2D<TSpace,TSequence>::
469halfSpace( 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//-----------------------------------------------------------------------------
490template <typename TSpace, typename TSequence>
491inline
492typename DGtal::LatticePolytope2D<TSpace,TSequence>::HalfSpace
493DGtal::LatticePolytope2D<TSpace,TSequence>::
494halfSpace( 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//-----------------------------------------------------------------------------
520template <typename TSpace, typename TSequence>
521template <typename DigitalSet>
522void
523DGtal::LatticePolytope2D<TSpace,TSequence>::
524getIncludedDigitalPoints( 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//-----------------------------------------------------------------------------
589template <typename TSpace, typename TSequence>
590bool
591DGtal::LatticePolytope2D<TSpace,TSequence>::
592getFirstPointsOfHull
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//-----------------------------------------------------------------------------
676template <typename TSpace, typename TSequence>
677void
678DGtal::LatticePolytope2D<TSpace,TSequence>::
679getAllPointsOfHull( 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//-----------------------------------------------------------------------------
800template <typename TSpace, typename TSequence>
801template <typename OutputIterator>
802OutputIterator
803DGtal::LatticePolytope2D<TSpace,TSequence>::
804computeConvexHullBorder
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 */
865template <typename TSpace, typename TSequence>
866inline
867void
868DGtal::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 */
877template <typename TSpace, typename TSequence>
878inline
879bool
880DGtal::LatticePolytope2D<TSpace,TSequence>::isValid() const
881{
882 return true;
883}
884//-----------------------------------------------------------------------------
885template <typename TSpace, typename TSequence>
886inline
887std::string
888DGtal::LatticePolytope2D<TSpace,TSequence>::className
889() const
890{
891 return "LatticePolytope2D";
892}
893
894
895
896///////////////////////////////////////////////////////////////////////////////
897// Implementation of inline functions //
898
899template <typename TSpace, typename TSequence>
900inline
901std::ostream&
902DGtal::operator<< ( std::ostream & out,
903 const LatticePolytope2D<TSpace,TSequence> & object )
904{
905 object.selfDisplay( out );
906 return out;
907}
908
909// //
910///////////////////////////////////////////////////////////////////////////////
Generated on Sun Nov 27 2022 15:15:34 for DGtal by doxygen 1.9.5