Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java @ 1119:7c4f81f74c47

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