OpendTect  6.6
polygon.h
Go to the documentation of this file.
1 #pragma once
2 
3 /*+
4 ________________________________________________________________________
5 
6  (C) dGB Beheer B.V.; (LICENSE) http://opendtect.org/OpendTect_license.txt
7  Author: J.C. Glas
8  Date: Dec 2006
9  RCS: $Id$
10 ________________________________________________________________________
11 
12 -*/
13 
14 #include "coord.h"
15 #include "typeset.h"
16 #include "iopar.h"
17 #include "bufstring.h"
18 #include "bendpointfinder.h"
19 #include <math.h>
20 
25 template <class T>
26 mClass(Algo) ODPolygon
27 {
28 public:
29  ODPolygon()
30  : closed_(true), udf_(Geom::Point2D<T>::udf())
31  , xrg_(mUdf(T),mUdf(T)), yrg_(mUdf(T),mUdf(T)) {}
32 
33  ODPolygon(const TypeSet<Geom::Point2D<T> >& plg)
34  : poly_(plg), closed_(true)
35  , udf_(Geom::Point2D<T>::udf())
36  , xrg_(mUdf(T),mUdf(T)), yrg_(mUdf(T),mUdf(T)) {}
37 
38  void erase();
39  bool isEmpty() const { return poly_.isEmpty(); }
40 
41  int size() const { return poly_.size(); }
42  bool validIdx(int idx) const { return poly_.validIdx(idx); }
43  void setEmpty() { erase(); }
44 
45  void add(const Geom::Point2D<T>& vtx);
46  void remove(int idx);
47  void insert(int idx,const Geom::Point2D<T>& vtx);
48 
49  bool isInside(const Geom::Point2D<T>&,
50  bool inclborder,T eps) const;
51 
52  int isInside(const ODPolygon& testpoly,T eps=0) const;
53  /* 0: testpoly fully outside (borders don't touch)
54  2: testpoly fully inside (borders don't touch)
55  1: all intermediate cases */
56 
57  bool segmentMeetsBorder(const Geom::Point2D<T>& pt1,
58  const Geom::Point2D<T>& pt2,T eps) const;
59  bool windowOverlaps(const Interval<T>& xrange,
60  const Interval<T>& yrange,T eps) const;
61 
62  // defined for closed polygon
63  const Geom::Point2D<T>& getVertex(int idx) const;
64  const Geom::Point2D<T>& nextVertex(int idx) const;
65  const Geom::Point2D<T>& prevVertex(int idx) const;
66 
67  void setClosed( bool yn ) { closed_ = yn; }
68  bool isClosed() const { return closed_; }
69  void setUdf( Geom::Point2D<T> pt ) { udf_ = pt; }
70  Geom::Point2D<T> getUdf() const { return udf_; }
71  const TypeSet<Geom::Point2D<T> >& data() const { return poly_; }
72 
73  Interval<T> getRange(bool of_x) const;
74  void getData(bool of_x,TypeSet<T>&) const;
75 
76  void removeZeroLengths();
77  bool isUTurn(int idx) const;
78  bool isSelfIntersecting() const;
79 
80  void convexHull();
81  void keepBendPoints(float eps);
82 
83  // not for self-intersecting polygons
84  double area() const { return fabs(sgnArea2<double>()/2.0); }
85  bool clockwise() const { return sgnArea2<T>()<0; }
86  bool anticlockwise() const { return sgnArea2<T>()>0; }
87 
88  void reverse();
89 
90  double distTo(const Geom::Point2D<T>& refpt,int* segmentidxptr=0,
91  double* fractionptr=0) const;
92 
93  double maxDistToBorderEstimate(double maxrelerr=0.001) const;
94 
95  bool operator==(const ODPolygon<T>&) const;
96  bool operator>( const ODPolygon<T>& plg ) const
97  { return poly_.size()>plg.size(); }
98 
99 protected:
100 
101  static int doSegmentsMeet( const Geom::Point2D<T>& p1,
102  const Geom::Point2D<T>& p2,
103  const Geom::Point2D<T>& q1,
104  const Geom::Point2D<T>& q2,
105  T eps );
106 
107  static bool isOnSegment( const Geom::Point2D<T>& pt,
108  const Geom::Point2D<T>& pt0,
109  const Geom::Point2D<T>& pt1,
110  T eps );
111  static bool isOnHalfLine( const Geom::Point2D<T>& point,
112  const Geom::Point2D<T>& dirvec,
113  const Geom::Point2D<T>& endvec,
114  T eps );
115 
116  static bool isEdgeCrossing( const Geom::Point2D<T>& raydir,
117  const Geom::Point2D<T>& raysrc,
118  const Geom::Point2D<T>& vtx1,
119  const Geom::Point2D<T>& vtx2 );
120 
121  static bool isOnLine( const Geom::Point2D<T>& point,
122  const Geom::Point2D<T>& dirvec,
123  const Geom::Point2D<T>& posvec,
124  T eps );
125 
126  static bool isRightOfLine( const Geom::Point2D<T>& point,
127  const Geom::Point2D<T>& dirvec,
128  const Geom::Point2D<T>& posvec );
129 
130  static bool doCoincide( const Geom::Point2D<T>& point1,
131  const Geom::Point2D<T>& point2,
132  T eps=mDefEps );
133 
134  static double sgnDistToLine( const Geom::Point2D<T>& point,
135  const Geom::Point2D<T>& dirvec,
136  const Geom::Point2D<T>& posvec );
137 
138  template <class ST> ST sgnArea2() const;
139 
140  static double distToSegment( const Geom::Point2D<T>& p1,
141  const Geom::Point2D<T>& p2,
142  const Geom::Point2D<T>& refpt,
143  double* fractionptr=0 );
144 
145  TypeSet<Geom::Point2D<T> > poly_;
146  bool closed_;
147  Geom::Point2D<T> udf_;
148 
149  mutable Interval<T> xrg_;
150  mutable Interval<T> yrg_;
151 };
152 
153 
154 template <class T> inline
155 bool ODPolygon<T>::operator==( const ODPolygon<T>& poly ) const
156 {
157  if ( data() != poly.data() ) return false;
158  return true;
159 }
160 
161 
162 template <class T> inline
163 void ODPolygon<T>::getData( bool forx, TypeSet<T>& ts ) const
164 {
165  for ( int idx=0; idx<size(); idx++ )
166  {
167  const Geom::Point2D<T>& vtx = poly_[idx];
168  ts += forx ? vtx.x : vtx.y;
169  }
170 }
171 
172 
173 template <class T> inline
174 void fillPar( IOPar& iop, const ODPolygon<T>& poly, const char* inpkey )
175 {
176  const BufferString keywd( inpkey ); const char* key = keywd.buf();
177  iop.setYN( IOPar::compKey(key,"Closed"), poly.isClosed() );
178  iop.set( IOPar::compKey(key,"Undef"), poly.getUdf().x, poly.getUdf().y );
179  TypeSet<T> ts; poly.getData( true, ts );
180  iop.set( IOPar::compKey(key,"Data.X"), ts );
181  ts.erase(); poly.getData( false, ts );
182  iop.set( IOPar::compKey(key,"Data.Y"), ts );
183 }
184 
185 
186 template <class T> inline
187 void usePar( const IOPar& iop, ODPolygon<T>& poly, const char* inpkey )
188 {
189  const BufferString keywd( inpkey ); const char* key = keywd.buf();
190  bool yn = false; iop.getYN( IOPar::compKey(key,"Closed"), yn );
191  poly.setClosed( yn );
192  Geom::Point2D<T> pt; iop.get( IOPar::compKey(key,"Undef"), pt.x, pt.y );
193  poly.setUdf( pt );
194 
195  if ( !iop.find( IOPar::compKey(key,"Data.X") )
196  || !iop.find( IOPar::compKey(key,"Data.Y") ) )
197  return;
198 
199  poly.setEmpty(); TypeSet<T> tsx, tsy;
200  iop.get( IOPar::compKey(key,"Data.X"), tsx );
201  iop.get( IOPar::compKey(key,"Data.Y"), tsy );
202  const int sz = tsx.size() > tsy.size() ? tsy.size() : tsx.size();
203  for ( int idx=0; idx<sz; idx++ )
204  poly.add( Geom::Point2D<T>(tsx[idx],tsy[idx]) );
205 }
206 
207 
208 template <class T> inline
209 void ODPolygon<T>::insert( int idx, const Geom::Point2D<T>& vtx )
210 {
211  if ( idx>=0 && idx<=size() )
212  poly_.insert( idx, vtx );
213  xrg_.set( mUdf(T), mUdf(T) );
214  yrg_.set( mUdf(T), mUdf(T) );
215 }
216 
217 
218 template <class T> inline
219 void ODPolygon<T>::erase()
220 {
221  poly_.erase();
222  xrg_.set( mUdf(T), mUdf(T) );
223  yrg_.set( mUdf(T), mUdf(T) );
224 }
225 
226 
227 template <class T> inline
228 void ODPolygon<T>::add( const Geom::Point2D<T>& vtx )
229 {
230  poly_+=vtx;
231  xrg_.set( mUdf(T), mUdf(T) );
232  yrg_.set( mUdf(T), mUdf(T) );
233 }
234 
235 
236 template <class T> inline
237 void ODPolygon<T>::remove( int idx )
238 {
239  if ( poly_.validIdx(idx) )
240  poly_.removeSingle( idx );
241  xrg_.set( mUdf(T), mUdf(T) );
242  yrg_.set( mUdf(T), mUdf(T) );
243 }
244 
245 
246 template <class T> inline
247 const Geom::Point2D<T>& ODPolygon<T>::getVertex( int idx ) const
248 { return poly_.validIdx(idx) ? poly_[idx] : udf_; }
249 
250 
251 template <class T> inline
252 const Geom::Point2D<T>& ODPolygon<T>::nextVertex( int idx ) const
253 { return getVertex( idx+1<size() ? idx+1 : 0 ); }
254 
255 
256 template <class T> inline
257 const Geom::Point2D<T>& ODPolygon<T>::prevVertex( int idx ) const
258 { return getVertex( idx ? idx-1 : size()-1 ); }
259 
260 
261 template <class T> inline
262 Interval<T> ODPolygon<T>::getRange( bool forx ) const
263 {
264  if ( poly_.isEmpty() ) return Interval<T>( udf_.x, udf_.y );
265  Geom::Point2D<T> vtx0 = poly_[0];
266  Interval<T> ret = forx ? xrg_ : yrg_;
267  if ( !mIsUdf(ret.start) && !mIsUdf(ret.stop) )
268  return ret;
269  ret.start = ret.stop = forx ? vtx0.x : vtx0.y;
270  for ( int idx=1; idx<size(); idx++ )
271  {
272  const Geom::Point2D<T>& vtx = poly_[idx];
273  const T val = forx ? vtx.x : vtx.y;
274  if ( val < ret.start ) ret.start = val;
275  else if ( val > ret.stop ) ret.stop = val;
276  }
277  if ( forx )
278  xrg_ = ret;
279  else
280  yrg_ = ret;
281  return ret;
282 }
283 
284 
285 template <class T> inline
286 bool ODPolygon<T>::isInside( const Geom::Point2D<T>& point,
287  bool inclborder, T eps ) const
288 {
289  const T abseps = eps<0 ? -eps : eps;
290  if ( (!mIsUdf(xrg_.start) && !mIsUdf(xrg_.stop) &&
291  (xrg_.start>point.x+abseps || xrg_.stop<point.x-abseps)) ||
292  (!mIsUdf(yrg_.start) && !mIsUdf(yrg_.stop) &&
293  (yrg_.start>point.y+abseps || yrg_.stop<point.y-abseps)) )
294  return false;
295 
296  const Geom::Point2D<T> arbitrarydir( 1, 0 );
297 
298  bool nrcrossingsodd = false;
299  for ( int idx=0; idx<size(); idx++ )
300  {
301  const Geom::Point2D<T>& vtxcurr = poly_[idx];
302  const Geom::Point2D<T>& vtxnext = nextVertex( idx );
303 
304  if ( isOnSegment(point, vtxcurr, vtxnext, eps) )
305  return inclborder;
306  if ( isEdgeCrossing(arbitrarydir, point, vtxcurr, vtxnext) )
307  nrcrossingsodd = !nrcrossingsodd;
308  }
309 
310  return nrcrossingsodd;
311 }
312 
313 
314 template <class T> inline
315 bool ODPolygon<T>::segmentMeetsBorder( const Geom::Point2D<T>& pt1,
316  const Geom::Point2D<T>& pt2,
317  T eps ) const
318 {
319  for ( int idx=0; idx<size(); idx++ )
320  {
321  const Geom::Point2D<T>& vtxcurr = poly_[idx];
322  const Geom::Point2D<T>& vtxnext = nextVertex( idx );
323 
324  if ( doSegmentsMeet(pt1, pt2, vtxcurr, vtxnext, eps) )
325  return true;
326  }
327 
328  return false;
329 }
330 
331 
332 template <class T> inline
333 bool ODPolygon<T>::windowOverlaps( const Interval<T>& xrange,
334  const Interval<T>& yrange,
335  T eps ) const
336 {
337  ODPolygon window;
338  window.add( Geom::Point2D<T>(xrange.start, yrange.start) );
339  window.add( Geom::Point2D<T>(xrange.stop, yrange.start) );
340  window.add( Geom::Point2D<T>(xrange.stop, yrange.stop) );
341  window.add( Geom::Point2D<T>( xrange.start, yrange.stop) );
342 
343  return isInside( window, eps );
344 }
345 
346 
347 template <class T> inline
348 int ODPolygon<T>::isInside( const ODPolygon& testpoly, T eps ) const
349 {
350  if ( isEmpty() || testpoly.isEmpty() )
351  return 0;
352 
353  for ( int idx=0; idx<size(); idx++ )
354  {
355  const Geom::Point2D<T>& vtxcurr = poly_[idx];
356  const Geom::Point2D<T>& vtxnext = nextVertex( idx );
357 
358  if ( testpoly.segmentMeetsBorder(vtxcurr, vtxnext, eps) )
359  return 1;
360  }
361 
362  if ( isInside(testpoly.poly_[0], false, eps) )
363  return 2;
364 
365  return testpoly.isInside(poly_[0], false, eps) ? 1 : 0;
366 }
367 
368 
369 template <class T> inline
370 void ODPolygon<T>::removeZeroLengths()
371 {
372  const int startidx = isClosed() ? size()-1 : size()-2;
373  for ( int idx=startidx; idx>=0; idx-- )
374  {
375  if ( poly_[idx]==nextVertex(idx) && size()>1 )
376  remove(idx);
377  }
378 }
379 
380 
381 template <class T> inline
382 bool ODPolygon<T>::isUTurn( int idx ) const
383 {
384  if ( !validIdx(idx) || ( !isClosed() && (idx==0 || idx==size()-1) ) )
385  return false;
386 
387  const Geom::Point2D<T>& vec1 = prevVertex(idx) - poly_[idx];
388  const Geom::Point2D<T>& vec2 = nextVertex(idx) - poly_[idx];
389 
390  return vec1.x*vec2.y-vec1.y*vec2.x==0 && vec1.x*vec2.x+vec1.y*vec2.y>0;
391 }
392 
393 
394 template <class T> inline
395 bool ODPolygon<T>::isSelfIntersecting() const
396 {
397  ODPolygon<T> plg = *this;
398  plg.removeZeroLengths();
399 
400  const int stopidx = plg.isClosed() ? plg.size() : plg.size()-1;
401  for ( int idx=0; idx<stopidx; idx++ )
402  {
403  if ( plg.isUTurn(idx) )
404  return true;
405 
406  const Geom::Point2D<T>& vtxcurr = plg.poly_[idx];
407  const Geom::Point2D<T>& vtxnext = plg.nextVertex( idx );
408 
409  for ( int idy=0; idy<stopidx; idy++ )
410  {
411  const int dif = abs( idx-idy );
412  if ( dif<=1 || dif>=plg.size()-1 )
413  continue;
414 
415  const Geom::Point2D<T>& pt1 = plg.poly_[idy];
416  const Geom::Point2D<T>& pt2 = plg.nextVertex( idy );
417 
418  if ( vtxcurr==pt1 || vtxcurr==pt2 )
419  return true;
420 
421  if ( isEdgeCrossing(vtxnext-vtxcurr, vtxcurr, pt1, pt2) &&
422  isEdgeCrossing(vtxcurr-vtxnext, vtxnext, pt1, pt2) )
423  return true;
424  }
425  }
426  return false;
427 }
428 
429 
430 template <class T> inline
431 int ODPolygon<T>::doSegmentsMeet( const Geom::Point2D<T>& p1,
432  const Geom::Point2D<T>& p2,
433  const Geom::Point2D<T>& q1,
434  const Geom::Point2D<T>& q2,
435  T eps )
436 {
437  if ( isOnSegment(p1, q1, q2, eps) || isOnSegment(p2, q1, q2, eps) ||
438  isOnSegment(q1, p1, p2, eps) || isOnSegment(q2, p1, p2, eps) )
439  return 1;
440 
441  if ( p1==p2 || q1==q2 || !isEdgeCrossing(p2-p1, p1, q1, q2) ||
442  !isEdgeCrossing(p1-p2, p2, q1, q2) )
443  return 0;
444 
445  return 2;
446 }
447 
448 
449 template <class T> inline
450 bool ODPolygon<T>::isOnSegment( const Geom::Point2D<T>& pt,
451  const Geom::Point2D<T>& pt0,
452  const Geom::Point2D<T>& pt1,
453  T eps )
454 {
455  return isOnHalfLine( pt, pt1-pt0, pt0, eps ) &&
456  isOnHalfLine( pt, pt0-pt1, pt1, eps );
457 }
458 
459 
460 template <class T> inline
461 bool ODPolygon<T>::isOnHalfLine( const Geom::Point2D<T>& point,
462  const Geom::Point2D<T>& dirvec,
463  const Geom::Point2D<T>& endvec,
464  T eps )
465 {
466  if ( doCoincide(point, endvec, eps) )
467  return true;
468  if ( !isOnLine(point, dirvec, endvec, eps) )
469  return false;
470  const Geom::Point2D<T> rot90dirvec( -dirvec.y, dirvec.x );
471  return isRightOfLine( point, rot90dirvec, endvec );
472 }
473 
474 
475 template <class T> inline
476 bool ODPolygon<T>::isEdgeCrossing( const Geom::Point2D<T>& raydir,
477  const Geom::Point2D<T>& raysrc,
478  const Geom::Point2D<T>& vtx1,
479  const Geom::Point2D<T>& vtx2 )
480 {
481  const bool vtx1right = isRightOfLine( vtx1, raydir, raysrc );
482  const bool vtx2right = isRightOfLine( vtx2, raydir, raysrc );
483 
484  if ( vtx1right && !vtx2right )
485  return !isRightOfLine( raysrc, vtx2-vtx1, vtx1 );
486  if ( !vtx1right && vtx2right )
487  return !isRightOfLine( raysrc, vtx1-vtx2, vtx2 );
488  return false;
489 }
490 
491 
492 template <class T> inline
493 bool ODPolygon<T>::isOnLine( const Geom::Point2D<T>& point,
494  const Geom::Point2D<T>& dirvec,
495  const Geom::Point2D<T>& posvec,
496  T eps )
497 {
498  const double signeddist = sgnDistToLine( point, dirvec, posvec );
499  return signeddist * signeddist <= eps * eps;
500 }
501 
502 template <class T> inline
503 bool ODPolygon<T>::isRightOfLine( const Geom::Point2D<T>& point,
504  const Geom::Point2D<T>& dirvec,
505  const Geom::Point2D<T>& posvec )
506 {
507  return sgnDistToLine( point, dirvec, posvec ) > 0;
508 }
509 
510 
511 template <class T> inline
512 bool ODPolygon<T>::doCoincide( const Geom::Point2D<T>& point1,
513  const Geom::Point2D<T>& point2,
514  T eps )
515 {
516  return point1.sqDistTo( point2 ) <= eps * eps;
517 }
518 
519 
520 template <class T> inline
521 double ODPolygon<T>::sgnDistToLine( const Geom::Point2D<T>& point,
522  const Geom::Point2D<T>& dirvec,
523  const Geom::Point2D<T>& posvec )
524 {
525  const double nolinedist = 0;
526 
527  const double dirveclen = dirvec.distTo( Geom::Point2D<T>(0,0) );
528  if ( mIsZero(dirveclen, mDefEps) )
529  return nolinedist;
530  const double substpointinlineeqn =
531  dirvec.y * ( point.x - posvec.x )-dirvec.x * ( point.y - posvec.y );
532  return substpointinlineeqn / dirveclen;
533 }
534 
535 
536 template <class T>
537 template <class ST> inline
538 ST ODPolygon<T>::sgnArea2() const
539 {
540  ST area2 = 0;
541 
542  if ( poly_.isEmpty() )
543  return area2;
544 
545  const Geom::Point2D<T>& pt0 = poly_[0];
546  for ( int idx=1; idx<size()-1; idx++ )
547  {
548  const Geom::Point2D<T>& pt1 = poly_[idx];
549  const Geom::Point2D<T>& pt2 = nextVertex( idx );
550  area2 += (ST) ( (pt1.x-pt0.x) * (pt2.y-pt0.y) -
551  (pt2.x-pt0.x) * (pt1.y-pt0.y) );
552  }
553 
554  return area2;
555 }
556 
557 
558 template <class T> inline
559 void ODPolygon<T>::convexHull()
560 {
561  // Code based on the Graham scan algorithm (1972)
562 
563  setClosed( true );
564  if ( size()<2 )
565  return;
566 
567  // Find guaranteed vertex of the convex hull to become pivot
568  Geom::Point2D<T> pivot = poly_[0];
569  for ( int idx=1; idx<size(); idx++ )
570  {
571  const Geom::Point2D<T>& vtx = poly_[idx];
572  if ( vtx.x<pivot.x || (vtx.x==pivot.x && vtx.y<pivot.y) )
573  pivot = vtx;
574  }
575 
576  // Remove all pivot copies
577  for ( int idx=size()-1; idx>=0; idx-- )
578  {
579  if ( pivot == poly_[idx] )
580  //poly_.removeSingle( idx, false );
581  poly_.removeSingle( idx );
582  }
583 
584  // Angular sort of all pivot-to-point segments
585  for ( int idx=size()-2; idx>=0; idx-- )
586  {
587  const Geom::Point2D<T>& vtx = poly_[idx];
588  for ( int idy=size()-1; idy>idx; idy-- )
589  {
590  const Geom::Point2D<T>& vty = poly_[idy];
591  const double dist = sgnDistToLine( vty, vtx-pivot, pivot );
592 
593  if ( dist > 0 )
594  continue;
595 
596  if ( dist < 0 )
597  poly_.insert( idy+1, vtx );
598  else if ( pivot.sqDistTo(vtx) > pivot.sqDistTo(vty) )
599  poly_[idy] = vtx;
600 
601  poly_.removeSingle( idx );
602  break;
603  }
604  }
605 
606  // Expand convex hull incrementally by backward removal of inner points
607  for ( int idx=size()-3; idx>=0; idx-- )
608  {
609  const Geom::Point2D<T>& vtx = poly_[idx];
610  while ( idx<size()-2 )
611  {
612  const Geom::Point2D<T>& vty = poly_[idx+1];
613  const Geom::Point2D<T>& vtz = poly_[idx+2];
614  if ( isRightOfLine(vtz, vty-vtx, vtx) )
615  break;
616 
617  poly_.removeSingle( idx+1 );
618  }
619  }
620 
621  poly_ += pivot;
622 
623  xrg_.set( mUdf(T), mUdf(T) );
624  yrg_.set( mUdf(T), mUdf(T) );
625 }
626 
627 
628 template <class T> inline
629 void ODPolygon<T>::reverse()
630 {
631  const int sz = poly_.size();
632  for ( int idx=0; idx<sz/2-1; idx++ )
633  {
634  Geom::Point2D<T> temp = poly_[idx];
635  poly_[idx] = poly_[sz-1-idx];
636  poly_[sz-1-idx] = temp;
637  }
638 
639  xrg_.set( mUdf(T), mUdf(T) );
640  yrg_.set( mUdf(T), mUdf(T) );
641 }
642 
643 
644 template <class T> inline
645 double ODPolygon<T>::distToSegment( const Geom::Point2D<T>& p1,
646  const Geom::Point2D<T>& p2,
647  const Geom::Point2D<T>& refpt,
648  double* fractionptr )
649 {
650  double frac = 0;
651 
652  if ( p1 != p2 )
653  {
654  const Geom::Point2D<T> dif = p2 - p1;
655  const double numerator = dif.x*(refpt.x-p1.x) + dif.y*(refpt.y-p1.y);
656  frac = numerator / (dif.x*dif.x + dif.y*dif.y);
657 
658  if ( frac < 0 )
659  frac = 0;
660  if ( frac > 1 )
661  frac = 1;
662  }
663 
664  if ( fractionptr )
665  *fractionptr = frac;
666 
667  const Geom::Point2D<T> pointonseg( (T)(p1.x * (1-frac) + p2.x * frac),
668  (T)(p1.y * (1-frac) + p2.y * frac) );
669  return refpt.distTo( pointonseg );
670 }
671 
672 
673 template <class T> inline
674 double ODPolygon<T>::distTo( const Geom::Point2D<T>& refpt,
675  int* segmentidxptr, double* fractionptr ) const
676 {
677  const int sz = size();
678  if ( !sz )
679  return mUdf(double);
680 
681  double mindist = MAXDOUBLE;
682 
683  for ( int idx=(isClosed() ? sz-1 : sz-2); idx>=0; idx-- )
684  {
685  const Geom::Point2D<T>& pt1 = getVertex( idx );
686  const Geom::Point2D<T>& pt2 = nextVertex( idx );
687  double frac;
688  const double dist = distToSegment( pt1, pt2, refpt, &frac );
689 
690  if ( mindist >= dist )
691  {
692  mindist = dist;
693  if ( segmentidxptr )
694  *segmentidxptr = idx;
695  if ( fractionptr )
696  *fractionptr = frac;
697  }
698  }
699 
700  return mindist;
701 }
702 
703 
704 template <class T> inline
705 double ODPolygon<T>::maxDistToBorderEstimate( double maxrelerr ) const
706 {
707  if ( maxrelerr <= 0.0 )
708  return mUdf(double);
709 
710  if ( size() < 3 )
711  return isEmpty() ? mUdf(double) : 0.0;
712 
713  const double upperbound = mMIN( 0.5 * getRange(true).width(),
714  0.5 * getRange(false).width() );
715  if ( !upperbound )
716  return 0.0;
717 
718  ODPolygon<double> poly;
719  for ( int idx=0; idx<size(); idx++ )
720  {
721  const Geom::Point2D<T>& pt = getVertex(idx);
722  poly.add( Geom::Point2D<double>(pt.x, pt.y) );
723  }
724 
725  double maxdist = 0.0;
726  for ( int idx=0; idx<poly.size(); idx++ )
727  {
728  Geom::Point2D<double> curpt =
729  (poly.prevVertex(idx)+poly.getVertex(idx)+poly.nextVertex(idx)) / 3;
730 
731  if ( !poly.isInside(curpt, false, mDefEps) )
732  continue;
733 
734  double curdist = poly.distTo( curpt );
735  double gamma = 0.1 * upperbound;
736 
737  for ( int step=0; step<100; step++ )
738  {
739  if ( curdist > maxdist )
740  maxdist = curdist;
741 
742  // Gradient ascent
743  Geom::Point2D<double> pt1 = curpt;
744  Geom::Point2D<double> pt2 = curpt;
745  const double delta = 0.0001 * upperbound;
746  pt1.x += delta;
747  pt2.y += delta;
748  const double dist1 = poly.distTo( pt1 );
749  const double dist2 = poly.distTo( pt2 );
750  Geom::Point2D<double> nextpt( dist1-curdist, dist2-curdist );
751  nextpt *= gamma/delta;
752  nextpt += curpt;
753 
754  if ( !poly.isInside(nextpt, false, mDefEps) )
755  {
756  gamma *= 0.5;
757  continue;
758  }
759 
760  double nextdist = poly.distTo( nextpt );
761 
762  if ( nextdist <= curdist )
763  {
764  gamma *= 0.5;
765 
766  if ( curpt.distTo(nextpt) <= 2*maxrelerr*curdist )
767  break;
768  }
769 
770  curpt = nextpt;
771  curdist = nextdist;
772  }
773  }
774  return maxdist;
775 }
776 
777 
778 template <class T> inline
779 void ODPolygon<T>::keepBendPoints( float eps )
780 {
781  removeZeroLengths();
782 
783  const int sz = size();
784  if ( sz < 3 )
785  return;
786 
787  const int extra = closed_ ? 1 : 0;
788 
789  TypeSet<Coord> coords;
790  for ( int idx=-extra; idx<sz+extra; idx++ )
791  {
792  const Geom::Point2D<T>& vtx = getVertex( (sz+idx)%sz );
793  coords += Coord( vtx.x, vtx.y );
794  }
795 
796  BendPointFinder2D finder( coords, eps );
797  finder.execute();
798  const TypeSet<int>& bendpoints = finder.bendPoints();
799 
800  int bpidx = bendpoints.size()-extra-1;
801  if ( bpidx<0 )
802  return;
803 
804  for ( int vtxidx=sz-1; vtxidx>=0; vtxidx-- )
805  {
806  if ( bpidx>=0 && vtxidx==bendpoints[bpidx]-extra )
807  bpidx--;
808  else
809  remove( vtxidx );
810  }
811 
812  if ( closed_ && poly_.size()>1 && poly_.first() != poly_.last() )
813  poly_ += poly_.first();
814 }
815 
ODPolygon::segmentMeetsBorder
bool segmentMeetsBorder(const Geom::Point2D< T > &pt1, const Geom::Point2D< T > &pt2, T eps) const
Definition: polygon.h:315
IOPar::set
void set(const char *ky, const char *val)
ODPolygon::setEmpty
void setEmpty()
Definition: polygon.h:43
ODPolygon
(Closed) sequence of connected 2-D coordinates.
Definition: polygon.h:27
ODPolygon::isOnHalfLine
static bool isOnHalfLine(const Geom::Point2D< T > &point, const Geom::Point2D< T > &dirvec, const Geom::Point2D< T > &endvec, T eps)
Definition: polygon.h:461
ODPolygon::anticlockwise
bool anticlockwise() const
Definition: polygon.h:86
ODPolygon::yrg_
Interval< T > yrg_
Definition: polygon.h:150
ODPolygon::isRightOfLine
static bool isRightOfLine(const Geom::Point2D< T > &point, const Geom::Point2D< T > &dirvec, const Geom::Point2D< T > &posvec)
Definition: polygon.h:503
ODPolygon::closed_
bool closed_
Definition: polygon.h:146
BendPointFinderBase::bendPoints
const TypeSet< int > & bendPoints() const
Definition: bendpointfinder.h:32
ODPolygon::removeZeroLengths
void removeZeroLengths()
Definition: polygon.h:370
IOPar::get
bool get(const char *, short &) const
BendPointFinder2D
Used to find bendpoints in set of XY Coordinates.
Definition: bendpointfinder.h:82
ODPolygon::keepBendPoints
void keepBendPoints(float eps)
Definition: polygon.h:779
ODPolygon::isInside
int isInside(const ODPolygon &testpoly, T eps=0) const
Definition: polygon.h:348
ODPolygon::poly_
TypeSet< Geom::Point2D< T > > poly_
Definition: polygon.h:145
ODPolygon::setUdf
void setUdf(Geom::Point2D< T > pt)
Definition: polygon.h:69
ParallelTask::execute
bool execute()
Definition: paralleltask.h:70
ODPolygon::remove
void remove(int idx)
Definition: polygon.h:237
ODPolygon::getUdf
Geom::Point2D< T > getUdf() const
Definition: polygon.h:70
mIsUdf
#define mIsUdf(val)
Use mIsUdf to check for undefinedness of simple types.
Definition: undefval.h:289
ODPolygon::udf_
Geom::Point2D< T > udf_
Definition: polygon.h:147
ODPolygon::insert
void insert(int idx, const Geom::Point2D< T > &vtx)
Definition: polygon.h:209
usePar
void usePar(const IOPar &iop, ODPolygon< T > &poly, const char *inpkey)
Definition: polygon.h:187
mDefEps
#define mDefEps
Definition: commondefs.h:71
ODPolygon::isOnLine
static bool isOnLine(const Geom::Point2D< T > &point, const Geom::Point2D< T > &dirvec, const Geom::Point2D< T > &posvec, T eps)
Definition: polygon.h:493
ODPolygon::size
int size() const
Definition: polygon.h:41
typeset.h
ODPolygon::xrg_
Interval< T > xrg_
Definition: polygon.h:149
bufstring.h
ODPolygon::isEdgeCrossing
static bool isEdgeCrossing(const Geom::Point2D< T > &raydir, const Geom::Point2D< T > &raysrc, const Geom::Point2D< T > &vtx1, const Geom::Point2D< T > &vtx2)
Definition: polygon.h:476
ODPolygon::data
const TypeSet< Geom::Point2D< T > > & data() const
Definition: polygon.h:71
ODPolygon::getVertex
const Geom::Point2D< T > & getVertex(int idx) const
Definition: polygon.h:247
isEmpty
bool isEmpty(const NLAModel *mdl)
ODPolygon::isInside
bool isInside(const Geom::Point2D< T > &, bool inclborder, T eps) const
Definition: polygon.h:286
ODPolygon::add
void add(const Geom::Point2D< T > &vtx)
Definition: polygon.h:228
ODPolygon::setClosed
void setClosed(bool yn)
Definition: polygon.h:67
Coord
A cartesian coordinate in 2D space.
Definition: coord.h:25
File::remove
bool remove(const char *)
ODPolygon::isUTurn
bool isUTurn(int idx) const
Definition: polygon.h:382
IOPar::compKey
static const char * compKey(const char *, const char *)
The composite key: (a,b) -> a.b.
mIsZero
#define mIsZero(x, eps)
Definition: commondefs.h:66
mClass
#define mClass(module)
Definition: commondefs.h:181
ODPolygon::erase
void erase()
Definition: polygon.h:219
Geom::Point2D::distTo
double distTo(const Point2D< T > &) const
Definition: geometry.h:386
ODPolygon::sgnArea2
ST sgnArea2() const
Definition: polygon.h:538
OD::String::buf
const char * buf() const
Definition: odstring.h:46
Geom::Point2D::sqDistTo
T sqDistTo(const Point2D< T > &) const
Definition: geometry.h:391
ODPolygon::validIdx
bool validIdx(int idx) const
Definition: polygon.h:42
ODPolygon::ODPolygon
ODPolygon(const TypeSet< Geom::Point2D< T > > &plg)
Definition: polygon.h:33
ODPolygon::prevVertex
const Geom::Point2D< T > & prevVertex(int idx) const
Definition: polygon.h:257
ODPolygon::isEmpty
bool isEmpty() const
Definition: polygon.h:39
fillPar
void fillPar(IOPar &iop, const ODPolygon< T > &poly, const char *inpkey)
Definition: polygon.h:174
ODPolygon::doSegmentsMeet
static int doSegmentsMeet(const Geom::Point2D< T > &p1, const Geom::Point2D< T > &p2, const Geom::Point2D< T > &q1, const Geom::Point2D< T > &q2, T eps)
Definition: polygon.h:431
ODPolygon::maxDistToBorderEstimate
double maxDistToBorderEstimate(double maxrelerr=0.001) const
Definition: polygon.h:705
MAXDOUBLE
#define MAXDOUBLE
Definition: commondefs.h:120
BufferString
OD::String with its own variable length buffer. The buffer has a guaranteed minimum size.
Definition: bufstring.h:40
ODPolygon::area
double area() const
Definition: polygon.h:84
ODPolygon::reverse
void reverse()
Definition: polygon.h:629
Geom
Definition: geometry.h:19
Geom::Point2D
Basic point class.
Definition: geometry.h:27
IOPar::find
const char * find(const char *) const
returns null if not found
ODPolygon::sgnDistToLine
static double sgnDistToLine(const Geom::Point2D< T > &point, const Geom::Point2D< T > &dirvec, const Geom::Point2D< T > &posvec)
Definition: polygon.h:521
ODPolygon::operator==
bool operator==(const ODPolygon< T > &) const
Definition: polygon.h:155
mMIN
#define mMIN(x, y)
Definition: commondefs.h:62
ODPolygon::distToSegment
static double distToSegment(const Geom::Point2D< T > &p1, const Geom::Point2D< T > &p2, const Geom::Point2D< T > &refpt, double *fractionptr=0)
Definition: polygon.h:645
ODPolygon::isOnSegment
static bool isOnSegment(const Geom::Point2D< T > &pt, const Geom::Point2D< T > &pt0, const Geom::Point2D< T > &pt1, T eps)
Definition: polygon.h:450
ODPolygon::operator>
bool operator>(const ODPolygon< T > &plg) const
Definition: polygon.h:96
Geom::Point2D::y
T y
Definition: geometry.h:68
IOPar::setYN
void setYN(const char *, bool)
Geom::Point2D::x
T x
Definition: geometry.h:67
ODPolygon::nextVertex
const Geom::Point2D< T > & nextVertex(int idx) const
Definition: polygon.h:252
ODPolygon::distTo
double distTo(const Geom::Point2D< T > &refpt, int *segmentidxptr=0, double *fractionptr=0) const
Definition: polygon.h:674
mUdf
#define mUdf(type)
Use this macro to get the undefined for simple types.
Definition: undefval.h:274
ODPolygon::clockwise
bool clockwise() const
Definition: polygon.h:85
ODPolygon::isClosed
bool isClosed() const
Definition: polygon.h:68
ODPolygon::getRange
Interval< T > getRange(bool of_x) const
Definition: polygon.h:262
iopar.h
ODPolygon::convexHull
void convexHull()
Definition: polygon.h:559
Interval
Interval of values.
Definition: commontypes.h:30
ODPolygon::windowOverlaps
bool windowOverlaps(const Interval< T > &xrange, const Interval< T > &yrange, T eps) const
Definition: polygon.h:333
IOPar::getYN
bool getYN(const char *, bool &) const
bendpointfinder.h
IOPar
Generalized set of parameters of the keyword-value type.
Definition: iopar.h:55
ODPolygon::ODPolygon
ODPolygon()
Definition: polygon.h:29
ODPolygon::doCoincide
static bool doCoincide(const Geom::Point2D< T > &point1, const Geom::Point2D< T > &point2, T eps=mDefEps)
Definition: polygon.h:512
ODPolygon::getData
void getData(bool of_x, TypeSet< T > &) const
Definition: polygon.h:163
TypeSet
Sets of (small) copyable elements.
Definition: commontypes.h:29
ODPolygon::isSelfIntersecting
bool isSelfIntersecting() const
Definition: polygon.h:395
coord.h

Generated at for the OpendTect seismic interpretation project. Copyright (C): dGB Beheer B.V. 1995-2021