Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java @ 875:5e9efdda6894

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged gnv-artifacts/1.0
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:13:56 +0200
parents 05bf8534a35a
children f953c9a559d8
comparison
equal deleted inserted replaced
722:bb3ffe7d719e 875:5e9efdda6894
1 package de.intevation.gnv.math;
2
3 import com.vividsolutions.jts.geom.Coordinate;
4 import com.vividsolutions.jts.geom.Envelope;
5
6 import com.vividsolutions.jts.index.strtree.STRtree;
7
8 import java.util.List;
9
10 import org.apache.log4j.Logger;
11
12 /**
13 * Interpolates along a given line string. This used to generate
14 * "horizontale Schnittprofile".
15 *
16 * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
17 */
18 public final class Interpolation2D
19 {
20 private static Logger log = Logger.getLogger(Interpolation2D.class);
21
22 /**
23 * If the number of input points are over this threshold
24 * some culling strategies are applied to reduce the amount of points.
25 */
26 public static final int CULL_POINT_THRESHOLD = Integer.getInteger(
27 "gnv.interpolation2d.cull.point.threshold", 1000);
28
29 /**
30 * Numerical stability epsilon.
31 */
32 public static final double EPS = 1e-6d;
33
34 /**
35 * Callback to send the interpolated values back to the
36 * caller without building large temporary strcutures.
37 */
38 public interface Consumer {
39 /**
40 * Sends an interpolated point back to the called.
41 * The interpolated parameter is stored in the z component.
42 * @param point The interpolated point.
43 * @param success true if interpolation at the point
44 * succeed, else false.
45 */
46 void interpolated(Coordinate point, boolean success);
47 } // interface Consumer
48
49 private Interpolation2D() {
50 }
51
52 /**
53 * Calculates the relevant area for a given line string.
54 * @param path The line string.
55 * @param points The sample points.
56 * @return The relevant area.
57 */
58 public static final Envelope relevantArea(
59 List<? extends Coordinate> path,
60 List<? extends Coordinate> points
61 ) {
62 return relevantArea(path, points, CULL_POINT_THRESHOLD);
63 }
64
65 /**
66 * Calculates the relevant area for a given bounding box.
67 * @param pathBBox The bounding box.
68 * @param points The sample points.
69 * @return The relevant area.
70 */
71 public static final Envelope relevantArea(
72 Envelope pathBBox,
73 List<? extends Coordinate> points
74 ) {
75 return relevantArea(pathBBox, points, CULL_POINT_THRESHOLD);
76 }
77
78 /**
79 * Calculates the relevant area for a given bounding box.
80 * @param pathBBox The bounding box.
81 * @param points The sample points.
82 * @param threshold If the number of sample points
83 * are below this threshold no relevant area is calculated.
84 * @return The relevant area.
85 */
86 public static final Envelope relevantArea(
87 Envelope pathBBox,
88 List<? extends Coordinate> points,
89 int threshold
90 ) {
91 return points.size() < threshold
92 ? null
93 : relevantArea(
94 pathBBox,
95 pointsBoundingBox(points));
96 }
97
98 /**
99 * Calculates the relevant area for a given line string.
100 * @param path The line string.
101 * @param points The sample points.
102 * @param threshold If the number of sample points
103 * are below this threshold no relevant area is calculated.
104 * @return The relevant area.
105 */
106 public static final Envelope relevantArea(
107 List<? extends Coordinate> path,
108 List<? extends Coordinate> points,
109 int threshold
110 ) {
111 return points.size() < threshold || path.isEmpty()
112 ? null
113 : relevantArea(
114 pointsBoundingBox(path),
115 pointsBoundingBox(points));
116 }
117
118 /**
119 * Heuristic function to define a 'relevant area'.
120 * @param pathBBox The bounding box of the line line string.
121 * @param pointsBBox The bounding box of the sample points.
122 * @return The relevant area.
123 */
124 public static final Envelope relevantArea(
125 Envelope pathBBox,
126 Envelope pointsBBox
127 ) {
128 double pathArea = pathBBox.getWidth()*pathBBox.getHeight();
129 double pointsArea = pointsBBox.getWidth()*pointsBBox.getHeight();
130
131 if (pathArea > 0.8d*pointsArea) { return null; }
132
133 double nArea = 1.44d * pathArea;
134 if (nArea < 0.1d*pointsArea) nArea = 0.1d*pointsArea;
135 double w = pathBBox.getWidth();
136 double h = pathBBox.getHeight();
137 double [] d = solveQuadratic(1d, w+h, pathArea - nArea);
138
139 if (d == null) { return null; }
140
141 double extra = pos(d);
142
143 pathBBox.expandBy(extra);
144
145 return pathBBox;
146 }
147
148 /**
149 * Solves quadratic equation a*x*x + b*x + c = 0.
150 * @param a a-coefficent.
151 * @param b b-coefficent
152 * @param c c-coefficent.
153 * @return the solution of the equation, or null if not solvable.
154 */
155 public static final double [] solveQuadratic(
156 double a, double b, double c
157 ) {
158 double d = b*b - 4d*a*c;
159 if (d < 0d) { return null; }
160
161 d = Math.sqrt(d);
162 a = 1d/(2d*a);
163 b = -b;
164
165 return new double [] { a*(b + d), a*(b - d) };
166 }
167
168 /**
169 * Return the element of a two element array which
170 * is greater or equal zero.
171 * @param x The two values.
172 * @return The value which is greater or equal zero.
173 */
174 public static final double pos(double [] x) {
175 return x[0] >= 0 ? x[0] : x[1];
176 }
177
178
179 /**
180 * Calculates the bounding box of a given line string.
181 * @param path The line string.
182 * @return The bounding box.
183 */
184 public static Envelope pointsBoundingBox(
185 List<? extends Coordinate> path
186 ) {
187 int N = path.size();
188 Envelope area = new Envelope(path.get(N-1));
189
190 for (int i = N-2; i >= 0; --i) {
191 area.expandToInclude(path.get(i));
192 }
193
194 return area;
195 }
196
197 /**
198 * Interpolates linearly a number of coordinates and parameter values along
199 * a given line string path. The results are issued to a consumer.
200 * @param path The line string path.
201 * @param points The sample points.
202 * @param from Start point as a scalar value linear
203 * referenced on the line string.
204 * @param to End point of as a scalar value linear
205 * referenced on the line string.
206 * @param steps Number of points to be interpolated.
207 * @param metrics The used metric.
208 * @param consumer The callback to retrieve the result points.
209 */
210 public static void interpolate(
211 List<? extends Coordinate> path,
212 List<? extends Point2d> points,
213 double from,
214 double to,
215 int steps,
216 Metrics metrics,
217 Consumer consumer
218 ) {
219 boolean debug = log.isDebugEnabled();
220
221 int N = path.size();
222 int M = points.size();
223
224 if (debug) {
225 log.debug("Size of path: " + N);
226 log.debug("Size of points: " + M);
227 }
228
229 if (M < 1 || N < 2) { // nothing to do
230 return;
231 }
232
233 List<GridCell> cells = GridCell.pointsToGridCells(
234 points, relevantArea(path, points));
235
236 if (cells.isEmpty()) {
237 log.warn("no cells found");
238 return;
239 }
240
241 // put into spatial index to speed up finding neighbors.
242 STRtree spatialIndex = new STRtree();
243
244 for (GridCell cell: cells) {
245 spatialIndex.insert(cell.getEnvelope(), cell);
246 }
247
248 LinearToMap linearToMap = new LinearToMap(
249 path, from, to, metrics);
250
251 double dP = (to - from)/steps;
252
253 Coordinate center = new Coordinate();
254 Envelope queryBuffer = new Envelope();
255 GridCell.CellFinder finder = new GridCell.CellFinder();
256
257 int missedInterpolations = 0;
258 int interpolations = 0;
259
260 for (double p = from; p <= to; p += dP) {
261 if (!linearToMap.locate(p, center)) {
262 continue;
263 }
264 queryBuffer.init(
265 center.x - EPS, center.x + EPS,
266 center.y - EPS, center.y + EPS);
267
268 finder.prepare(center);
269 spatialIndex.query(queryBuffer, finder);
270
271 GridCell found = finder.found;
272
273 if (found == null) {
274 consumer.interpolated(center, false);
275 ++missedInterpolations;
276 continue;
277 }
278
279 Point2d n0 = found.p1;
280 Point2d n1 = found.p2;
281 Point2d n2 = found.p3;
282 Point2d n3 = found.p4;
283
284 double z1 = interpolate(
285 n0.x, n0.z,
286 n1.x, n1.z,
287 center.x);
288 double z2 = interpolate(
289 n2.x, n2.z,
290 n3.x, n3.z,
291 center.x);
292 double y1 = interpolate(
293 n0.x, n0.y,
294 n1.x, n1.y,
295 center.x);
296 double y2 = interpolate(
297 n2.x, n2.y,
298 n3.x, n3.y,
299 center.x);
300 center.z = interpolate(
301 y1, z1,
302 y2, z2,
303 center.y);
304 consumer.interpolated(center, true);
305 ++interpolations;
306 }
307
308 if (debug) {
309 log.debug("interpolations: " +
310 interpolations + " / " + missedInterpolations);
311 }
312 }
313
314 /**
315 * Linear interpolate a value between (x1, y1) and (x2, y2) at
316 * a given x-value.
317 * @param x1 x component of first point.
318 * @param y1 y component of first point.
319 * @param x2 x component of second point.
320 * @param y2 y component of second point.
321 * @param x The x value.
322 * @return The intepolated result.
323 */
324 public static final double interpolate(
325 double x1, double y1,
326 double x2, double y2,
327 double x
328 ) {
329 if (x2 == x1) {
330 return (y1 + y2)*0.5d;
331 }
332 double m = (y2-y1)/(x2-x1);
333 double b = y1 - m*x1;
334 return m*x + b;
335 }
336 }
337 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org