OGS 6.2.1-97-g73d1aeda3
Polyline.cpp
Go to the documentation of this file.
1 
15 #include "Polyline.h"
16 
17 #include <algorithm>
18 #include <logog/include/logog.hpp>
19 #include "BaseLib/Error.h"
20 #include "AnalyticalGeometry.h"
22 
23 namespace GeoLib
24 {
25 Polyline::Polyline(const std::vector<Point*>& pnt_vec) : _ply_pnts(pnt_vec)
26 {
27  _length.push_back (0.0);
28 }
29 
31  : _ply_pnts(ply._ply_pnts),
33  _length(ply._length)
34 {}
35 
36 void Polyline::write(std::ostream &os) const
37 {
38  std::size_t size(_ply_pnt_ids.size());
39  for (std::size_t k(0); k < size; k++)
40  {
41  os << *(_ply_pnts[_ply_pnt_ids[k]]) << "\n";
42  }
43 }
44 
45 bool Polyline::addPoint(std::size_t pnt_id)
46 {
47  assert(pnt_id < _ply_pnts.size());
48  std::size_t const n_pnts(_ply_pnt_ids.size());
49 
50  // don't insert point if this would result in identical IDs for two adjacent points
51  if (n_pnts > 0 && _ply_pnt_ids.back() == pnt_id)
52  {
53  return false;
54  }
55 
56  _ply_pnt_ids.push_back(pnt_id);
57 
58  if (n_pnts > 0)
59  {
60  double const act_dist(std::sqrt(MathLib::sqrDist(
61  *_ply_pnts[_ply_pnt_ids[n_pnts-1]], *_ply_pnts[pnt_id])));
62  double dist_until_now(0.0);
63  if (n_pnts > 1)
64  {
65  dist_until_now = _length[n_pnts - 1];
66  }
67 
68  _length.push_back(dist_until_now + act_dist);
69  }
70  return true;
71 }
72 
73 bool Polyline::insertPoint(std::size_t pos, std::size_t pnt_id)
74 {
75  assert(pnt_id < _ply_pnts.size());
76  assert(pos <= _ply_pnt_ids.size());
77 
78  if (pos == _ply_pnt_ids.size()) {
79  return addPoint(pnt_id);
80  }
81 
82  // check if inserting pnt_id would result in two identical IDs for adjacent points
83  if (pos == 0 && pnt_id == _ply_pnt_ids[0]) {
84  return false;
85  }
86  if (pos != 0)
87  {
88  if (pos == (_ply_pnt_ids.size() - 1) && pnt_id == _ply_pnt_ids[pos])
89  {
90  return false;
91  }
92  if (pnt_id == _ply_pnt_ids[pos - 1] || pnt_id == _ply_pnt_ids[pos])
93  {
94  return false;
95  }
96  }
97 
98  auto const pos_dt(
99  static_cast<std::vector<std::size_t>::difference_type>(pos));
100  auto it(_ply_pnt_ids.begin() + pos_dt);
101  _ply_pnt_ids.insert(it, pnt_id);
102 
103  if (_ply_pnt_ids.size() > 1) {
104  // update the _length vector
105  if (pos == 0) {
106  // insert at first position
107  double const act_dist(std::sqrt(MathLib::sqrDist(
108  *_ply_pnts[_ply_pnt_ids[1]], *_ply_pnts[pnt_id])));
109  _length.insert(_length.begin() + 1, act_dist);
110  const std::size_t s(_length.size());
111  for (std::size_t k(2); k < s; k++)
112  {
113  _length[k] += _length[1];
114  }
115  } else {
116  if (pos == _ply_pnt_ids.size() - 1) {
117  // insert at last position
118  double const act_dist(std::sqrt(MathLib::sqrDist(
119  *_ply_pnts[_ply_pnt_ids[_ply_pnt_ids.size() - 2]],
120  *_ply_pnts[pnt_id])));
121  double dist_until_now (0.0);
122  if (_ply_pnt_ids.size() > 2)
123  {
124  dist_until_now = _length[_ply_pnt_ids.size() - 2];
125  }
126 
127  _length.insert(_length.begin() + pos_dt,
128  dist_until_now + act_dist);
129  } else {
130  // insert at arbitrary position within the vector
131  double dist_until_now (0.0);
132  if (pos > 1)
133  {
134  dist_until_now = _length[pos - 1];
135  }
136  double len_seg0(std::sqrt(MathLib::sqrDist(
137  *_ply_pnts[_ply_pnt_ids[pos - 1]],
138  *_ply_pnts[pnt_id])));
139  double len_seg1(std::sqrt(MathLib::sqrDist(
140  *_ply_pnts[_ply_pnt_ids[pos + 1]],
141  *_ply_pnts[pnt_id])));
142  double update_dist(
143  len_seg0 + len_seg1 - (_length[pos] - dist_until_now));
144  _length[pos] = dist_until_now + len_seg0;
145  auto it1(_length.begin() + pos_dt + 1);
146  _length.insert(it1, _length[pos] + len_seg1);
147  for (it1 = _length.begin() + pos_dt + 2; it1 != _length.end();
148  ++it1)
149  {
150  *it1 += update_dist;
151  }
152  }
153  }
154  }
155  return true;
156 }
157 
158 void Polyline::removePoint(std::size_t pos)
159 {
160  if (pos >= _ply_pnt_ids.size())
161  {
162  return;
163  }
164 
165  auto const pos_dt(
166  static_cast<std::vector<std::size_t>::difference_type>(pos));
167  _ply_pnt_ids.erase(_ply_pnt_ids.begin() + pos_dt);
168 
169  if (pos == _ply_pnt_ids.size())
170  {
171  _length.erase(_length.begin() + pos_dt);
172  return;
173  }
174 
175  const std::size_t n_ply_pnt_ids(_ply_pnt_ids.size());
176  if (pos == 0) {
177  double seg_length(_length[0]);
178  for (std::size_t k(0); k < n_ply_pnt_ids; k++)
179  {
180  _length[k] = _length[k + 1] - seg_length;
181  }
182  _length.pop_back();
183  } else {
184  const double len_seg0(_length[pos] - _length[pos - 1]);
185  const double len_seg1(_length[pos + 1] - _length[pos]);
186  _length.erase(_length.begin() + pos_dt);
187  const double len_new_seg(std::sqrt(MathLib::sqrDist(*_ply_pnts[_ply_pnt_ids[pos - 1]],
188  *_ply_pnts[_ply_pnt_ids[pos]])));
189  double seg_length_diff(len_new_seg - len_seg0 - len_seg1);
190 
191  for (std::size_t k(pos); k < n_ply_pnt_ids; k++)
192  {
193  _length[k] += seg_length_diff;
194  }
195  }
196 }
197 
198 std::size_t Polyline::getNumberOfPoints() const
199 {
200  return _ply_pnt_ids.size();
201 }
202 
203 std::size_t Polyline::getNumberOfSegments() const
204 {
205  return _ply_pnt_ids.empty() ? 0 : _ply_pnt_ids.size()-1;
206 }
207 
208 bool Polyline::isClosed() const
209 {
210  if (_ply_pnt_ids.size() < 3)
211  {
212  return false;
213  }
214 
215  return _ply_pnt_ids.front() == _ply_pnt_ids.back();
216 }
217 
219 {
220  std::size_t const n_points (_ply_pnt_ids.size());
221  if (n_points < 4)
222  {
223  return true;
224  }
225 
226  GeoLib::Point const& p0 (*this->getPoint(0));
227  GeoLib::Point const& p1 (*this->getPoint(1));
228  GeoLib::Point const& p2 (*this->getPoint(2));
229  for (std::size_t i=3; i<n_points; ++i)
230  {
231  if (!MathLib::isCoplanar(p0, p1, p2, *this->getPoint(i)))
232  {
233  DBUG ("Point %d is not coplanar to the first three points of the line.", i);
234  return false;
235  }
236  }
237  return true;
238 }
239 
240 bool Polyline::isPointIDInPolyline(std::size_t pnt_id) const
241 {
242  return std::find(_ply_pnt_ids.begin(), _ply_pnt_ids.end(), pnt_id) != _ply_pnt_ids.end();
243 }
244 
245 std::size_t Polyline::getPointID(std::size_t i) const
246 {
247  assert(i < _ply_pnt_ids.size());
248  return _ply_pnt_ids[i];
249 }
250 
251 LineSegment Polyline::getSegment(std::size_t i) const
252 {
253  assert(i < getNumberOfSegments());
255  _ply_pnts[_ply_pnt_ids[i + 1]], false);
256 }
257 
259 {
260  assert(i < getNumberOfSegments());
262  _ply_pnts[_ply_pnt_ids[i + 1]], false);
263 }
264 
265 void Polyline::setPointID(std::size_t idx, std::size_t id)
266 {
267  assert(idx < _ply_pnt_ids.size());
268  _ply_pnt_ids[idx] = id;
269 }
270 
271 const Point* Polyline::getPoint(std::size_t i) const
272 {
273  assert(i < _ply_pnt_ids.size());
274  return _ply_pnts[_ply_pnt_ids[i]];
275 }
276 
277 std::vector<Point*> const& Polyline::getPointsVec () const
278 {
279  return _ply_pnts;
280 }
281 
282 double Polyline::getLength (std::size_t k) const
283 {
284  assert(k < _length.size());
285  return _length[k];
286 }
287 
288 Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> &ply_vec,
289  double prox)
290 {
291  std::size_t nLines = ply_vec.size();
292 
293  auto* new_ply = new Polyline(*ply_vec[0]);
294  std::vector<GeoLib::Point*> pnt_vec(new_ply->getPointsVec());
295 
296  std::vector<Polyline*> local_ply_vec;
297  for (std::size_t i = 1; i < nLines; i++)
298  {
299  local_ply_vec.push_back(ply_vec[i]);
300  }
301 
302  while (!local_ply_vec.empty())
303  {
304  bool ply_found(false);
305  prox *= prox; // square distance once to save time later
306  for (auto it = local_ply_vec.begin(); it != local_ply_vec.end(); ++it)
307  {
308  if (pnt_vec == (*it)->getPointsVec())
309  {
310  std::size_t nPoints((*it)->getNumberOfPoints());
311 
312  //if (new_ply->getPointID(0) == (*it)->getPointID(0))
313  if (pointsAreIdentical(pnt_vec, new_ply->getPointID(0),
314  (*it)->getPointID(0), prox))
315  {
316  auto* tmp = new Polyline((*it)->getPointsVec());
317  for (std::size_t k = 0; k < nPoints; k++)
318  {
319  tmp->addPoint((*it)->getPointID(nPoints - k - 1));
320  }
321 
322  std::size_t new_ply_size(new_ply->getNumberOfPoints());
323  for (std::size_t k = 1; k < new_ply_size; k++)
324  {
325  tmp->addPoint(new_ply->getPointID(k));
326  }
327  delete new_ply;
328  new_ply = tmp;
329  ply_found = true;
330  }
331  //else if (new_ply->getPointID(0) == (*it)->getPointID(nPoints-1))
332  else if (pointsAreIdentical(pnt_vec, new_ply->getPointID(0),
333  (*it)->getPointID(nPoints - 1), prox))
334  {
335  auto* tmp = new Polyline(**it);
336  std::size_t new_ply_size(new_ply->getNumberOfPoints());
337  for (std::size_t k = 1; k < new_ply_size; k++)
338  {
339  tmp->addPoint(new_ply->getPointID(k));
340  }
341  delete new_ply;
342  new_ply = tmp;
343  ply_found = true;
344  }
345  //else if (new_ply->getPointID(new_ply->getNumberOfPoints()-1) == (*it)->getPointID(0))
346  else if (pointsAreIdentical(pnt_vec,
347  new_ply->getPointID(new_ply->
349  - 1),
350  (*it)->getPointID(0), prox))
351  {
352  for (std::size_t k = 1; k < nPoints; k++)
353  {
354  new_ply->addPoint((*it)->getPointID(k));
355  }
356  ply_found = true;
357  }
358  //else if (new_ply->getPointID(new_ply->getNumberOfPoints()-1) == (*it)->getPointID(nPoints-1))
359  else if (pointsAreIdentical(pnt_vec,
360  new_ply->getPointID(new_ply->
362  - 1),
363  (*it)->getPointID(nPoints - 1), prox))
364  {
365  for (std::size_t k = 1; k < nPoints; k++)
366  {
367  new_ply->addPoint((*it)->getPointID(nPoints - k - 1));
368  }
369  ply_found = true;
370  }
371  if (ply_found)
372  {
373  local_ply_vec.erase(it);
374  break;
375  }
376  }
377  else
378  ERR("Error in Polyline::contructPolylineFromSegments() - Line segments use different point vectors.");
379  }
380 
381  if (!ply_found)
382  {
383  ERR("Error in Polyline::contructPolylineFromSegments() - Not all segments are connected.");
384  delete new_ply;
385  new_ply = nullptr;
386  break;
387  }
388  }
389  return new_ply;
390 }
391 
393 {
394  if (getNumberOfPoints() < 2) {
395  ERR("Polyline::closePolyline(): Input polyline needs to be composed of at least three points.");
396  }
397  if (! isClosed()) {
398  addPoint(getPointID(0));
399  }
400 }
401 
402 Location Polyline::getLocationOfPoint (std::size_t k, GeoLib::Point const & pnt) const
403 {
404  assert (k < _ply_pnt_ids.size() - 1);
405 
406  GeoLib::Point const& source (*(_ply_pnts[_ply_pnt_ids[k]]));
407  GeoLib::Point const& dest (*(_ply_pnts[_ply_pnt_ids[k + 1]]));
408  long double a[2] = {dest[0] - source[0], dest[1] - source[1]}; // vector
409  long double b[2] = {pnt[0] - source[0], pnt[1] - source[1]}; // vector
410 
411  long double det_2x2 (a[0] * b[1] - a[1] * b[0]);
412 
413  if (det_2x2 > std::numeric_limits<double>::epsilon())
414  {
415  return Location::LEFT;
416  }
417  if (std::numeric_limits<double>::epsilon() < std::abs(det_2x2))
418  {
419  return Location::RIGHT;
420  }
421  if (a[0] * b[0] < 0.0 || a[1] * b[1] < 0.0)
422  {
423  return Location::BEHIND;
424  }
425  if (a[0] * a[0] + a[1] * a[1] < b[0] * b[0] + b[1] * b[1])
426  {
427  return Location::BEYOND;
428  }
429  if (MathLib::sqrDist(pnt, *_ply_pnts[_ply_pnt_ids[k]]) <
430  pow(std::numeric_limits<double>::epsilon(), 2))
431  {
432  return Location::SOURCE;
433  }
434  if (MathLib::sqrDist(pnt, *_ply_pnts[_ply_pnt_ids[k + 1]]) <
435  std::sqrt(std::numeric_limits<double>::epsilon()))
436  {
437  return Location::DESTINATION;
438  }
439  return Location::BETWEEN;
440 }
441 
443  const double epsilon_radius) const
444 {
445  double dist(-1.0);
446  double lambda;
447  bool found = false;
448  // loop over all line segments of the polyline
449  for (std::size_t k = 0; k < getNumberOfSegments(); k++) {
450  // is the orthogonal projection of the j-th node to the
451  // line g(lambda) = _ply->getPoint(k) + lambda * (_ply->getPoint(k+1) - _ply->getPoint(k))
452  // at the k-th line segment of the polyline, i.e. 0 <= lambda <= 1?
454  (getPoint(k))->getCoords(), (getPoint(k + 1))->getCoords(),
455  lambda, dist) <= epsilon_radius) {
456 
457  double act_length_of_ply(getLength(k));
458  double seg_length (getLength(k+1)-getLength(k));
459  double lower_lambda (- epsilon_radius / seg_length);
460  double upper_lambda (1 + epsilon_radius / seg_length);
461 
462  if (lower_lambda <= lambda && lambda <= upper_lambda) {
463  found = true;
464  dist = act_length_of_ply + dist;
465  break;
466  } // end if lambda
467  }
468  } // end line segment loop
469 
470  if (!found)
471  {
472  dist = -1.0;
473  }
474  return dist;
475 }
476 
478  std::size_t segment_number)
479  : _polyline(&polyline),
480  _segment_number(
481  static_cast<std::vector<GeoLib::Point*>::size_type>(segment_number))
482 {}
483 
486 {}
487 
489  SegmentIterator const& rhs)
490 {
491  if (&rhs == this)
492  {
493  return *this;
494  }
495 
496  _polyline = rhs._polyline;
498  return *this;
499 }
500 
502 {
503  return static_cast<std::size_t>(_segment_number);
504 }
505 
507 {
508  ++_segment_number;
509  return *this;
510 }
511 
513 {
515 }
516 
518 {
520 }
521 
523 {
524  return !(*this != other);
525 }
526 
528 {
529  return other._segment_number != _segment_number ||
530  other._polyline != _polyline;
531 }
532 
534  std::vector<GeoLib::Point>::difference_type n)
535 {
536  if (n < 0) {
537  _segment_number -=
538  static_cast<std::vector<GeoLib::Point>::size_type>(-n);
539  } else {
540  _segment_number +=
541  static_cast<std::vector<GeoLib::Point>::size_type>(n);
542  }
544  {
545  OGS_FATAL("");
546  }
547  return *this;
548 }
549 
551  std::vector<GeoLib::Point>::difference_type n)
552 {
553  SegmentIterator t(*this);
554  t += n;
555  return t;
556 }
557 
559  std::vector<GeoLib::Point>::difference_type n)
560 {
561  if (n >= 0) {
562  _segment_number -=
563  static_cast<std::vector<GeoLib::Point>::size_type>(n);
564  } else {
565  _segment_number +=
566  static_cast<std::vector<GeoLib::Point>::size_type>(-n);
567  }
569  {
570  OGS_FATAL("");
571  }
572  return *this;
573 }
574 
576  std::vector<GeoLib::Point>::difference_type n)
577 {
578  Polyline::SegmentIterator t(*this);
579  t -= n;
580  return t;
581 }
582 
583 std::ostream& operator<< (std::ostream &os, const Polyline &pl)
584 {
585  pl.write (os);
586  return os;
587 }
588 
589 bool containsEdge (const Polyline& ply, std::size_t id0, std::size_t id1)
590 {
591  if (id0 == id1)
592  {
593  ERR("no valid edge id0 == id1 == %d.", id0);
594  return false;
595  }
596  if (id0 > id1)
597  {
598  std::swap(id0, id1);
599  }
600  const std::size_t n (ply.getNumberOfPoints() - 1);
601  for (std::size_t k(0); k < n; k++)
602  {
603  std::size_t ply_pnt0 (ply.getPointID (k));
604  std::size_t ply_pnt1 (ply.getPointID (k + 1));
605  if (ply_pnt0 > ply_pnt1)
606  {
607  std::swap(ply_pnt0, ply_pnt1);
608  }
609  if (ply_pnt0 == id0 && ply_pnt1 == id1)
610  {
611  return true;
612  }
613  }
614  return false;
615 }
616 
617 bool operator==(Polyline const& lhs, Polyline const& rhs)
618 {
619  if (lhs.getNumberOfPoints() != rhs.getNumberOfPoints())
620  {
621  return false;
622  }
623 
624  const std::size_t n(lhs.getNumberOfPoints());
625  for (std::size_t k(0); k < n; k++)
626  {
627  if (lhs.getPointID(k) != rhs.getPointID(k))
628  {
629  return false;
630  }
631  }
632 
633  return true;
634 }
635 
636 bool pointsAreIdentical(const std::vector<Point*> &pnt_vec,
637  std::size_t i,
638  std::size_t j,
639  double prox)
640 {
641  if (i == j)
642  {
643  return true;
644  }
645  return MathLib::sqrDist(*pnt_vec[i], *pnt_vec[j]) < prox;
646 }
647 } // end namespace GeoLib
double getLength(std::size_t k) const
Definition: Polyline.cpp:282
LineSegment operator*() const
Definition: Polyline.cpp:512
Definition of the PolyLine class.
double sqrDist(const double *p0, const double *p1)
Definition: MathTools.h:83
std::vector< std::size_t > _ply_pnt_ids
Definition: Polyline.h:236
std::size_t getNumberOfSegments() const
Definition: Polyline.cpp:203
GeoLib::Polyline const * _polyline
Definition: Polyline.h:88
bool isClosed() const
Definition: Polyline.cpp:208
void setPointID(std::size_t idx, std::size_t id)
Definition: Polyline.cpp:265
double getDistanceAlongPolyline(const MathLib::Point3d &pnt, const double epsilon_radius) const
Definition: Polyline.cpp:442
const Point * getPoint(std::size_t i) const
returns the i-th point contained in the polyline
Definition: Polyline.cpp:271
Definition of analytical geometry functions.
virtual void removePoint(std::size_t pos)
Definition: Polyline.cpp:158
Location
Definition: Polyline.h:30
std::vector< double > _length
Definition: Polyline.h:240
static Polyline * constructPolylineFromSegments(const std::vector< Polyline *> &ply_vec, double prox=0.0)
Definition: Polyline.cpp:288
std::size_t getNumberOfPoints() const
Definition: Polyline.cpp:198
std::vector< Point * > const & getPointsVec() const
Definition: Polyline.cpp:277
Definition of the GEOObjects class.
Definition: BaseItem.h:20
void write(std::ostream &os) const
Definition: Polyline.cpp:36
bool isCoplanar(const MathLib::Point3d &a, const MathLib::Point3d &b, const MathLib::Point3d &c, const MathLib::Point3d &d)
Checks if the four given points are located on a plane.
std::size_t getPointID(std::size_t i) const
Definition: Polyline.cpp:245
std::ostream & operator<<(std::ostream &os, LineSegment const &s)
Definition: LineSegment.cpp:82
const std::vector< Point * > & _ply_pnts
Definition: Polyline.h:234
const T * getCoords() const
Definition: TemplatePoint.h:77
virtual bool insertPoint(std::size_t pos, std::size_t pnt_id)
Definition: Polyline.cpp:73
Polyline(const std::vector< Point *> &pnt_vec)
Definition: Polyline.cpp:25
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition: Polyline.h:50
SegmentIterator operator-(std::vector< GeoLib::Point >::difference_type n)
Definition: Polyline.cpp:575
double calcProjPntToLineAndDists(const double p[3], const double a[3], const double b[3], double &lambda, double &d0)
Definition: MathTools.cpp:18
Location getLocationOfPoint(std::size_t k, GeoLib::Point const &pnt) const
Definition: Polyline.cpp:402
bool isPointIDInPolyline(std::size_t pnt_id) const
Definition: Polyline.cpp:240
std::vector< GeoLib::Point * >::size_type _segment_number
Definition: Polyline.h:89
virtual bool addPoint(std::size_t pnt_id)
Definition: Polyline.cpp:45
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
bool containsEdge(const Polyline &ply, std::size_t id0, std::size_t id1)
Definition: Polyline.cpp:589
SegmentIterator operator+(std::vector< GeoLib::Point >::difference_type n)
Definition: Polyline.cpp:550
bool pointsAreIdentical(const std::vector< Point *> &pnt_vec, std::size_t i, std::size_t j, double prox)
Definition: Polyline.cpp:636
SegmentIterator & operator++()
Definition: Polyline.cpp:506
std::size_t getSegmentNumber() const
Definition: Polyline.cpp:501
SegmentIterator & operator+=(std::vector< GeoLib::Point >::difference_type n)
Definition: Polyline.cpp:533
bool operator==(SegmentIterator const &other)
Definition: Polyline.cpp:522
bool isCoplanar() const
Definition: Polyline.cpp:218
bool operator!=(SegmentIterator const &other)
Definition: Polyline.cpp:527
SegmentIterator & operator-=(std::vector< GeoLib::Point >::difference_type n)
Definition: Polyline.cpp:558
SegmentIterator & operator=(SegmentIterator const &rhs)
Definition: Polyline.cpp:488
LineSegment getSegment(std::size_t i) const
Definition: Polyline.cpp:251
void closePolyline()
Definition: Polyline.cpp:392