DGtal 1.4.0
Loading...
Searching...
No Matches
COBANaivePlaneComputer.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 COBANaivePlaneComputer.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 *
22 * @date 2012/09/20
23 *
24 * Implementation of inline methods defined in COBANaivePlaneComputer.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29
30//////////////////////////////////////////////////////////////////////////////
31#include <cstdlib>
32//////////////////////////////////////////////////////////////////////////////
33
34///////////////////////////////////////////////////////////////////////////////
35// IMPLEMENTATION of inline methods.
36///////////////////////////////////////////////////////////////////////////////
37
38///////////////////////////////////////////////////////////////////////////////
39// ----------------------- Standard services ------------------------------
40
41//-----------------------------------------------------------------------------
42template <typename TSpace, typename TInternalInteger>
43inline
44DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
45~COBANaivePlaneComputer()
46{ // Nothing to do.
47}
48//-----------------------------------------------------------------------------
49template <typename TSpace, typename TInternalInteger>
50inline
51DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
52COBANaivePlaneComputer()
53 : myG( NumberTraits<TInternalInteger>::ZERO )
54{ // Object is invalid
55}
56//-----------------------------------------------------------------------------
57template <typename TSpace, typename TInternalInteger>
58inline
59DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
60COBANaivePlaneComputer( const COBANaivePlaneComputer & other )
61 : myAxis( other.myAxis ),
62 myG( other.myG ),
63 myWidth( other.myWidth ),
64 myPointSet( other.myPointSet ),
65 myState( other.myState ),
66 myCst1( other.myCst1 ),
67 myCst2( other.myCst2 )
68{
69}
70//-----------------------------------------------------------------------------
71template <typename TSpace, typename TInternalInteger>
72inline
73DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger> &
74DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
75operator=( const COBANaivePlaneComputer & other )
76{
77 if ( this != &other )
78 {
79 myAxis = other.myAxis;
80 myG = other.myG;
81 myWidth = other.myWidth;
82 myPointSet = other.myPointSet;
83 myState = other.myState;
84 myCst1 = other.myCst1;
85 myCst2 = other.myCst2;
86 }
87 return *this;
88}
89//-----------------------------------------------------------------------------
90template <typename TSpace, typename TInternalInteger>
91inline
92typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::MyIntegerComputer &
93DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
94ic() const
95{
96 return myState.cip.ic();
97}
98//-----------------------------------------------------------------------------
99template <typename TSpace, typename TInternalInteger>
100inline
101void
102DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
103clear()
104{
105 myPointSet.clear();
106 myState.cip.clear();
107 // initialize the search space as a square.
108 myState.cip.pushBack( InternalPoint2( -myG, -myG ) );
109 myState.cip.pushBack( InternalPoint2( myG, -myG ) );
110 myState.cip.pushBack( InternalPoint2( myG, myG ) );
111 myState.cip.pushBack( InternalPoint2( -myG, myG ) );
112 computeCentroidAndNormal( myState );
113}
114//-----------------------------------------------------------------------------
115template <typename TSpace, typename TInternalInteger>
116inline
117void
118DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
119init( Dimension axis, InternalInteger diameter,
120 InternalInteger widthNumerator,
121 InternalInteger widthDenominator )
122{
123 myAxis = axis;
124 myWidth[ 0 ] = widthNumerator;
125 myWidth[ 1 ] = widthDenominator;
126 // initialize the grid step.
127 myG = 2*diameter; myG *= diameter; myG *= diameter;
128 // Initializes some constants
129 // _cst1 = ( (int) ceil( get_si( myG ) * myWidth ) + 1 );
130 // _cst2 = ( (int) floor( get_si( myG ) * myWidth ) - 1 );
131 myCst1 = ( ( myG * myWidth[ 0 ] - 1 ) / myWidth[ 1 ] ) + 2;
132 myCst2 = ( ( myG * myWidth[ 0 ] ) / myWidth[ 1 ] ) - 1;
133 clear();
134}
135//-----------------------------------------------------------------------------
136template <typename TSpace, typename TInternalInteger>
137inline
138typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::ConstIterator
139DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
140begin() const
141{
142 return myPointSet.begin();
143}
144//-----------------------------------------------------------------------------
145template <typename TSpace, typename TInternalInteger>
146inline
147typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::ConstIterator
148DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
149end() const
150{
151 return myPointSet.end();
152}
153
154//-----------------------------------------------------------------------------
155template <typename TSpace, typename TInternalInteger>
156inline
157typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Size
158DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
159size() const
160{
161 return myPointSet.size();
162}
163//-----------------------------------------------------------------------------
164template <typename TSpace, typename TInternalInteger>
165inline
166bool
167DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
168empty() const
169{
170 return myPointSet.empty();
171}
172//-----------------------------------------------------------------------------
173template <typename TSpace, typename TInternalInteger>
174inline
175typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Size
176DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
177max_size() const
178{
179 return myPointSet.max_size();
180}
181//-----------------------------------------------------------------------------
182template <typename TSpace, typename TInternalInteger>
183inline
184typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Size
185DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
186maxSize() const
187{
188 return max_size();
189}
190//-----------------------------------------------------------------------------
191template <typename TSpace, typename TInternalInteger>
192inline
193typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Size
194DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
195complexity() const
196{
197 return myState.cip.size();
198}
199//-----------------------------------------------------------------------------
200template <typename TSpace, typename TInternalInteger>
201inline
202bool
203DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
204operator()( const Point & p ) const
205{
206 ic().getDotProduct( _v, myState.N, p );
207 return ( _v >= myState.min ) && ( _v <= myState.max );
208}
209//-----------------------------------------------------------------------------
210template <typename TSpace, typename TInternalInteger>
211inline
212bool
213DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
214extendAsIs( const Point & p )
215{
216 ASSERT( isValid() && ! empty() );
217 bool ok = this->operator()( p );
218 if ( ok ) myPointSet.insert( p );
219 return ok;
220}
221
222//-----------------------------------------------------------------------------
223template <typename TSpace, typename TInternalInteger>
224bool
225DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
226extend( const Point & p )
227{
228 ASSERT( isValid() );
229 // Checks if first point.
230 if ( empty() )
231 {
232 myPointSet.insert( p );
233 ic().getDotProduct( myState.max, myState.N, p );
234 myState.min = myState.max;
235 myState.ptMax = myState.ptMin = p;
236 return true;
237 }
238
239 // Check first if p is already a point of the plane.
240 if ( myPointSet.find( p ) != myPointSet.end() ) // already in set
241 return true;
242 // Check if p lies within the current bounds of the plane.
243 _state.N = myState.N;
244 _state.min = myState.min;
245 _state.max = myState.max;
246 _state.ptMin = myState.ptMin;
247 _state.ptMax = myState.ptMax;
248 bool changed = updateMinMax( _state, &p, (&p)+1 );
249 // Check if point is already within bounds.
250 if ( ! changed )
251 {
252 myPointSet.insert( p );
253 return true;
254 }
255 // Check if width is still ok
256 if ( checkPlaneWidth( _state ) )
257 {
258 myState.min = _state.min;
259 myState.max = _state.max;
260 myState.ptMin = _state.ptMin;
261 myState.ptMax = _state.ptMax;
262 myPointSet.insert( p );
263 return true;
264 }
265 // We have to find a new normal. First, update gradient.
266 computeGradient( _grad, _state );
267
268 // Checks if we can change the normal so as to find another digital plane.
269 if( ( ( _grad[ 0 ] == NumberTraits<InternalInteger>::ZERO )
270 && ( _grad[ 1 ] == NumberTraits<InternalInteger>::ZERO ) ) )
271 {
272 // Unable to update solution.
273 return false;
274 }
275
276 // There is a gradient, tries to optimize.
277 _state.cip = myState.cip;
278 doubleCut( _grad, _state );
279
280 // While at least 1 point left on the search space
281 while ( ! _state.cip.empty() )
282 {
283 computeCentroidAndNormal( _state );
284 // Calls oracle
285 computeMinMax( _state, myPointSet.begin(), myPointSet.end() );
286 updateMinMax( _state, &p, (&p)+1 );
287 // Check if width is now ok
288 if ( checkPlaneWidth( _state ) )
289 { // Found a plane.
290 myState.min = _state.min;
291 myState.max = _state.max;
292 myState.ptMin = _state.ptMin;
293 myState.ptMax = _state.ptMax;
294 myState.cip.swap( _state.cip );
295 myState.centroid = _state.centroid;
296 myState.N = _state.N;
297 myPointSet.insert( p );
298 return true;
299 }
300
301 // We have to find a new normal. First, update gradient.
302 computeGradient( _grad, _state );
303
304 // Checks if we can change the normal so as to find another digital plane.
305 if ( ( ( _grad[ 0 ] == NumberTraits<InternalInteger>::ZERO )
306 && ( _grad[ 1 ] == NumberTraits<InternalInteger>::ZERO ) ) )
307 {
308 // Unable to update solution. Removes point from set.
309 break;
310 }
311
312 // There is a gradient, tries to optimize.
313 doubleCut( _grad, _state );
314 }
315 // was unable to find a correct plane.
316 return false;
317}
318//-----------------------------------------------------------------------------
319template <typename TSpace, typename TInternalInteger>
320bool
321DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
322isExtendable( const Point & p ) const
323{
324 ASSERT( isValid() );
325 // Checks if first point.
326 if ( empty() ) return true;
327
328 // Check first if p is already a point of the plane.
329 if ( myPointSet.find( p ) != myPointSet.end() ) // already in set
330 return true;
331 // Check if p lies within the current bounds of the plane.
332 _state.N = myState.N;
333 _state.min = myState.min;
334 _state.max = myState.max;
335 _state.ptMin = myState.ptMin;
336 _state.ptMax = myState.ptMax;
337 bool changed = updateMinMax( _state, (&p), (&p)+1 );
338 // Check if point is already within bounds.
339 if ( ! changed ) return true;
340 // Check if width is still ok
341 if ( checkPlaneWidth( _state ) )
342 return true;
343 // We have to find a new normal. First, update gradient.
344 computeGradient( _grad, _state );
345
346 // Checks if we can change the normal so as to find another digital plane.
347 if( ( ( _grad[ 0 ] == NumberTraits<InternalInteger>::ZERO )
348 && ( _grad[ 1 ] == NumberTraits<InternalInteger>::ZERO ) ) )
349 // Unable to update solution.
350 return false;
351
352 // There is a gradient, tries to optimize.
353 _state.cip = myState.cip;
354 doubleCut( _grad, _state );
355
356 // While at least 1 point left on the search space
357 while ( ! _state.cip.empty() )
358 {
359 computeCentroidAndNormal( _state );
360 // Calls oracle
361 computeMinMax( _state, myPointSet.begin(), myPointSet.end() );
362 updateMinMax( _state, (&p), (&p)+1 );
363 // Check if width is now ok
364 if ( checkPlaneWidth( _state ) )
365 // Found a plane.
366 return true;
367
368 // We have to find a new normal. First, update gradient.
369 computeGradient( _grad, _state );
370
371 // Checks if we can change the normal so as to find another digital plane.
372 if ( ( ( _grad[ 0 ] == NumberTraits<InternalInteger>::ZERO )
373 && ( _grad[ 1 ] == NumberTraits<InternalInteger>::ZERO ) ) )
374 // Unable to update solution.
375 return false;
376
377 // There is a gradient, tries to optimize.
378 doubleCut( _grad, _state );
379 }
380 // was unable to find a correct plane.
381 return false;
382}
383//-----------------------------------------------------------------------------
384template <typename TSpace, typename TInternalInteger>
385template <typename TInputIterator>
386bool
387DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
388extend( TInputIterator it, TInputIterator itE )
389{
390 BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
391
392 ASSERT( isValid() );
393 if ( it == itE ) return true;
394 // Check if points lies within the current bounds of the plane.
395 bool changed;
396 _state.N = myState.N;
397 if ( empty() )
398 {
399 changed = true;
400 computeMinMax( _state, it, itE );
401 }
402 else
403 {
404 _state.min = myState.min;
405 _state.max = myState.max;
406 _state.ptMin = myState.ptMin;
407 _state.ptMax = myState.ptMax;
408 changed = updateMinMax( _state, it, itE );
409 }
410 // Check if points are already within bounds.
411 if ( ! changed )
412 { // All points are within bounds. Put them in pointset.
413 for ( TInputIterator tmpIt = it; tmpIt != itE; ++tmpIt )
414 myPointSet.insert( *tmpIt );
415 return true;
416 }
417 // Check if width is still ok
418 if ( checkPlaneWidth( _state ) )
419 {
420 myState.min = _state.min;
421 myState.max = _state.max;
422 myState.ptMin = _state.ptMin;
423 myState.ptMax = _state.ptMax;
424 for ( TInputIterator tmpIt = it; tmpIt != itE; ++tmpIt )
425 myPointSet.insert( *tmpIt );
426 return true;
427 }
428 // We have to find a new normal. First, update gradient.
429 computeGradient( _grad, _state );
430
431 // Checks if we can change the normal so as to find another digital plane.
432 if( ( ( _grad[ 0 ] == NumberTraits<InternalInteger>::ZERO )
433 && ( _grad[ 1 ] == NumberTraits<InternalInteger>::ZERO ) ) )
434 {
435 // Unable to update solution.
436 return false;
437 }
438
439 // There is a gradient, tries to optimize.
440 _state.cip = myState.cip;
441 doubleCut( _grad, _state );
442
443 // While at least 1 point left on the search space
444 while ( ! _state.cip.empty() )
445 {
446 computeCentroidAndNormal( _state );
447 // Calls oracle
448 if ( ! myPointSet.empty() ) {
449 computeMinMax( _state, myPointSet.begin(), myPointSet.end() );
450 updateMinMax( _state, it, itE );
451 }
452 else computeMinMax( _state, it, itE );
453 // Check if width is now ok
454 if ( checkPlaneWidth( _state ) )
455 { // Found a plane.
456 myState.min = _state.min;
457 myState.max = _state.max;
458 myState.ptMin = _state.ptMin;
459 myState.ptMax = _state.ptMax;
460 myState.cip.swap( _state.cip );
461 myState.centroid = _state.centroid;
462 myState.N = _state.N;
463 for ( TInputIterator tmpIt = it; tmpIt != itE; ++tmpIt )
464 myPointSet.insert( *tmpIt );
465 return true;
466 }
467
468 // We have to find a new normal. First, update gradient.
469 computeGradient( _grad, _state );
470
471 // Checks if we can change the normal so as to find another digital plane.
472 if ( ( ( _grad[ 0 ] == NumberTraits<InternalInteger>::ZERO )
473 && ( _grad[ 1 ] == NumberTraits<InternalInteger>::ZERO ) ) )
474 {
475 // Unable to update solution. Removes point from set.
476 break;
477 }
478
479 // There is a gradient, tries to optimize.
480 doubleCut( _grad, _state );
481 }
482 // was unable to find a correct plane.
483 return false;
484}
485//-----------------------------------------------------------------------------
486template <typename TSpace, typename TInternalInteger>
487template <typename TInputIterator>
488bool
489DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
490isExtendable( TInputIterator it, TInputIterator itE ) const
491{
492 BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
493
494 ASSERT( isValid() );
495 if ( it == itE ) return true;
496 // Check if points lies within the current bounds of the plane.
497 bool changed;
498 _state.N = myState.N;
499 if ( empty() )
500 {
501 changed = true;
502 computeMinMax( _state, it, itE );
503 }
504 else
505 {
506 _state.min = myState.min;
507 _state.max = myState.max;
508 _state.ptMin = myState.ptMin;
509 _state.ptMax = myState.ptMax;
510 changed = updateMinMax( _state, it, itE );
511 }
512
513 // Check if point is already within bounds.
514 if ( ! changed ) return true;
515 // Check if width is still ok
516 if ( checkPlaneWidth( _state ) )
517 return true;
518 // We have to find a new normal. First, update gradient.
519 computeGradient( _grad, _state );
520
521 // Checks if we can change the normal so as to find another digital plane.
522 if( ( ( _grad[ 0 ] == NumberTraits<InternalInteger>::ZERO )
523 && ( _grad[ 1 ] == NumberTraits<InternalInteger>::ZERO ) ) )
524 // Unable to update solution.
525 return false;
526
527 // There is a gradient, tries to optimize.
528 _state.cip = myState.cip;
529 doubleCut( _grad, _state );
530
531 // While at least 1 point left on the search space
532 while ( ! _state.cip.empty() )
533 {
534 computeCentroidAndNormal( _state );
535 // Calls oracle
536 if ( ! myPointSet.empty() ) {
537 computeMinMax( _state, myPointSet.begin(), myPointSet.end() );
538 updateMinMax( _state, it, itE );
539 }
540 else computeMinMax( _state, it, itE );
541 // Check if width is now ok
542 if ( checkPlaneWidth( _state ) )
543 // Found a plane.
544 return true;
545
546 // We have to find a new normal. First, update gradient.
547 computeGradient( _grad, _state );
548
549 // Checks if we can change the normal so as to find another digital plane.
550 if ( ( ( _grad[ 0 ] == NumberTraits<InternalInteger>::ZERO )
551 && ( _grad[ 1 ] == NumberTraits<InternalInteger>::ZERO ) ) )
552 // Unable to update solution.
553 return false;
554
555 // There is a gradient, tries to optimize.
556 doubleCut( _grad, _state );
557 }
558 // was unable to find a correct plane.
559 return false;
560}
561
562//-----------------------------------------------------------------------------
563template <typename TSpace, typename TInternalInteger>
564inline
565typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Primitive
566DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
567primitive() const
568{
569 typedef typename Space::RealVector RealVector;
570 typedef typename RealVector::Component Scalar;
571 RealVector N;
572 N[ 0 ] = NumberTraits<InternalInteger>::castToDouble( myState.N[ 0 ] );
573 N[ 1 ] = NumberTraits<InternalInteger>::castToDouble( myState.N[ 1 ] );
574 N[ 2 ] = NumberTraits<InternalInteger>::castToDouble( myState.N[ 2 ] );
575 Scalar min = NumberTraits<InternalInteger>::castToDouble( myState.min );
576 Scalar max = NumberTraits<InternalInteger>::castToDouble( myState.max );
577 return Primitive( min, N, max - min );
578}
579
580//-----------------------------------------------------------------------------
581template <typename TSpace, typename TInternalInteger>
582template <typename Vector3D>
583inline
584void
585DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
586getNormal( Vector3D & normal ) const
587{
588 switch( myAxis ) {
589 case 0 :
590 normal[0] = 1.0;
591 normal[1] = NumberTraits<InternalInteger>::castToDouble( myState.N[ 1 ]) / NumberTraits<InternalInteger>::castToDouble( myState.N[ 0 ] );
592 normal[2] = NumberTraits<InternalInteger>::castToDouble( myState.N[ 2 ]) / NumberTraits<InternalInteger>::castToDouble( myState.N[ 0 ] );
593 break;
594 case 1 :
595 normal[0] = NumberTraits<InternalInteger>::castToDouble( myState.N[ 0 ]) / NumberTraits<InternalInteger>::castToDouble( myState.N[ 1 ] );
596 normal[1] = 1.0;
597 normal[2] = NumberTraits<InternalInteger>::castToDouble( myState.N[ 2 ]) / NumberTraits<InternalInteger>::castToDouble( myState.N[ 1 ] );
598 break;
599 case 2 :
600 normal[0] = NumberTraits<InternalInteger>::castToDouble( myState.N[ 0 ]) / NumberTraits<InternalInteger>::castToDouble( myState.N[ 2 ] );
601 normal[1] = NumberTraits<InternalInteger>::castToDouble( myState.N[ 1 ]) / NumberTraits<InternalInteger>::castToDouble( myState.N[ 2 ] );
602 normal[2] = 1.0;
603 break;
604 }
605}
606//-----------------------------------------------------------------------------
607//-----------------------------------------------------------------------------
608template <typename TSpace, typename TInternalInteger>
609inline
610const typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::IntegerVector3 &
611DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
612exactNormal() const
613{
614 return myState.N;
615}
616//-----------------------------------------------------------------------------
617//-----------------------------------------------------------------------------
618template <typename TSpace, typename TInternalInteger>
619template <typename Vector3D>
620inline
621void
622DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
623getUnitNormal( Vector3D & normal ) const
624{
625 getNormal( normal );
626 double l = sqrt( normal[ 0 ] * normal[ 0 ]
627 + normal[ 1 ] * normal[ 1 ]
628 + normal[ 2 ] * normal[ 2 ] );
629 normal[ 0 ] /= l;
630 normal[ 1 ] /= l;
631 normal[ 2 ] /= l;
632}
633//-----------------------------------------------------------------------------
634template <typename TSpace, typename TInternalInteger>
635inline
636void
637DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
638getBounds( double & min, double & max ) const
639{
640 double nx = NumberTraits<InternalInteger>::castToDouble( myState.N[ 0 ] );
641 double ny = NumberTraits<InternalInteger>::castToDouble( myState.N[ 1 ] );
642 double nz = NumberTraits<InternalInteger>::castToDouble( myState.N[ 2 ] );
643 double l = sqrt( nx*nx + ny*ny + nz*nz );
644 min = NumberTraits<InternalInteger>::castToDouble( myState.min ) / l;
645 max = NumberTraits<InternalInteger>::castToDouble( myState.max ) / l;
646}
647//-----------------------------------------------------------------------------
648template <typename TSpace, typename TInternalInteger>
649inline
650const typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Point &
651DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
652minimalPoint() const
653{
654 ASSERT( ! this->empty() );
655 return myState.ptMin;
656}
657//-----------------------------------------------------------------------------
658template <typename TSpace, typename TInternalInteger>
659inline
660const typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Point &
661DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
662maximalPoint() const
663{
664 ASSERT( ! this->empty() );
665 return myState.ptMax;
666}
667
668
669
670///////////////////////////////////////////////////////////////////////////////
671// Interface - public :
672
673/**
674 * Writes/Displays the object on an output stream.
675 * @param out the output stream where the object is written.
676 */
677template <typename TSpace, typename TInternalInteger>
678inline
679void
680DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::selfDisplay ( std::ostream & out ) const
681{
682 double min, max;
683 double N[] = {0., 0., 0.};
684 out << "[COBANaivePlaneComputer"
685 << " axis=" << myAxis << " w=" << myWidth[ 0 ] << "/" << myWidth[ 1 ]
686 << " size=" << size() << " complexity=" << complexity() << " N=" << myState.N << ": ";
687 this->getUnitNormal( N );
688 this->getBounds( min, max );
689 out << min << " <= "
690 << N[ 0 ] << " * x + "
691 << N[ 1 ] << " * y + "
692 << N[ 2 ] << " * z "
693 << " <= " << max << " ]";
694}
695
696/**
697 * Checks the validity/consistency of the object.
698 * @return 'true' if the object is valid, 'false' otherwise.
699 */
700template <typename TSpace, typename TInternalInteger>
701inline
702bool
703DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::isValid() const
704{
705 return myG != NumberTraits< InternalInteger >::ZERO;
706}
707
708
709///////////////////////////////////////////////////////////////////////////////
710// Internals
711///////////////////////////////////////////////////////////////////////////////
712//-----------------------------------------------------------------------------
713template <typename TSpace, typename TInternalInteger>
714inline
715void
716DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
717computeCentroidAndNormal( State & state ) const
718{
719 if ( state.cip.empty() ) return;
720 state.centroid = state.cip.centroid();
721 ic().reduce( state.centroid );
722 switch( myAxis ){
723 case 0 :
724 state.N[ 0 ] = state.centroid[ 2 ] * myG;
725 state.N[ 1 ] = state.centroid[ 0 ];
726 state.N[ 2 ] = state.centroid[ 1 ];
727 break;
728 case 1 :
729 state.N[ 0 ] = state.centroid[ 0 ];
730 state.N[ 1 ] = state.centroid[ 2 ] * myG;
731 state.N[ 2 ] = state.centroid[ 1 ];
732 break;
733 case 2 :
734 state.N[ 0 ] = state.centroid[ 0 ];
735 state.N[ 1 ] = state.centroid[ 1 ];
736 state.N[ 2 ] = state.centroid[ 2 ] * myG;
737 break;
738 }
739
740}
741//-----------------------------------------------------------------------------
742template <typename TSpace, typename TInternalInteger>
743inline
744void
745DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
746doubleCut( InternalPoint2 & grad, State & state ) const
747{
748 // 2 cuts on the search space:
749 // Gradient.p <= cst1 - _v
750 // -Gradient.p <= cst2 + _v
751 _v = myG * ( state.ptMin[ myAxis ]
752 - state.ptMax[ myAxis ] );
753 state.cip.cut( HalfSpace( grad, myCst1 - _v ) );
754 grad.negate();
755 state.cip.cut( HalfSpace( grad, myCst2 + _v ) );
756 grad.negate();
757}
758
759//-----------------------------------------------------------------------------
760template <typename TSpace, typename TInternalInteger>
761template <typename TInputIterator>
762void
763DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
764computeMinMax( State & state, TInputIterator itB, TInputIterator itE ) const
765{
766 BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
767
768 ASSERT( itB != itE );
769 ic().getDotProduct( state.min, state.N, *itB );
770 state.max = state.min;
771 state.ptMax = state.ptMin = *itB;
772 ++itB;
773 // look for the points defining the min dot product and the max dot product
774 for ( ; itB != itE; ++itB )
775 {
776 ic().getDotProduct( _v, state.N, *itB );
777 if ( _v > state.max )
778 {
779 state.max = _v;
780 state.ptMax = *itB;
781 }
782 else if ( _v < state.min )
783 {
784 state.min = _v;
785 state.ptMin = *itB;
786 }
787 }
788}
789//-----------------------------------------------------------------------------
790template <typename TSpace, typename TInternalInteger>
791template <typename TInputIterator>
792bool
793DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
794updateMinMax( State & state, TInputIterator itB, TInputIterator itE ) const
795
796{
797 BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
798
799 bool changed = false;
800 // look for the points defining the min dot product and the max dot product
801 for ( ; itB != itE; ++itB )
802 {
803 ic().getDotProduct( _v, state.N, *itB );
804 if ( _v > state.max )
805 {
806 state.max = _v;
807 state.ptMax = *itB;
808 changed = true;
809 }
810 else if ( _v < state.min )
811 {
812 state.min = _v;
813 state.ptMin = *itB;
814 changed = true;
815 }
816 }
817 return changed;
818}
819//-----------------------------------------------------------------------------
820template <typename TSpace, typename TInternalInteger>
821bool
822DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
823checkPlaneWidth( const State & state ) const
824{
825 _v = ic().abs( state.N[ myAxis ] );
826 return ( ( state.max - state.min )
827 < ( _v * myWidth[ 0 ] / myWidth[ 1 ] ) );
828}
829//-----------------------------------------------------------------------------
830template <typename TSpace, typename TInternalInteger>
831void
832DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
833computeGradient( InternalPoint2 & grad, const State & state ) const
834{
835 // computation of the gradient
836 switch( myAxis ){
837 case 0 :
838 grad[ 0 ] = state.ptMin[ 1 ] - state.ptMax[ 1 ];
839 grad[ 1 ] = state.ptMin[ 2 ] - state.ptMax[ 2 ];
840 break;
841 case 1:
842 grad[ 0 ] = state.ptMin[ 0 ] - state.ptMax[ 0 ];
843 grad[ 1 ] = state.ptMin[ 2 ] - state.ptMax[ 2 ];
844 break;
845 case 2:
846 grad[ 0 ] = state.ptMin[ 0 ] - state.ptMax[ 0 ];
847 grad[ 1 ] = state.ptMin[ 1 ] - state.ptMax[ 1 ];
848 break;
849 }
850}
851
852
853///////////////////////////////////////////////////////////////////////////////
854// Implementation of inline functions //
855
856template <typename TSpace, typename TInternalInteger>
857inline
858std::ostream&
859DGtal::operator<< ( std::ostream & out,
860 const COBANaivePlaneComputer<TSpace, TInternalInteger> & object )
861{
862 object.selfDisplay( out );
863 return out;
864}
865
866// //
867///////////////////////////////////////////////////////////////////////////////
868
869