Doc Classic

Class CGAlgorithms

  1 /*
  2  * Copyright (c) 2016 Vivid Solutions.
  3  *
  4  * All rights reserved. This program and the accompanying materials
  5  * are made available under the terms of the Eclipse Public License 2.0
  6  * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
  7  * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
  8  * and the Eclipse Distribution License is available at
  9  *
 10  * http://www.eclipse.org/org/documents/edl-v10.php.
 11  */
 12 package org.locationtech.jts.algorithm;
 13  
 14 import org.locationtech.jts.geom.Coordinate;
 15 import org.locationtech.jts.geom.CoordinateSequence;
 16 import org.locationtech.jts.geom.Envelope;
 17 import org.locationtech.jts.geom.Location;
 18 import org.locationtech.jts.math.MathUtil;
 19  
 20 /**
 21  * Specifies and implements various fundamental Computational Geometric
 22  * algorithms. The algorithms supplied in this class are robust for
 23  * double-precision floating point.
 24  * 
 25  * @version 1.7
 26  * @deprecated See {@link Length}, {@link Area}, {@link Distance},
 27  *             {@link Orientation}, {@link PointLocation}
 28  */
 29 public class CGAlgorithms
 30 {
 31  
 32   /**
 33    * A value that indicates an orientation of clockwise, or a right turn.
 34    * 
 35    * @deprecated Use {@link Orientation#CLOCKWISE} instead.
 36    */
 37   public static final int CLOCKWISE = -1;
 38  
 39   /**
 40    * A value that indicates an orientation of clockwise, or a right turn.
 41    * 
 42    * @deprecated Use {@link Orientation#RIGHT} instead.
 43    */
 44   public static final int RIGHT = CLOCKWISE;
 45  
 46   /**
 47    * A value that indicates an orientation of counterclockwise, or a left turn.
 48    * 
 49    * @deprecated Use {@link Orientation#COUNTERCLOCKWISE} instead.
 50    */
 51   public static final int COUNTERCLOCKWISE = 1;
 52  
 53   /**
 54    * A value that indicates an orientation of counterclockwise, or a left turn.
 55    * 
 56    * @deprecated Use {@link Orientation#LEFT} instead.
 57    */
 58   public static final int LEFT = COUNTERCLOCKWISE;
 59  
 60   /**
 61    * A value that indicates an orientation of collinear, or no turn (straight).
 62    * 
 63    * @deprecated Use {@link Orientation#COLLINEAR} instead.
 64    */
 65   public static final int COLLINEAR = 0;
 66  
 67   /**
 68    * A value that indicates an orientation of collinear, or no turn (straight).
 69    * 
 70    * @deprecated Use {@link Orientation#STRAIGHT} instead.
 71    */
 72   public static final int STRAIGHT = COLLINEAR;
 73  
 74   /**
 75    * Returns the index of the direction of the point <code>q</code> relative to
 76    * a vector specified by <code>p1-p2</code>.
 77    * 
 78    * @param p1 the origin point of the vector
 79    * @param p2 the final point of the vector
 80    * @param q the point to compute the direction to
 81    * 
 82    * @return 1 if q is counter-clockwise (left) from p1-p2
 83    * @return -1 if q is clockwise (right) from p1-p2
 84    * @return 0 if q is collinear with p1-p2
 85    * @deprecated Use {@link Orientation#index(Coordinate, Coordinate, Coordinate)}
 86    *             instead.
 87    */
 88   public static int orientationIndex(Coordinate p1, Coordinate p2, Coordinate q)
 89   {
 90     /**
 91      * MD - 9 Aug 2010 It seems that the basic algorithm is slightly orientation
 92      * dependent, when computing the orientation of a point very close to a
 93      * line. This is possibly due to the arithmetic in the translation to the
 94      * origin.
 95      * 
 96      * For instance, the following situation produces identical results in spite
 97      * of the inverse orientation of the line segment:
 98      * 
 99      * Coordinate p0 = new Coordinate(219.3649559090992, 140.84159161824724);
100      * Coordinate p1 = new Coordinate(168.9018919682399, -5.713787599646864);
101      * 
102      * Coordinate p = new Coordinate(186.80814046338352, 46.28973405831556); int
103      * orient = orientationIndex(p0, p1, p); int orientInv =
104      * orientationIndex(p1, p0, p);
105      * 
106      * A way to force consistent results is to normalize the orientation of the
107      * vector using the following code. However, this may make the results of
108      * orientationIndex inconsistent through the triangle of points, so it's not
109      * clear this is an appropriate patch.
110      * 
111      */
112     return CGAlgorithmsDD.orientationIndex(p1, p2, q);
113     // testing only
114     //return ShewchuksDeterminant.orientationIndex(p1, p2, q);
115     // previous implementation - not quite fully robust
116     //return RobustDeterminant.orientationIndex(p1, p2, q);
117     
118   }
119  
120   public CGAlgorithms()
121   {
122   }
123  
124   /**
125    * Tests whether a point lies inside or on a ring. The ring may be oriented in
126    * either direction. A point lying exactly on the ring boundary is considered
127    * to be inside the ring.
128    * <p>
129    * This method does <i>not</i> first check the point against the envelope of
130    * the ring.
131    * 
132    * @param p
133    *          point to check for ring inclusion
134    * @param ring
135    *          an array of coordinates representing the ring (which must have
136    *          first point identical to last point)
137    * @return true if p is inside ring
138    * 
139    * @see locatePointInRing
140    * @deprecated Use {@link PointLocation#isInRing(Coordinate, Coordinate[])}
141    *             instead.
142    */
143   public static boolean isPointInRing(Coordinate p, Coordinate[] ring)
144   {
145     return locatePointInRing(p, ring) != Location.EXTERIOR;
146   }
147  
148   /**
149    * Determines whether a point lies in the interior, on the boundary, or in the
150    * exterior of a ring. The ring may be oriented in either direction.
151    * <p>
152    * This method does <i>not</i> first check the point against the envelope of
153    * the ring.
154    * 
155    * @param p
156    *          point to check for ring inclusion
157    * @param ring
158    *          an array of coordinates representing the ring (which must have
159    *          first point identical to last point)
160    * @return the {@link Location} of p relative to the ring
161    * @deprecated Use
162    *             {@link PointLocation#locateInRing(Coordinate, Coordinate[])}
163    *             instead.
164    */
165   public static int locatePointInRing(Coordinate p, Coordinate[] ring)
166   {
167     return RayCrossingCounter.locatePointInRing(p, ring);
168   }
169  
170   /**
171    * Tests whether a point lies on the line segments defined by a list of
172    * coordinates.
173    * 
174    * @return true if the point is a vertex of the line or lies in the interior
175    *         of a line segment in the linestring
176    * @deprecated Use {@link PointLocation#isOnLine(Coordinate, Coordinate[])}
177    *             instead.
178    */
179   public static boolean isOnLine(Coordinate p, Coordinate[] pt)
180   {
181     LineIntersector lineIntersector = new RobustLineIntersector();
182     for (int i = 1; i < pt.length; i++) {
183       Coordinate p0 = pt[i - 1];
184       Coordinate p1 = pt[i];
185       lineIntersector.computeIntersection(p, p0, p1);
186       if (lineIntersector.hasIntersection()) {
187         return true;
188       }
189     }
190     return false;
191   }
192  
193   /**
194    * Computes whether a ring defined by an array of {@link Coordinate}s is
195    * oriented counter-clockwise.
196    * <ul>
197    * <li>The list of points is assumed to have the first and last points equal.
198    * <li>This will handle coordinate lists which contain repeated points.
199    * </ul>
200    * This algorithm is <b>only</b> guaranteed to work with valid rings. If the
201    * ring is invalid (e.g. self-crosses or touches), the computed result may not
202    * be correct.
203    * 
204    * @param ring
205    *          an array of Coordinates forming a ring
206    * @return true if the ring is oriented counter-clockwise.
207    * @throws IllegalArgumentException
208    *           if there are too few points to determine orientation (< 4)
209    * @deprecated Use {@link Orientation#isCCW(Coordinate[])} instead.
210    */
211   public static boolean isCCW(Coordinate[] ring)
212   {
213     // # of points without closing endpoint
214     int nPts = ring.length - 1;
215     // sanity check
216     if (nPts < 3)
217       throw new IllegalArgumentException(
218           "Ring has fewer than 4 points, so orientation cannot be determined");
219  
220     // find highest point
221     Coordinate hiPt = ring[0];
222     int hiIndex = 0;
223     for (int i = 1; i <= nPts; i++) {
224       Coordinate p = ring[i];
225       if (p.y > hiPt.y) {
226         hiPt = p;
227         hiIndex = i;
228       }
229     }
230  
231     // find distinct point before highest point
232     int iPrev = hiIndex;
233     do {
234       iPrev = iPrev - 1;
235       if (iPrev < 0)
236         iPrev = nPts;
237     } while (ring[iPrev].equals2D(hiPt) && iPrev != hiIndex);
238  
239     // find distinct point after highest point
240     int iNext = hiIndex;
241     do {
242       iNext = (iNext + 1) % nPts;
243     } while (ring[iNext].equals2D(hiPt) && iNext != hiIndex);
244  
245     Coordinate prev = ring[iPrev];
246     Coordinate next = ring[iNext];
247  
248     /**
249      * This check catches cases where the ring contains an A-B-A configuration
250      * of points. This can happen if the ring does not contain 3 distinct points
251      * (including the case where the input array has fewer than 4 elements), or
252      * it contains coincident line segments.
253      */
254     if (prev.equals2D(hiPt) || next.equals2D(hiPt) || prev.equals2D(next))
255       return false;
256  
257     int disc = computeOrientation(prev, hiPt, next);
258  
259     /**
260      * If disc is exactly 0, lines are collinear. There are two possible cases:
261      * (1) the lines lie along the x axis in opposite directions (2) the lines
262      * lie on top of one another
263      * 
264      * (1) is handled by checking if next is left of prev ==> CCW (2) will never
265      * happen if the ring is valid, so don't check for it (Might want to assert
266      * this)
267      */
268     boolean isCCW;
269     if (disc == 0) {
270       // poly is CCW if prev x is right of next x
271       isCCW = (prev.x > next.x);
272     }
273     else {
274       // if area is positive, points are ordered CCW
275       isCCW = (disc > 0);
276     }
277     return isCCW;
278   }
279  
280   /**
281    * Computes the orientation of a point q to the directed line segment p1-p2.
282    * The orientation of a point relative to a directed line segment indicates
283    * which way you turn to get to q after travelling from p1 to p2.
284    * 
285    * @param p1 the first vertex of the line segment
286    * @param p2 the second vertex of the line segment
287    * @param q the point to compute the relative orientation of
288    * @return 1 if q is counter-clockwise from p1-p2,
289    * or -1 if q is clockwise from p1-p2,
290    * or 0 if q is collinear with p1-p2
291    * @deprecated Use {@link Orientation#index(Coordinate, Coordinate, Coordinate)}
292    *             instead.
293    */
294   public static int computeOrientation(Coordinate p1, Coordinate p2,
295       Coordinate q)
296   {
297     return orientationIndex(p1, p2, q);
298   }
299  
300   /**
301    * Computes the distance from a point p to a line segment AB
302    * 
303    * Note: NON-ROBUST!
304    * 
305    * @param p
306    *          the point to compute the distance for
307    * @param A
308    *          one point of the line
309    * @param B
310    *          another point of the line (must be different to A)
311    * @return the distance from p to line segment AB
312    * @deprecated Use
313    *             {@link Distance#pointToSegment(Coordinate, Coordinate, Coordinate)}
314    *             instead.
315    */
316   public static double distancePointLine(Coordinate p, Coordinate A,
317       Coordinate B)
318   {
319     // if start = end, then just compute distance to one of the endpoints
320     if (A.x == B.x && A.y == B.y)
321       return p.distance(A);
322  
323     // otherwise use comp.graphics.algorithms Frequently Asked Questions method
324     /*
325      * (1) r = AC dot AB 
326      *         --------- 
327      *         ||AB||^2 
328      *         
329      * r has the following meaning: 
330      *   r=0 P = A 
331      *   r=1 P = B 
332      *   r<0 P is on the backward extension of AB 
333      *   r>1 P is on the forward extension of AB 
334      *   0<r<1 P is interior to AB
335      */
336  
337     double len2 = (B.x - A.x) * (B.x - A.x) + (B.y - A.y) * (B.y - A.y);
338     double r = ((p.x - A.x) * (B.x - A.x) + (p.y - A.y) * (B.y - A.y))
339         / len2;
340  
341     if (r <= 0.0)
342       return p.distance(A);
343     if (r >= 1.0)
344       return p.distance(B);
345  
346     /*
347      * (2) s = (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) 
348      *         ----------------------------- 
349      *                    L^2
350      * 
351      * Then the distance from C to P = |s|*L.
352      * 
353      * This is the same calculation as {@link #distancePointLinePerpendicular}.
354      * Unrolled here for performance.
355      */
356     double s = ((A.y - p.y) * (B.x - A.x) - (A.x - p.x) * (B.y - A.y))
357         / len2;
358     return Math.abs(s) * Math.sqrt(len2);
359   }
360  
361   /**
362    * Computes the perpendicular distance from a point p to the (infinite) line
363    * containing the points AB
364    * 
365    * @param p
366    *          the point to compute the distance for
367    * @param A
368    *          one point of the line
369    * @param B
370    *          another point of the line (must be different to A)
371    * @return the distance from p to line AB
372    * @deprecated Use
373    *             {@link Distance#pointToLinePerpendicular(Coordinate, Coordinate, Coordinate)}
374    *             instead.
375    */
376   public static double distancePointLinePerpendicular(Coordinate p,
377       Coordinate A, Coordinate B)
378   {
379     // use comp.graphics.algorithms Frequently Asked Questions method
380     /*
381      * (2) s = (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) 
382      *         ----------------------------- 
383      *                    L^2
384      * 
385      * Then the distance from C to P = |s|*L.
386      */
387     double len2 = (B.x - A.x) * (B.x - A.x) + (B.y - A.y) * (B.y - A.y);
388     double s = ((A.y - p.y) * (B.x - A.x) - (A.x - p.x) * (B.y - A.y))
389         / len2;
390  
391     return Math.abs(s) * Math.sqrt(len2);
392   }
393  
394   /**
395    * Computes the distance from a point to a sequence of line segments.
396    * 
397    * @param p
398    *          a point
399    * @param line
400    *          a sequence of contiguous line segments defined by their vertices
401    * @return the minimum distance between the point and the line segments
402    * @deprecated Use
403    *             {@link Distance#pointToSegmentString(Coordinate, Coordinate[])}
404    *             instead.
405    */
406   public static double distancePointLine(Coordinate p, Coordinate[] line)
407   {
408     if (line.length == 0)
409       throw new IllegalArgumentException(
410           "Line array must contain at least one vertex");
411     // this handles the case of length = 1
412     double minDistance = p.distance(line[0]);
413     for (int i = 0; i < line.length - 1; i++) {
414       double dist = distancePointLine(p, line[i], line[i + 1]);
415       if (dist < minDistance) {
416         minDistance = dist;
417       }
418     }
419     return minDistance;
420   }
421  
422   /**
423    * Computes the distance from a line segment AB to a line segment CD
424    * 
425    * Note: NON-ROBUST!
426    * 
427    * @param A
428    *          a point of one line
429    * @param B
430    *          the second point of (must be different to A)
431    * @param C
432    *          one point of the line
433    * @param D
434    *          another point of the line (must be different to A)
435    * @deprecated Use
436    *             {@link Distance#segmentToSegment(Coordinate, Coordinate, Coordinate, Coordinate)}
437    *             instead.
438    */
439   public static double distanceLineLine(Coordinate A, Coordinate B,
440       Coordinate C, Coordinate D)
441   {
442     // check for zero-length segments
443     if (A.equals(B))
444       return distancePointLine(A, C, D);
445     if (C.equals(D))
446       return distancePointLine(D, A, B);
447  
448     // AB and CD are line segments
449     /*
450      * from comp.graphics.algo
451      * 
452      * Solving the above for r and s yields 
453      * 
454      *     (Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy) 
455      * r = ----------------------------- (eqn 1) 
456      *     (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
457      * 
458      *     (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)  
459      * s = ----------------------------- (eqn 2)
460      *     (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx) 
461      *     
462      * Let P be the position vector of the
463      * intersection point, then 
464      *   P=A+r(B-A) or 
465      *   Px=Ax+r(Bx-Ax) 
466      *   Py=Ay+r(By-Ay) 
467      * By examining the values of r & s, you can also determine some other limiting
468      * conditions: 
469      *   If 0<=r<=1 & 0<=s<=1, intersection exists 
470      *      r<0 or r>1 or s<0 or s>1 line segments do not intersect 
471      *   If the denominator in eqn 1 is zero, AB & CD are parallel 
472      *   If the numerator in eqn 1 is also zero, AB & CD are collinear.
473      */
474  
475     boolean noIntersection = false;
476     if (! Envelope.intersects(A, B, C, D)) {
477       noIntersection = true;
478     }
479     else {
480       double denom = (B.x - A.x) * (D.y - C.y) - (B.y - A.y) * (D.x - C.x);
481       
482       if (denom == 0) {
483         noIntersection = true;
484       }
485       else {
486         double r_num = (A.y - C.y) * (D.x - C.x) - (A.x - C.x) * (D.y - C.y);
487         double s_num = (A.y - C.y) * (B.x - A.x) - (A.x - C.x) * (B.y - A.y);
488         
489         double s = s_num / denom;
490         double r = r_num / denom;
491   
492         if ((r < 0) || (r > 1) || (s < 0) || (s > 1)) {
493           noIntersection = true;
494         }
495       }
496     }
497     if (noIntersection) {
498       return MathUtil.min(
499             distancePointLine(A, C, D),
500             distancePointLine(B, C, D),
501             distancePointLine(C, A, B),
502             distancePointLine(D, A, B));
503     }
504     // segments intersect
505     return 0.0; 
506   }
507  
508   /**
509    * Computes the signed area for a ring. The signed area is positive if the
510    * ring is oriented CW, negative if the ring is oriented CCW, and zero if the
511    * ring is degenerate or flat.
512    * 
513    * @param ring
514    *          the coordinates forming the ring
515    * @return the signed area of the ring
516    * @deprecated Use {@link Area#ofRing(Coordinate[])} or
517    *             {@link Area#ofRingSigned(Coordinate[])} instead.
518    */
519   public static double signedArea(Coordinate[] ring)
520   {
521     if (ring.length < 3)
522       return 0.0;
523     double sum = 0.0;
524     /**
525      * Based on the Shoelace formula.
526      * http://en.wikipedia.org/wiki/Shoelace_formula
527      */
528     double x0 = ring[0].x;
529     for (int i = 1; i < ring.length - 1; i++) {
530       double x = ring[i].x - x0;
531       double y1 = ring[i + 1].y;
532       double y2 = ring[i - 1].y;
533       sum += x * (y2 - y1);
534     }
535     return sum / 2.0;
536   }
537  
538   /**
539    * Computes the signed area for a ring. The signed area is:
540    * <ul>
541    * <li>positive if the ring is oriented CW
542    * <li>negative if the ring is oriented CCW
543    * <li>zero if the ring is degenerate or flat
544    * </ul>
545    * 
546    * @param ring
547    *          the coordinates forming the ring
548    * @return the signed area of the ring
549    * @deprecated Use {@link Area#ofRing(CoordinateSequence)} or
550    *             {@link Area#ofRingSigned(CoordinateSequence)} instead.
551    */
552   public static double signedArea(CoordinateSequence ring)
553   {
554     int n = ring.size();
555     if (n < 3)
556       return 0.0;
557     /**
558      * Based on the Shoelace formula.
559      * http://en.wikipedia.org/wiki/Shoelace_formula
560      */
561     Coordinate p0 = new Coordinate();
562     Coordinate p1 = new Coordinate();
563     Coordinate p2 = new Coordinate();
564     ring.getCoordinate(0, p1);
565     ring.getCoordinate(1, p2);
566     double x0 = p1.x;
567     p2.x -= x0;
568     double sum = 0.0;
569     for (int i = 1; i < n - 1; i++) {
570       p0.y = p1.y;
571       p1.x = p2.x;
572       p1.y = p2.y;
573       ring.getCoordinate(i + 1, p2);
574       p2.x -= x0;
575       sum += p1.x * (p0.y - p2.y);
576     }
577     return sum / 2.0;
578   }
579  
580   /**
581    * Computes the length of a linestring specified by a sequence of points.
582    * 
583    * @param pts
584    *          the points specifying the linestring
585    * @return the length of the linestring
586    * @deprecated Use {@link Length#ofLine(CoordinateSequence)} instead.
587    */
588   public static double length(CoordinateSequence pts)
589   {
590     // optimized for processing CoordinateSequences
591     int n = pts.size();
592     if (n <= 1)
593       return 0.0;
594  
595     double len = 0.0;
596  
597     Coordinate p = new Coordinate();
598     pts.getCoordinate(0, p);
599     double x0 = p.x;
600     double y0 = p.y;
601  
602     for (int i = 1; i < n; i++) {
603       pts.getCoordinate(i, p);
604       double x1 = p.x;
605       double y1 = p.y;
606       double dx = x1 - x0;
607       double dy = y1 - y0;
608  
609       len += Math.sqrt(dx * dx + dy * dy);
610  
611       x0 = x1;
612       y0 = y1;
613     }
614     return len;
615   }
616  
617 }
618