DGtal  1.5.beta
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 //-----------------------------------------------------------------------------
42 template <typename TSpace, typename TInternalInteger>
43 inline
44 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
45 ~COBANaivePlaneComputer()
46 { // Nothing to do.
47 }
48 //-----------------------------------------------------------------------------
49 template <typename TSpace, typename TInternalInteger>
50 inline
51 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
52 COBANaivePlaneComputer()
53  : myG( NumberTraits<TInternalInteger>::ZERO )
54 { // Object is invalid
55 }
56 //-----------------------------------------------------------------------------
57 template <typename TSpace, typename TInternalInteger>
58 inline
59 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
60 COBANaivePlaneComputer( 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 //-----------------------------------------------------------------------------
71 template <typename TSpace, typename TInternalInteger>
72 inline
73 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger> &
74 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
75 operator=( 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 //-----------------------------------------------------------------------------
90 template <typename TSpace, typename TInternalInteger>
91 inline
92 typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::MyIntegerComputer &
93 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
94 ic() const
95 {
96  return myState.cip.ic();
97 }
98 //-----------------------------------------------------------------------------
99 template <typename TSpace, typename TInternalInteger>
100 inline
101 void
102 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
103 clear()
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 //-----------------------------------------------------------------------------
115 template <typename TSpace, typename TInternalInteger>
116 inline
117 void
118 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
119 init( 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 //-----------------------------------------------------------------------------
136 template <typename TSpace, typename TInternalInteger>
137 inline
138 typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::ConstIterator
139 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
140 begin() const
141 {
142  return myPointSet.begin();
143 }
144 //-----------------------------------------------------------------------------
145 template <typename TSpace, typename TInternalInteger>
146 inline
147 typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::ConstIterator
148 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
149 end() const
150 {
151  return myPointSet.end();
152 }
153 
154 //-----------------------------------------------------------------------------
155 template <typename TSpace, typename TInternalInteger>
156 inline
157 typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Size
158 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
159 size() const
160 {
161  return myPointSet.size();
162 }
163 //-----------------------------------------------------------------------------
164 template <typename TSpace, typename TInternalInteger>
165 inline
166 bool
167 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
168 empty() const
169 {
170  return myPointSet.empty();
171 }
172 //-----------------------------------------------------------------------------
173 template <typename TSpace, typename TInternalInteger>
174 inline
175 typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Size
176 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
177 max_size() const
178 {
179  return myPointSet.max_size();
180 }
181 //-----------------------------------------------------------------------------
182 template <typename TSpace, typename TInternalInteger>
183 inline
184 typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Size
185 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
186 maxSize() const
187 {
188  return max_size();
189 }
190 //-----------------------------------------------------------------------------
191 template <typename TSpace, typename TInternalInteger>
192 inline
193 typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Size
194 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
195 complexity() const
196 {
197  return myState.cip.size();
198 }
199 //-----------------------------------------------------------------------------
200 template <typename TSpace, typename TInternalInteger>
201 inline
202 bool
203 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
204 operator()( const Point & p ) const
205 {
206  ic().getDotProduct( _v, myState.N, p );
207  return ( _v >= myState.min ) && ( _v <= myState.max );
208 }
209 //-----------------------------------------------------------------------------
210 template <typename TSpace, typename TInternalInteger>
211 inline
212 bool
213 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
214 extendAsIs( 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 //-----------------------------------------------------------------------------
223 template <typename TSpace, typename TInternalInteger>
224 bool
225 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
226 extend( 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 //-----------------------------------------------------------------------------
319 template <typename TSpace, typename TInternalInteger>
320 bool
321 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
322 isExtendable( 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 //-----------------------------------------------------------------------------
384 template <typename TSpace, typename TInternalInteger>
385 template <typename TInputIterator>
386 bool
387 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
388 extend( 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 //-----------------------------------------------------------------------------
486 template <typename TSpace, typename TInternalInteger>
487 template <typename TInputIterator>
488 bool
489 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
490 isExtendable( 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 //-----------------------------------------------------------------------------
563 template <typename TSpace, typename TInternalInteger>
564 inline
565 typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Primitive
566 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
567 primitive() 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 //-----------------------------------------------------------------------------
581 template <typename TSpace, typename TInternalInteger>
582 template <typename Vector3D>
583 inline
584 void
585 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
586 getNormal( 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 //-----------------------------------------------------------------------------
608 template <typename TSpace, typename TInternalInteger>
609 inline
610 const typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::IntegerVector3 &
611 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
612 exactNormal() const
613 {
614  return myState.N;
615 }
616 //-----------------------------------------------------------------------------
617 //-----------------------------------------------------------------------------
618 template <typename TSpace, typename TInternalInteger>
619 template <typename Vector3D>
620 inline
621 void
622 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
623 getUnitNormal( 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 //-----------------------------------------------------------------------------
634 template <typename TSpace, typename TInternalInteger>
635 inline
636 void
637 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
638 getBounds( 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 //-----------------------------------------------------------------------------
648 template <typename TSpace, typename TInternalInteger>
649 inline
650 const typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Point &
651 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
652 minimalPoint() const
653 {
654  ASSERT( ! this->empty() );
655  return myState.ptMin;
656 }
657 //-----------------------------------------------------------------------------
658 template <typename TSpace, typename TInternalInteger>
659 inline
660 const typename DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::Point &
661 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
662 maximalPoint() 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  */
677 template <typename TSpace, typename TInternalInteger>
678 inline
679 void
680 DGtal::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  */
700 template <typename TSpace, typename TInternalInteger>
701 inline
702 bool
703 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::isValid() const
704 {
705  return myG != NumberTraits< InternalInteger >::ZERO;
706 }
707 
708 
709 ///////////////////////////////////////////////////////////////////////////////
710 // Internals
711 ///////////////////////////////////////////////////////////////////////////////
712 //-----------------------------------------------------------------------------
713 template <typename TSpace, typename TInternalInteger>
714 inline
715 void
716 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
717 computeCentroidAndNormal( 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 //-----------------------------------------------------------------------------
742 template <typename TSpace, typename TInternalInteger>
743 inline
744 void
745 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
746 doubleCut( 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 //-----------------------------------------------------------------------------
760 template <typename TSpace, typename TInternalInteger>
761 template <typename TInputIterator>
762 void
763 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
764 computeMinMax( 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 //-----------------------------------------------------------------------------
790 template <typename TSpace, typename TInternalInteger>
791 template <typename TInputIterator>
792 bool
793 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
794 updateMinMax( 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 //-----------------------------------------------------------------------------
820 template <typename TSpace, typename TInternalInteger>
821 bool
822 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
823 checkPlaneWidth( 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 //-----------------------------------------------------------------------------
830 template <typename TSpace, typename TInternalInteger>
831 void
832 DGtal::COBANaivePlaneComputer<TSpace, TInternalInteger>::
833 computeGradient( 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 
856 template <typename TSpace, typename TInternalInteger>
857 inline
858 std::ostream&
859 DGtal::operator<< ( std::ostream & out,
860  const COBANaivePlaneComputer<TSpace, TInternalInteger> & object )
861 {
862  object.selfDisplay( out );
863  return out;
864 }
865 
866 // //
867 ///////////////////////////////////////////////////////////////////////////////
868 
869