DGtal  1.5.beta
ChordNaivePlaneComputer.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 ChordNaivePlaneComputer.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 Yan Gérard
22  * @author Isabelle Debled-Rennesson
23  * @author Paul Zimmermann
24  *
25  * @date 2012/09/20
26  *
27  * Implementation of inline methods defined in ChordNaivePlaneComputer.h
28  *
29  * This file is part of the DGtal library.
30  */
31 
32 
33 //////////////////////////////////////////////////////////////////////////////
34 #include <cstdlib>
35 //////////////////////////////////////////////////////////////////////////////
36 
37 ///////////////////////////////////////////////////////////////////////////////
38 // IMPLEMENTATION of inline methods.
39 ///////////////////////////////////////////////////////////////////////////////
40 
41 ///////////////////////////////////////////////////////////////////////////////
42 // ----------------------- Standard services ------------------------------
43 
44 //-----------------------------------------------------------------------------
45 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
46 inline
47 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
48 ~ChordNaivePlaneComputer()
49 { // Nothing to do.
50 }
51 //-----------------------------------------------------------------------------
52 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
53 inline
54 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
55 ChordNaivePlaneComputer() : z( -1 )
56 { // Object is invalid
57 }
58 //-----------------------------------------------------------------------------
59 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
60 inline
61 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
62 ChordNaivePlaneComputer( const ChordNaivePlaneComputer & other )
63  : z( other.z ), x( other.x ), y( other.y ),
64  myWidth0( other.myWidth0 ),
65  myWidth1( other.myWidth1 ),
66  myPointSet( other.myPointSet ),
67  myState( other.myState )
68 {
69 }
70 //-----------------------------------------------------------------------------
71 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
72 inline
73 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar> &
74 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
75 operator=( const ChordNaivePlaneComputer & other )
76 {
77  if ( this != &other )
78  {
79  z = other.z;
80  x = other.x;
81  y = other.y;
82  myWidth0 = other.myWidth0;
83  myWidth1 = other.myWidth1;
84  myPointSet = other.myPointSet;
85  myState = other.myState;
86  }
87  return *this;
88 }
89 //-----------------------------------------------------------------------------
90 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
91 inline
92 void
93 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
94 clear()
95 {
96  myPointSet.clear();
97  myState.nbValid = 0;
98 }
99 //-----------------------------------------------------------------------------
100 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
101 inline
102 void
103 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
104 init( Dimension axis,
105  InternalScalar widthNumerator,
106  InternalScalar widthDenominator )
107 {
108  ASSERT( ( axis < 3 ) );
109  myWidth0 = widthNumerator;
110  myWidth1 = widthDenominator;
111  switch ( axis ) {
112  case 0: x = 1; y = 2; z = 0; break;
113  case 1: x = 2; y = 0; z = 1; break;
114  case 2: x = 0; y = 1; z = 2; break;
115  }
116  clear();
117 }
118 //-----------------------------------------------------------------------------
119 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
120 inline
121 typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::Size
122 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
123 size() const
124 {
125  return myPointSet.size();
126 }
127 //-----------------------------------------------------------------------------
128 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
129 inline
130 bool
131 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
132 empty() const
133 {
134  return myPointSet.empty();
135 }
136 //-----------------------------------------------------------------------------
137 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
138 inline
139 typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::ConstIterator
140 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
141 begin() const
142 {
143  return myPointSet.begin();
144 }
145 //-----------------------------------------------------------------------------
146 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
147 inline
148 typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::ConstIterator
149 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
150 end() const
151 {
152  return myPointSet.end();
153 }
154 //-----------------------------------------------------------------------------
155 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
156 inline
157 typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::Size
158 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
159 max_size() const
160 {
161  return myPointSet.max_size();
162 }
163 //-----------------------------------------------------------------------------
164 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
165 inline
166 typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::Size
167 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
168 maxSize() const
169 {
170  return max_size();
171 }
172 //-----------------------------------------------------------------------------
173 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
174 inline
175 bool
176 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
177 operator()( const Point & p ) const
178 {
179  _d = internalDot( myState.N, p );
180  return ( _d >= myState.min ) && ( _d <= myState.max );
181 }
182 //-----------------------------------------------------------------------------
183 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
184 inline
185 bool
186 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
187 extendAsIs( const InputPoint & p )
188 {
189  ASSERT( isValid() );
190  if ( empty() ) {
191  myPointSet.insert( p );
192  return setUp1( p );
193  }
194  bool ok = this->operator()( p );
195  if ( ok ) myPointSet.insert( p );
196  return ok;
197 }
198 //-----------------------------------------------------------------------------
199 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
200 bool
201 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
202 extend( const InputPoint & p )
203 {
204  unsigned int loop;
205  ASSERT( isValid() );
206  // Checks if algorithm is initialized.
207  if ( extendAsIs( p ) )
208  {
209  // std::cout << "extended as is: " << p << std::endl;
210  return true;
211  }
212  std::pair<Iterator,bool> ins = myPointSet.insert( p );
213  Iterator itP = ins.first;
214  ASSERT( ins.second == true );
215  if ( myState.nbValid < 3 )
216  { // initial case
217  _state.nbValid = findTriangle( _state, myPointSet.begin(), myPointSet.end() );
218  // std::cout << "findTriangle #=" << _state.nbValid << std::endl;
219  setUpNormal( _state );
220  }
221  else
222  _state = myState;
223  unsigned int result = 0; // 0: computing, 1: too large, 2: found.
224  for ( loop = 0; result == 0; ++loop )
225  {
226  computeHeight( _state );
227  if ( ( _state.height * myWidth1 ) >= ( _state.N[ z ] * myWidth0 ) ) // because width is strict
228  result = 1;
229  else
230  {
231  computeMinMax( _state, myPointSet.begin(), myPointSet.end() );
232  if ( checkWidth( _state ) )
233  { // Ok. Found it.
234  result = 2;
235  }
236  else if ( _state.nbValid >= 2 )
237  { // in case 2 and 3, we have a starting triangle.
238  InputPoint M = _state.ptMax - _state.ptMin;
239  if ( ( M[ x ] == 0 ) && ( M[ y ] == 0 ) )
240  { // Case where M is aligned with axis, width is M[ z ].
241  result = M[ z ] * myWidth1 >= myWidth0 ? 1 : 2;
242  }
243  else
244  {
245  newCurrentTriangle( _state, M );
246  computeNormal( _state );
247  }
248  }
249  else
250  { // in case 1, normal is already aligned with main axis.
251  _state.A = _state.ptMax - _state.ptMin;
252  }
253  }
254  if (loop>=1000) {
255  trace.warning() << "[ChordNaivePlaneComputer::extend()]"
256  << " more than 1000 loops, computing error suspected" << std::endl;
257  trace.warning() << "- Former state: " << *this << std::endl;
258  trace.warning() << "- Current state: ";
259  selfDisplay( trace.warning(), _state );
260  trace.warning() << std::endl;
261  trace.warning() << "- Points: ";
262  for ( ConstIterator it = this->begin(), itE = this->end();
263  it != itE; ++it )
264  trace.warning() << " " << *it;
265  trace.warning() << std::endl;
266  trace.warning() << "- Added Points: " << p;
267  trace.warning() << std::endl;
268  result = 1;
269  }
270  }
271  // was unable to find a correct plane.
272  if ( result == 2 )
273  myState = _state;
274  else
275  myPointSet.erase( itP );
276  return result == 2;
277 }
278 //-----------------------------------------------------------------------------
279 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
280 bool
281 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
282 isExtendable( const InputPoint & p ) const
283 {
284  unsigned int loop;
285  ASSERT( isValid() );
286  if ( empty() || this->operator()( p ) ) return true;
287  std::pair<Iterator,bool> ins = myPointSet.insert( p );
288  Iterator itP = ins.first;
289  ASSERT( ins.second == true );
290  if ( myState.nbValid < 3 )
291  { // initial case
292  _state.nbValid = findTriangle( _state, myPointSet.begin(), myPointSet.end() );
293  // std::cout << "findTriangle #=" << _state.nbValid << std::endl;
294  setUpNormal( _state );
295  }
296  else
297  _state = myState;
298  unsigned int result = 0; // 0: computing, 1: too large, 2: found.
299  for ( loop = 0; result == 0; ++loop )
300  {
301  computeHeight( _state );
302  if ( ( _state.height * myWidth1 ) >= ( _state.N[ z ] * myWidth0 ) ) // check >=
303  result = 1;
304  else
305  {
306  computeMinMax( _state, myPointSet.begin(), myPointSet.end() );
307  if ( checkWidth( _state ) )
308  { // Ok. Found it.
309  result = 2;
310  }
311  else if ( _state.nbValid >= 2 )
312  { // in case 2 and 3, we have a starting triangle.
313  InputPoint M = _state.ptMax - _state.ptMin;
314  if ( ( M[ x ] == 0 ) && ( M[ y ] == 0 ) )
315  { // Case where M is aligned with axis, width is M[ z ].
316  result = M[ z ] * myWidth1 >= myWidth0 ? 1 : 2;
317  }
318  else
319  {
320  newCurrentTriangle( _state, M );
321  computeNormal( _state );
322  }
323  }
324  else
325  { // in case 1, normal is already aligned with main axis.
326  _state.A = _state.ptMax - _state.ptMin;
327  }
328  }
329  if (loop>=1000) {
330  trace.warning() << "[ChordNaivePlaneComputer::isExtendable()]"
331  << " more than 1000 loops, computing error suspected" << std::endl;
332  trace.warning() << "- Former state: " << *this << std::endl;
333  trace.warning() << "- Current state: ";
334  selfDisplay( trace.warning(), _state );
335  trace.warning() << std::endl;
336  trace.warning() << "- Points: ";
337  for ( ConstIterator it = this->begin(), itE = this->end();
338  it != itE; ++it )
339  trace.warning() << " " << *it;
340  trace.warning() << std::endl;
341  trace.warning() << "- Added Points: " << p;
342  trace.warning() << std::endl;
343  result = 1;
344  }
345  }
346  // Goes back to starting state.
347  myPointSet.erase( itP );
348  return result == 2;
349 }
350 //-----------------------------------------------------------------------------
351 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
352 template <typename TInputIterator>
353 bool
354 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
355 satisfies( State & state, TInputIterator itB, TInputIterator itE ) const
356 {
357  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
358  ASSERT( isValid() );
359  state.nbValid = 0;
360  if ( itB == itE ) return true;
361  unsigned int loop;
362  state.nbValid = findTriangle( state, itB, itE );
363  setUpNormal( state );
364  unsigned int result = 0; // 0: computing, 1: too large, 2: found.
365  for ( loop = 0; result == 0; ++loop )
366  {
367  computeHeight( state );
368  if ( ( state.height * myWidth1 ) >= ( state.N[ z ] * myWidth0 ) ) // because width is strict
369  result = 1;
370  else
371  {
372  computeMinMax( state, itB, itE );
373  if ( checkWidth( state ) )
374  result = 2;
375  else if ( state.nbValid >= 2 )
376  { // in case 2 and 3, we have a starting triangle.
377  InputPoint M = _state.ptMax - _state.ptMin;
378  if ( ( M[ x ] == 0 ) && ( M[ y ] == 0 ) )
379  { // Case where M is aligned with axis, width is M[ z ].
380  result = M[ z ] * myWidth1 >= myWidth0 ? 1 : 2;
381  }
382  else
383  {
384  newCurrentTriangle( _state, M );
385  computeNormal( _state );
386  }
387  }
388  else
389  { // in case 1, normal is already aligned with main axis.
390  state.A = state.ptMax - state.ptMin;
391  }
392  }
393  if (loop>=1000) {
394  trace.warning() << "[ChordNaivePlaneComputer::satisfies()]"
395  << " more than 1000 loops, computing error suspected" << std::endl;
396  trace.warning() << "- Former state: " << *this << std::endl;
397  trace.warning() << "- Current state: ";
398  selfDisplay( trace.warning(), state );
399  trace.warning() << std::endl;
400  trace.warning() << "- Added Points: ";
401  for ( TInputIterator it = itB; it != itE; ++it )
402  trace.warning() << " " << *it;
403  trace.warning() << std::endl;
404  result = 1;
405  }
406  }
407  return result == 2;
408 }
409 //-----------------------------------------------------------------------------
410 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
411 template <typename TInputIterator>
412 inline
413 bool
414 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
415 satisfies( TInputIterator itB, TInputIterator itE ) const
416 {
417  return satisfies( _state, itB, itE );
418 }
419 
420 //-----------------------------------------------------------------------------
421 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
422 template <typename TInputIterator>
423 inline
424 std::pair<TInternalScalar, TInternalScalar>
425 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
426 computeAxisWidth( Dimension axis, TInputIterator itB, TInputIterator itE )
427 {
428  ChordNaivePlaneComputer plane;
429  plane.init( axis );
430  return plane.axisWidth( itB, itE );
431 }
432 
433 //-----------------------------------------------------------------------------
434 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
435 template <typename TInputIterator>
436 std::pair<TInternalScalar, TInternalScalar>
437 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
438 axisWidth( State & state, TInputIterator itB, TInputIterator itE ) const
439 {
440  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
441  ASSERT( isValid() );
442  state.nbValid = 0;
443  if ( itB == itE ) return std::make_pair( (TInternalScalar) 0, (TInternalScalar) 0 );
444  state.nbValid = findTriangle( state, itB, itE );
445  setUpNormal( state );
446  computeMinMax( state, itB, itE );
447  InternalVector formerN; // memorizes former normal vector.
448  formerN[ 0 ] = state.N[ 0 ]; formerN[ 1 ] = state.N[ 1 ]; formerN[ 2 ] = state.N[ 2 ];
449  std::pair<TInternalScalar, TInternalScalar> width( state.max - state.min, state.N[ z ] );
450  bool optimum = false;
451  for ( unsigned int loop = 0; ! optimum; ++loop )
452  {
453  computeHeight( state );
454  if ( state.nbValid >= 2 )
455  { // in case 2 and 3, we have a starting triangle.
456  InputPoint M = state.ptMax - state.ptMin;
457  if ( ( M[ x ] == 0 ) && ( M[ y ] == 0 ) )
458  { // Case where M is aligned with axis, width is M[ z ].
459  width.first = M[ z ];
460  width.second = 1;
461  optimum = true;
462  break;
463  }
464  else
465  {
466  newCurrentTriangle( state, M );
467  computeNormal( state );
468  }
469  }
470  else
471  { // in case 1, normal is already aligned with main axis.
472  state.A = state.ptMax - state.ptMin;
473  }
474  computeMinMax( state, itB, itE );
475  // To determine the optimum position:
476  // - it is not enough to look at the current width and check for the smallest one.
477  // - it is not enough to look at the current vector ptMax - ptMin to see if it is changing
478  // - it suffices to look at the normal vector and check that it has changed.
479  if ( ( formerN[ 0 ] == state.N[ 0 ] )
480  && ( formerN[ 1 ] == state.N[ 1 ] )
481  && ( formerN[ 2 ] == state.N[ 2 ] ) )
482  { // optimum
483  optimum = true;
484  }
485  else
486  { // better width found.
487  width.first = state.max - state.min;
488  width.second = state.N[ z ];
489  formerN[ 0 ] = state.N[ 0 ]; formerN[ 1 ] = state.N[ 1 ]; formerN[ 2 ] = state.N[ 2 ];
490  }
491  if (loop>=1000) {
492  trace.warning() << "[ChordNaivePlaneComputer::axisWidth()]"
493  << " more than 1000 loops, computing error suspected" << std::endl;
494  trace.warning() << "- Former state: " << *this << std::endl;
495  trace.warning() << "- Current state: ";
496  selfDisplay( trace.warning(), state );
497  trace.warning() << std::endl;
498  trace.warning() << "- Added Points: ";
499  for ( TInputIterator it = itB; it != itE; ++it )
500  trace.warning() << " " << *it;
501  trace.warning() << std::endl;
502  optimum = true;
503  }
504  }
505  return width;
506 }
507 //-----------------------------------------------------------------------------
508 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
509 template <typename TInputIterator>
510 inline
511 std::pair<TInternalScalar, TInternalScalar>
512 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
513 axisWidth( TInputIterator itB, TInputIterator itE ) const
514 {
515  return axisWidth( _state, itB, itE );
516 }
517 
518 
519 //-----------------------------------------------------------------------------
520 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
521 template <typename TInputIterator>
522 bool
523 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
524 extend( TInputIterator itB, TInputIterator itE )
525 {
526  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
527  if ( itB == itE ) return true;
528  unsigned int loop;
529  ASSERT( isValid() );
530  if ( empty() ) { // case where object is empty.
531  bool ok = satisfies( _state, itB, itE );
532  if ( ok ) // ideal case: all points false within current plane.
533  {
534  myState = _state;
535  myPointSet.insert( itB, itE );
536  // std::cout << "- (Group extend) satisfies for " << *this << std::endl;
537  }
538  return ok;
539  }
540  // Object contains already some points.
541  _state = myState;
542  bool changed = updateMinMax( _state, itB, itE );
543  if ( ! changed ) { // cases where added points did not change the bounds.
544  myPointSet.insert( itB, itE );
545  // std::cout << "- (Group extend) unchanged for " << *this << std::endl;
546  return true;
547  }
548  if ( myState.nbValid < 3 )
549  { // initial case
550  _state.nbValid = findMixedTriangle( _state,
551  myPointSet.begin(), myPointSet.end(),
552  itB, itE );
553  setUpNormal( _state );
554  // std::cout << "- (Group extend) findTriangle #=" << _state.nbValid << std::endl;
555  }
556  else
557  {
558  _state = myState;
559  // std::cout << "- (Group extend) Getting state " << std::endl;
560  }
561  unsigned int result = 0; // 0: computing, 1: too large, 2: found.
562  for ( loop = 0; result == 0; ++loop )
563  {
564  computeHeight( _state );
565  // std::cout << " - " << loop << "/h=" << _state.height
566  // << "/N=" << _state.N[ z ]
567  // << " : ";
568  // selfDisplay( std::cout, _state );
569  // std::cout << std::endl;
570  if ( ( _state.height * myWidth1 ) >= ( _state.N[ z ] * myWidth0 ) ) // because width is strict
571  result = 1;
572  else
573  {
574  computeMinMax( _state, myPointSet.begin(), myPointSet.end() );
575  updateMinMax( _state, itB, itE );
576  if ( checkWidth( _state ) )
577  { // ok, plane is fine.
578  result = 2;
579  }
580  else if ( _state.nbValid >= 2 )
581  { // in case 2 and 3, we have a starting triangle.
582  InputPoint M = _state.ptMax - _state.ptMin;
583  if ( ( M[ x ] == 0 ) && ( M[ y ] == 0 ) )
584  { // Case where M is aligned with axis, width is M[ z ].
585  result = M[ z ] * myWidth1 >= myWidth0 ? 1 : 2;
586  }
587  else
588  {
589  newCurrentTriangle( _state, M );
590  computeNormal( _state );
591  }
592  }
593  else
594  { // in case 1, normal is already aligned with main axis.
595  _state.A = _state.ptMax - _state.ptMin;
596  }
597  }
598  if (loop>=1000) {
599  trace.warning() << "[ChordNaivePlaneComputer::extend(Iterator,Iterator)]"
600  << " more than 1000 loops, computing error suspected" << std::endl;
601  trace.warning() << "- Former state: " << *this << std::endl;
602  trace.warning() << "- Current state: ";
603  selfDisplay( trace.warning(), _state );
604  trace.warning() << std::endl;
605  trace.warning() << "- Points: ";
606  for ( ConstIterator it = this->begin(), itEnd = this->end();
607  it != itEnd; ++it )
608  trace.warning() << " " << *it;
609  trace.warning() << std::endl;
610  trace.warning() << "- Added Points: ";
611  for ( TInputIterator it2 = itB; it2 != itE; ++it2 )
612  trace.warning() << " " << *it2;
613  trace.warning() << std::endl;
614  result = 1;
615  }
616  }
617  // find a correct plane.
618  if ( result == 2 )
619  {
620  myState = _state;
621  myPointSet.insert( itB, itE );
622  return true;
623  }
624  // was unable to find a correct plane.
625  return false;
626 }
627 //-----------------------------------------------------------------------------
628 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
629 template <typename TInputIterator>
630 bool
631 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
632 isExtendable( TInputIterator itB, TInputIterator itE ) const
633 {
634  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
635  if ( itB == itE ) return true;
636  unsigned int loop;
637  ASSERT( isValid() );
638  if ( empty() )
639  { // case where object is empty.
640  return satisfies( _state, itB, itE );
641  }
642  // Object contains already some points.
643  _state = myState;
644  bool changed = updateMinMax( _state, itB, itE );
645  if ( ! changed )
646  { // cases where added points did not change the bounds.
647  return true;
648  }
649  if ( myState.nbValid < 3 )
650  { // initial case
651  _state.nbValid = findMixedTriangle( _state,
652  myPointSet.begin(), myPointSet.end(),
653  itB, itE );
654  setUpNormal( _state );
655  }
656  else
657  {
658  _state = myState;
659  }
660  unsigned int result = 0; // 0: computing, 1: too large, 2: found.
661  for ( loop = 0; result == 0; ++loop )
662  {
663  computeHeight( _state );
664  if ( ( _state.height * myWidth1 ) >= ( _state.N[ z ] * myWidth0 ) ) // because width is strict
665  result = 1;
666  else
667  {
668  computeMinMax( _state, myPointSet.begin(), myPointSet.end() );
669  updateMinMax( _state, itB, itE );
670  if ( checkWidth( _state ) )
671  { // ok, plane is fine.
672  result = 2;
673  }
674  else if ( _state.nbValid >= 2 )
675  { // in case 2 and 3, we have a starting triangle.
676  InputPoint M = _state.ptMax - _state.ptMin;
677  if ( ( M[ x ] == 0 ) && ( M[ y ] == 0 ) )
678  { // Case where M is aligned with axis, width is M[ z ].
679  result = M[ z ] * myWidth1 >= myWidth0 ? 1 : 2;
680  }
681  else
682  {
683  newCurrentTriangle( _state, M );
684  computeNormal( _state );
685  }
686  }
687  else
688  { // in case 1, normal is already aligned with main axis.
689  _state.A = _state.ptMax - _state.ptMin;
690  }
691  }
692  if (loop>=1000) {
693  trace.warning() << "[ChordNaivePlaneComputer::isExtendable(Iterator,Iterator)]"
694  << " more than 1000 loops, computing error suspected" << std::endl;
695  trace.warning() << "- Former state: " << *this << std::endl;
696  trace.warning() << "- Current state: ";
697  selfDisplay( trace.warning(), _state );
698  trace.warning() << std::endl;
699  trace.warning() << "- Points: ";
700  for ( ConstIterator it = this->begin(), itEnd = this->end();
701  it != itEnd; ++it )
702  trace.warning() << " " << *it;
703  trace.warning() << std::endl;
704  trace.warning() << "- Added Points: ";
705  for ( TInputIterator it2 = itB; it2 != itE; ++it2 )
706  trace.warning() << " " << *it2;
707  trace.warning() << std::endl;
708  result = 1;
709  }
710  }
711  // find a correct plane.
712  return result == 2;
713 }
714 
715 //-----------------------------------------------------------------------------
716 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
717 inline
718 typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::Primitive
719 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
720 primitive() const
721 {
722  typedef typename Space::RealVector RealVector;
723  typedef typename RealVector::Component Scalar;
724  RealVector N;
725  for ( Dimension i = 0; i < 3; ++i )
726  N[ i ] = NumberTraits<InternalScalar>::castToDouble( myState.N[ i ] );
727  Scalar min = NumberTraits<InternalScalar>::castToDouble( myState.min );
728  Scalar max = NumberTraits<InternalScalar>::castToDouble( myState.max );
729  return Primitive( min, N, max - min );
730 }
731 
732 //-----------------------------------------------------------------------------
733 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
734 template <typename Vector3D>
735 inline
736 void
737 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
738 getNormal( Vector3D & normal ) const
739 {
740  for ( Dimension i = 0; i < 3; ++i )
741  normal[ i ] = NumberTraits<InternalScalar>::castToDouble( myState.N[ i ] );
742 }
743 //-----------------------------------------------------------------------------
744 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
745 inline
746 const typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::InternalVector&
747 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
748 exactNormal() const
749 {
750  return myState.N;
751 }
752 //-----------------------------------------------------------------------------
753 //-----------------------------------------------------------------------------
754 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
755 template <typename Vector3D>
756 inline
757 void
758 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
759 getUnitNormal( Vector3D & normal ) const
760 {
761  getNormal( normal );
762  double l = sqrt( normal[ 0 ] * normal[ 0 ]
763  + normal[ 1 ] * normal[ 1 ]
764  + normal[ 2 ] * normal[ 2 ] );
765  normal[ 0 ] /= l;
766  normal[ 1 ] /= l;
767  normal[ 2 ] /= l;
768 }
769 //-----------------------------------------------------------------------------
770 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
771 inline
772 void
773 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
774 getBounds( double & min, double & max ) const
775 {
776  double nx = NumberTraits<InternalScalar>::castToDouble( myState.N[ 0 ] );
777  double ny = NumberTraits<InternalScalar>::castToDouble( myState.N[ 1 ] );
778  double nz = NumberTraits<InternalScalar>::castToDouble( myState.N[ 2 ] );
779  double l = sqrt( nx*nx + ny*ny + nz*nz );
780  min = NumberTraits<InternalScalar>::castToDouble( myState.min ) / l;
781  max = NumberTraits<InternalScalar>::castToDouble( myState.max ) / l;
782 }
783 //-----------------------------------------------------------------------------
784 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
785 inline
786 const typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::InputPoint &
787 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
788 minimalPoint() const
789 {
790  ASSERT( ! this->empty() );
791  return myState.ptMin;
792 }
793 //-----------------------------------------------------------------------------
794 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
795 inline
796 const typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::InputPoint &
797 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
798 maximalPoint() const
799 {
800  ASSERT( ! this->empty() );
801  return myState.ptMax;
802 }
803 
804 //-----------------------------------------------------------------------------
805 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
806 template <typename TVector1, typename TVector2>
807 inline
808 typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::InternalScalar
809 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
810 internalDot( const TVector1 & u, const TVector2 & v )
811 {
812  return (InternalScalar) u[ 0 ] * (InternalScalar) v[ 0 ]
813  + (InternalScalar) u[ 1 ] * (InternalScalar) v[ 1 ]
814  + (InternalScalar) u[ 2 ] * (InternalScalar) v[ 2 ];
815 }
816 //-----------------------------------------------------------------------------
817 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
818 template <typename TVector1, typename TVector2>
819 inline
820 void
821 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
822 internalCross( InternalVector & n, const TVector1 & u, const TVector2 & v )
823 {
824  n[ 0 ] = (InternalScalar) u[ 1 ] * (InternalScalar) v[ 2 ] - (InternalScalar) u[ 2 ] * (InternalScalar) v[ 1 ];
825  n[ 1 ] = (InternalScalar) u[ 2 ] * (InternalScalar) v[ 0 ] - (InternalScalar) u[ 0 ] * (InternalScalar) v[ 2 ];
826  n[ 2 ] = (InternalScalar) u[ 0 ] * (InternalScalar) v[ 1 ] - (InternalScalar) u[ 1 ] * (InternalScalar) v[ 0 ];
827 }
828 
829 /**
830  * Writes/Displays the object on an output stream.
831  * @param out the output stream where the object is written.
832  */
833 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
834 inline
835 void
836 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::selfDisplay ( std::ostream & out ) const
837 {
838  selfDisplay( out, myState );
839 }
840 /**
841  * Writes/Displays the object on an output stream.
842  * @param out the output stream where the object is written.
843  */
844 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
845 inline
846 void
847 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::selfDisplay ( std::ostream & out, const State & state ) const
848 {
849  double min, max;
850  double N[ 3 ];
851  out << "[ChordNaivePlaneComputer"
852  << " axis=" << z << " w=" << myWidth0 << "/" << myWidth1
853  << " size=" << size()
854  << " N=(" << state.N[ 0 ] << "," << state.N[ 1 ] << "," << state.N[ 2 ] << ")"
855  << " A=(" << state.A[ 0 ] << "," << state.A[ 1 ] << "," << state.A[ 2 ] << ")"
856  << " B=(" << state.B[ 0 ] << "," << state.B[ 1 ] << "," << state.B[ 2 ] << ")"
857  << " C=(" << state.C[ 0 ] << "," << state.C[ 1 ] << "," << state.C[ 2 ] << ")"
858  << " #=" << state.nbValid
859  << " min=" << state.min
860  << " max=" << state.max
861  << ": ";
862  this->getUnitNormal( N );
863  this->getBounds( min, max );
864  out << min << " <= "
865  << N[ 0 ] << " * x + "
866  << N[ 1 ] << " * y + "
867  << N[ 2 ] << " * z "
868  << " <= " << max << " ]";
869 }
870 
871 /**
872  * Checks the validity/consistency of the object.
873  * @return 'true' if the object is valid, 'false' otherwise.
874  */
875 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
876 inline
877 bool
878 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::isValid() const
879 {
880  return z < 3;
881 }
882 
883 
884 ///////////////////////////////////////////////////////////////////////////////
885 // Internals
886 ///////////////////////////////////////////////////////////////////////////////
887 //-----------------------------------------------------------------------------
888 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
889 void
890 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
891 setUpNormal( State & state ) const
892 {
893  switch ( state.nbValid ) {
894  case 1: setUpNormal1( state ); break;
895  case 2: setUpNormal2( state ); break;
896  case 3: setUpNormal3( state ); break;
897  }
898 }
899 //-----------------------------------------------------------------------------
900 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
901 void
902 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
903 setUpNormal1( State & state ) const
904 {
905  state.A.reset();
906  state.B.reset();
907  state.C.reset();
908  state.N[ x ] = state.N[ y ] = NumberTraits<InternalScalar>::ZERO;
909  state.N[ z ] = NumberTraits<InternalScalar>::ONE;
910 }
911 //-----------------------------------------------------------------------------
912 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
913 void
914 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
915 setUpNormal2( State & state ) const
916 { // In this case, we create an imaginary point in the plane 0xy, orthogonal to direction AB.
917  state.N[ x ] = state.N[ y ] = NumberTraits<InternalScalar>::ZERO;
918  state.N[ z ] = NumberTraits<InternalScalar>::ONE;
919  state.A = state.B - state.A;
920  // 3 next lines compute internalCross( state.C, state.A, state.N );
921  state.C[ x ] = state.A[ y ];
922  state.C[ y ] = -state.A[ x ];
923  state.C[ z ] = 0;
924  if ( ( state.C[ x ] != 0 ) || ( state.C[ y ] != 0 ) )
925  // zero iff _state.A is aligned with main axis
926  internalCross( state.N, state.C, state.A );
927  state.B = -state.A - state.C;
928  if ( state.N[ z ] < 0 ) {
929  state.N[ 0 ] = -state.N[ 0 ];
930  state.N[ 1 ] = -state.N[ 1 ];
931  state.N[ 2 ] = -state.N[ 2 ];
932  InputPoint M = state.A;
933  state.A = state.B;
934  state.B = M;
935  }
936 }
937 //-----------------------------------------------------------------------------
938 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
939 inline
940 void
941 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
942 setUpNormal3( State & state ) const
943 {
944  state.A = state.B - state.A; // p2 - p1;
945  state.B = state.C - state.B; // p3 - p2;
946  state.C = -state.A - state.B;// p1 - p3;
947  computeNormal( state );
948 }
949 //-----------------------------------------------------------------------------
950 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
951 inline
952 bool
953 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
954 setUp1( const InputPoint & p1 )
955 {
956  myState.nbValid = 1;
957  myState.height = NumberTraits<InternalScalar>::ZERO;
958  setUpNormal1( myState );
959  myState.ptMax = p1;
960  myState.ptMin = p1;
961  myState.max = (InternalScalar) p1[ z ];
962  myState.min = myState.max;
963  return true;
964 }
965 //-----------------------------------------------------------------------------
966 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
967 inline
968 int
969 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
970 signDelta( const InputPoint & A, const InputPoint & B, const InputPoint & C ) const
971 {
972  InternalScalar res =
973  ( (InternalScalar) (B[ x ]-A[ x ]) ) * ( (InternalScalar) (C[ y ]-A[ y ]) )
974  - ( (InternalScalar) (B[ y ]-A[ y ]) ) * ( (InternalScalar) (C[ x ]-A[ x ]) );
975  return ( res > (InternalScalar) 0 )
976  ? 1 : ( ( res < (InternalScalar) 0 ) ? -1 : 0 );
977 }
978 //-----------------------------------------------------------------------------
979 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
980 inline
981 int
982 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
983 signDelta( const InputPoint & A, const InputPoint & C ) const
984 { // B is zero
985  InternalScalar res =
986  ( (InternalScalar) (-A[ x ]) ) * ( (InternalScalar) (C[ y ]-A[ y ]) )
987  + ( (InternalScalar) (A[ y ]) ) * ( (InternalScalar) (C[ x ]-A[ x ]) );
988  return ( res > (InternalScalar) 0 )
989  ? 1 : ( ( res < (InternalScalar) 0 ) ? -1 : 0 );
990 }
991 
992 //-----------------------------------------------------------------------------
993 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
994 template <typename TInputIterator>
995 bool
996 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
997 updateMinMax( State & state, TInputIterator itB, TInputIterator itE ) const
998 {
999  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
1000  bool changed = false;
1001  for ( ; itB != itE; ++itB )
1002  {
1003  _d = internalDot( state.N, *itB );
1004  if ( _d > state.max ) {
1005  state.max = _d;
1006  state.ptMax = *itB;
1007  changed = true;
1008  } else if ( _d < state.min ) {
1009  state.min = _d;
1010  state.ptMin = *itB;
1011  changed = true;
1012  }
1013  }
1014  return changed;
1015 }
1016 //-----------------------------------------------------------------------------
1017 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1018 template <typename TInputIterator>
1019 void
1020 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1021 computeMinMax( State & state, TInputIterator itB, TInputIterator itE ) const
1022 {
1023  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
1024  ASSERT( itB != itE );
1025  _d = internalDot( state.N, *itB );
1026  state.min = state.max = _d;
1027  state.ptMin = state.ptMax = *itB;
1028  ++itB;
1029  for ( ; itB != itE; ++itB )
1030  {
1031  _d = internalDot( state.N, *itB );
1032  if ( _d > state.max ) {
1033  state.max = _d;
1034  state.ptMax = *itB;
1035  } else if ( _d < state.min ) {
1036  state.min = _d;
1037  state.ptMin = *itB;
1038  }
1039  }
1040 }
1041 //-----------------------------------------------------------------------------
1042 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1043 template <typename TInputIterator>
1044 unsigned int
1045 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1046 findTriangle( State & state, TInputIterator itB, TInputIterator itE ) const
1047 {
1048  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
1049  if ( itB == itE ) return 0;
1050  state.A = *itB++; // first point
1051  for ( ; ( itB != itE ) && alignedAlongAxis( state.A, *itB ); ++itB )
1052  ;
1053  if ( itB == itE ) return 1;
1054  state.B = *itB++; // second point
1055  for ( ; ( itB != itE ) && ( signDelta( state.A, state.B, *itB ) == 0 ); ++itB )
1056  ;
1057  if ( itB == itE ) return 2;
1058  state.C = *itB; // third point
1059  return 3;
1060 }
1061 //-----------------------------------------------------------------------------
1062 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1063 template <typename TInputIterator>
1064 unsigned int
1065 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1066 findTriangle1( State & state, TInputIterator itB, TInputIterator itE ) const
1067 {
1068  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
1069  for ( ; ( itB != itE ) && alignedAlongAxis( state.A, *itB ); ++itB )
1070  ;
1071  if ( itB == itE ) return 1;
1072  state.B = *itB++; // second point
1073  for ( ; ( itB != itE ) && ( signDelta( state.A, state.B, *itB ) == 0 ); ++itB )
1074  ;
1075  if ( itB == itE ) return 2;
1076  state.C = *itB; // third point
1077  return 3;
1078 }
1079 //-----------------------------------------------------------------------------
1080 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1081 template <typename TInputIterator>
1082 unsigned int
1083 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1084 findTriangle2( State & state, TInputIterator itB, TInputIterator itE ) const
1085 {
1086  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator> ));
1087  for ( ; ( itB != itE ) && ( signDelta( state.A, state.B, *itB ) == 0 ); ++itB )
1088  ;
1089  if ( itB == itE ) return 2;
1090  state.C = *itB; // third point
1091  return 3;
1092 }
1093 //-----------------------------------------------------------------------------
1094 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1095 template <typename TInputIterator1, typename TInputIterator2>
1096 unsigned int
1097 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1098 findMixedTriangle( State & state,
1099  TInputIterator1 itB1, TInputIterator1 itE1,
1100  TInputIterator2 itB2, TInputIterator2 itE2 ) const
1101 {
1102  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator1> ));
1103  BOOST_CONCEPT_ASSERT(( boost::InputIterator<TInputIterator2> ));
1104  if ( itB1 == itE1 ) return findTriangle( state, itB2, itE2 ) ;
1105  state.A = *itB1++; // first point
1106  for ( ; ( itB1 != itE1 ) && alignedAlongAxis( state.A, *itB1 ); ++itB1 )
1107  ;
1108  if ( itB1 == itE1 ) return findTriangle1( state, itB2, itE2 );
1109 
1110  state.B = *itB1++; // second point
1111  for ( ; ( itB1 != itE1 ) && ( signDelta( state.A, state.B, *itB1 ) == 0 ); ++itB1 )
1112  ;
1113  if ( itB1 == itE1 ) return findTriangle2( state, itB2, itE2 );
1114  state.C = *itB1; // third point
1115  return 3;
1116 }
1117 //-----------------------------------------------------------------------------
1118 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1119 inline
1120 bool
1121 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1122 alignedAlongAxis( const InputPoint & p1, const InputPoint & p2 ) const
1123 {
1124  return ( p1[ x ] == p2[ x ] ) && ( p1[ y ] == p2[ y ] );
1125 }
1126 //-----------------------------------------------------------------------------
1127 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1128 inline
1129 void
1130 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1131 computeHeight( State & state ) const
1132 {
1133  state.height = internalDot( state.A, state.N );
1134 }
1135 //-----------------------------------------------------------------------------
1136 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1137 inline
1138 void
1139 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1140 computeNormal( State & state ) const
1141 {
1142  ASSERT( state.nbValid >= 2 );
1143  internalCross( state.N, state.B - state.A, state.C - state.A );
1144  if ( state.N[ z ] < 0 ) {
1145  state.N[ 0 ] = -state.N[ 0 ];
1146  state.N[ 1 ] = -state.N[ 1 ];
1147  state.N[ 2 ] = -state.N[ 2 ];
1148  InputPoint M = state.A;
1149  state.A = state.B;
1150  state.B = M;
1151  }
1152 }
1153 //-----------------------------------------------------------------------------
1154 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1155 inline
1156 bool
1157 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1158 checkWidth( const State & state ) const
1159 {
1160  return ( state.max - state.min ) <= state.height;
1161 }
1162 //-----------------------------------------------------------------------------
1163 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1164 bool
1165 DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1166 newCurrentTriangle( State & state, const InputPoint & M ) const
1167 {
1168  ASSERT( state.nbValid >= 2 );
1169  int da = signDelta( state.A, M ); // A, O, M
1170  int db = signDelta( state.B, M ); // B, O, M
1171  int dc = signDelta( state.C, M ); // C, O, M
1172  if ( ( da >= 0 ) && ( db <= 0 ) ) {
1173  if ( ( signDelta( state.A, state.B ) == 0 ) && ( da == 0 ) ) {
1174  if ( ( state.A[ x ] * M[ x ] > 0 ) || ( state.A[ y ] * M[ y ] > 0 ) )
1175  state.A = M;
1176  else
1177  state.B = M;
1178  } else
1179  state.C = M;
1180  } else if ( ( db >= 0 ) && ( dc <= 0 ) ) {
1181  if ( ( signDelta( state.B, state.C ) ==0 ) && ( db == 0 ) ) {
1182  if ( ( state.B[ x ] * M[ x ] > 0 ) || ( state.B[ y ] * M[ y ] > 0 ) )
1183  state.B = M;
1184  else
1185  state.C = M;
1186  } else
1187  state.A = M;
1188  } else if ( ( dc >= 0 ) && ( da <= 0 ) ) {
1189  if ( ( signDelta( state.C, state.A ) ==0 ) && ( dc == 0 ) ) {
1190  if ( ( state.C[ x ] * M[ x ] > 0 ) || ( state.C[ y ] * M[ y ] > 0 ) )
1191  state.C = M;
1192  else
1193  state.A = M;
1194  } else
1195  state.B = M;
1196  } else {
1197  trace.warning() << "[ChordNaivePlaneComputer::newCurrentTriangle]"
1198  << " cannot find a triangle cutting the main axis" << std::endl;
1199  return false;
1200  }
1201  return true;
1202 }
1203 
1204 ///////////////////////////////////////////////////////////////////////////////
1205 // Implementation of inline functions //
1206 
1207 template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1208 inline
1209 std::ostream&
1210 DGtal::operator<< ( std::ostream & out,
1211  const ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar> & object )
1212 {
1213  object.selfDisplay( out );
1214  return out;
1215 }
1216 
1217 // //
1218 ///////////////////////////////////////////////////////////////////////////////