InteriorPointArea.java
/*
* Copyright (c) 2016 Vivid Solutions.
*
* All rights reserved. This program and the accompanying materials
* are made available under the terms of the Eclipse Public License 2.0
* and Eclipse Distribution License v. 1.0 which accompanies this distribution.
* The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
* and the Eclipse Distribution License is available at
*
* http://www.eclipse.org/org/documents/edl-v10.php.
*/
package org.locationtech.jts.algorithm;
import java.util.ArrayList;
import java.util.List;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.CoordinateSequence;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryCollection;
import org.locationtech.jts.geom.LineString;
import org.locationtech.jts.geom.LinearRing;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.util.Assert;
/**
* Computes a point in the interior of an areal geometry.
* The point will lie in the geometry interior
* in all except certain pathological cases.
*
* <h2>Algorithm</h2>
* For each input polygon:
* <ul>
* <li>Determine a horizontal scan line on which the interior
* point will be located.
* To increase the chance of the scan line
* having non-zero-width intersection with the polygon
* the scan line Y ordinate is chosen to be near the centre of the polygon's
* Y extent but distinct from all of vertex Y ordinates.
* <li>Compute the sections of the scan line
* which lie in the interior of the polygon.
* <li>Choose the widest interior section
* and take its midpoint as the interior point.
* </ul>
* The final interior point is chosen as
* the one occurring in the widest interior section.
* <p>
* This algorithm is a tradeoff between performance
* and point quality (where points further from the geometry
* boundary are considered to be higher quality)
* Priority is given to performance.
* This means that the computed interior point
* may not be suitable for some uses
* (such as label positioning).
* <p>
* The algorithm handles some kinds of invalid/degenerate geometry,
* including zero-area and self-intersecting polygons.
* <p>
* Empty geometry is handled by returning a <code>null</code> point.
*
* <h3>KNOWN BUGS</h3>
* <ul>
* <li>If a fixed precision model is used, in some cases this method may return
* a point which does not lie in the interior.
* <li>If the input polygon is <i>extremely</i> narrow the computed point
* may not lie in the interior of the polygon.
* </ul>
*
* @version 1.17
*/
public class InteriorPointArea {
/**
* Computes an interior point for the
* polygonal components of a Geometry.
*
* @param geom the geometry to compute
* @return the computed interior point,
* or <code>null</code> if the geometry has no polygonal components
*/
public static Coordinate getInteriorPoint(Geometry geom) {
InteriorPointArea intPt = new InteriorPointArea(geom);
return intPt.getInteriorPoint();
}
private static double avg(double a, double b) {
return (a + b) / 2.0;
}
private Coordinate interiorPoint = null;
private double maxWidth = -1;
/**
* Creates a new interior point finder for an areal geometry.
*
* @param g an areal geometry
*/
public InteriorPointArea(Geometry g) {
process(g);
}
/**
* Gets the computed interior point.
*
* @return the coordinate of an interior point
* or <code>null</code> if the input geometry is empty
*/
public Coordinate getInteriorPoint() {
return interiorPoint;
}
/**
* Processes a geometry to determine
* the best interior point for
* all component polygons.
*
* @param geom the geometry to process
*/
private void process(Geometry geom) {
if ( geom.isEmpty() )
return;
if ( geom instanceof Polygon ) {
processPolygon((Polygon) geom);
} else if ( geom instanceof GeometryCollection ) {
GeometryCollection gc = (GeometryCollection) geom;
for (int i = 0; i < gc.getNumGeometries(); i++) {
process(gc.getGeometryN(i));
}
}
}
/**
* Computes an interior point of a component Polygon
* and updates current best interior point
* if appropriate.
*
* @param polygon the polygon to process
*/
private void processPolygon(Polygon polygon) {
InteriorPointPolygon intPtPoly = new InteriorPointPolygon(polygon);
intPtPoly.process();
double width = intPtPoly.getWidth();
if ( width > maxWidth ) {
maxWidth = width;
interiorPoint = intPtPoly.getInteriorPoint();
}
}
/**
* Computes an interior point in a single {@link Polygon},
* as well as the width of the scan-line section it occurs in
* to allow choosing the widest section occurrence.
*
* @author mdavis
*
*/
private static class InteriorPointPolygon {
private Polygon polygon;
private double interiorPointY;
private double interiorSectionWidth = 0.0;
private Coordinate interiorPoint = null;
/**
* Creates a new InteriorPointPolygon instance.
*
* @param polygon the polygon to test
*/
public InteriorPointPolygon(Polygon polygon) {
this.polygon = polygon;
interiorPointY = ScanLineYOrdinateFinder.getScanLineY(polygon);
}
/**
* Gets the computed interior point.
*
* @return the interior point coordinate,
* or <code>null</code> if the input geometry is empty
*/
public Coordinate getInteriorPoint() {
return interiorPoint;
}
/**
* Gets the width of the scanline section containing the interior point.
* Used to determine the best point to use.
*
* @return the width
*/
public double getWidth() {
return interiorSectionWidth;
}
/**
* Compute the interior point.
*
*/
public void process() {
/**
* This results in returning a null Coordinate
*/
if (polygon.isEmpty()) return;
/**
* set default interior point in case polygon has zero area
*/
interiorPoint = new Coordinate(polygon.getCoordinate());
List<Double> crossings = new ArrayList<Double>();
scanRing((LinearRing) polygon.getExteriorRing(), crossings);
for (int i = 0; i < polygon.getNumInteriorRing(); i++) {
scanRing((LinearRing) polygon.getInteriorRingN(i), crossings);
}
findBestMidpoint(crossings);
}
private void scanRing(LinearRing ring, List<Double> crossings) {
// skip rings which don't cross scan line
if ( !intersectsHorizontalLine(ring.getEnvelopeInternal(), interiorPointY) )
return;
CoordinateSequence seq = ring.getCoordinateSequence();
for (int i = 1; i < seq.size(); i++) {
Coordinate ptPrev = seq.getCoordinate(i - 1);
Coordinate pt = seq.getCoordinate(i);
addEdgeCrossing(ptPrev, pt, interiorPointY, crossings);
}
}
private void addEdgeCrossing(Coordinate p0, Coordinate p1, double scanY, List<Double> crossings) {
// skip non-crossing segments
if ( !intersectsHorizontalLine(p0, p1, scanY) )
return;
if (! isEdgeCrossingCounted(p0, p1, scanY) )
return;
// edge intersects scan line, so add a crossing
double xInt = intersection(p0, p1, scanY);
crossings.add(xInt);
//checkIntersectionDD(p0, p1, scanY, xInt);
}
/**
* Finds the midpoint of the widest interior section.
* Sets the {@link #interiorPoint} location
* and the {@link #interiorSectionWidth}
*
* @param crossings the list of scan-line crossing X ordinates
*/
private void findBestMidpoint(List<Double> crossings) {
// zero-area polygons will have no crossings
if (crossings.size() == 0) return;
// TODO: is there a better way to verify the crossings are correct?
Assert.isTrue(0 == crossings.size() % 2, "Interior Point robustness failure: odd number of scanline crossings");
crossings.sort(Double::compare);
/*
* Entries in crossings list are expected to occur in pairs representing a
* section of the scan line interior to the polygon (which may be zero-length)
*/
for (int i = 0; i < crossings.size(); i += 2) {
double x1 = crossings.get(i);
// crossings count must be even so this should be safe
double x2 = crossings.get(i + 1);
double width = x2 - x1;
if ( width > interiorSectionWidth ) {
interiorSectionWidth = width;
double interiorPointX = avg(x1, x2);
interiorPoint = new Coordinate(interiorPointX, interiorPointY);
}
}
}
/**
* Tests if an edge intersection contributes to the crossing count.
* Some crossing situations are not counted,
* to ensure that the list of crossings
* captures strict inside/outside topology.
*
* @param p0 an endpoint of the segment
* @param p1 an endpoint of the segment
* @param scanY the Y-ordinate of the horizontal line
* @return true if the edge crossing is counted
*/
private static boolean isEdgeCrossingCounted(Coordinate p0, Coordinate p1, double scanY) {
double y0 = p0.getY();
double y1 = p1.getY();
// skip horizontal lines
if ( y0 == y1 )
return false;
// handle cases where vertices lie on scan-line
// downward segment does not include start point
if ( y0 == scanY && y1 < scanY )
return false;
// upward segment does not include endpoint
if ( y1 == scanY && y0 < scanY )
return false;
return true;
}
/**
* Computes the intersection of a segment with a horizontal line.
* The segment is expected to cross the horizontal line
* - this condition is not checked.
* Computation uses regular double-precision arithmetic.
* Test seems to indicate this is as good as using DD arithmetic.
*
* @param p0 an endpoint of the segment
* @param p1 an endpoint of the segment
* @param Y the Y-ordinate of the horizontal line
* @return
*/
private static double intersection(Coordinate p0, Coordinate p1, double Y) {
double x0 = p0.getX();
double x1 = p1.getX();
if ( x0 == x1 )
return x0;
// Assert: segDX is non-zero, due to previous equality test
double segDX = x1 - x0;
double segDY = p1.getY() - p0.getY();
double m = segDY / segDX;
double x = x0 + ((Y - p0.getY()) / m);
return x;
}
/**
* Tests if an envelope intersects a horizontal line.
*
* @param env the envelope to test
* @param y the Y-ordinate of the horizontal line
* @return true if the envelope and line intersect
*/
private static boolean intersectsHorizontalLine(Envelope env, double y) {
if ( y < env.getMinY() )
return false;
if ( y > env.getMaxY() )
return false;
return true;
}
/**
* Tests if a line segment intersects a horizontal line.
*
* @param p0 a segment endpoint
* @param p1 a segment endpoint
* @param y the Y-ordinate of the horizontal line
* @return true if the segment and line intersect
*/
private static boolean intersectsHorizontalLine(Coordinate p0, Coordinate p1, double y) {
// both ends above?
if ( p0.getY() > y && p1.getY() > y )
return false;
// both ends below?
if ( p0.getY() < y && p1.getY() < y )
return false;
// segment must intersect line
return true;
}
/*
// for testing only
private static void checkIntersectionDD(Coordinate p0, Coordinate p1, double scanY, double xInt) {
double xIntDD = intersectionDD(p0, p1, scanY);
System.out.println(
((xInt != xIntDD) ? ">>" : "")
+ "IntPt x - DP: " + xInt + ", DD: " + xIntDD
+ " y: " + scanY + " " + WKTWriter.toLineString(p0, p1) );
}
private static double intersectionDD(Coordinate p0, Coordinate p1, double Y) {
double x0 = p0.getX();
double x1 = p1.getX();
if ( x0 == x1 )
return x0;
DD segDX = DD.valueOf(x1).selfSubtract(x0);
// Assert: segDX is non-zero, due to previous equality test
DD segDY = DD.valueOf(p1.getY()).selfSubtract(p0.getY());
DD m = segDY.divide(segDX);
DD dy = DD.valueOf(Y).selfSubtract(p0.getY());
DD dx = dy.divide(m);
DD xInt = DD.valueOf(x0).selfAdd(dx);
return xInt.doubleValue();
}
*/
}
/**
* Finds a safe scan line Y ordinate by projecting
* the polygon segments
* to the Y axis and finding the
* Y-axis interval which contains the centre of the Y extent.
* The centre of
* this interval is returned as the scan line Y-ordinate.
* <p>
* Note that in the case of (degenerate, invalid)
* zero-area polygons the computed Y value
* may be equal to a vertex Y-ordinate.
*
* @author mdavis
*
*/
private static class ScanLineYOrdinateFinder {
public static double getScanLineY(Polygon poly) {
ScanLineYOrdinateFinder finder = new ScanLineYOrdinateFinder(poly);
return finder.getScanLineY();
}
private Polygon poly;
private double centreY;
private double hiY = Double.MAX_VALUE;
private double loY = -Double.MAX_VALUE;
public ScanLineYOrdinateFinder(Polygon poly) {
this.poly = poly;
// initialize using extremal values
hiY = poly.getEnvelopeInternal().getMaxY();
loY = poly.getEnvelopeInternal().getMinY();
centreY = avg(loY, hiY);
}
public double getScanLineY() {
process(poly.getExteriorRing());
for (int i = 0; i < poly.getNumInteriorRing(); i++) {
process(poly.getInteriorRingN(i));
}
double scanLineY = avg(hiY, loY);
return scanLineY;
}
private void process(LineString line) {
CoordinateSequence seq = line.getCoordinateSequence();
for (int i = 0; i < seq.size(); i++) {
double y = seq.getY(i);
updateInterval(y);
}
}
private void updateInterval(double y) {
if ( y <= centreY ) {
if ( y > loY )
loY = y;
} else if ( y > centreY ) {
if ( y < hiY ) {
hiY = y;
}
}
}
}
}