Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java @ 434:0eed5749fd63

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implemented the raster interpolation for the 'Profilschnitt'. gnv-artifacts/trunk@482 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 23 Dec 2009 15:28:40 +0000
parents
children f42ed4f10b79
comparison
equal deleted inserted replaced
433:828df3ddb758 434:0eed5749fd63
1 package de.intevation.gnv.math;
2
3 import java.util.List;
4 import java.util.Arrays;
5 import java.util.Collections;
6
7 import com.vividsolutions.jts.geom.Coordinate;
8 import com.vividsolutions.jts.geom.Envelope;
9
10 import com.vividsolutions.jts.index.quadtree.Quadtree;
11
12 import org.apache.log4j.Logger;
13
14 public class Interpolation3D
15 {
16 private static Logger log = Logger.getLogger(Interpolation3D.class);
17
18 public static final int DEFAULT_WIDTH = 1024;
19 public static final int DEFAULT_HEIGHT = 768;
20
21 protected int width;
22 protected int height;
23
24 protected double [] raster;
25 protected double [] depths;
26
27 public Interpolation3D() {
28 this(DEFAULT_WIDTH, DEFAULT_HEIGHT);
29 }
30
31 public Interpolation3D(int width, int height) {
32 this.width = width;
33 this.height = height;
34 }
35
36 public int getHeight() {
37 return height;
38 }
39
40 public int getWidth() {
41 return width;
42 }
43
44 public double [] getRaster() {
45 return raster;
46 }
47
48 public double [] getDepths() {
49 return depths;
50 }
51
52 public boolean interpolate(
53 List<? extends Coordinate> path,
54 List<? extends XYColumn> points,
55 double from,
56 double to,
57 Metrics metrics,
58 XYDepth xyDepth
59 ) {
60 boolean debug = log.isDebugEnabled();
61
62 int N = path.size();
63 int M = points.size();
64
65 if (debug) {
66 log.debug("Size of path: " + N);
67 log.debug("Size of points: " + M);
68 }
69
70 if (M < 1 || N < 2) { // nothing to do
71 return false;
72 }
73
74 double [] buffer = Interpolation2D.calculateBuffer(points);
75 double dxMax = buffer[0];
76 double dyMax = buffer[1];
77
78 if (debug) {
79 log.debug("buffer size: " + dxMax + " / " + dyMax);
80 }
81
82 // put into spatial index to speed up finding neighbors.
83 Quadtree spatialIndex = new Quadtree();
84
85 for (int i = 0; i < M; ++i) {
86 Point2d p = points.get(i);
87 spatialIndex.insert(p.envelope(), p);
88 }
89
90 LinearToMap linearToMap = new LinearToMap(
91 path, from, to, metrics);
92
93 double [] depths = new double[width];
94 Arrays.fill(depths, Double.NaN);
95
96 double cellWidth = (to - from)/width;
97
98 double maxDepth = Double.MAX_VALUE;
99
100 int i = 0;
101 Coordinate center = new Coordinate();
102 for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
103 if (linearToMap.locate(p, center)) {
104 double depth = xyDepth.depth(center);
105 depths[i] = depth;
106 if (depth < maxDepth) {
107 maxDepth = depth;
108 }
109 }
110 }
111
112 if (maxDepth == Double.MAX_VALUE) {
113 log.warn("no depth found -> no interpolation");
114 return false;
115 }
116
117 if (debug) {
118 log.debug("max depth found: " + maxDepth);
119 }
120
121 double cellHeight = Math.abs(maxDepth)/height;
122
123 if (debug) {
124 log.debug("cell size: " + cellWidth + " x " + cellHeight);
125 }
126
127 double [] raster = new double[width*height];
128 Arrays.fill(raster, Double.NaN);
129
130 Envelope queryBuffer = new Envelope();
131 XYColumn [] neighbors = new XYColumn[4];
132
133 i = 0;
134 for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
135 double depth = depths[i];
136 if (Double.isNaN(depth)) {
137 continue;
138 }
139 linearToMap.locate(p, center);
140
141 queryBuffer.init(
142 center.x - dxMax, center.x + dxMax,
143 center.y - dyMax, center.y + dyMax);
144
145 List potential = spatialIndex.query(queryBuffer);
146
147 L1Comparator invL1 = new L1Comparator(center);
148 Collections.sort(potential, invL1);
149
150 neighbors[0] = neighbors[1] =
151 neighbors[2] = neighbors[3] = null;
152
153 /* bit code of neighbors
154 0---1
155 | x |
156 2---3
157 */
158
159 int mask = 0;
160 // reversed order is essential here!
161 for (int j = potential.size()-1; j >= 0; --j) {
162 XYColumn n = (XYColumn)potential.get(j);
163 int code = n.x > center.x ? 1 : 0;
164 if (n.y > center.y) code |= 2;
165 neighbors[code] = n;
166 mask |= 1 << code;
167 }
168
169 int numNeighbors = Integer.bitCount(mask);
170
171 // only interpolate if we have all four neighbors
172 // and we do not have any gaps.
173 if (numNeighbors == 4
174 && !neighbors[0].hasIGap(neighbors[1])
175 && !neighbors[1].hasJGap(neighbors[3])
176 && !neighbors[3].hasIGap(neighbors[2])
177 && !neighbors[2].hasJGap(neighbors[0])
178 && neighbors[0].prepare(xyDepth)
179 && neighbors[1].prepare(xyDepth)
180 && neighbors[2].prepare(xyDepth)
181 && neighbors[3].prepare(xyDepth)
182 ) {
183 double y1 = Interpolation2D.interpolate(
184 neighbors[0].x, neighbors[0].y,
185 neighbors[1].x, neighbors[1].y,
186 center.x);
187 double y2 = Interpolation2D.interpolate(
188 neighbors[2].x, neighbors[2].y,
189 neighbors[3].x, neighbors[3].y,
190 center.x);
191 int j = i;
192 for (double z = -cellHeight*0.5;
193 j < raster.length; z -= cellHeight, j += width) {
194
195 double v0 = neighbors[0].value(z);
196 double v1 = neighbors[1].value(z);
197 double v2 = neighbors[2].value(z);
198 double v3 = neighbors[3].value(z);
199
200 double z1 = Interpolation2D.interpolate(
201 neighbors[0].x, v0,
202 neighbors[1].x, v1,
203 center.x);
204 double z2 = Interpolation2D.interpolate(
205 neighbors[2].x, v2,
206 neighbors[3].x, v3,
207 center.x);
208 double value = Interpolation2D.interpolate(
209 y1, z1,
210 y2, z2,
211 center.y);
212 raster[j] = value;
213 }
214 } // down the x/y column
215 } // all along the track
216
217 this.depths = depths;
218 this.raster = raster;
219
220 return true;
221 }
222 }
223 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:
http://dive4elements.wald.intevation.org