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