GEOS 3.12.0
LineIntersector.h
1/**********************************************************************
2 *
3 * GEOS - Geometry Engine Open Source
4 * http://geos.osgeo.org
5 *
6 * Copyright (C) 2005-2006 Refractions Research Inc.
7 * Copyright (C) 2001-2002 Vivid Solutions Inc.
8 *
9 * This is free software; you can redistribute and/or modify it under
10 * the terms of the GNU Lesser General Public Licence as published
11 * by the Free Software Foundation.
12 * See the COPYING file for more information.
13 *
14 **********************************************************************
15 *
16 * Last port: algorithm/RobustLineIntersector.java r785 (JTS-1.13+)
17 *
18 **********************************************************************/
19
20#pragma once
21
22#include <geos/export.h>
23#include <geos/algorithm/Intersection.h>
24#include <geos/algorithm/Interpolate.h>
25#include <geos/algorithm/Orientation.h>
26#include <geos/geom/Coordinate.h>
27#include <geos/geom/Envelope.h>
28#include <geos/geom/PrecisionModel.h>
29
30#include <string>
31
32// Forward declarations
33namespace geos {
34namespace geom {
35class PrecisionModel;
36}
37}
38
39namespace geos {
40namespace algorithm { // geos::algorithm
41
53class GEOS_DLL LineIntersector {
54public:
55
74 static double computeEdgeDistance(const geom::CoordinateXY& p, const geom::CoordinateXY& p0, const geom::CoordinateXY& p1);
75
76 static double nonRobustComputeEdgeDistance(const geom::Coordinate& p, const geom::Coordinate& p1,
77 const geom::Coordinate& p2);
78
79 explicit LineIntersector(const geom::PrecisionModel* initialPrecisionModel = nullptr)
80 :
81 precisionModel(initialPrecisionModel),
82 result(0),
83 inputLines(),
84 isProperVar(false)
85 {}
86
87 ~LineIntersector() = default;
88
97 {
98 if(isInteriorIntersection(0)) {
99 return true;
100 }
101 if(isInteriorIntersection(1)) {
102 return true;
103 }
104 return false;
105 };
106
114 bool isInteriorIntersection(std::size_t inputLineIndex)
115 {
116 for(std::size_t i = 0; i < result; ++i) {
117 if(!(intPt[i].equals2D(*inputLines[inputLineIndex][0])
118 || intPt[i].equals2D(*inputLines[inputLineIndex][1]))) {
119 return true;
120 }
121 }
122 return false;
123 };
124
131 void
133 {
134 precisionModel = newPM;
135 }
136
143 void computeIntersection(const geom::CoordinateXY& p, const geom::CoordinateXY& p1, const geom::CoordinateXY& p2);
144
146 static bool hasIntersection(const geom::CoordinateXY& p, const geom::CoordinateXY& p1, const geom::CoordinateXY& p2);
147
148 enum intersection_type : uint8_t {
150 NO_INTERSECTION = 0,
151
153 POINT_INTERSECTION = 1,
154
156 COLLINEAR_INTERSECTION = 2
157 };
158
160 template<typename C1, typename C2>
161 void computeIntersection(const C1& p1, const C1& p2,
162 const C2& p3, const C2& p4)
163 {
164 inputLines[0][0] = &p1;
165 inputLines[0][1] = &p2;
166 inputLines[1][0] = &p3;
167 inputLines[1][1] = &p4;
168 result = computeIntersect(p1, p2, p3, p4);
169 }
170
172 void computeIntersection(const geom::CoordinateSequence& p, std::size_t p0,
173 const geom::CoordinateSequence& q, std::size_t q0);
174
175 std::string toString() const;
176
182 bool
184 {
185 return result != NO_INTERSECTION;
186 }
187
188
196 const geom::CoordinateXY*
197 getEndpoint(std::size_t segmentIndex, std::size_t ptIndex) const
198 {
199 return inputLines[segmentIndex][ptIndex];
200 }
201
206 size_t
208 {
209 return result;
210 }
211
212
219 const geom::CoordinateXYZM&
220 getIntersection(std::size_t intIndex) const
221 {
222 return intPt[intIndex];
223 }
224
229 static bool isSameSignAndNonZero(double a, double b);
230
241 bool isIntersection(const geom::Coordinate& pt) const
242 {
243 for(std::size_t i = 0; i < result; ++i) {
244 if(intPt[i].equals2D(pt)) {
245 return true;
246 }
247 }
248 return false;
249 };
250
265 bool
266 isProper() const
267 {
268 return hasIntersection() && isProperVar;
269 }
270
281 const geom::Coordinate& getIntersectionAlongSegment(std::size_t segmentIndex, std::size_t intIndex);
282
292 std::size_t getIndexAlongSegment(std::size_t segmentIndex, std::size_t intIndex);
293
303 double getEdgeDistance(std::size_t geomIndex, std::size_t intIndex) const;
304
305private:
306
311 const geom::PrecisionModel* precisionModel;
312
313 std::size_t result;
314
315 const geom::CoordinateXY* inputLines[2][2];
316
321 geom::CoordinateXYZM intPt[2];
322
327 std::size_t intLineIndex[2][2];
328
329 bool isProperVar;
330 //Coordinate &pa;
331 //Coordinate &pb;
332
333 bool
334 isCollinear() const
335 {
336 return result == COLLINEAR_INTERSECTION;
337 }
338
339 template<typename C1, typename C2>
340 uint8_t computeIntersect(const C1& p1, const C1& p2, const C2& q1, const C2& q2)
341 {
342 isProperVar = false;
343
344 // first try a fast test to see if the envelopes of the lines intersect
345 if(!geom::Envelope::intersects(p1, p2, q1, q2)) {
346 return NO_INTERSECTION;
347 }
348
349 // for each endpoint, compute which side of the other segment it lies
350 // if both endpoints lie on the same side of the other segment,
351 // the segments do not intersect
352 int Pq1 = Orientation::index(p1, p2, q1);
353 int Pq2 = Orientation::index(p1, p2, q2);
354
355 if((Pq1 > 0 && Pq2 > 0) || (Pq1 < 0 && Pq2 < 0)) {
356 return NO_INTERSECTION;
357 }
358
359 int Qp1 = Orientation::index(q1, q2, p1);
360 int Qp2 = Orientation::index(q1, q2, p2);
361
362 if((Qp1 > 0 && Qp2 > 0) || (Qp1 < 0 && Qp2 < 0)) {
363 return NO_INTERSECTION;
364 }
365
369 bool collinear = Pq1 == 0 && Pq2 == 0 && Qp1 == 0 && Qp2 == 0;
370 if(collinear) {
371 return computeCollinearIntersection(p1, p2, q1, q2);
372 }
373
374 /*
375 * At this point we know that there is a single intersection point
376 * (since the lines are not collinear).
377 */
378
379 /*
380 * Check if the intersection is an endpoint.
381 * If it is, copy the endpoint as
382 * the intersection point. Copying the point rather than
383 * computing it ensures the point has the exact value,
384 * which is important for robustness. It is sufficient to
385 * simply check for an endpoint which is on the other line,
386 * since at this point we know that the inputLines must
387 * intersect.
388 */
389 geom::CoordinateXYZM p;
390 double z = DoubleNotANumber;
391 double m = DoubleNotANumber;
392
393 if(Pq1 == 0 || Pq2 == 0 || Qp1 == 0 || Qp2 == 0) {
394
395 isProperVar = false;
396
397 /* Check for two equal endpoints.
398 * This is done explicitly rather than by the orientation tests
399 * below in order to improve robustness.
400 *
401 * (A example where the orientation tests fail
402 * to be consistent is:
403 *
404 * LINESTRING ( 19.850257749638203 46.29709338043669,
405 * 20.31970698357233 46.76654261437082 )
406 * and
407 * LINESTRING ( -48.51001596420236 -22.063180333403878,
408 * 19.850257749638203 46.29709338043669 )
409 *
410 * which used to produce the INCORRECT result:
411 * (20.31970698357233, 46.76654261437082, NaN)
412 */
413
414 if (p1.equals2D(q1)) {
415 p = p1;
416 z = Interpolate::zGet(p1, q1);
417 m = Interpolate::mGet(p1, q1);
418 }
419 else if (p1.equals2D(q2)) {
420 p = p1;
421 z = Interpolate::zGet(p1, q2);
422 m = Interpolate::mGet(p1, q2);
423 }
424 else if (p2.equals2D(q1)) {
425 p = p2;
426 z = Interpolate::zGet(p2, q1);
427 m = Interpolate::mGet(p2, q1);
428 }
429 else if (p2.equals2D(q2)) {
430 p = p2;
431 z = Interpolate::zGet(p2, q2);
432 m = Interpolate::mGet(p2, q2);
433 }
434 /*
435 * Now check to see if any endpoint lies on the interior of the other segment.
436 */
437 else if(Pq1 == 0) {
438 p = q1;
439 z = Interpolate::zGetOrInterpolate(q1, p1, p2);
440 m = Interpolate::mGetOrInterpolate(q1, p1, p2);
441 }
442 else if(Pq2 == 0) {
443 p = q2;
444 z = Interpolate::zGetOrInterpolate(q2, p1, p2);
445 m = Interpolate::mGetOrInterpolate(q2, p1, p2);
446 }
447 else if(Qp1 == 0) {
448 p = p1;
449 z = Interpolate::zGetOrInterpolate(p1, q1, q2);
450 m = Interpolate::mGetOrInterpolate(p1, q1, q2);
451 }
452 else if(Qp2 == 0) {
453 p = p2;
454 z = Interpolate::zGetOrInterpolate(p2, q1, q2);
455 m = Interpolate::mGetOrInterpolate(p2, q1, q2);
456 }
457 } else {
458 isProperVar = true;
459 p = intersection(p1, p2, q1, q2);
460 z = Interpolate::zInterpolate(p, p1, p2, q1, q2);
461 m = Interpolate::mInterpolate(p, p1, p2, q1, q2);
462 }
463 intPt[0] = geom::CoordinateXYZM(p.x, p.y, z, m);
464 #if GEOS_DEBUG
465 std::cerr << " POINT_INTERSECTION; intPt[0]:" << intPt[0].toString() << std::endl;
466 #endif // GEOS_DEBUG
467 return POINT_INTERSECTION;
468 }
469
470 bool
471 isEndPoint() const
472 {
473 return hasIntersection() && !isProperVar;
474 }
475
476 void computeIntLineIndex();
477
478 void computeIntLineIndex(std::size_t segmentIndex);
479
480 template<typename C1, typename C2>
481 uint8_t computeCollinearIntersection(const C1& p1, const C1& p2, const C2& q1, const C2& q2)
482 {
483 bool q1inP = geom::Envelope::intersects(p1, p2, q1);
484 bool q2inP = geom::Envelope::intersects(p1, p2, q2);
485 bool p1inQ = geom::Envelope::intersects(q1, q2, p1);
486 bool p2inQ = geom::Envelope::intersects(q1, q2, p2);
487
488 if(q1inP && q2inP) {
489 intPt[0] = zmGetOrInterpolateCopy(q1, p1, p2);
490 intPt[1] = zmGetOrInterpolateCopy(q2, p1, p2);
491 return COLLINEAR_INTERSECTION;
492 }
493 if(p1inQ && p2inQ) {
494 intPt[0] = zmGetOrInterpolateCopy(p1, q1, q2);
495 intPt[1] = zmGetOrInterpolateCopy(p2, q1, q2);
496 return COLLINEAR_INTERSECTION;
497 }
498 if(q1inP && p1inQ) {
499 // if pts are equal Z is chosen arbitrarily
500 intPt[0] = zmGetOrInterpolateCopy(q1, p1, p2);
501 intPt[1] = zmGetOrInterpolateCopy(p1, q1, q2);
502
503 return (q1 == p1) && !q2inP && !p2inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
504 }
505 if(q1inP && p2inQ) {
506 // if pts are equal Z is chosen arbitrarily
507 intPt[0] = zmGetOrInterpolateCopy(q1, p1, p2);
508 intPt[1] = zmGetOrInterpolateCopy(p2, q1, q2);
509
510 return (q1 == p2) && !q2inP && !p1inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
511 }
512 if(q2inP && p1inQ) {
513 // if pts are equal Z is chosen arbitrarily
514 intPt[0] = zmGetOrInterpolateCopy(q2, p1, p2);
515 intPt[1] = zmGetOrInterpolateCopy(p1, q1, q2);
516
517 return (q2 == p1) && !q1inP && !p2inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
518 }
519 if(q2inP && p2inQ) {
520 // if pts are equal Z is chosen arbitrarily
521 intPt[0] = zmGetOrInterpolateCopy(q2, p1, p2);
522 intPt[1] = zmGetOrInterpolateCopy(p2, q1, q2);
523 return (q2 == p2) && !q1inP && !p1inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
524 }
525 return NO_INTERSECTION;
526 }
527
537 template<typename C1, typename C2>
538 geom::CoordinateXYZM intersection (const C1& p1, const C1& p2, const C2& q1, const C2& q2) const {
539 auto intPtOut = intersectionSafe(p1, p2, q1, q2);
540
541 /*
542 * Due to rounding it can happen that the computed intersection is
543 * outside the envelopes of the input segments. Clearly this
544 * is inconsistent.
545 * This code checks this condition and forces a more reasonable answer
546 *
547 * MD - May 4 2005 - This is still a problem. Here is a failure case:
548 *
549 * LINESTRING (2089426.5233462777 1180182.3877339689,
550 * 2085646.6891757075 1195618.7333999649)
551 * LINESTRING (1889281.8148903656 1997547.0560044837,
552 * 2259977.3672235999 483675.17050843034)
553 * int point = (2097408.2633752143,1144595.8008114607)
554 */
555
556 if(! isInSegmentEnvelopes(intPtOut)) {
557 //intPt = CentralEndpointIntersector::getIntersection(p1, p2, q1, q2);
558 intPtOut = nearestEndpoint(p1, p2, q1, q2);
559 }
560
561 if(precisionModel != nullptr) {
562 precisionModel->makePrecise(intPtOut);
563 }
564
565 return intPtOut;
566 }
567
578 bool isInSegmentEnvelopes(const geom::CoordinateXY& pt) const
579 {
580 geom::Envelope env0(*inputLines[0][0], *inputLines[0][1]);
581 geom::Envelope env1(*inputLines[1][0], *inputLines[1][1]);
582 return env0.contains(pt) && env1.contains(pt);
583 };
584
597 template<typename C1, typename C2>
598 geom::CoordinateXYZM intersectionSafe(const C1& p1, const C1& p2,
599 const C2& q1, const C2& q2) const
600 {
601 geom::CoordinateXYZM ptInt(Intersection::intersection(p1, p2, q1, q2));
602 if (ptInt.isNull()) {
603 // FIXME need to cast to correct type in mixed-dimensionality case
604 ptInt = static_cast<const C1&>(nearestEndpoint(p1, p2, q1, q2));
605 }
606 return ptInt;
607 }
608
628 static const geom::CoordinateXY& nearestEndpoint(const geom::CoordinateXY& p1,
629 const geom::CoordinateXY& p2,
630 const geom::CoordinateXY& q1,
631 const geom::CoordinateXY& q2);
632
633
634 template<typename C1, typename C2>
635 static geom::CoordinateXYZM zmGetOrInterpolateCopy(
636 const C1& p,
637 const C2& p1,
638 const C2& p2)
639 {
640 geom::CoordinateXYZM pCopy(p);
641 pCopy.z = Interpolate::zGetOrInterpolate(p, p1, p2);
642 pCopy.m = Interpolate::mGetOrInterpolate(p, p1, p2);
643 return pCopy;
644 }
645
646};
647
648
649} // namespace geos::algorithm
650} // namespace geos
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
A LineIntersector is an algorithm that can both test whether two line segments intersect and compute ...
Definition LineIntersector.h:53
intersection_type
Definition LineIntersector.h:148
void computeIntersection(const geom::CoordinateSequence &p, std::size_t p0, const geom::CoordinateSequence &q, std::size_t q0)
Compute the intersection between two segments, given a sequence and starting index of each.
void computeIntersection(const geom::CoordinateXY &p, const geom::CoordinateXY &p1, const geom::CoordinateXY &p2)
bool isInteriorIntersection(std::size_t inputLineIndex)
Tests whether either intersection point is an interior point of the specified input segment.
Definition LineIntersector.h:114
const geom::Coordinate & getIntersectionAlongSegment(std::size_t segmentIndex, std::size_t intIndex)
Computes the intIndex'th intersection point in the direction of a specified input line segment.
void setPrecisionModel(const geom::PrecisionModel *newPM)
Definition LineIntersector.h:132
bool isInteriorIntersection()
Tests whether either intersection point is an interior point of one of the input segments.
Definition LineIntersector.h:96
const geom::CoordinateXYZM & getIntersection(std::size_t intIndex) const
Definition LineIntersector.h:220
static double computeEdgeDistance(const geom::CoordinateXY &p, const geom::CoordinateXY &p0, const geom::CoordinateXY &p1)
bool hasIntersection() const
Definition LineIntersector.h:183
static bool isSameSignAndNonZero(double a, double b)
static bool hasIntersection(const geom::CoordinateXY &p, const geom::CoordinateXY &p1, const geom::CoordinateXY &p2)
Same as above but doesn't compute intersection point. Faster.
void computeIntersection(const C1 &p1, const C1 &p2, const C2 &p3, const C2 &p4)
Computes the intersection of the lines p1-p2 and p3-p4.
Definition LineIntersector.h:161
const geom::CoordinateXY * getEndpoint(std::size_t segmentIndex, std::size_t ptIndex) const
Definition LineIntersector.h:197
bool isProper() const
Tests whether an intersection is proper.
Definition LineIntersector.h:266
size_t getIntersectionNum() const
Definition LineIntersector.h:207
double getEdgeDistance(std::size_t geomIndex, std::size_t intIndex) const
Computes the "edge distance" of an intersection point along the specified input line segment.
std::size_t getIndexAlongSegment(std::size_t segmentIndex, std::size_t intIndex)
Computes the index of the intIndex'th intersection point in the direction of a specified input line s...
bool isIntersection(const geom::Coordinate &pt) const
Test whether a point is a intersection point of two line segments.
Definition LineIntersector.h:241
The internal representation of a list of coordinates inside a Geometry.
Definition CoordinateSequence.h:56
Coordinate is the lightweight class used to store coordinates.
Definition Coordinate.h:216
Specifies the precision model of the Coordinate in a Geometry.
Definition PrecisionModel.h:90
double makePrecise(double val) const
Rounds a numeric value to the PrecisionModel grid.
Basic namespace for all GEOS functionalities.
Definition geos.h:39