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