Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.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.awt.Dimension;
17
18 import java.io.Serializable;
19
20 import java.util.Arrays;
21 import java.util.List;
22
23 import org.apache.log4j.Logger;
24
25 /**
26 * Interpolates parameter values along a given line string from surface
27 * to the sea ground to generate "Profilschnitte".
28 *
29 * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
30 */
31 public class Interpolation3D
32 implements Serializable
33 {
34 private static Logger log = Logger.getLogger(Interpolation3D.class);
35
36 /**
37 * The default width of the interpolation: {@value}
38 */
39 public static final int DEFAULT_WIDTH = 1024;
40
41 /**
42 * The default height of the interpolation: {@value}
43 */
44 public static final int DEFAULT_HEIGHT = 768;
45
46 /**
47 * Epsilon for numerical stability.
48 */
49 public static final double EPS = 1e-6d;
50
51 /**
52 * The width of the interpolation.
53 */
54 protected int width;
55
56 /**
57 * The height of the interpolation.
58 */
59 protected int height;
60
61 /**
62 * The cell width of the interpolation in world units.
63 */
64 protected double cellWidth;
65
66 /**
67 * The cell height of the interpolation in world units.
68 */
69 protected double cellHeight;
70
71 /**
72 * The coordinates of the interpolation.
73 */
74 protected Coordinate[] coordinates;
75
76 /**
77 * The generated raster.
78 */
79 protected double [] raster;
80
81 /**
82 * The sea ground depth along the line string.
83 */
84 protected double [] depths;
85
86 /**
87 * Default constructor.
88 */
89 public Interpolation3D() {
90 this(DEFAULT_WIDTH, DEFAULT_HEIGHT);
91 }
92
93 /**
94 * Constructor to create a Interpolation3D with a given size.
95 * @param size The size of the interpolation.
96 */
97 public Interpolation3D(Dimension size) {
98 this(size.width, size.height);
99 }
100
101 /**
102 * Constructor to create a Interpolation3D with a given size.
103 * @param width The width of the interpolation.
104 * @param height the height of the interpolation.
105 */
106 public Interpolation3D(int width, int height) {
107 this.width = width;
108 this.height = height;
109 }
110
111 /**
112 * Returns the raster height of the interpolation.
113 * @return The raster height of the interpolation.
114 */
115 public int getHeight() {
116 return height;
117 }
118
119 /**
120 * Returns the raster width of the interpolation.
121 * @return The raster width of the interpolation.
122 */
123 public int getWidth() {
124 return width;
125 }
126
127 /**
128 * Returns the cell width of the interpolation in world units.
129 * @return The cell width of the interpolation in world units.
130 */
131 public double getCellWidth() {
132 return cellWidth;
133 }
134
135 /**
136 * Returns the cell height of the interpolation in world units.
137 * @return The cell height of the interpolation in world units.
138 */
139 public double getCellHeight() {
140 return cellHeight;
141 }
142
143 /**
144 * Returns the coordinates used for the interpolation.
145 * @return the coordinates used for the interpolation.
146 */
147 public Coordinate[] getCoordinates() {
148 return coordinates;
149 }
150
151 /**
152 * Returns the generated raster.
153 * @return The generated raster.
154 */
155 public double [] getRaster() {
156 return raster;
157 }
158
159 /**
160 * Returns the depths along the line string path.
161 * @return The depth along the line string path.
162 */
163 public double [] getDepths() {
164 return depths;
165 }
166
167 /**
168 * Returns the deepest depth along the line string path.
169 * @return The deepest depth along the line string path.
170 */
171 public double getMaxDepth() {
172 double maxDepth = Double.MAX_VALUE;
173 for (int i = depths!=null?depths.length-1:0; i >= 0; --i) {
174 double d = depths[i];
175 if (!Double.isNaN(d) && d < maxDepth) {
176 maxDepth = d;
177 }
178 }
179 return maxDepth;
180 }
181
182 /**
183 * Interpolates parameters along a given line string path from the surface
184 * to the sea ground.
185 * @param path The line string path.
186 * @param points The sample points.
187 * @param from Start point in scalar terms of the line string.
188 * @param to End point in scalar terms of the line string.
189 * @param metrics The used metric.
190 * @param xyDepth The callback to query the depth at a given point.
191 * @return true if the interpolation succeeds, else false.
192 */
193 public boolean interpolate(
194 List<? extends Coordinate> path,
195 List<? extends XYColumn> points,
196 double from,
197 double to,
198 Metrics metrics,
199 XYDepth xyDepth
200 ) {
201 boolean debug = log.isDebugEnabled();
202
203 int N = path.size();
204 int M = points.size();
205
206 if (debug) {
207 log.debug("Size of path: " + N);
208 log.debug("Size of points: " + M);
209 }
210
211 if (M < 1 || N < 2) { // nothing to do
212 return false;
213 }
214
215 List<GridCell> cells = GridCell.pointsToGridCells(
216 points, Interpolation2D.relevantArea(path, points));
217
218 if (cells.isEmpty()) {
219 log.warn("no cells found");
220 return false;
221 }
222
223 // put into spatial index to speed up finding neighbors.
224 STRtree spatialIndex = new STRtree();
225
226 for (GridCell cell: cells) {
227 spatialIndex.insert(cell.getEnvelope(), cell);
228 }
229
230 LinearToMap linearToMap = new LinearToMap(
231 path, from, to, metrics);
232
233 double [] depths = new double[width];
234 Arrays.fill(depths, Double.NaN);
235
236 Coordinate[] coordinates = new Coordinate[width];
237
238 double cellWidth = (to - from)/width;
239
240 double maxDepth = Double.MAX_VALUE;
241
242 int i = 0;
243 Coordinate center = new Coordinate();
244 for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
245 if (linearToMap.locate(p, center)) {
246 coordinates[i] = (Coordinate) center.clone();
247 double depth = xyDepth.depth(center);
248 depths[i] = depth;
249 if (depth < maxDepth) {
250 maxDepth = depth;
251 }
252 }
253 }
254
255 if (maxDepth == Double.MAX_VALUE) {
256 log.warn("no depth found -> no interpolation");
257 return false;
258 }
259
260 double cellHeight = Math.abs(maxDepth)/height;
261
262 if (debug) {
263 log.debug("max depth found: " + maxDepth);
264 log.debug("cell size: " + cellWidth + " x " + cellHeight);
265 }
266
267 double [] raster = new double[width*height];
268 Arrays.fill(raster, Double.NaN);
269
270 Envelope queryBuffer = new Envelope();
271 GridCell.CellFinder finder = new GridCell.CellFinder();
272
273 i = 0;
274 for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
275 double depth = depths[i];
276 if (Double.isNaN(depth) || depth >= 0d) {
277 continue;
278 }
279 linearToMap.locate(p, center);
280
281 queryBuffer.init(
282 center.x - EPS, center.x + EPS,
283 center.y - EPS, center.y + EPS);
284
285 finder.prepare(center);
286 spatialIndex.query(queryBuffer, finder);
287
288 GridCell found = finder.found;
289
290 if (found == null) {
291 continue;
292 }
293
294 XYColumn n0 = (XYColumn)found.p1;
295 XYColumn n1 = (XYColumn)found.p2;
296 XYColumn n2 = (XYColumn)found.p3;
297 XYColumn n3 = (XYColumn)found.p4;
298
299 if (n0.prepare(xyDepth)
300 && n1.prepare(xyDepth)
301 && n2.prepare(xyDepth)
302 && n3.prepare(xyDepth)
303 ) {
304 double y1 = Interpolation2D.interpolate(
305 n0.x, n0.y,
306 n1.x, n1.y,
307 center.x);
308 double y2 = Interpolation2D.interpolate(
309 n2.x, n2.y,
310 n3.x, n3.y,
311 center.x);
312 double z = -cellHeight*0.5;
313 for (int j = i;
314 j < raster.length && z >= depth;
315 z -= cellHeight, j += width) {
316
317 double v0 = n0.value(z);
318 double v1 = n1.value(z);
319 double v2 = n2.value(z);
320 double v3 = n3.value(z);
321
322 double z1 = Interpolation2D.interpolate(
323 n0.x, v0,
324 n1.x, v1,
325 center.x);
326 double z2 = Interpolation2D.interpolate(
327 n2.x, v2,
328 n3.x, v3,
329 center.x);
330 double value = Interpolation2D.interpolate(
331 y1, z1,
332 y2, z2,
333 center.y);
334 raster[j] = value;
335 }
336 // XXX: Adjusted depth to create no gap
337 // between last value and ground.
338 depths[i] = z+0.5d*cellHeight;
339 } // down the x/y column
340 } // all along the track
341
342 this.coordinates = coordinates;
343 this.depths = depths;
344 this.raster = raster;
345 this.cellWidth = cellWidth;
346 this.cellHeight = cellHeight;
347
348 return true;
349 }
350 }
351 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org