Doc Classic

Class RobustLineIntersector

  1  
  2  
  3 /*
  4  * Copyright (c) 2016 Vivid Solutions.
  5  *
  6  * All rights reserved. This program and the accompanying materials
  7  * are made available under the terms of the Eclipse Public License 2.0
  8  * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
  9  * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
 10  * and the Eclipse Distribution License is available at
 11  *
 12  * http://www.eclipse.org/org/documents/edl-v10.php.
 13  */
 14 package org.locationtech.jts.algorithm;
 15  
 16 /**
 17  *@version 1.7
 18  */
 19 import org.locationtech.jts.geom.Coordinate;
 20 import org.locationtech.jts.geom.Envelope;
 21  
 22 /**
 23  * A robust version of {@link LineIntersector}.
 24  *
 25  * @version 1.7
 26  */
 27 public class RobustLineIntersector
 28     extends LineIntersector
 29 {
 30  
 31   public RobustLineIntersector() {
 32   }
 33  
 34   public void computeIntersection(Coordinate p, Coordinate p1, Coordinate p2) {
 35     isProper = false;
 36     // do between check first, since it is faster than the orientation test
 37     if (Envelope.intersects(p1, p2, p)) {
 38       if ((Orientation.index(p1, p2, p) == 0)
 39           && (Orientation.index(p2, p1, p) == 0)) {
 40         isProper = true;
 41         if (p.equals(p1) || p.equals(p2)) {
 42           isProper = false;
 43         }
 44         result = POINT_INTERSECTION;
 45         return;
 46       }
 47     }
 48     result = NO_INTERSECTION;
 49   }
 50  
 51   protected int computeIntersect(
 52                 Coordinate p1, Coordinate p2,
 53                 Coordinate q1, Coordinate q2  ) {
 54     isProper = false;
 55  
 56     // first try a fast test to see if the envelopes of the lines intersect
 57     if (! Envelope.intersects(p1, p2, q1, q2))
 58       return NO_INTERSECTION;
 59  
 60     // for each endpoint, compute which side of the other segment it lies
 61     // if both endpoints lie on the same side of the other segment,
 62     // the segments do not intersect
 63     int Pq1 = Orientation.index(p1, p2, q1);
 64     int Pq2 = Orientation.index(p1, p2, q2);
 65  
 66     if ((Pq1>0 && Pq2>0) || (Pq1<0 && Pq2<0)) {
 67       return NO_INTERSECTION;
 68     }
 69  
 70     int Qp1 = Orientation.index(q1, q2, p1);
 71     int Qp2 = Orientation.index(q1, q2, p2);
 72  
 73     if ((Qp1>0 && Qp2>0) || (Qp1<0 && Qp2<0)) {
 74         return NO_INTERSECTION;
 75     }
 76  
 77     boolean collinear = Pq1 == 0
 78          && Pq2 == 0
 79          && Qp1 == 0
 80          && Qp2 == 0;
 81     if (collinear) {
 82       return computeCollinearIntersection(p1, p2, q1, q2);
 83     }
 84     
 85     /**
 86      * At this point we know that there is a single intersection point
 87      * (since the lines are not collinear).
 88      */
 89     
 90     /**
 91      *  Check if the intersection is an endpoint. If it is, copy the endpoint as
 92      *  the intersection point. Copying the point rather than computing it
 93      *  ensures the point has the exact value, which is important for
 94      *  robustness. It is sufficient to simply check for an endpoint which is on
 95      *  the other line, since at this point we know that the inputLines must
 96      *  intersect.
 97      */
 98     if (Pq1 == 0 || Pq2 == 0 || Qp1 == 0 || Qp2 == 0) {
 99       isProper = false;
100       
101       /**
102        * Check for two equal endpoints.  
103        * This is done explicitly rather than by the orientation tests
104        * below in order to improve robustness.
105        * 
106        * [An example where the orientation tests fail to be consistent is
107        * the following (where the true intersection is at the shared endpoint
108        * POINT (19.850257749638203 46.29709338043669)
109        * 
110        * LINESTRING ( 19.850257749638203 46.29709338043669, 20.31970698357233 46.76654261437082 ) 
111        * and 
112        * LINESTRING ( -48.51001596420236 -22.063180333403878, 19.850257749638203 46.29709338043669 )
113        * 
114        * which used to produce the INCORRECT result: (20.31970698357233, 46.76654261437082, NaN)
115        * 
116        */
117       if (p1.equals2D(q1) 
118               || p1.equals2D(q2)) {
119           intPt[0] = p1;
120       }
121       else if (p2.equals2D(q1) 
122               || p2.equals2D(q2)) {
123           intPt[0] = p2;
124       }
125       
126       /**
127        * Now check to see if any endpoint lies on the interior of the other segment.
128        */
129       else if (Pq1 == 0) {
130         intPt[0] = new Coordinate(q1);
131       }
132       else if (Pq2 == 0) {
133         intPt[0] = new Coordinate(q2);
134       }
135       else if (Qp1 == 0) {
136         intPt[0] = new Coordinate(p1);
137       }
138       else if (Qp2 == 0) {
139         intPt[0] = new Coordinate(p2);
140       }
141     }
142     else {
143       isProper = true;
144       intPt[0] = intersection(p1, p2, q1, q2);
145     }
146     return POINT_INTERSECTION;
147   }
148  
149   private int computeCollinearIntersection(Coordinate p1, Coordinate p2,
150       Coordinate q1, Coordinate q2) {
151     boolean p1q1p2 = Envelope.intersects(p1, p2, q1);
152     boolean p1q2p2 = Envelope.intersects(p1, p2, q2);
153     boolean q1p1q2 = Envelope.intersects(q1, q2, p1);
154     boolean q1p2q2 = Envelope.intersects(q1, q2, p2);
155  
156     if (p1q1p2 && p1q2p2) {
157       intPt[0] = q1;
158       intPt[1] = q2;
159       return COLLINEAR_INTERSECTION;
160     }
161     if (q1p1q2 && q1p2q2) {
162       intPt[0] = p1;
163       intPt[1] = p2;
164       return COLLINEAR_INTERSECTION;
165     }
166     if (p1q1p2 && q1p1q2) {
167       intPt[0] = q1;
168       intPt[1] = p1;
169       return q1.equals(p1) && !p1q2p2 && !q1p2q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
170     }
171     if (p1q1p2 && q1p2q2) {
172       intPt[0] = q1;
173       intPt[1] = p2;
174       return q1.equals(p2) && !p1q2p2 && !q1p1q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
175     }
176     if (p1q2p2 && q1p1q2) {
177       intPt[0] = q2;
178       intPt[1] = p1;
179       return q2.equals(p1) && !p1q1p2 && !q1p2q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
180     }
181     if (p1q2p2 && q1p2q2) {
182       intPt[0] = q2;
183       intPt[1] = p2;
184       return q2.equals(p2) && !p1q1p2 && !q1p1q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
185     }
186     return NO_INTERSECTION;
187   }
188  
189   /**
190    * This method computes the actual value of the intersection point.
191    * To obtain the maximum precision from the intersection calculation,
192    * the coordinates are normalized by subtracting the minimum
193    * ordinate values (in absolute value).  This has the effect of
194    * removing common significant digits from the calculation to
195    * maintain more bits of precision.
196    */
197   private Coordinate intersection(
198     Coordinate p1, Coordinate p2, Coordinate q1, Coordinate q2)
199   {
200     Coordinate intPt = intersectionSafe(p1, p2, q1, q2);
201     
202     /*
203     // TESTING ONLY
204     Coordinate intPtDD = CGAlgorithmsDD.intersection(p1, p2, q1, q2);
205     double dist = intPt.distance(intPtDD);
206     System.out.println(intPt + " - " + intPtDD + " dist = " + dist);
207     //intPt = intPtDD;
208     */
209     
210     /**
211      * Due to rounding it can happen that the computed intersection is
212      * outside the envelopes of the input segments.  Clearly this
213      * is inconsistent. 
214      * This code checks this condition and forces a more reasonable answer
215      * 
216      * MD - May 4 2005 - This is still a problem.  Here is a failure case:
217      *
218      * LINESTRING (2089426.5233462777 1180182.3877339689, 2085646.6891757075 1195618.7333999649)
219      * LINESTRING (1889281.8148903656 1997547.0560044837, 2259977.3672235999 483675.17050843034)
220      * int point = (2097408.2633752143,1144595.8008114607)
221      * 
222      * MD - Dec 14 2006 - This does not seem to be a failure case any longer
223      */
224     if (! isInSegmentEnvelopes(intPt)) {
225 //      System.out.println("Intersection outside segment envelopes: " + intPt);
226       
227       // compute a safer result
228       // copy the coordinate, since it may be rounded later
229       intPt = new Coordinate(nearestEndpoint(p1, p2, q1, q2));
230 //    intPt = CentralEndpointIntersector.getIntersection(p1, p2, q1, q2);
231       
232 //      System.out.println("Segments: " + this);
233 //      System.out.println("Snapped to " + intPt);
234 //      checkDD(p1, p2, q1, q2, intPt);
235     }
236     if (precisionModel != null) {
237       precisionModel.makePrecise(intPt);
238     }
239     return intPt;
240   }
241  
242   private void checkDD(Coordinate p1, Coordinate p2, Coordinate q1,
243       Coordinate q2, Coordinate intPt)
244   {
245     Coordinate intPtDD = CGAlgorithmsDD.intersection(p1, p2, q1, q2);
246     boolean isIn = isInSegmentEnvelopes(intPtDD);
247     System.out.println(   "DD in env = " + isIn + "  --------------------- " + intPtDD);
248     if (intPt.distance(intPtDD) > 0.0001) {
249       System.out.println("Distance = " + intPt.distance(intPtDD));
250     }
251   }
252   
253   /**
254    * Computes a segment intersection using homogeneous coordinates.
255    * Round-off error can cause the raw computation to fail, 
256    * (usually due to the segments being approximately parallel).
257    * If this happens, a reasonable approximation is computed instead.
258    * 
259    * @param p1 a segment endpoint
260    * @param p2 a segment endpoint
261    * @param q1 a segment endpoint
262    * @param q2 a segment endpoint
263    * @return the computed intersection point
264    */
265   private Coordinate intersectionSafe(Coordinate p1, Coordinate p2, Coordinate q1, Coordinate q2)
266   {
267     Coordinate intPt = Intersection.intersection(p1, p2, q1, q2);
268     if (intPt == null)
269       intPt = nearestEndpoint(p1, p2, q1, q2);
270  //     System.out.println("Snapped to " + intPt);
271     return intPt;
272   }
273  
274   /**
275    * Tests whether a point lies in the envelopes of both input segments.
276    * A correctly computed intersection point should return <code>true</code>
277    * for this test.
278    * Since this test is for debugging purposes only, no attempt is
279    * made to optimize the envelope test.
280    *
281    * @return <code>true</code> if the input point lies within both input segment envelopes
282    */
283   private boolean isInSegmentEnvelopes(Coordinate intPt)
284   {
285     Envelope env0 = new Envelope(inputLines[0][0], inputLines[0][1]);
286     Envelope env1 = new Envelope(inputLines[1][0], inputLines[1][1]);
287     return env0.contains(intPt) && env1.contains(intPt);
288   }
289  
290   /**
291    * Finds the endpoint of the segments P and Q which 
292    * is closest to the other segment.
293    * This is a reasonable surrogate for the true 
294    * intersection points in ill-conditioned cases
295    * (e.g. where two segments are nearly coincident,
296    * or where the endpoint of one segment lies almost on the other segment).
297    * <p>
298    * This replaces the older CentralEndpoint heuristic,
299    * which chose the wrong endpoint in some cases
300    * where the segments had very distinct slopes 
301    * and one endpoint lay almost on the other segment.
302    * 
303    * @param p1 an endpoint of segment P
304    * @param p2 an endpoint of segment P
305    * @param q1 an endpoint of segment Q
306    * @param q2 an endpoint of segment Q
307    * @return the nearest endpoint to the other segment
308    */
309   private static Coordinate nearestEndpoint(Coordinate p1, Coordinate p2,
310       Coordinate q1, Coordinate q2)
311   {
312     Coordinate nearestPt = p1;
313     double minDist = Distance.pointToSegment(p1, q1, q2);
314     
315     double dist = Distance.pointToSegment(p2, q1, q2);
316     if (dist < minDist) {
317       minDist = dist;
318       nearestPt = p2;
319     }
320     dist = Distance.pointToSegment(q1, p1, p2);
321     if (dist < minDist) {
322       minDist = dist;
323       nearestPt = q1;
324     }
325     dist = Distance.pointToSegment(q2, p1, p2);
326     if (dist < minDist) {
327       minDist = dist;
328       nearestPt = q2;
329     }
330     return nearestPt;
331   }
332  
333  
334 }
335