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