Class Distance3DOp

1	/*
2	* Copyright (c) 2016 Martin Davis.
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.operation.distance3d;
13
14	import org.locationtech.jts.algorithm.CGAlgorithms3D;
15	import org.locationtech.jts.geom.Coordinate;
16	import org.locationtech.jts.geom.CoordinateSequence;
17	import org.locationtech.jts.geom.Geometry;
18	import org.locationtech.jts.geom.GeometryCollection;
19	import org.locationtech.jts.geom.LineSegment;
20	import org.locationtech.jts.geom.LineString;
21	import org.locationtech.jts.geom.Point;
22	import org.locationtech.jts.geom.Polygon;
23	import org.locationtech.jts.operation.distance.GeometryLocation;
24
25	/**
26	* Find two points on two 3D {@link Geometry}s which lie within a given distance,
27	* or else are the nearest points on the geometries (in which case this also
28	* provides the distance between the geometries).
29	* <p>
30	* 3D geometries have vertex Z ordinates defined.
31	* 3D {@link Polygon}s are assumed to lie in a single plane (which is enforced if not actually the case).
32	* 3D {@link LineString}s and {@link Point}s may have any configuration.
33	* <p>
34	* The distance computation also finds a pair of points in the input geometries
35	* which have the minimum distance between them. If a point lies in the interior
36	* of a line segment, the coordinate computed is a close approximation to the
37	* exact point for X and Y ordinates. Z ordinate is not interpolated.
38	* <p>
39	* The algorithms used are straightforward O(n^2) comparisons. This worst-case
40	* performance could be improved on by using Voronoi techniques or spatial
41	* indexes.
42	*
43	* @version 1.7
44	*/
45	public class Distance3DOp {
46	/**
47	* Compute the distance between the nearest points of two geometries.
48	*
49	* @param g0
50	* a {@link Geometry}
51	* @param g1
52	* another {@link Geometry}
53	* @return the distance between the geometries
54	*/
55	public static double distance(Geometry g0, Geometry g1) {
56	Distance3DOp distOp = new Distance3DOp(g0, g1);
57	return distOp.distance();
58	}
59
60	/**
61	* Test whether two geometries lie within a given distance of each other.
62	*
63	* @param g0
64	* a {@link Geometry}
65	* @param g1
66	* another {@link Geometry}
67	* @param distance
68	* the distance to test
69	* @return true if g0.distance(g1) <= distance
70	*/
71	public static boolean isWithinDistance(Geometry g0, Geometry g1,
72	double distance) {
73	Distance3DOp distOp = new Distance3DOp(g0, g1, distance);
74	return distOp.distance() <= distance;
75	}
76
77	/**
78	* Compute the the nearest points of two geometries. The points are
79	* presented in the same order as the input Geometries.
80	*
81	* @param g0
82	* a {@link Geometry}
83	* @param g1
84	* another {@link Geometry}
85	* @return the nearest points in the geometries
86	*/
87	public static Coordinate[] nearestPoints(Geometry g0, Geometry g1) {
88	Distance3DOp distOp = new Distance3DOp(g0, g1);
89	return distOp.nearestPoints();
90	}
91
92	// input
93	private Geometry[] geom;
94	private double terminateDistance = 0.0;
95	// working
96	private GeometryLocation[] minDistanceLocation;
97	private double minDistance = Double.MAX_VALUE;
98	private boolean isDone = false;
99
100	/**
101	* Constructs a DistanceOp that computes the distance and nearest points
102	* between the two specified geometries.
103	*
104	* @param g0
105	* a Geometry
106	* @param g1
107	* a Geometry
108	*/
109	public Distance3DOp(Geometry g0, Geometry g1) {
110	this(g0, g1, 0.0);
111	}
112
113	/**
114	* Constructs a DistanceOp that computes the distance and nearest points
115	* between the two specified geometries.
116	*
117	* @param g0
118	* a Geometry
119	* @param g1
120	* a Geometry
121	* @param terminateDistance
122	* the distance on which to terminate the search
123	*/
124	public Distance3DOp(Geometry g0, Geometry g1, double terminateDistance) {
125	this.geom = new Geometry[2];
126	geom[0] = g0;
127	geom[1] = g1;
128	this.terminateDistance = terminateDistance;
129	}
130
131	/**
132	* Report the distance between the nearest points on the input geometries.
133	*
134	* @return the distance between the geometries, or 0 if either input geometry is empty
135	* @throws IllegalArgumentException
136	* if either input geometry is null
137	*/
138	public double distance() {
139	if (geom[0] == null \|\| geom[1] == null)
140	throw new IllegalArgumentException(
141	"null geometries are not supported");
142	if (geom[0].isEmpty() \|\| geom[1].isEmpty())
143	return 0.0;
144
145	computeMinDistance();
146	return minDistance;
147	}
148
149	/**
150	* Report the coordinates of the nearest points in the input geometries. The
151	* points are presented in the same order as the input Geometries.
152	*
153	* @return a pair of {@link Coordinate}s of the nearest points
154	*/
155	public Coordinate[] nearestPoints() {
156	computeMinDistance();
157	Coordinate[] nearestPts = new Coordinate[] {
158	minDistanceLocation[0].getCoordinate(),
159	minDistanceLocation[1].getCoordinate() };
160	return nearestPts;
161	}
162
163	/**
164	* Report the locations of the nearest points in the input geometries. The
165	* locations are presented in the same order as the input Geometries.
166	*
167	* @return a pair of {@link GeometryLocation}s for the nearest points
168	*/
169	public GeometryLocation[] nearestLocations() {
170	computeMinDistance();
171	return minDistanceLocation;
172	}
173
174	private void updateDistance(double dist,
175	GeometryLocation loc0, GeometryLocation loc1,
176	boolean flip) {
177	this.minDistance = dist;
178	int index = flip ? 1 : 0;
179	minDistanceLocation[index] = loc0;
180	minDistanceLocation[1-index] = loc1;
181	if (minDistance < terminateDistance)
182	isDone = true;
183	}
184
185	private void computeMinDistance() {
186	// only compute once
187	if (minDistanceLocation != null)
188	return;
189	minDistanceLocation = new GeometryLocation[2];
190
191	int geomIndex = mostPolygonalIndex();
192	boolean flip = geomIndex == 1;
193	computeMinDistanceMultiMulti(geom[geomIndex], geom[1-geomIndex], flip);
194	}
195
196	/**
197	* Finds the index of the "most polygonal" input geometry.
198	* This optimizes the computation of the best-fit plane,
199	* since it is cached only for the left-hand geometry.
200	*
201	* @return the index of the most polygonal geometry
202	*/
203	private int mostPolygonalIndex() {
204	int dim0 = geom[0].getDimension();
205	int dim1 = geom[1].getDimension();
206	if (dim0 >= 2 && dim1 >= 2) {
207	if (geom[0].getNumPoints() > geom[1].getNumPoints())
208	return 0;
209	return 1;
210	}
211	// no more than one is dim 2
212	if (dim0 >= 2) return 0;
213	if (dim1 >= 2) return 1;
214	// both dim <= 1 - don't flip
215	return 0;
216	}
217
218	private void computeMinDistanceMultiMulti(Geometry g0, Geometry g1, boolean flip) {
219	if (g0 instanceof GeometryCollection) {
220	int n = g0.getNumGeometries();
221	for (int i = 0; i < n; i++) {
222	Geometry g = g0.getGeometryN(i);
223	computeMinDistanceMultiMulti(g, g1, flip);
224	if (isDone) return;
225	}
226	}
227	else {
228	// handle case of multigeom component being empty
229	if (g0.isEmpty())
230	return;
231
232	// compute planar polygon only once for efficiency
233	if (g0 instanceof Polygon) {
234	computeMinDistanceOneMulti(polyPlane(g0), g1, flip);
235	}
236	else
237	computeMinDistanceOneMulti(g0, g1, flip);
238	}
239	}
240
241	private void computeMinDistanceOneMulti(Geometry g0, Geometry g1, boolean flip) {
242	if (g1 instanceof GeometryCollection) {
243	int n = g1.getNumGeometries();
244	for (int i = 0; i < n; i++) {
245	Geometry g = g1.getGeometryN(i);
246	computeMinDistanceOneMulti(g0, g, flip);
247	if (isDone) return;
248	}
249	}
250	else {
251	computeMinDistance(g0, g1, flip);
252	}
253	}
254
255	private void computeMinDistanceOneMulti(PlanarPolygon3D poly, Geometry geom, boolean flip) {
256	if (geom instanceof GeometryCollection) {
257	int n = geom.getNumGeometries();
258	for (int i = 0; i < n; i++) {
259	Geometry g = geom.getGeometryN(i);
260	computeMinDistanceOneMulti(poly, g, flip);
261	if (isDone) return;
262	}
263	}
264	else {
265	if (geom instanceof Point) {
266	computeMinDistancePolygonPoint(poly, (Point) geom, flip);
267	return;
268	}
269	if (geom instanceof LineString) {
270	computeMinDistancePolygonLine(poly, (LineString) geom, flip);
271	return;
272	}
273	if (geom instanceof Polygon) {
274	computeMinDistancePolygonPolygon(poly, (Polygon) geom, flip);
275	return;
276	}
277	}
278	}
279
280	/**
281	* Convenience method to create a Plane3DPolygon
282	* @param poly
283	* @return
284	*/
285	private static PlanarPolygon3D polyPlane(Geometry poly)
286	{
287	return new PlanarPolygon3D((Polygon) poly);
288	}
289
290	private void computeMinDistance(Geometry g0, Geometry g1, boolean flip) {
291	if (g0 instanceof Point) {
292	if (g1 instanceof Point) {
293	computeMinDistancePointPoint((Point) g0, (Point) g1, flip);
294	return;
295	}
296	if (g1 instanceof LineString) {
297	computeMinDistanceLinePoint((LineString) g1, (Point) g0, ! flip);
298	return;
299	}
300	if (g1 instanceof Polygon) {
301	computeMinDistancePolygonPoint(polyPlane(g1), (Point) g0, ! flip);
302	return;
303	}
304	}
305	if (g0 instanceof LineString) {
306	if (g1 instanceof Point) {
307	computeMinDistanceLinePoint((LineString) g0, (Point) g1, flip);
308	return;
309	}
310	if (g1 instanceof LineString) {
311	computeMinDistanceLineLine((LineString) g0, (LineString) g1, flip);
312	return;
313	}
314	if (g1 instanceof Polygon) {
315	computeMinDistancePolygonLine(polyPlane(g1), (LineString) g0, ! flip);
316	return;
317	}
318	}
319	if (g0 instanceof Polygon) {
320	if (g1 instanceof Point) {
321	computeMinDistancePolygonPoint(polyPlane(g0), (Point) g1, flip);
322	return;
323	}
324	if (g1 instanceof LineString) {
325	computeMinDistancePolygonLine(polyPlane(g0), (LineString) g1, flip);
326	return;
327	}
328	if (g1 instanceof Polygon) {
329	computeMinDistancePolygonPolygon(polyPlane(g0), (Polygon) g1, flip);
330	return;
331	}
332	}
333	}
334
335	/**
336	* Computes distance between two polygons.
337	*
338	* To compute the distance, compute the distance
339	* between the rings of one polygon and the other polygon,
340	* and vice-versa.
341	* If the polygons intersect, then at least one ring must
342	* intersect the other polygon.
343	* Note that it is NOT sufficient to test only the shell rings.
344	* A counter-example is a "figure-8" polygon A
345	* and a simple polygon B at right angles to A, with the ring of B
346	* passing through the holes of A.
347	* The polygons intersect,
348	* but A's shell does not intersect B, and B's shell does not intersect A.
349	*
350	* @param poly0
351	* @param poly1
352	* @param geomIndex
353	*/
354	private void computeMinDistancePolygonPolygon(PlanarPolygon3D poly0, Polygon poly1,
355	boolean flip) {
356	computeMinDistancePolygonRings(poly0, poly1, flip);
357	if (isDone) return;
358	PlanarPolygon3D polyPlane1 = new PlanarPolygon3D(poly1);
359	computeMinDistancePolygonRings(polyPlane1, poly0.getPolygon(), flip);
360	}
361
362	/**
363	* Compute distance between a polygon and the rings of another.
364	*
365	* @param poly
366	* @param ringPoly
367	* @param geomIndex
368	*/
369	private void computeMinDistancePolygonRings(PlanarPolygon3D poly, Polygon ringPoly,
370	boolean flip) {
371	// compute shell ring
372	computeMinDistancePolygonLine(poly, ringPoly.getExteriorRing(), flip);
373	if (isDone) return;
374	// compute hole rings
375	int nHole = ringPoly.getNumInteriorRing();
376	for (int i = 0; i < nHole; i++) {
377	computeMinDistancePolygonLine(poly, ringPoly.getInteriorRingN(i), flip);
378	if (isDone) return;
379	}
380	}
381
382	private void computeMinDistancePolygonLine(PlanarPolygon3D poly,LineString line,
383	boolean flip) {
384
385	// first test if line intersects polygon
386	Coordinate intPt = intersection(poly, line);
387	if (intPt != null) {
388	updateDistance(0,
389	new GeometryLocation(poly.getPolygon(), 0, intPt),
390	new GeometryLocation(line, 0, intPt),
391	flip
392	);
393	return;
394	}
395
396	// if no intersection, then compute line distance to polygon rings
397	computeMinDistanceLineLine(poly.getPolygon().getExteriorRing(), line, flip);
398	if (isDone) return;
399	int nHole = poly.getPolygon().getNumInteriorRing();
400	for (int i = 0; i < nHole; i++) {
401	computeMinDistanceLineLine(poly.getPolygon().getInteriorRingN(i), line, flip);
402	if (isDone) return;
403	}
404	}
405
406	private Coordinate intersection(PlanarPolygon3D poly,LineString line) {
407	CoordinateSequence seq = line.getCoordinateSequence();
408	if (seq.size() == 0)
409	return null;
410
411	// start point of line
412	Coordinate p0 = new Coordinate();
413	seq.getCoordinate(0, p0);
414	double d0 = poly.getPlane().orientedDistance(p0);
415
416	// for each segment in the line
417	Coordinate p1 = new Coordinate();
418	for (int i = 0; i < seq.size() - 1; i++) {
419	seq.getCoordinate(i, p0);
420	seq.getCoordinate(i + 1, p1);
421	double d1 = poly.getPlane().orientedDistance(p1);
422
423	/**
424	* If the oriented distances of the segment endpoints have the same sign,
425	* the segment does not cross the plane, and is skipped.
426	*/
427	if (d0 * d1 > 0)
428	continue;
429
430	/**
431	* Compute segment-plane intersection point
432	* which is then used for a point-in-polygon test.
433	* The endpoint distances to the plane d0 and d1
434	* give the proportional distance of the intersection point
435	* along the segment.
436	*/
437	Coordinate intPt = segmentPoint(p0, p1, d0, d1);
438	// Coordinate intPt = polyPlane.intersection(p0, p1, s0, s1);
439	if (poly.intersects(intPt)) {
440	return intPt;
441	}
442
443	// shift to next segment
444	d0 = d1;
445	}
446	return null;
447	}
448
449	private void computeMinDistancePolygonPoint(PlanarPolygon3D polyPlane, Point point,
450	boolean flip) {
451	Coordinate pt = point.getCoordinate();
452
453	LineString shell = polyPlane.getPolygon().getExteriorRing();
454	if (polyPlane.intersects(pt, shell)) {
455	// point is either inside or in a hole
456
457	int nHole = polyPlane.getPolygon().getNumInteriorRing();
458	for (int i = 0; i < nHole; i++) {
459	LineString hole = polyPlane.getPolygon().getInteriorRingN(i);
460	if (polyPlane.intersects(pt, hole)) {
461	computeMinDistanceLinePoint(hole, point, flip);
462	return;
463	}
464	}
465	// point is in interior of polygon
466	// distance is distance to polygon plane
467	double dist = Math.abs(polyPlane.getPlane().orientedDistance(pt));
468	updateDistance(dist,
469	new GeometryLocation(polyPlane.getPolygon(), 0, pt),
470	new GeometryLocation(point, 0, pt),
471	flip
472	);
473	}
474	// point is outside polygon, so compute distance to shell linework
475	computeMinDistanceLinePoint(shell, point, flip);
476	}
477
478	private void computeMinDistanceLineLine(LineString line0, LineString line1,
479	boolean flip) {
480	Coordinate[] coord0 = line0.getCoordinates();
481	Coordinate[] coord1 = line1.getCoordinates();
482	// brute force approach!
483	for (int i = 0; i < coord0.length - 1; i++) {
484	for (int j = 0; j < coord1.length - 1; j++) {
485	double dist = CGAlgorithms3D.distanceSegmentSegment(coord0[i],
486	coord0[i + 1], coord1[j], coord1[j + 1]);
487	if (dist < minDistance) {
488	minDistance = dist;
489	// TODO: compute closest pts in 3D
490	LineSegment seg0 = new LineSegment(coord0[i], coord0[i + 1]);
491	LineSegment seg1 = new LineSegment(coord1[j], coord1[j + 1]);
492	Coordinate[] closestPt = seg0.closestPoints(seg1);
493	updateDistance(dist,
494	new GeometryLocation(line0, i, closestPt[0]),
495	new GeometryLocation(line1, j, closestPt[1]),
496	flip
497	);
498	}
499	if (isDone) return;
500	}
501	}
502	}
503
504	private void computeMinDistanceLinePoint(LineString line,Point point,
505	boolean flip) {
506	Coordinate[] lineCoord = line.getCoordinates();
507	Coordinate coord = point.getCoordinate();
508	// brute force approach!
509	for (int i = 0; i < lineCoord.length - 1; i++) {
510	double dist = CGAlgorithms3D.distancePointSegment(coord, lineCoord[i],
511	lineCoord[i + 1]);
512	if (dist < minDistance) {
513	LineSegment seg = new LineSegment(lineCoord[i], lineCoord[i + 1]);
514	Coordinate segClosestPoint = seg.closestPoint(coord);
515	updateDistance(dist,
516	new GeometryLocation(line, i, segClosestPoint),
517	new GeometryLocation(point, 0, coord),
518	flip);
519	}
520	if (isDone) return;
521	}
522	}
523
524	private void computeMinDistancePointPoint(Point point0, Point point1, boolean flip) {
525	double dist = CGAlgorithms3D.distance(
526	point0.getCoordinate(),
527	point1.getCoordinate());
528	if (dist < minDistance) {
529	updateDistance(dist,
530	new GeometryLocation(point0, 0, point0.getCoordinate()),
531	new GeometryLocation(point1, 0, point1.getCoordinate()),
532	flip);
533	}
534	}
535
536	/**
537	* Computes a point at a distance along a segment
538	* specified by two relatively proportional values.
539	* The fractional distance along the segment is d0/(d0+d1).
540	*
541	* @param p0
542	* start point of the segment
543	* @param p1
544	* end point of the segment
545	* @param d0
546	* proportional distance from start point to computed point
547	* @param d1
548	* proportional distance from computed point to end point
549	* @return the computed point
550	*/
551	private static Coordinate segmentPoint(Coordinate p0, Coordinate p1, double d0,
552	double d1) {
553	if (d0 <= 0) return new Coordinate(p0);
554	if (d1 <= 0) return new Coordinate(p1);
555
556	double f = Math.abs(d0) / (Math.abs(d0) + Math.abs(d1));
557	double intx = p0.x + f * (p1.x - p0.x);
558	double inty = p0.y + f * (p1.y - p0.y);
559	double intz = p0.getZ() + f * (p1.getZ() - p0.getZ());
560	return new Coordinate(intx, inty, intz);
561	}
562
563
564
565
566	}
567