OTB  6.7.0
Orfeo Toolbox
otbPolygon.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbPolygon_hxx
22 #define otbPolygon_hxx
23 
24 #include "otbPolygon.h"
25 
26 namespace otb
27 {
28 
29 template<class TValue>
30 void
33 {
34  Superclass::AddVertex(vertex);
35  m_AreaIsValid = false;
36 }
37 
44 template<class TValue>
45 bool
47 ::IsInside(VertexType point) const
48 {
49  double x = point[0];
50  double y = point[1];
51  unsigned int crossingCount = 0;
52  VertexListConstIteratorType it = this->GetVertexList()->Begin();
53  double xa = it.Value()[0];
54  double ya = it.Value()[1];
55  ++it;
56  while (it != this->GetVertexList()->End())
57  {
58  double xb = it.Value()[0];
59  double yb = it.Value()[1];
60  if (std::abs(xb - xa) < m_Epsilon)
61  {
62  if (ya > yb && xa > x && y >= yb && y < ya)
63  {
64  ++crossingCount;
65  }
66  else if (ya < yb && xa > x && y >= ya && y < yb)
67  {
68  ++crossingCount;
69  }
70  }
71  else if (std::abs(yb - ya) >= m_Epsilon)
72  {
73  double xcross = xa + (xb - xa) * (y - ya) / (yb - ya);
74 
75  if (ya > yb && xcross > x && y >= yb && y < ya)
76  {
77  ++crossingCount;
78  }
79  else if (ya < yb && xcross > x && y >= ya && y < yb)
80  {
81  ++crossingCount;
82  }
83  }
84  xa = xb;
85  ya = yb;
86  ++it;
87  }
88  double xb = this->GetVertexList()->Begin().Value()[0];
89  double yb = this->GetVertexList()->Begin().Value()[1];
90  if (std::abs(xb - xa) < m_Epsilon)
91  {
92  if (ya > yb && xa > x && y >= yb && y < ya)
93  {
94  ++crossingCount;
95  }
96  else if (ya < yb && xa > x && y >= ya && y < yb)
97  {
98  ++crossingCount;
99  }
100  }
101  else if (std::abs(yb - ya) >= m_Epsilon)
102  {
103  double xcross = xa + (xb - xa) * (y - ya) / (yb - ya);
104 
105  if (ya > yb && xcross > x && y >= yb && y < ya)
106  {
107  ++crossingCount;
108  }
109  else if (ya < yb && xcross > x && y >= ya && y < yb)
110  {
111  ++crossingCount;
112  }
113  }
114  //std::cout<<"Crossing count: "<<crossingCount<<std::endl;
115  return (crossingCount % 2 == 1);
116 }
117 
123 template<class TValue>
124 bool
126 ::IsOnEdge(VertexType point) const
127 {
128  //std::cout<<"Checking point "<<point<<std::endl;
129  bool resp = false;
130  double x = point[0];
131  double y = point[1];
132  double xb, yb;
133  VertexListConstIteratorType it = this->GetVertexList()->Begin();
134  double xa = it.Value()[0];
135  double ya = it.Value()[1];
136  double xbegin = xa;
137  double ybegin = ya;
138  ++it;
139  while (!resp && it != this->GetVertexList()->End())
140  {
141  xb = it.Value()[0];
142  yb = it.Value()[1];
143  if (std::abs(xb - xa) >= m_Epsilon)
144  {
145  double cd = (yb - ya) / (xb - xa);
146  double oo = (ya - cd * xa);
147  double xmin = std::min(xa, xb);
148  double xmax = std::max(xa, xb);
149  if ((std::abs(y - cd * x - oo) < m_Epsilon)
150  && (x <= xmax)
151  && (x >= xmin))
152  {
153  //std::cout<<"Case 1: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
154  resp = true;
155  }
156  }
157  else
158  {
159  double ymin = std::min(ya, yb);
160  double ymax = std::max(ya, yb);
161  if ((std::abs(x - xa) < m_Epsilon)
162  && (y <= ymax)
163  && (y >= ymin))
164  {
165  resp = true;
166  //std::cout<<"Case 2: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
167  }
168  }
169  xa = xb;
170  ya = yb;
171  ++it;
172  }
173  xb = xbegin;
174  yb = ybegin;
175  if (std::abs(xb - xa) >= m_Epsilon)
176  {
177  double cd = (yb - ya) / (xb - xa);
178  double oo = (ya - cd * xa);
179  double xmin = std::min(xa, xb);
180  double xmax = std::max(xa, xb);
182 
183  if ((std::abs(y - cd * x - oo) < m_Epsilon)
184  && (x <= xmax)
185  && (x >= xmin))
186  {
187  resp = true;
188  //std::cout<<"Case 1: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
189  }
190  }
191  else
192  {
193  double ymin = std::min(ya, yb);
194  double ymax = std::max(ya, yb);
195  if ((std::abs(x - xa) <= m_Epsilon)
196  && (y <= ymax)
197  && (y >= ymin))
198  {
199  resp = true;
200  //std::cout<<"Case 2: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
201  }
202  }
203  //std::cout<<"Result: "<<resp<<std::endl;
204  return resp;
205 }
212 template<class TValue>
213 unsigned int
216 {
217  unsigned int resp = 0;
218  VertexListConstIteratorType it = this->GetVertexList()->Begin();
219  VertexListConstIteratorType it_end = this->GetVertexList()->End();
220  VertexType current = it.Value();
221  VertexType first = current;
222  ++it;
223  while (it != it_end)
224  {
225  //std::cout<<"Testing Is crossing "<<a<<" "<<b<<current<<it.Value()<<" = "<<IsCrossing(a, b, current, it.Value())<<std::endl;
226  if (IsCrossing(a, b, current, it.Value()))
227  {
228  ++resp;
229  }
230  current = it.Value();
231  ++it;
232  }
233  //std::cout<<"Testing Is crossing "<<a<<" "<<b<<current<<first<<" = "<<IsCrossing(a, b, current, first)<<std::endl;
234  if (IsCrossing(a, b, current, first))
235  {
236  ++resp;
237  }
238  return resp;
239 }
240 
247 template<class TValue>
248 unsigned int
251 {
252  unsigned int resp = 0;
253  VertexListConstIteratorType it = this->GetVertexList()->Begin();
254  VertexListConstIteratorType it_end = this->GetVertexList()->End();
255  VertexType current = it.Value();
256  VertexType first = current;
257  ++it;
258  while (it != it_end)
259  {
260  //std::cout<<"IsTouching "<<a<<" "<<b<<", "<<current<<" "<<it.Value()<<" -> "<<IsTouching(a, b, current, it.Value())<<std::endl;
261  if (IsTouching(a, b, current, it.Value()))
262  {
263  ++resp;
264  }
265  current = it.Value();
266  ++it;
267  }
268  //std::cout<<"IsTouching "<<a<<" "<<b<<", "<<current<<" "<<first<<" -> "<<IsTouching(a, b, current, first)<<std::endl;
269  if (IsTouching(a, b, current, first))
270  {
271  ++resp;
272  }
273  return resp;
274 }
275 
284 template<class TValue>
285 bool
288 {
289  bool resp = false;
290  double xbmin = std::min(b1[0], b2[0]);
291  double xbmax = std::max(b1[0], b2[0]);
292  double ybmin = std::min(b1[1], b2[1]);
293  double ybmax = std::max(b1[1], b2[1]);
294  double xamin = std::min(a1[0], a2[0]);
295  double xamax = std::max(a1[0], a2[0]);
296  double yamin = std::min(a1[1], a2[1]);
297  double yamax = std::max(a1[1], a2[1]);
298  if (std::abs(a1[0] - a2[0]) < m_Epsilon && std::abs(b1[0] - b2[0]) < m_Epsilon)
299  {
300  resp = false;
301  }
302  else if (std::abs(a1[0] - a2[0]) < m_Epsilon)
303  {
304  double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
305  double oo2 = b1[1] - cd2 * b1[0];
306  double ycross = cd2 * a1[0] + oo2;
307  resp = (xbmin <a1[0] && xbmax> a1[0]
308  && yamin <ycross && yamax> ycross);
309  }
310  else if (std::abs(b1[0] - b2[0]) < m_Epsilon)
311  {
312  double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
313  double oo1 = a1[1] - cd1 * a1[0];
314  double ycross = cd1 * b1[0] + oo1;
315  resp = (xamin <b1[0] && xamax> b1[0]
316  && ybmin <ycross && ybmax> ycross);
317  }
318  else
319  {
320  double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
321  double oo1 = a1[1] - cd1 * a1[0];
322  double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
323  double oo2 = b1[1] - cd2 * b1[0];
324  if (cd1 != cd2)
325  {
326  double xcross = (oo2 - oo1) / (cd1 - cd2);
327  resp = (xamin <xcross && xbmin <xcross
328  && xamax> xcross && xbmax> xcross);
329  }
330  }
331  return resp;
332 }
333 
342 template<class TValue>
343 bool
346 {
347  bool resp = false;
348  double xbmin = std::min(b1[0], b2[0]);
349  double xbmax = std::max(b1[0], b2[0]);
350  double ybmin = std::min(b1[1], b2[1]);
351  double ybmax = std::max(b1[1], b2[1]);
352  double xamin = std::min(a1[0], a2[0]);
353  double xamax = std::max(a1[0], a2[0]);
354  double yamin = std::min(a1[1], a2[1]);
355  double yamax = std::max(a1[1], a2[1]);
356  if (std::abs(a1[0] - a2[0]) < m_Epsilon && std::abs(b1[0] - b2[0]) < m_Epsilon)
357  {
358  resp = (std::abs(a1[0] - b1[0]) < m_Epsilon)
359  && ((a1[1] <= ybmax && a1[1] >= ybmin)
360  || (a2[1] <= ybmax && a2[1] >= ybmin)
361  || (b1[1] <= yamax && b1[1] >= yamin)
362  || (b2[1] <= yamax && b2[1] >= yamin));
363  }
364  else if (std::abs(a1[0] - a2[0]) < m_Epsilon)
365  {
366  double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
367  double oo2 = b1[1] - cd2 * b1[0];
369 
370  if (std::abs(a1[1] - cd2 * a1[0] - oo2) < m_Epsilon)
371  {
372  resp = (a1[0] >= xbmin && a1[0] <= xbmax);
373  }
374  else if (std::abs(a2[1] - cd2 * a2[0] - oo2) < m_Epsilon)
375  {
376  resp = (a2[0] >= xbmin && a2[0] <= xbmax);
377  }
378  }
379  else if (std::abs(b1[0] - b2[0]) < m_Epsilon)
380  {
381  double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
382  double oo1 = a1[1] - cd1 * a1[0];
383 
384  if (std::abs(b1[1] - cd1 * b1[0] - oo1) < m_Epsilon)
385  {
386  resp = (b1[0] >= xamin && b1[0] <= xamax);
387  }
388  else if (std::abs(b2[1] - cd1 * b2[0] - oo1) < m_Epsilon)
389  {
390  resp = (b2[0] >= xamin && b2[0] <= xamax);
391  }
392  }
393  else
394  {
395  double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
396  double oo1 = a1[1] - cd1 * a1[0];
397  double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
398  double oo2 = b1[1] - cd2 * b1[0];
399  if (std::abs(cd1 - cd2) < m_Epsilon && std::abs(oo1 - oo2) < m_Epsilon)
400  {
401  resp = ((xamin <= xbmax && xamin >= xbmin)
402  || (xamax <= xbmax && xamax >= xbmin)
403  || (xbmin <= xamax && xbmin >= xamin)
404  || (xbmax <= xamax && xbmax >= xamin));
405  }
406  else
407  {
408  if (std::abs(a1[1] - cd2 * a1[0] - oo2) < m_Epsilon)
409  {
410  resp = (a1[0] >= xbmin && a1[0] <= xbmax);
411  }
412  else if (std::abs(a2[1] - cd2 * a2[0] - oo2) < m_Epsilon)
413  {
414  resp = (a2[0] >= xbmin && a2[0] <= xbmax);
415  }
416  if (std::abs(b1[1] - cd1 * b1[0] - oo1) < m_Epsilon)
417  {
418  resp = (b1[0] >= xamin && b1[0] <= xamax);
419  }
420  else if (std::abs(b2[1] - cd1 * b2[0] - oo1) < m_Epsilon)
421  {
422  resp = (b2[0] >= xamin && b2[0] <= xamax);
423  }
424  }
425  }
426  return resp;
427 }
428 
432 template<class TValue>
433 void
436 {
437  VertexListConstIteratorType it = this->GetVertexList()->Begin();
438 
439  if (this->GetVertexList()->Size() > 2)
440  {
441  double area = 0.0;
442  VertexType origin = it.Value();
443  ++it;
444  VertexType pt1 = it.Value();
445  VertexType pt2 = it.Value();
446 
447  while (it != this->GetVertexList()->End())
448  {
449  pt1 = pt2;
450  pt2 = it.Value();
451 
452  double vector1x = pt1[0] - origin[0];
453  double vector1y = pt1[1] - origin[1];
454  double vector2x = pt2[0] - origin[0];
455  double vector2y = pt2[1] - origin[1];
456  double crossProdduct = vector1x * vector2y - vector2x * vector1y;
457  area += crossProdduct;
458  ++it;
459  }
460 
461  m_Area = fabs(area / 2.0);
462 
463  }
464  else //if there is strictly less than 3 points, surface is 0
465  {
466  m_Area = 0.0;
467  }
468 
469  m_AreaIsValid = true;
470 }
471 
475 template<class TValue>
476 double
478 ::GetArea() const
479 {
480  if (!m_AreaIsValid)
481  {
482  ComputeArea();
483  }
484  return m_Area;
485 }
487 
491 template <class TValue>
492 double Polygon<TValue>
493 ::GetLength() const
494 {
495  double length = 0.0;
496  VertexListConstIteratorType it = this->GetVertexList()->Begin();
498 
499  VertexType origin = it.Value();
500  if (this->GetVertexList()->Size() > 1)
501  {
502 
503  VertexType pt1 = it.Value(); //just init, won't be used like that
504  VertexType pt2 = it.Value();
505 
506  ++it;
507  while (it != this->GetVertexList()->End())
508  {
509  pt1 = pt2;
510  pt2 = it.Value();
511  double accum = 0.0;
512  for (int i = 0; i < 2; ++i)
513  {
514  accum += (pt1[i] - pt2[i]) * (pt1[i] - pt2[i]);
515  }
516  length += std::sqrt(accum);
517  ++it;
518  }
519 
520  //Adding the last segment (between first and last point)
521  double accum = 0.0;
522  for (int i = 0; i < 2; ++i)
523  {
524  accum += (origin[i] - pt2[i]) * (origin[i] - pt2[i]);
525  }
526  length += std::sqrt(accum);
527 
528  }
529  else //if there is strictly less than 2 points, length is 0
530  {
531  length = 0.0;
532  }
533 
534  return length;
535 }
536 
537 template <class TValue>
538 void
540 ::Modified() const
541 {
542  Superclass::Modified();
543  m_AreaIsValid = false;
544 }
545 
549 template<class TValue>
550 void
552 ::PrintSelf(std::ostream& os, itk::Indent indent) const
553 {
554  Superclass::PrintSelf(os, indent);
555 }
556 
557 } // End namespace otb
558 
559 #endif
virtual double GetArea() const
Definition: otbPolygon.hxx:478
bool IsCrossing(VertexType a1, VertexType a2, VertexType b1, VertexType b2) const
Definition: otbPolygon.hxx:287
void AddVertex(const ContinuousIndexType &vertex) override
Definition: otbPolygon.hxx:32
void Modified() const override
Definition: otbPolygon.hxx:540
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbPolygon.hxx:552
bool IsOnEdge(VertexType point) const
Definition: otbPolygon.hxx:126
unsigned int NbCrossing(VertexType a, VertexType b) const
Definition: otbPolygon.hxx:215
bool IsInside(VertexType point) const
Definition: otbPolygon.hxx:47
Iterator Begin()
unsigned int NbTouching(VertexType a, VertexType b) const
Definition: otbPolygon.hxx:250
double GetLength() const override
Definition: otbPolygon.hxx:493
Superclass::VertexListConstIteratorType VertexListConstIteratorType
Definition: otbPolygon.h:65
bool IsTouching(VertexType a1, VertexType a2, VertexType b1, VertexType b2) const
Definition: otbPolygon.hxx:345
virtual void ComputeArea() const
Definition: otbPolygon.hxx:435