OTB  9.0.0
Orfeo Toolbox
otbPolygon.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 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>
31 {
32  Superclass::AddVertex(vertex);
33  m_AreaIsValid = false;
34 }
35 
42 template <class TValue>
44 {
45  double x = point[0];
46  double y = point[1];
47  unsigned int crossingCount = 0;
48  VertexListConstIteratorType it = this->GetVertexList()->Begin();
49  double xa = it.Value()[0];
50  double ya = it.Value()[1];
51  ++it;
52  while (it != this->GetVertexList()->End())
53  {
54  double xb = it.Value()[0];
55  double yb = it.Value()[1];
56  if (std::abs(xb - xa) < m_Epsilon)
57  {
58  if (ya > yb && xa > x && y >= yb && y < ya)
59  {
60  ++crossingCount;
61  }
62  else if (ya < yb && xa > x && y >= ya && y < yb)
63  {
64  ++crossingCount;
65  }
66  }
67  else if (std::abs(yb - ya) >= m_Epsilon)
68  {
69  double xcross = xa + (xb - xa) * (y - ya) / (yb - ya);
70 
71  if (ya > yb && xcross > x && y >= yb && y < ya)
72  {
73  ++crossingCount;
74  }
75  else if (ya < yb && xcross > x && y >= ya && y < yb)
76  {
77  ++crossingCount;
78  }
79  }
80  xa = xb;
81  ya = yb;
82  ++it;
83  }
84  double xb = this->GetVertexList()->Begin().Value()[0];
85  double yb = this->GetVertexList()->Begin().Value()[1];
86  if (std::abs(xb - xa) < m_Epsilon)
87  {
88  if (ya > yb && xa > x && y >= yb && y < ya)
89  {
90  ++crossingCount;
91  }
92  else if (ya < yb && xa > x && y >= ya && y < yb)
93  {
94  ++crossingCount;
95  }
96  }
97  else if (std::abs(yb - ya) >= m_Epsilon)
98  {
99  double xcross = xa + (xb - xa) * (y - ya) / (yb - ya);
100 
101  if (ya > yb && xcross > x && y >= yb && y < ya)
102  {
103  ++crossingCount;
104  }
105  else if (ya < yb && xcross > x && y >= ya && y < yb)
106  {
107  ++crossingCount;
108  }
109  }
110  // std::cout<<"Crossing count: "<<crossingCount<<std::endl;
111  return (crossingCount % 2 == 1);
112 }
113 
119 template <class TValue>
121 {
122  // std::cout<<"Checking point "<<point<<std::endl;
123  bool resp = false;
124  double x = point[0];
125  double y = point[1];
126  double xb, yb;
127  VertexListConstIteratorType it = this->GetVertexList()->Begin();
128  double xa = it.Value()[0];
129  double ya = it.Value()[1];
130  double xbegin = xa;
131  double ybegin = ya;
132  ++it;
133  while (!resp && it != this->GetVertexList()->End())
134  {
135  xb = it.Value()[0];
136  yb = it.Value()[1];
137  if (std::abs(xb - xa) >= m_Epsilon)
138  {
139  double cd = (yb - ya) / (xb - xa);
140  double oo = (ya - cd * xa);
141  double xmin = std::min(xa, xb);
142  double xmax = std::max(xa, xb);
143  if ((std::abs(y - cd * x - oo) < m_Epsilon) && (x <= xmax) && (x >= xmin))
144  {
145  // std::cout<<"Case 1: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
146  resp = true;
147  }
148  }
149  else
150  {
151  double ymin = std::min(ya, yb);
152  double ymax = std::max(ya, yb);
153  if ((std::abs(x - xa) < m_Epsilon) && (y <= ymax) && (y >= ymin))
154  {
155  resp = true;
156  // std::cout<<"Case 2: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
157  }
158  }
159  xa = xb;
160  ya = yb;
161  ++it;
162  }
163  xb = xbegin;
164  yb = ybegin;
165  if (std::abs(xb - xa) >= m_Epsilon)
166  {
167  double cd = (yb - ya) / (xb - xa);
168  double oo = (ya - cd * xa);
169  double xmin = std::min(xa, xb);
170  double xmax = std::max(xa, xb);
172 
173  if ((std::abs(y - cd * x - oo) < m_Epsilon) && (x <= xmax) && (x >= xmin))
174  {
175  resp = true;
176  // std::cout<<"Case 1: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
177  }
178  }
179  else
180  {
181  double ymin = std::min(ya, yb);
182  double ymax = std::max(ya, yb);
183  if ((std::abs(x - xa) <= m_Epsilon) && (y <= ymax) && (y >= ymin))
184  {
185  resp = true;
186  // std::cout<<"Case 2: Point is on segment "<<xa<<", "<<ya<<" <-> "<<xb<<", "<<yb<<std::endl;
187  }
188  }
189  // std::cout<<"Result: "<<resp<<std::endl;
190  return resp;
191 }
198 template <class TValue>
200 {
201  unsigned int resp = 0;
202  VertexListConstIteratorType it = this->GetVertexList()->Begin();
203  VertexListConstIteratorType it_end = this->GetVertexList()->End();
204  VertexType current = it.Value();
205  VertexType first = current;
206  ++it;
207  while (it != it_end)
208  {
209  // std::cout<<"Testing Is crossing "<<a<<" "<<b<<current<<it.Value()<<" = "<<IsCrossing(a, b, current, it.Value())<<std::endl;
210  if (IsCrossing(a, b, current, it.Value()))
211  {
212  ++resp;
213  }
214  current = it.Value();
215  ++it;
216  }
217  // std::cout<<"Testing Is crossing "<<a<<" "<<b<<current<<first<<" = "<<IsCrossing(a, b, current, first)<<std::endl;
218  if (IsCrossing(a, b, current, first))
219  {
220  ++resp;
221  }
222  return resp;
223 }
224 
231 template <class TValue>
233 {
234  unsigned int resp = 0;
235  VertexListConstIteratorType it = this->GetVertexList()->Begin();
236  VertexListConstIteratorType it_end = this->GetVertexList()->End();
237  VertexType current = it.Value();
238  VertexType first = current;
239  ++it;
240  while (it != it_end)
241  {
242  // std::cout<<"IsTouching "<<a<<" "<<b<<", "<<current<<" "<<it.Value()<<" -> "<<IsTouching(a, b, current, it.Value())<<std::endl;
243  if (IsTouching(a, b, current, it.Value()))
244  {
245  ++resp;
246  }
247  current = it.Value();
248  ++it;
249  }
250  // std::cout<<"IsTouching "<<a<<" "<<b<<", "<<current<<" "<<first<<" -> "<<IsTouching(a, b, current, first)<<std::endl;
251  if (IsTouching(a, b, current, first))
252  {
253  ++resp;
254  }
255  return resp;
256 }
257 
266 template <class TValue>
268 {
269  bool resp = false;
270  double xbmin = std::min(b1[0], b2[0]);
271  double xbmax = std::max(b1[0], b2[0]);
272  double ybmin = std::min(b1[1], b2[1]);
273  double ybmax = std::max(b1[1], b2[1]);
274  double xamin = std::min(a1[0], a2[0]);
275  double xamax = std::max(a1[0], a2[0]);
276  double yamin = std::min(a1[1], a2[1]);
277  double yamax = std::max(a1[1], a2[1]);
278  if (std::abs(a1[0] - a2[0]) < m_Epsilon && std::abs(b1[0] - b2[0]) < m_Epsilon)
279  {
280  resp = false;
281  }
282  else if (std::abs(a1[0] - a2[0]) < m_Epsilon)
283  {
284  double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
285  double oo2 = b1[1] - cd2 * b1[0];
286  double ycross = cd2 * a1[0] + oo2;
287  resp = (xbmin < a1[0] && xbmax > a1[0] && yamin < ycross && yamax > ycross);
288  }
289  else if (std::abs(b1[0] - b2[0]) < m_Epsilon)
290  {
291  double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
292  double oo1 = a1[1] - cd1 * a1[0];
293  double ycross = cd1 * b1[0] + oo1;
294  resp = (xamin < b1[0] && xamax > b1[0] && ybmin < ycross && ybmax > ycross);
295  }
296  else
297  {
298  double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
299  double oo1 = a1[1] - cd1 * a1[0];
300  double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
301  double oo2 = b1[1] - cd2 * b1[0];
302  if (cd1 != cd2)
303  {
304  double xcross = (oo2 - oo1) / (cd1 - cd2);
305  resp = (xamin < xcross && xbmin < xcross && xamax > xcross && xbmax > xcross);
306  }
307  }
308  return resp;
309 }
310 
319 template <class TValue>
321 {
322  bool resp = false;
323  double xbmin = std::min(b1[0], b2[0]);
324  double xbmax = std::max(b1[0], b2[0]);
325  double ybmin = std::min(b1[1], b2[1]);
326  double ybmax = std::max(b1[1], b2[1]);
327  double xamin = std::min(a1[0], a2[0]);
328  double xamax = std::max(a1[0], a2[0]);
329  double yamin = std::min(a1[1], a2[1]);
330  double yamax = std::max(a1[1], a2[1]);
331  if (std::abs(a1[0] - a2[0]) < m_Epsilon && std::abs(b1[0] - b2[0]) < m_Epsilon)
332  {
333  resp = (std::abs(a1[0] - b1[0]) < m_Epsilon) && ((a1[1] <= ybmax && a1[1] >= ybmin) || (a2[1] <= ybmax && a2[1] >= ybmin) ||
334  (b1[1] <= yamax && b1[1] >= yamin) || (b2[1] <= yamax && b2[1] >= yamin));
335  }
336  else if (std::abs(a1[0] - a2[0]) < m_Epsilon)
337  {
338  double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
339  double oo2 = b1[1] - cd2 * b1[0];
341 
342  if (std::abs(a1[1] - cd2 * a1[0] - oo2) < m_Epsilon)
343  {
344  resp = (a1[0] >= xbmin && a1[0] <= xbmax);
345  }
346  else if (std::abs(a2[1] - cd2 * a2[0] - oo2) < m_Epsilon)
347  {
348  resp = (a2[0] >= xbmin && a2[0] <= xbmax);
349  }
350  }
351  else if (std::abs(b1[0] - b2[0]) < m_Epsilon)
352  {
353  double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
354  double oo1 = a1[1] - cd1 * a1[0];
355 
356  if (std::abs(b1[1] - cd1 * b1[0] - oo1) < m_Epsilon)
357  {
358  resp = (b1[0] >= xamin && b1[0] <= xamax);
359  }
360  else if (std::abs(b2[1] - cd1 * b2[0] - oo1) < m_Epsilon)
361  {
362  resp = (b2[0] >= xamin && b2[0] <= xamax);
363  }
364  }
365  else
366  {
367  double cd1 = (a2[1] - a1[1]) / (a2[0] - a1[0]);
368  double oo1 = a1[1] - cd1 * a1[0];
369  double cd2 = (b2[1] - b1[1]) / (b2[0] - b1[0]);
370  double oo2 = b1[1] - cd2 * b1[0];
371  if (std::abs(cd1 - cd2) < m_Epsilon && std::abs(oo1 - oo2) < m_Epsilon)
372  {
373  resp = ((xamin <= xbmax && xamin >= xbmin) || (xamax <= xbmax && xamax >= xbmin) || (xbmin <= xamax && xbmin >= xamin) ||
374  (xbmax <= xamax && xbmax >= xamin));
375  }
376  else
377  {
378  if (std::abs(a1[1] - cd2 * a1[0] - oo2) < m_Epsilon)
379  {
380  resp = (a1[0] >= xbmin && a1[0] <= xbmax);
381  }
382  else if (std::abs(a2[1] - cd2 * a2[0] - oo2) < m_Epsilon)
383  {
384  resp = (a2[0] >= xbmin && a2[0] <= xbmax);
385  }
386  if (std::abs(b1[1] - cd1 * b1[0] - oo1) < m_Epsilon)
387  {
388  resp = (b1[0] >= xamin && b1[0] <= xamax);
389  }
390  else if (std::abs(b2[1] - cd1 * b2[0] - oo1) < m_Epsilon)
391  {
392  resp = (b2[0] >= xamin && b2[0] <= xamax);
393  }
394  }
395  }
396  return resp;
397 }
398 
402 template <class TValue>
404 {
405  VertexListConstIteratorType it = this->GetVertexList()->Begin();
406 
407  if (this->GetVertexList()->Size() > 2)
408  {
409  double area = 0.0;
410  VertexType origin = it.Value();
411  ++it;
412  VertexType pt1 = it.Value();
413  VertexType pt2 = it.Value();
414 
415  while (it != this->GetVertexList()->End())
416  {
417  pt1 = pt2;
418  pt2 = it.Value();
419 
420  double vector1x = pt1[0] - origin[0];
421  double vector1y = pt1[1] - origin[1];
422  double vector2x = pt2[0] - origin[0];
423  double vector2y = pt2[1] - origin[1];
424  double crossProdduct = vector1x * vector2y - vector2x * vector1y;
425  area += crossProdduct;
426  ++it;
427  }
428 
429  m_Area = fabs(area / 2.0);
430  }
431  else // if there is strictly less than 3 points, surface is 0
432  {
433  m_Area = 0.0;
434  }
435 
436  m_AreaIsValid = true;
437 }
438 
442 template <class TValue>
444 {
445  if (!m_AreaIsValid)
446  {
447  ComputeArea();
448  }
449  return m_Area;
450 }
452 
456 template <class TValue>
458 {
459  double length = 0.0;
460  VertexListConstIteratorType it = this->GetVertexList()->Begin();
462 
463  VertexType origin = it.Value();
464  if (this->GetVertexList()->Size() > 1)
465  {
466 
467  VertexType pt1 = it.Value(); // just init, won't be used like that
468  VertexType pt2 = it.Value();
469 
470  ++it;
471  while (it != this->GetVertexList()->End())
472  {
473  pt1 = pt2;
474  pt2 = it.Value();
475  double accum = 0.0;
476  for (int i = 0; i < 2; ++i)
477  {
478  accum += (pt1[i] - pt2[i]) * (pt1[i] - pt2[i]);
479  }
480  length += std::sqrt(accum);
481  ++it;
482  }
483 
484  // Adding the last segment (between first and last point)
485  double accum = 0.0;
486  for (int i = 0; i < 2; ++i)
487  {
488  accum += (origin[i] - pt2[i]) * (origin[i] - pt2[i]);
489  }
490  length += std::sqrt(accum);
491  }
492  else // if there is strictly less than 2 points, length is 0
493  {
494  length = 0.0;
495  }
496 
497  return length;
498 }
499 
500 template <class TValue>
502 {
503  Superclass::Modified();
504  m_AreaIsValid = false;
505 }
506 
510 template <class TValue>
511 void Polygon<TValue>::PrintSelf(std::ostream& os, itk::Indent indent) const
512 {
513  Superclass::PrintSelf(os, indent);
514 }
515 
516 } // End namespace otb
517 
518 #endif
otb::Polygon::ComputeArea
virtual void ComputeArea() const
Definition: otbPolygon.hxx:403
otbPolygon.h
otb::Polygon::Modified
void Modified() const override
Definition: otbPolygon.hxx:501
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::Polygon::IsInside
bool IsInside(VertexType point) const
Definition: otbPolygon.hxx:43
otb::Polygon::NbCrossing
unsigned int NbCrossing(VertexType a, VertexType b) const
Definition: otbPolygon.hxx:199
otb::Polygon::IsTouching
bool IsTouching(VertexType a1, VertexType a2, VertexType b1, VertexType b2) const
Definition: otbPolygon.hxx:320
otb::Polygon::GetLength
double GetLength() const override
Definition: otbPolygon.hxx:457
otb::Polygon::AddVertex
void AddVertex(const ContinuousIndexType &vertex) override
Definition: otbPolygon.hxx:30
otb::Polygon::ContinuousIndexType
Superclass::ContinuousIndexType ContinuousIndexType
Definition: otbPolygon.h:63
otb::Polygon::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbPolygon.hxx:511
otb::Polygon::VertexType
Superclass::VertexType VertexType
Definition: otbPolygon.h:58
otb::Polygon::IsOnEdge
bool IsOnEdge(VertexType point) const
Definition: otbPolygon.hxx:120
otb::Polygon::GetArea
virtual double GetArea() const
Definition: otbPolygon.hxx:443
otb::Polygon::IsCrossing
bool IsCrossing(VertexType a1, VertexType a2, VertexType b1, VertexType b2) const
Definition: otbPolygon.hxx:267
otb::Polygon::NbTouching
unsigned int NbTouching(VertexType a, VertexType b) const
Definition: otbPolygon.hxx:232
otb::Polygon::VertexListConstIteratorType
Superclass::VertexListConstIteratorType VertexListConstIteratorType
Definition: otbPolygon.h:64