Skip to content

Commit 44e7b41

Browse files
fmeurissefmeurisse
committed
S2PolylineSimplifier helper class.
1 parent 5bd8781 commit 44e7b41

File tree

3 files changed

+443
-0
lines changed

3 files changed

+443
-0
lines changed

‎s2-geometry.iml

+23
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
<?xml version="1.0" encoding="UTF-8"?>
2+
<module org.jetbrains.idea.maven.project.MavenProjectsManager.isMavenModule="true" type="JAVA_MODULE" version="4">
3+
<component name="NewModuleRootManager" LANGUAGE_LEVEL="JDK_1_8">
4+
<output url="file://$MODULE_DIR$/target/classes" />
5+
<output-test url="file://$MODULE_DIR$/target/test-classes" />
6+
<content url="file://$MODULE_DIR$">
7+
<sourceFolder url="file://$MODULE_DIR$/src" isTestSource="false" />
8+
<sourceFolder url="file://$MODULE_DIR$/tests" isTestSource="true" />
9+
<excludeFolder url="file://$MODULE_DIR$/target" />
10+
</content>
11+
<orderEntry type="inheritedJdk" />
12+
<orderEntry type="sourceFolder" forTests="false" />
13+
<orderEntry type="library" name="Maven: com.google.guava:guava:29.0-jre" level="project" />
14+
<orderEntry type="library" name="Maven: com.google.guava:failureaccess:1.0.1" level="project" />
15+
<orderEntry type="library" name="Maven: com.google.guava:listenablefuture:9999.0-empty-to-avoid-conflict-with-guava" level="project" />
16+
<orderEntry type="library" name="Maven: org.checkerframework:checker-qual:2.11.1" level="project" />
17+
<orderEntry type="library" name="Maven: com.google.errorprone:error_prone_annotations:2.3.4" level="project" />
18+
<orderEntry type="library" name="Maven: com.google.j2objc:j2objc-annotations:1.3" level="project" />
19+
<orderEntry type="library" name="Maven: com.google.code.findbugs:jsr305:3.0.1" level="project" />
20+
<orderEntry type="library" scope="TEST" name="Maven: junit:junit:4.13.1" level="project" />
21+
<orderEntry type="library" scope="TEST" name="Maven: org.hamcrest:hamcrest-core:1.3" level="project" />
22+
</component>
23+
</module>

‎src/com/google/common/geometry/S2PolylineSimplifier.java

