Class InteriorPointArea

1
2	/*
3	* Copyright (c) 2016 Vivid Solutions.
4	*
5	* All rights reserved. This program and the accompanying materials
6	* are made available under the terms of the Eclipse Public License 2.0
7	* and Eclipse Distribution License v. 1.0 which accompanies this distribution.
8	* The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
9	* and the Eclipse Distribution License is available at
10	*
11	* http://www.eclipse.org/org/documents/edl-v10.php.
12	*/
13	package org.locationtech.jts.algorithm;
14
15	import java.util.ArrayList;
16	import java.util.List;
17
18	import org.locationtech.jts.geom.Coordinate;
19	import org.locationtech.jts.geom.CoordinateSequence;
20	import org.locationtech.jts.geom.Envelope;
21	import org.locationtech.jts.geom.Geometry;
22	import org.locationtech.jts.geom.GeometryCollection;
23	import org.locationtech.jts.geom.LineString;
24	import org.locationtech.jts.geom.LinearRing;
25	import org.locationtech.jts.geom.Polygon;
26	import org.locationtech.jts.util.Assert;
27
28	/**
29	* Computes a point in the interior of an areal geometry.
30	* The point will lie in the geometry interior
31	* in all except certain pathological cases.
32	*
33	* <h2>Algorithm</h2>
34	* For each input polygon:
35	* <ul>
36	* <li>Determine a horizontal scan line on which the interior
37	* point will be located.
38	* To increase the chance of the scan line
39	* having non-zero-width intersection with the polygon
40	* the scan line Y ordinate is chosen to be near the centre of the polygon's
41	* Y extent but distinct from all of vertex Y ordinates.
42	* <li>Compute the sections of the scan line
43	* which lie in the interior of the polygon.
44	* <li>Choose the widest interior section
45	* and take its midpoint as the interior point.
46	* </ul>
47	* The final interior point is chosen as
48	* the one occurring in the widest interior section.
49	* <p>
50	* This algorithm is a tradeoff between performance
51	* and point quality (where points further from the geometry
52	* boundary are considered to be higher quality)
53	* Priority is given to performance.
54	* This means that the computed interior point
55	* may not be suitable for some uses
56	* (such as label positioning).
57	* <p>
58	* The algorithm handles some kinds of invalid/degenerate geometry,
59	* including zero-area and self-intersecting polygons.
60	* <p>
61	* Empty geometry is handled by returning a <code>null</code> point.
62	*
63	* <h3>KNOWN BUGS</h3>
64	* <ul>
65	* <li>If a fixed precision model is used, in some cases this method may return
66	* a point which does not lie in the interior.
67	* <li>If the input polygon is <i>extremely</i> narrow the computed point
68	* may not lie in the interior of the polygon.
69	* </ul>
70	*
71	* @version 1.17
72	*/
73	public class InteriorPointArea {
74
75	/**
76	* Computes an interior point for the
77	* polygonal components of a Geometry.
78	*
79	* @param geom the geometry to compute
80	* @return the computed interior point,
81	* or <code>null</code> if the geometry has no polygonal components
82	*/
83	public static Coordinate getInteriorPoint(Geometry geom) {
84	InteriorPointArea intPt = new InteriorPointArea(geom);
85	return intPt.getInteriorPoint();
86	}
87
88	private static double avg(double a, double b) {
89	return (a + b) / 2.0;
90	}
91
92	private Coordinate interiorPoint = null;
93	private double maxWidth = -1;
94
95	/**
96	* Creates a new interior point finder for an areal geometry.
97	*
98	* @param g an areal geometry
99	*/
100	public InteriorPointArea(Geometry g) {
101	process(g);
102	}
103
104	/**
105	* Gets the computed interior point.
106	*
107	* @return the coordinate of an interior point
108	* or <code>null</code> if the input geometry is empty
109	*/
110	public Coordinate getInteriorPoint() {
111	return interiorPoint;
112	}
113
114	/**
115	* Processes a geometry to determine
116	* the best interior point for
117	* all component polygons.
118	*
119	* @param geom the geometry to process
120	*/
121	private void process(Geometry geom) {
122	if ( geom.isEmpty() )
123	return;
124
125	if ( geom instanceof Polygon ) {
126	processPolygon((Polygon) geom);
127	} else if ( geom instanceof GeometryCollection ) {
128	GeometryCollection gc = (GeometryCollection) geom;
129	for (int i = 0; i < gc.getNumGeometries(); i++) {
130	process(gc.getGeometryN(i));
131	}
132	}
133	}
134
135	/**
136	* Computes an interior point of a component Polygon
137	* and updates current best interior point
138	* if appropriate.
139	*
140	* @param polygon the polygon to process
141	*/
142	private void processPolygon(Polygon polygon) {
143	InteriorPointPolygon intPtPoly = new InteriorPointPolygon(polygon);
144	intPtPoly.process();
145	double width = intPtPoly.getWidth();
146	if ( width > maxWidth ) {
147	maxWidth = width;
148	interiorPoint = intPtPoly.getInteriorPoint();
149	}
150	}
151
152	/**
153	* Computes an interior point in a single {@link Polygon},
154	* as well as the width of the scan-line section it occurs in
155	* to allow choosing the widest section occurrence.
156	*
157	* @author mdavis
158	*
159	*/
160	private static class InteriorPointPolygon {
161	private Polygon polygon;
162	private double interiorPointY;
163	private double interiorSectionWidth = 0.0;
164	private Coordinate interiorPoint = null;
165
166	/**
167	* Creates a new InteriorPointPolygon instance.
168	*
169	* @param polygon the polygon to test
170	*/
171	public InteriorPointPolygon(Polygon polygon) {
172	this.polygon = polygon;
173	interiorPointY = ScanLineYOrdinateFinder.getScanLineY(polygon);
174	}
175
176	/**
177	* Gets the computed interior point.
178	*
179	* @return the interior point coordinate,
180	* or <code>null</code> if the input geometry is empty
181	*/
182	public Coordinate getInteriorPoint() {
183	return interiorPoint;
184	}
185
186	/**
187	* Gets the width of the scanline section containing the interior point.
188	* Used to determine the best point to use.
189	*
190	* @return the width
191	*/
192	public double getWidth() {
193	return interiorSectionWidth;
194	}
195
196	/**
197	* Compute the interior point.
198	*
199	*/
200	public void process() {
201	/**
202	* This results in returning a null Coordinate
203	*/
204	if (polygon.isEmpty()) return;
205
206	/**
207	* set default interior point in case polygon has zero area
208	*/
209	interiorPoint = new Coordinate(polygon.getCoordinate());
210
211	List<Double> crossings = new ArrayList<Double>();
212	scanRing((LinearRing) polygon.getExteriorRing(), crossings);
213	for (int i = 0; i < polygon.getNumInteriorRing(); i++) {
214	scanRing((LinearRing) polygon.getInteriorRingN(i), crossings);
215	}
216	findBestMidpoint(crossings);
217	}
218
219	private void scanRing(LinearRing ring, List<Double> crossings) {
220	// skip rings which don't cross scan line
221	if ( !intersectsHorizontalLine(ring.getEnvelopeInternal(), interiorPointY) )
222	return;
223
224	CoordinateSequence seq = ring.getCoordinateSequence();
225	for (int i = 1; i < seq.size(); i++) {
226	Coordinate ptPrev = seq.getCoordinate(i - 1);
227	Coordinate pt = seq.getCoordinate(i);
228	addEdgeCrossing(ptPrev, pt, interiorPointY, crossings);
229	}
230	}
231
232	private void addEdgeCrossing(Coordinate p0, Coordinate p1, double scanY, List<Double> crossings) {
233	// skip non-crossing segments
234	if ( !intersectsHorizontalLine(p0, p1, scanY) )
235	return;
236	if (! isEdgeCrossingCounted(p0, p1, scanY) )
237	return;
238
239	// edge intersects scan line, so add a crossing
240	double xInt = intersection(p0, p1, scanY);
241	crossings.add(xInt);
242	//checkIntersectionDD(p0, p1, scanY, xInt);
243	}
244
245	/**
246	* Finds the midpoint of the widest interior section.
247	* Sets the {@link #interiorPoint} location
248	* and the {@link #interiorSectionWidth}
249	*
250	* @param crossings the list of scan-line crossing X ordinates
251	*/
252	private void findBestMidpoint(List<Double> crossings) {
253	// zero-area polygons will have no crossings
254	if (crossings.size() == 0) return;
255
256	// TODO: is there a better way to verify the crossings are correct?
257	Assert.isTrue(0 == crossings.size() % 2, "Interior Point robustness failure: odd number of scanline crossings");
258
259	crossings.sort(Double::compare);
260	/*
261	* Entries in crossings list are expected to occur in pairs representing a
262	* section of the scan line interior to the polygon (which may be zero-length)
263	*/
264	for (int i = 0; i < crossings.size(); i += 2) {
265	double x1 = crossings.get(i);
266	// crossings count must be even so this should be safe
267	double x2 = crossings.get(i + 1);
268
269	double width = x2 - x1;
270	if ( width > interiorSectionWidth ) {
271	interiorSectionWidth = width;
272	double interiorPointX = avg(x1, x2);
273	interiorPoint = new Coordinate(interiorPointX, interiorPointY);
274	}
275	}
276	}
277
278	/**
279	* Tests if an edge intersection contributes to the crossing count.
280	* Some crossing situations are not counted,
281	* to ensure that the list of crossings
282	* captures strict inside/outside topology.
283	*
284	* @param p0 an endpoint of the segment
285	* @param p1 an endpoint of the segment
286	* @param scanY the Y-ordinate of the horizontal line
287	* @return true if the edge crossing is counted
288	*/
289	private static boolean isEdgeCrossingCounted(Coordinate p0, Coordinate p1, double scanY) {
290	double y0 = p0.getY();
291	double y1 = p1.getY();
292	// skip horizontal lines
293	if ( y0 == y1 )
294	return false;
295	// handle cases where vertices lie on scan-line
296	// downward segment does not include start point
297	if ( y0 == scanY && y1 < scanY )
298	return false;
299	// upward segment does not include endpoint
300	if ( y1 == scanY && y0 < scanY )
301	return false;
302	return true;
303	}
304
305	/**
306	* Computes the intersection of a segment with a horizontal line.
307	* The segment is expected to cross the horizontal line
308	* - this condition is not checked.
309	* Computation uses regular double-precision arithmetic.
310	* Test seems to indicate this is as good as using DD arithmetic.
311	*
312	* @param p0 an endpoint of the segment
313	* @param p1 an endpoint of the segment
314	* @param Y the Y-ordinate of the horizontal line
315	* @return
316	*/
317	private static double intersection(Coordinate p0, Coordinate p1, double Y) {
318	double x0 = p0.getX();
319	double x1 = p1.getX();
320
321	if ( x0 == x1 )
322	return x0;
323
324	// Assert: segDX is non-zero, due to previous equality test
325	double segDX = x1 - x0;
326	double segDY = p1.getY() - p0.getY();
327	double m = segDY / segDX;
328	double x = x0 + ((Y - p0.getY()) / m);
329	return x;
330	}
331
332	/**
333	* Tests if an envelope intersects a horizontal line.
334	*
335	* @param env the envelope to test
336	* @param y the Y-ordinate of the horizontal line
337	* @return true if the envelope and line intersect
338	*/
339	private static boolean intersectsHorizontalLine(Envelope env, double y) {
340	if ( y < env.getMinY() )
341	return false;
342	if ( y > env.getMaxY() )
343	return false;
344	return true;
345	}
346
347	/**
348	* Tests if a line segment intersects a horizontal line.
349	*
350	* @param p0 a segment endpoint
351	* @param p1 a segment endpoint
352	* @param y the Y-ordinate of the horizontal line
353	* @return true if the segment and line intersect
354	*/
355	private static boolean intersectsHorizontalLine(Coordinate p0, Coordinate p1, double y) {
356	// both ends above?
357	if ( p0.getY() > y && p1.getY() > y )
358	return false;
359	// both ends below?
360	if ( p0.getY() < y && p1.getY() < y )
361	return false;
362	// segment must intersect line
363	return true;
364	}
365
366	/*
367	// for testing only
368	private static void checkIntersectionDD(Coordinate p0, Coordinate p1, double scanY, double xInt) {
369	double xIntDD = intersectionDD(p0, p1, scanY);
370	System.out.println(
371	((xInt != xIntDD) ? ">>" : "")
372	+ "IntPt x - DP: " + xInt + ", DD: " + xIntDD
373	+ " y: " + scanY + " " + WKTWriter.toLineString(p0, p1) );
374	}
375
376	private static double intersectionDD(Coordinate p0, Coordinate p1, double Y) {
377	double x0 = p0.getX();
378	double x1 = p1.getX();
379
380	if ( x0 == x1 )
381	return x0;
382
383	DD segDX = DD.valueOf(x1).selfSubtract(x0);
384	// Assert: segDX is non-zero, due to previous equality test
385	DD segDY = DD.valueOf(p1.getY()).selfSubtract(p0.getY());
386	DD m = segDY.divide(segDX);
387	DD dy = DD.valueOf(Y).selfSubtract(p0.getY());
388	DD dx = dy.divide(m);
389	DD xInt = DD.valueOf(x0).selfAdd(dx);
390	return xInt.doubleValue();
391	}
392	*/
393	}
394
395	/**
396	* Finds a safe scan line Y ordinate by projecting
397	* the polygon segments
398	* to the Y axis and finding the
399	* Y-axis interval which contains the centre of the Y extent.
400	* The centre of
401	* this interval is returned as the scan line Y-ordinate.
402	* <p>
403	* Note that in the case of (degenerate, invalid)
404	* zero-area polygons the computed Y value
405	* may be equal to a vertex Y-ordinate.
406	*
407	* @author mdavis
408	*
409	*/
410	private static class ScanLineYOrdinateFinder {
411	public static double getScanLineY(Polygon poly) {
412	ScanLineYOrdinateFinder finder = new ScanLineYOrdinateFinder(poly);
413	return finder.getScanLineY();
414	}
415
416	private Polygon poly;
417
418	private double centreY;
419	private double hiY = Double.MAX_VALUE;
420	private double loY = -Double.MAX_VALUE;
421
422	public ScanLineYOrdinateFinder(Polygon poly) {
423	this.poly = poly;
424
425	// initialize using extremal values
426	hiY = poly.getEnvelopeInternal().getMaxY();
427	loY = poly.getEnvelopeInternal().getMinY();
428	centreY = avg(loY, hiY);
429	}
430
431	public double getScanLineY() {
432	process(poly.getExteriorRing());
433	for (int i = 0; i < poly.getNumInteriorRing(); i++) {
434	process(poly.getInteriorRingN(i));
435	}
436	double scanLineY = avg(hiY, loY);
437	return scanLineY;
438	}
439
440	private void process(LineString line) {
441	CoordinateSequence seq = line.getCoordinateSequence();
442	for (int i = 0; i < seq.size(); i++) {
443	double y = seq.getY(i);
444	updateInterval(y);
445	}
446	}
447
448	private void updateInterval(double y) {
449	if ( y <= centreY ) {
450	if ( y > loY )
451	loY = y;
452	} else if ( y > centreY ) {
453	if ( y < hiY ) {
454	hiY = y;
455	}
456	}
457	}
458	}
459	}
460