Mercurial > dive4elements > gnv-client

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