Class DiscreteHausdorffDistance

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
13	package org.locationtech.jts.algorithm.distance;
14
15	import org.locationtech.jts.geom.Coordinate;
16	import org.locationtech.jts.geom.CoordinateFilter;
17	import org.locationtech.jts.geom.CoordinateSequence;
18	import org.locationtech.jts.geom.CoordinateSequenceFilter;
19	import org.locationtech.jts.geom.Geometry;
20
21	/**
22	* An algorithm for computing a distance metric
23	* which is an approximation to the Hausdorff Distance
24	* based on a discretization of the input {@link Geometry}.
25	* The algorithm computes the Hausdorff distance restricted to discrete points
26	* for one of the geometries.
27	* The points can be either the vertices of the geometries (the default),
28	* or the geometries with line segments densified by a given fraction.
29	* Also determines two points of the Geometries which are separated by the computed distance.
30	* <p>
31	* This algorithm is an approximation to the standard Hausdorff distance.
32	* Specifically,
33	* <pre>
34	* for all geometries a, b: DHD(a, b) <= HD(a, b)
35	* </pre>
36	* The approximation can be made as close as needed by densifying the input geometries.
37	* In the limit, this value will approach the true Hausdorff distance:
38	* <pre>
39	* DHD(A, B, densifyFactor) -> HD(A, B) as densifyFactor -> 0.0
40	* </pre>
41	* The default approximation is exact or close enough for a large subset of useful cases.
42	* Examples of these are:
43	* <ul>
44	* <li>computing distance between Linestrings that are roughly parallel to each other,
45	* and roughly equal in length. This occurs in matching linear networks.
46	* <li>Testing similarity of geometries.
47	* </ul>
48	* An example where the default approximation is not close is:
49	* <pre>
50	* A = LINESTRING (0 0, 100 0, 10 100, 10 100)
51	* B = LINESTRING (0 100, 0 10, 80 10)
52	*
53	* DHD(A, B) = 22.360679774997898
54	* HD(A, B) ~= 47.8
55	* </pre>
56	*/
57	public class DiscreteHausdorffDistance
58	{
59	public static double distance(Geometry g0, Geometry g1)
60	{
61	DiscreteHausdorffDistance dist = new DiscreteHausdorffDistance(g0, g1);
62	return dist.distance();
63	}
64
65	public static double distance(Geometry g0, Geometry g1, double densifyFrac)
66	{
67	DiscreteHausdorffDistance dist = new DiscreteHausdorffDistance(g0, g1);
68	dist.setDensifyFraction(densifyFrac);
69	return dist.distance();
70	}
71
72	private Geometry g0;
73	private Geometry g1;
74	private PointPairDistance ptDist = new PointPairDistance();
75
76	/**
77	* Value of 0.0 indicates that no densification should take place
78	*/
79	private double densifyFrac = 0.0;
80
81	public DiscreteHausdorffDistance(Geometry g0, Geometry g1)
82	{
83	this.g0 = g0;
84	this.g1 = g1;
85	}
86
87	/**
88	* Sets the fraction by which to densify each segment.
89	* Each segment will be (virtually) split into a number of equal-length
90	* subsegments, whose fraction of the total length is closest
91	* to the given fraction.
92	*
93	* @param densifyFrac
94	*/
95	public void setDensifyFraction(double densifyFrac)
96	{
97	if (densifyFrac > 1.0
98	\|\| densifyFrac <= 0.0)
99	throw new IllegalArgumentException("Fraction is not in range (0.0 - 1.0]");
100
101	this.densifyFrac = densifyFrac;
102	}
103
104	public double distance()
105	{
106	compute(g0, g1);
107	return ptDist.getDistance();
108	}
109
110	public double orientedDistance()
111	{
112	computeOrientedDistance(g0, g1, ptDist);
113	return ptDist.getDistance();
114	}
115
116	public Coordinate[] getCoordinates() { return ptDist.getCoordinates(); }
117
118	private void compute(Geometry g0, Geometry g1)
119	{
120	computeOrientedDistance(g0, g1, ptDist);
121	computeOrientedDistance(g1, g0, ptDist);
122	}
123
124	private void computeOrientedDistance(Geometry discreteGeom, Geometry geom, PointPairDistance ptDist)
125	{
126	MaxPointDistanceFilter distFilter = new MaxPointDistanceFilter(geom);
127	discreteGeom.apply(distFilter);
128	ptDist.setMaximum(distFilter.getMaxPointDistance());
129
130	if (densifyFrac > 0) {
131	MaxDensifiedByFractionDistanceFilter fracFilter = new MaxDensifiedByFractionDistanceFilter(geom, densifyFrac);
132	discreteGeom.apply(fracFilter);
133	ptDist.setMaximum(fracFilter.getMaxPointDistance());
134
135	}
136	}
137
138	public static class MaxPointDistanceFilter
139	implements CoordinateFilter
140	{
141	private PointPairDistance maxPtDist = new PointPairDistance();
142	private PointPairDistance minPtDist = new PointPairDistance();
143	private DistanceToPoint euclideanDist = new DistanceToPoint();
144	private Geometry geom;
145
146	public MaxPointDistanceFilter(Geometry geom)
147	{
148	this.geom = geom;
149	}
150
151	public void filter(Coordinate pt)
152	{
153	minPtDist.initialize();
154	DistanceToPoint.computeDistance(geom, pt, minPtDist);
155	maxPtDist.setMaximum(minPtDist);
156	}
157
158	public PointPairDistance getMaxPointDistance() { return maxPtDist; }
159	}
160
161	public static class MaxDensifiedByFractionDistanceFilter
162	implements CoordinateSequenceFilter
163	{
164	private PointPairDistance maxPtDist = new PointPairDistance();
165	private PointPairDistance minPtDist = new PointPairDistance();
166	private Geometry geom;
167	private int numSubSegs = 0;
168
169	public MaxDensifiedByFractionDistanceFilter(Geometry geom, double fraction) {
170	this.geom = geom;
171	numSubSegs = (int) Math.rint(1.0/fraction);
172	}
173
174	public void filter(CoordinateSequence seq, int index)
175	{
176	/**
177	* This logic also handles skipping Point geometries
178	*/
179	if (index == 0)
180	return;
181
182	Coordinate p0 = seq.getCoordinate(index - 1);
183	Coordinate p1 = seq.getCoordinate(index);
184
185	double delx = (p1.x - p0.x)/numSubSegs;
186	double dely = (p1.y - p0.y)/numSubSegs;
187
188	for (int i = 0; i < numSubSegs; i++) {
189	double x = p0.x + i*delx;
190	double y = p0.y + i*dely;
191	Coordinate pt = new Coordinate(x, y);
192	minPtDist.initialize();
193	DistanceToPoint.computeDistance(geom, pt, minPtDist);
194	maxPtDist.setMaximum(minPtDist);
195	}
196
197
198	}
199
200	public boolean isGeometryChanged() { return false; }
201
202	public boolean isDone() { return false; }
203
204	public PointPairDistance getMaxPointDistance() {
205	return maxPtDist;
206	}
207	}
208
209	}
210