Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java @ 657:af3f56758f59

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged gnv-artifacts/0.5
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:13:53 +0200
parents 44415ae01ddb
children 9a828e5a2390
comparison
equal deleted inserted replaced
590:5f5f273c8566 657:af3f56758f59
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 * @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
19 */
20 public class Interpolation3D
21 implements Serializable
22 {
23 private static Logger log = Logger.getLogger(Interpolation3D.class);
24
25 public static final int DEFAULT_WIDTH = 1024;
26 public static final int DEFAULT_HEIGHT = 768;
27
28 public static final double EPS = 1e-6d;
29
30 protected int width;
31 protected int height;
32
33 protected double cellWidth;
34 protected double cellHeight;
35
36 protected double [] raster;
37 protected double [] depths;
38
39 public Interpolation3D() {
40 this(DEFAULT_WIDTH, DEFAULT_HEIGHT);
41 }
42
43 public Interpolation3D(Dimension size) {
44 this(size.width, size.height);
45 }
46
47 public Interpolation3D(int width, int height) {
48 this.width = width;
49 this.height = height;
50 }
51
52 public int getHeight() {
53 return height;
54 }
55
56 public int getWidth() {
57 return width;
58 }
59
60 public double getCellWidth() {
61 return cellWidth;
62 }
63
64 public double getCellHeight() {
65 return cellHeight;
66 }
67
68 public double [] getRaster() {
69 return raster;
70 }
71
72 public double [] getDepths() {
73 return depths;
74 }
75
76 public double getMaxDepth() {
77 double maxDepth = Double.MAX_VALUE;
78 for (int i = depths!=null?depths.length-1:0; i >= 0; --i) {
79 double d = depths[i];
80 if (!Double.isNaN(d) && d < maxDepth) {
81 maxDepth = d;
82 }
83 }
84 return maxDepth;
85 }
86
87 public boolean interpolate(
88 List<? extends Coordinate> path,
89 List<? extends XYColumn> points,
90 double from,
91 double to,
92 Metrics metrics,
93 XYDepth xyDepth
94 ) {
95 boolean debug = log.isDebugEnabled();
96
97 int N = path.size();
98 int M = points.size();
99
100 if (debug) {
101 log.debug("Size of path: " + N);
102 log.debug("Size of points: " + M);
103 }
104
105 if (M < 1 || N < 2) { // nothing to do
106 return false;
107 }
108
109 List<GridCell> cells = GridCell.pointsToGridCells(
110 points, Interpolation2D.relevantArea(path, points));
111
112 if (cells.isEmpty()) {
113 log.warn("no cells found");
114 return false;
115 }
116
117 // put into spatial index to speed up finding neighbors.
118 STRtree spatialIndex = new STRtree();
119
120 for (GridCell cell: cells) {
121 spatialIndex.insert(cell.getEnvelope(), cell);
122 }
123
124 LinearToMap linearToMap = new LinearToMap(
125 path, from, to, metrics);
126
127 double [] depths = new double[width];
128 Arrays.fill(depths, Double.NaN);
129
130 double cellWidth = (to - from)/width;
131
132 double maxDepth = Double.MAX_VALUE;
133
134 int i = 0;
135 Coordinate center = new Coordinate();
136 for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
137 if (linearToMap.locate(p, center)) {
138 double depth = xyDepth.depth(center);
139 depths[i] = depth;
140 if (depth < maxDepth) {
141 maxDepth = depth;
142 }
143 }
144 }
145
146 if (maxDepth == Double.MAX_VALUE) {
147 log.warn("no depth found -> no interpolation");
148 return false;
149 }
150
151 double cellHeight = Math.abs(maxDepth)/height;
152
153 if (debug) {
154 log.debug("max depth found: " + maxDepth);
155 log.debug("cell size: " + cellWidth + " x " + cellHeight);
156 }
157
158 double [] raster = new double[width*height];
159 Arrays.fill(raster, Double.NaN);
160
161 Envelope queryBuffer = new Envelope();
162 GridCell.CellFinder finder = new GridCell.CellFinder();
163
164 i = 0;
165 for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
166 double depth = depths[i];
167 if (Double.isNaN(depth) || depth >= 0d) {
168 continue;
169 }
170 linearToMap.locate(p, center);
171
172 queryBuffer.init(
173 center.x - EPS, center.x + EPS,
174 center.y - EPS, center.y + EPS);
175
176 finder.prepare(center);
177 spatialIndex.query(queryBuffer, finder);
178
179 GridCell found = finder.found;
180
181 if (found == null) {
182 continue;
183 }
184
185 XYColumn n0 = (XYColumn)found.p1;
186 XYColumn n1 = (XYColumn)found.p2;
187 XYColumn n2 = (XYColumn)found.p3;
188 XYColumn n3 = (XYColumn)found.p4;
189
190 if (n0.prepare(xyDepth)
191 && n1.prepare(xyDepth)
192 && n2.prepare(xyDepth)
193 && n3.prepare(xyDepth)
194 ) {
195 double y1 = Interpolation2D.interpolate(
196 n0.x, n0.y,
197 n1.x, n1.y,
198 center.x);
199 double y2 = Interpolation2D.interpolate(
200 n2.x, n2.y,
201 n3.x, n3.y,
202 center.x);
203 double z = -cellHeight*0.5;
204 for (int j = i;
205 j < raster.length && z >= depth;
206 z -= cellHeight, j += width) {
207
208 double v0 = n0.value(z);
209 double v1 = n1.value(z);
210 double v2 = n2.value(z);
211 double v3 = n3.value(z);
212
213 double z1 = Interpolation2D.interpolate(
214 n0.x, v0,
215 n1.x, v1,
216 center.x);
217 double z2 = Interpolation2D.interpolate(
218 n2.x, v2,
219 n3.x, v3,
220 center.x);
221 double value = Interpolation2D.interpolate(
222 y1, z1,
223 y2, z2,
224 center.y);
225 raster[j] = value;
226 }
227 // XXX: Adjusted depth to create no gap
228 // between last value and ground.
229 depths[i] = z+0.5d*cellHeight;
230 } // down the x/y column
231 } // all along the track
232
233 this.depths = depths;
234 this.raster = raster;
235 this.cellWidth = cellWidth;
236 this.cellHeight = cellHeight;
237
238 return true;
239 }
240 }
241 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:
http://dive4elements.wald.intevation.org