DGtal 1.3.0
Loading...
Searching...
No Matches
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//-----------------------------------------------------------------------------
45template <typename TSpace, typename TInputPoint, typename TInternalScalar>
46inline
47DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
48~ChordNaivePlaneComputer()
49{ // Nothing to do.
50}
51//-----------------------------------------------------------------------------
52template <typename TSpace, typename TInputPoint, typename TInternalScalar>
53inline
54DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
55ChordNaivePlaneComputer() : z( -1 )
56{ // Object is invalid
57}
58//-----------------------------------------------------------------------------
59template <typename TSpace, typename TInputPoint, typename TInternalScalar>
60inline
61DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
62ChordNaivePlaneComputer( 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//-----------------------------------------------------------------------------
71template <typename TSpace, typename TInputPoint, typename TInternalScalar>
72inline
73DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar> &
74DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
75operator=( 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//-----------------------------------------------------------------------------
90template <typename TSpace, typename TInputPoint, typename TInternalScalar>
91inline
92void
93DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
94clear()
95{
96 myPointSet.clear();
97 myState.nbValid = 0;
98}
99//-----------------------------------------------------------------------------
100template <typename TSpace, typename TInputPoint, typename TInternalScalar>
101inline
102void
103DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
104init( 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//-----------------------------------------------------------------------------
119template <typename TSpace, typename TInputPoint, typename TInternalScalar>
120inline
121typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::Size
122DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
123size() const
124{
125 return myPointSet.size();
126}
127//-----------------------------------------------------------------------------
128template <typename TSpace, typename TInputPoint, typename TInternalScalar>
129inline
130bool
131DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
132empty() const
133{
134 return myPointSet.empty();
135}
136//-----------------------------------------------------------------------------
137template <typename TSpace, typename TInputPoint, typename TInternalScalar>
138inline
139typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::ConstIterator
140DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
141begin() const
142{
143 return myPointSet.begin();
144}
145//-----------------------------------------------------------------------------
146template <typename TSpace, typename TInputPoint, typename TInternalScalar>
147inline
148typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::ConstIterator
149DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
150end() const
151{
152 return myPointSet.end();
153}
154//-----------------------------------------------------------------------------
155template <typename TSpace, typename TInputPoint, typename TInternalScalar>
156inline
157typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::Size
158DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
159max_size() const
160{
161 return myPointSet.max_size();
162}
163//-----------------------------------------------------------------------------
164template <typename TSpace, typename TInputPoint, typename TInternalScalar>
165inline
166typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::Size
167DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
168maxSize() const
169{
170 return max_size();
171}
172//-----------------------------------------------------------------------------
173template <typename TSpace, typename TInputPoint, typename TInternalScalar>
174inline
175bool
176DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
177operator()( const Point & p ) const
178{
179 _d = internalDot( myState.N, p );
180 return ( _d >= myState.min ) && ( _d <= myState.max );
181}
182//-----------------------------------------------------------------------------
183template <typename TSpace, typename TInputPoint, typename TInternalScalar>
184inline
185bool
186DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
187extendAsIs( 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//-----------------------------------------------------------------------------
199template <typename TSpace, typename TInputPoint, typename TInternalScalar>
200bool
201DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
202extend( 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//-----------------------------------------------------------------------------
279template <typename TSpace, typename TInputPoint, typename TInternalScalar>
280bool
281DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
282isExtendable( 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//-----------------------------------------------------------------------------
351template <typename TSpace, typename TInputPoint, typename TInternalScalar>
352template <typename TInputIterator>
353bool
354DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
355satisfies( 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//-----------------------------------------------------------------------------
410template <typename TSpace, typename TInputPoint, typename TInternalScalar>
411template <typename TInputIterator>
412inline
413bool
414DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
415satisfies( TInputIterator itB, TInputIterator itE ) const
416{
417 return satisfies( _state, itB, itE );
418}
419
420//-----------------------------------------------------------------------------
421template <typename TSpace, typename TInputPoint, typename TInternalScalar>
422template <typename TInputIterator>
423inline
424std::pair<TInternalScalar, TInternalScalar>
425DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
426computeAxisWidth( Dimension axis, TInputIterator itB, TInputIterator itE )
427{
428 ChordNaivePlaneComputer plane;
429 plane.init( axis );
430 return plane.axisWidth( itB, itE );
431}
432
433//-----------------------------------------------------------------------------
434template <typename TSpace, typename TInputPoint, typename TInternalScalar>
435template <typename TInputIterator>
436std::pair<TInternalScalar, TInternalScalar>
437DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
438axisWidth( 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//-----------------------------------------------------------------------------
508template <typename TSpace, typename TInputPoint, typename TInternalScalar>
509template <typename TInputIterator>
510inline
511std::pair<TInternalScalar, TInternalScalar>
512DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
513axisWidth( TInputIterator itB, TInputIterator itE ) const
514{
515 return axisWidth( _state, itB, itE );
516}
517
518
519//-----------------------------------------------------------------------------
520template <typename TSpace, typename TInputPoint, typename TInternalScalar>
521template <typename TInputIterator>
522bool
523DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
524extend( 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//-----------------------------------------------------------------------------
628template <typename TSpace, typename TInputPoint, typename TInternalScalar>
629template <typename TInputIterator>
630bool
631DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
632isExtendable( 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//-----------------------------------------------------------------------------
716template <typename TSpace, typename TInputPoint, typename TInternalScalar>
717inline
718typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::Primitive
719DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
720primitive() 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//-----------------------------------------------------------------------------
733template <typename TSpace, typename TInputPoint, typename TInternalScalar>
734template <typename Vector3D>
735inline
736void
737DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
738getNormal( Vector3D & normal ) const
739{
740 for ( Dimension i = 0; i < 3; ++i )
741 normal[ i ] = NumberTraits<InternalScalar>::castToDouble( myState.N[ i ] );
742}
743//-----------------------------------------------------------------------------
744template <typename TSpace, typename TInputPoint, typename TInternalScalar>
745inline
746const typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::InternalVector&
747DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
748exactNormal() const
749{
750 return myState.N;
751}
752//-----------------------------------------------------------------------------
753//-----------------------------------------------------------------------------
754template <typename TSpace, typename TInputPoint, typename TInternalScalar>
755template <typename Vector3D>
756inline
757void
758DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
759getUnitNormal( 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//-----------------------------------------------------------------------------
770template <typename TSpace, typename TInputPoint, typename TInternalScalar>
771inline
772void
773DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
774getBounds( 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//-----------------------------------------------------------------------------
784template <typename TSpace, typename TInputPoint, typename TInternalScalar>
785inline
786const typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::InputPoint &
787DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
788minimalPoint() const
789{
790 ASSERT( ! this->empty() );
791 return myState.ptMin;
792}
793//-----------------------------------------------------------------------------
794template <typename TSpace, typename TInputPoint, typename TInternalScalar>
795inline
796const typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::InputPoint &
797DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
798maximalPoint() const
799{
800 ASSERT( ! this->empty() );
801 return myState.ptMax;
802}
803
804//-----------------------------------------------------------------------------
805template <typename TSpace, typename TInputPoint, typename TInternalScalar>
806template <typename TVector1, typename TVector2>
807inline
808typename DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::InternalScalar
809DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
810internalDot( 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//-----------------------------------------------------------------------------
817template <typename TSpace, typename TInputPoint, typename TInternalScalar>
818template <typename TVector1, typename TVector2>
819inline
820void
821DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
822internalCross( 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 */
833template <typename TSpace, typename TInputPoint, typename TInternalScalar>
834inline
835void
836DGtal::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 */
844template <typename TSpace, typename TInputPoint, typename TInternalScalar>
845inline
846void
847DGtal::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 */
875template <typename TSpace, typename TInputPoint, typename TInternalScalar>
876inline
877bool
878DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::isValid() const
879{
880 return z < 3;
881}
882
883
884///////////////////////////////////////////////////////////////////////////////
885// Internals
886///////////////////////////////////////////////////////////////////////////////
887//-----------------------------------------------------------------------------
888template <typename TSpace, typename TInputPoint, typename TInternalScalar>
889void
890DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
891setUpNormal( 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//-----------------------------------------------------------------------------
900template <typename TSpace, typename TInputPoint, typename TInternalScalar>
901void
902DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
903setUpNormal1( 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//-----------------------------------------------------------------------------
912template <typename TSpace, typename TInputPoint, typename TInternalScalar>
913void
914DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
915setUpNormal2( 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//-----------------------------------------------------------------------------
938template <typename TSpace, typename TInputPoint, typename TInternalScalar>
939inline
940void
941DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
942setUpNormal3( 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//-----------------------------------------------------------------------------
950template <typename TSpace, typename TInputPoint, typename TInternalScalar>
951inline
952bool
953DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
954setUp1( 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//-----------------------------------------------------------------------------
966template <typename TSpace, typename TInputPoint, typename TInternalScalar>
967inline
968int
969DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
970signDelta( 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//-----------------------------------------------------------------------------
979template <typename TSpace, typename TInputPoint, typename TInternalScalar>
980inline
981int
982DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
983signDelta( 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//-----------------------------------------------------------------------------
993template <typename TSpace, typename TInputPoint, typename TInternalScalar>
994template <typename TInputIterator>
995bool
996DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
997updateMinMax( 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//-----------------------------------------------------------------------------
1017template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1018template <typename TInputIterator>
1019void
1020DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1021computeMinMax( 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//-----------------------------------------------------------------------------
1042template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1043template <typename TInputIterator>
1044unsigned int
1045DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1046findTriangle( 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//-----------------------------------------------------------------------------
1062template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1063template <typename TInputIterator>
1064unsigned int
1065DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1066findTriangle1( 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//-----------------------------------------------------------------------------
1080template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1081template <typename TInputIterator>
1082unsigned int
1083DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1084findTriangle2( 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//-----------------------------------------------------------------------------
1094template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1095template <typename TInputIterator1, typename TInputIterator2>
1096unsigned int
1097DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1098findMixedTriangle( 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//-----------------------------------------------------------------------------
1118template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1119inline
1120bool
1121DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1122alignedAlongAxis( const InputPoint & p1, const InputPoint & p2 ) const
1123{
1124 return ( p1[ x ] == p2[ x ] ) && ( p1[ y ] == p2[ y ] );
1125}
1126//-----------------------------------------------------------------------------
1127template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1128inline
1129void
1130DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1131computeHeight( State & state ) const
1132{
1133 state.height = internalDot( state.A, state.N );
1134}
1135//-----------------------------------------------------------------------------
1136template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1137inline
1138void
1139DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1140computeNormal( 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//-----------------------------------------------------------------------------
1154template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1155inline
1156bool
1157DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1158checkWidth( const State & state ) const
1159{
1160 return ( state.max - state.min ) <= state.height;
1161}
1162//-----------------------------------------------------------------------------
1163template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1164bool
1165DGtal::ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar>::
1166newCurrentTriangle( 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
1207template <typename TSpace, typename TInputPoint, typename TInternalScalar>
1208inline
1209std::ostream&
1210DGtal::operator<< ( std::ostream & out,
1211 const ChordNaivePlaneComputer<TSpace, TInputPoint, TInternalScalar> & object )
1212{
1213 object.selfDisplay( out );
1214 return out;
1215}
1216
1217// //
1218///////////////////////////////////////////////////////////////////////////////
Generated on Sun Nov 27 2022 15:15:37 for DGtal by doxygen 1.9.5