+258
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,258 @@
1+
package com.google.common.geometry;
2+
3+
import static com.google.common.geometry.S2.M_PI;
4+
5+
/**
6+
* This is a helper class for simplifying polylines. It allows you to compute a maximal edge that intersects a sequence
7+
* of discs, and that optionally avoids a different sequence of discs. The results are conservative in that the edge
8+
* is guaranteed to intersect or avoid the specified discs using exact arithmetic (see S2Predicates).
9+
*
10+
* Note that S2Builder can also simplify polylines and supports more features (e.g., snapping to S2CellId centers),
11+
* so it is only recommended to use this class if S2Builder does not meet your needs.
12+
*
13+
* Here is a simple example showing how to simplify a polyline into a sequence of edges that stay within "max_error"
14+
* of the original edges:
15+
*
16+
* <pre>
17+
* List<S2Point> v = ...
18+
* List<S2Point> simplifiedPolyline = new ArrayList<S2Point>();
19+
* S2PolylineSimplifier simplifier = new S2PolylineSimplifier();
20+
* simplifier.init(v.get(0));
21+
* simplifiedPolyline.add(v.get(0));
22+
* for (int i = 1 ; i < v.size() ; ++i) {
23+
* if (simplifier.extend(v[i])) {
24+
* simplifiedPolyline.add(v[i]);
25+
* simplifier.init(v[i])
26+
* simplifier.targetDisc(v[i], max_error)
27+
* }
28+
* }
29+
* </pre>
30+
*
31+
* Note that the points targeted by TargetDisc do not need to be the same as the candidate endpoints passed to extend.
32+
* So for example, you could target the original vertices of a polyline, but only consider endpoints that are snapped to
33+
* E7 coordinates or S2CellId centers.
34+
*
35+
* Please be aware that this class works by maintaining a range of acceptable angles (bearings) from the start vertex
36+
* to the hypothetical destination vertex. It does not keep track of distances to any of the discs to be targeted or
37+
* avoided. Therefore to use this class correctly, constraints should be added in increasing order of distance.
38+
* (The actual requirement is slightly weaker than this, which is why it is not enforced, but basically you should only
39+
* call targetDisc() and avoidDisc() with arguments that you want to constrain the immediately following call to
40+
* extend().)
41+
*
42+
*/
43+
class S2PolylineSimplifier {
44+
45+
private S2Point src;
46+
private S2Point xDir = new S2Point();
47+
private S2Point yDir = new S2Point();
48+
private S1Interval window;
49+
50+
/**
51+
* Starts a new simplified edge at "src".
52+
*
53+
* @param src The source point.
54+
*/
55+
public void init(S2Point src) {
56+
this.src = src;
57+
this.window = S1Interval.full();
58+
59+
// Precompute basis vectors for the tangent space at "src". This is similar
60+
// to GetFrame() except that we don't normalize the vectors. As it turns
61+
// out, the two basis vectors below have the same magnitude (up to the
62+
// length error in S2Point::Normalize).
63+
64+
int i, j, k;
65+
// Find the index of the component whose magnitude is smallest.
66+
S2Point tmp = src.fabs();
67+
if (tmp.x < tmp.y) {
68+
if (tmp.x < tmp.z) i = 0;
69+
else i = 2;
70+
} else if (tmp.y < tmp.z) i = 1;
71+
else i = 2;
72+
73+
// We define the "y" basis vector as the cross product of "src" and the
74+
// basis vector for axis "i". Let "j" and "k" be the indices of the other
75+
// two components in cyclic order.
76+
if (i == 2) j = 0;
77+
else j = (i + 1);
78+
if (i == 0) k = 2;
79+
else k = (i - 1);
80+
81+
double[] yCoords = new double[3];
82+
yCoords[i] = 0.0;
83+
yCoords[j] = src.get(k);
84+
yCoords[k] = -src.get(j);
85+
yDir = new S2Point(yCoords[0], yCoords[1], yCoords[2]);
86+
87+
// Compute the cross product of "y_dir" and "src". We write out the cross
88+
// product here mainly for documentation purposes; it also happens to save a
89+
// few multiplies because unfortunately the optimizer does *not* get rid of
90+
// multiplies by zero (since these multiplies propagate NaN, for example).
91+
double[] xCoords = new double[3];
92+
xCoords[i] = src.get(j) * src.get(j) + src.get(k) * src.get(k);
93+
xCoords[j] = -src.get(j) * src.get(i);
94+
xCoords[k] = -src.get(k) * src.get(i);
95+
xDir = new S2Point(xCoords[0], xCoords[1], xCoords[2]);
96+
}
97+
98+
/** Returns the source vertex of the output edge. */
99+
private S2Point getSrc() {
100+
return src;
101+
}
102+
103+
/**
104+
* Returns true if the edge (src, dst) satisfies all of the targeting requirements so far.
105+
* Returns false if the edge would be longer than 90 degrees (such edges are not supported).
106+
*/
107+
public boolean extend(S2Point dst) {
108+
// We limit the maximum edge length to 90 degrees in order to simplify the
109+
// error bounds. (The error gets arbitrarily large as the edge length
110+
// approaches 180 degrees.)
111+
if (new S1ChordAngle(src, dst).getLength2() > S1ChordAngle.RIGHT.getLength2()) return false;
112+
113+
// Otherwise check whether this vertex is in the acceptable angle range.
114+
return window.contains(getAngle(dst));
115+
}
116+
117+
/** Requires that the output edge must pass through the given disc. */
118+
public boolean targetDisc(S2Point point, S1ChordAngle radius) {
119+
// Shrink the target interval by the maximum error from all sources. This
120+
// guarantees that the output edge will intersect the given disc.
121+
double semiwidth = getSemiwidth(point, radius, -1 /*round down*/);
122+
if (semiwidth >= M_PI) {
123+
// The target disc contains "src", so there is nothing to do.
124+
return true;
125+
}
126+
if (semiwidth < 0) {
127+
window = S1Interval.empty();
128+
return false;
129+
}
130+
// Otherwise compute the angle interval corresponding to the target disc and
131+
// intersect it with the current window.
132+
double center = getAngle(point);
133+
S1Interval target = S1Interval.fromPoint(center).expanded(semiwidth);
134+
window = window.intersection(target);
135+
return !window.isEmpty();
136+
}
137+
138+
/**
139+
* Requires that the output edge must avoid the given disc. "disc_on_left"
140+
* specifies whether the disc must be to the left or right of the edge.
141+
* (This feature allows the simplified edge to preserve the topology of the
142+
* original polyline with respect to other nearby points.)
143+
*
144+
* If your input is a polyline, you can compute "disc_on_left" as follows.
145+
* Let the polyline be ABCDE and assume that it already avoids a set of
146+
* points X_i. Suppose that you have already added ABC to the simplifer, and
147+
* now want to extend the edge chain to D. First find the X_i that are near
148+
* the edge CD, then discard the ones such that AX_i <= AC or AX_i >= AD
149+
* (since these points have either already been considered or aren't
150+
* relevant yet). Now X_i is to the left of the polyline if and only if
151+
* s2pred::OrderedCCW(A, D, X, C) (in other words, if X_i is to the left of
152+
* the angle wedge ACD).
153+
*/
154+
public boolean avoidDisc(S2Point point, S1ChordAngle radius, boolean discOnLeft) {
155+
// Expand the interval by the maximum error from all sources. This
156+
// guarantees that the final output edge will avoid the given disc.
157+
double semiwidth = getSemiwidth(point, radius, 1 /*round up*/);
158+
if (semiwidth >= M_PI) {
159+
// The avoidance disc contains "src", so it is impossible to avoid.
160+
window = S1Interval.empty();
161+
return false;
162+
}
163+
double center = getAngle(point);
164+
double opposite;
165+
if (center > 0) opposite = center - M_PI;
166+
else opposite = center + M_PI;
167+
S1Interval target;
168+
if (discOnLeft) target = new S1Interval(opposite, center);
169+
else target = new S1Interval(center, opposite);
170+
window = window.intersection(target.expanded(-semiwidth));
171+
return !window.isEmpty();
172+
}
173+
174+
private double getAngle(S2Point p) {
175+
return Math.atan2(p.dotProd(yDir), p.dotProd(xDir));
176+
}
177+
178+
private double getSemiwidth(S2Point p, S1ChordAngle r, int roundDirection) {
179+
double err = 0.5 * S2.DBL_EPSILON;
180+
181+
// Using spherical trigonometry,
182+
//
183+
// sin(semiwidth) = sin(r) / sin(a)
184+
//
185+
// where "a" is the angle between "src" and "p". Rather than measuring
186+
// these angles, instead we measure the squared chord lengths through the
187+
// interior of the sphere (i.e., Cartersian distance). Letting "r2" be the
188+
// squared chord distance corresponding to "r", and "a2" be the squared
189+
// chord distance corresponding to "a", we use the relationships
190+
//
191+
// sin^2(r) = r2 (1 - r2 / 4)
192+
// sin^2(a) = d2 (1 - d2 / 4)
193+
//
194+
// which follow from the fact that r2 = (2 * sin(r / 2)) ^ 2, etc.
195+
196+
// "a2" has a relative error up to 5 * DBL_ERR, plus an absolute error of up
197+
// to 64 * DBL_ERR^2 (because "src" and "p" may differ from unit length by
198+
// up to 4 * DBL_ERR). We can correct for the relative error later, but for
199+
// the absolute error we use "round_direction" to account for it now.
200+
double r2 = r.getLength2();
201+
double a2 = new S1ChordAngle(src, p).getLength2();
202+
a2 -= 64 * err * err * roundDirection;
203+
if (a2 <= r2) return M_PI; // The given disc contains "src".
204+
205+
double sin2_r = r2 * (1 - 0.25 * r2);
206+
double sin2_a = a2 * (1 - 0.25 * a2);
207+
double semiwidth = Math.asin(Math.sqrt(sin2_r / sin2_a));
208+
209+
// We compute bounds on the errors from all sources:
210+
//
211+
// - The call to GetSemiwidth (this call).
212+
// - The call to GetAngle that computes the center of the interval.
213+
// - The call to GetAngle in Extend that tests whether a given point
214+
// is an acceptable destination vertex.
215+
//
216+
// Summary of the errors in GetAngle:
217+
//
218+
// y_dir_ has no error.
219+
//
220+
// x_dir_ has a relative error of DBL_ERR in two components, a relative
221+
// error of 2 * DBL_ERR in the other component, plus an overall relative
222+
// length error of 4 * DBL_ERR (compared to y_dir_) because "src" is assumed
223+
// to be normalized only to within the tolerances of S2Point::Normalize().
224+
//
225+
// p.DotProd(y_dir_) has a relative error of 1.5 * DBL_ERR and an
226+
// absolute error of 1.5 * DBL_ERR * y_dir_.Norm().
227+
//
228+
// p.DotProd(x_dir_) has a relative error of 5.5 * DBL_ERR and an absolute
229+
// error of 3.5 * DBL_ERR * y_dir_.Norm() (noting that x_dir_ and y_dir_
230+
// have the same length to within a relative error of 4 * DBL_ERR).
231+
//
232+
// It's possible to show by taking derivatives that these errors can affect
233+
// the angle atan2(y, x) by up 7.093 * DBL_ERR radians. Rounding up and
234+
// including the call to atan2 gives a final error bound of 10 * DBL_ERR.
235+
//
236+
// Summary of the errors in GetSemiwidth:
237+
//
238+
// The distance a2 has a relative error of 5 * DBL_ERR plus an absolute
239+
// error of 64 * DBL_ERR^2 because the points "src" and "p" may differ from
240+
// unit length (by up to 4 * DBL_ERR). We have already accounted for the
241+
// absolute error above, leaving only the relative error.
242+
//
243+
// sin2_r has a relative error of 2 * DBL_ERR.
244+
//
245+
// sin2_a has a relative error of 12 * DBL_ERR assuming that a2 <= 2,
246+
// i.e. distance(src, p) <= 90 degrees. (The relative error gets
247+
// arbitrarily larger as this distance approaches 180 degrees.)
248+
//
249+
// semiwidth has a relative error of 17 * DBL_ERR.
250+
//
251+
// Finally, (center +/- semiwidth) has a rounding error of up to 4 * DBL_ERR
252+
// because in theory, the result magnitude may be as large as 1.5 * M_PI
253+
// which is larger than 4.0. This gives a total error of:
254+
double error = (2 * 10 + 4) * err + 17 * err * semiwidth;
255+
return semiwidth + roundDirection * error;
256+
}
257+
258+
}

0 commit comments

Comments
 (0)