Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java @ 515:234d9892e497

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed "Profilschnitte" for gnv/issue153. gnv-artifacts/trunk@609 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 22 Jan 2010 18:51:47 +0000
parents 537e663d6c0c
children 464e03bf786b
comparison
equal deleted inserted replaced
514:d9d933e06875 515:234d9892e497
1 package de.intevation.gnv.math; 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 java.io.Serializable;
8
9 import java.awt.Dimension;
10 2
11 import com.vividsolutions.jts.geom.Coordinate; 3 import com.vividsolutions.jts.geom.Coordinate;
12 import com.vividsolutions.jts.geom.Envelope; 4 import com.vividsolutions.jts.geom.Envelope;
13 5
14 import com.vividsolutions.jts.index.quadtree.Quadtree; 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;
15 14
16 import org.apache.log4j.Logger; 15 import org.apache.log4j.Logger;
17 16
18 /** 17 /**
19 * @author Sascha L. Teichmann (sascha.teichmann@intevation.de) 18 * @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
23 { 22 {
24 private static Logger log = Logger.getLogger(Interpolation3D.class); 23 private static Logger log = Logger.getLogger(Interpolation3D.class);
25 24
26 public static final int DEFAULT_WIDTH = 1024; 25 public static final int DEFAULT_WIDTH = 1024;
27 public static final int DEFAULT_HEIGHT = 768; 26 public static final int DEFAULT_HEIGHT = 768;
27
28 public static final double EPS = 1e-6d;
28 29
29 protected int width; 30 protected int width;
30 protected int height; 31 protected int height;
31 32
32 protected double [] raster; 33 protected double [] raster;
92 93
93 if (M < 1 || N < 2) { // nothing to do 94 if (M < 1 || N < 2) { // nothing to do
94 return false; 95 return false;
95 } 96 }
96 97
97 double [] buffer = Interpolation2D.calculateBuffer(points); 98 List<GridCell> cells = GridCell.pointsToGridCells(points);
98 double dxMax = buffer[0]; 99
99 double dyMax = buffer[1]; 100 if (cells.isEmpty()) {
100 101 log.warn("no cells found");
101 if (debug) { 102 return false;
102 log.debug("buffer size: " + dxMax + " / " + dyMax);
103 } 103 }
104 104
105 // put into spatial index to speed up finding neighbors. 105 // put into spatial index to speed up finding neighbors.
106 Quadtree spatialIndex = new Quadtree(); 106 STRtree spatialIndex = new STRtree();
107 107
108 for (int i = 0; i < M; ++i) { 108 for (GridCell cell: cells) {
109 Point2d p = points.get(i); 109 spatialIndex.insert(cell.getEnvelope(), cell);
110 spatialIndex.insert(p.envelope(), p);
111 } 110 }
112 111
113 LinearToMap linearToMap = new LinearToMap( 112 LinearToMap linearToMap = new LinearToMap(
114 path, from, to, metrics); 113 path, from, to, metrics);
115 114
145 } 144 }
146 145
147 double [] raster = new double[width*height]; 146 double [] raster = new double[width*height];
148 Arrays.fill(raster, Double.NaN); 147 Arrays.fill(raster, Double.NaN);
149 148
150 Envelope queryBuffer = new Envelope(); 149 Envelope queryBuffer = new Envelope();
151 XYColumn [] neighbors = new XYColumn[4]; 150 GridCell.CellFinder finder = new GridCell.CellFinder();
152 151
153 i = 0; 152 i = 0;
154 for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) { 153 for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
155 double depth = depths[i]; 154 double depth = depths[i];
156 if (Double.isNaN(depth)) { 155 if (Double.isNaN(depth)) {
157 continue; 156 continue;
158 } 157 }
159 linearToMap.locate(p, center); 158 linearToMap.locate(p, center);
160 159
161 queryBuffer.init( 160 queryBuffer.init(
162 center.x - dxMax, center.x + dxMax, 161 center.x - EPS, center.x + EPS,
163 center.y - dyMax, center.y + dyMax); 162 center.y - EPS, center.y + EPS);
164 163
165 List potential = spatialIndex.query(queryBuffer); 164 finder.prepare(center);
166 165 spatialIndex.query(queryBuffer, finder);
167 L1Comparator invL1 = new L1Comparator(center); 166
168 Collections.sort(potential, invL1); 167 GridCell found = finder.found;
169 168
170 neighbors[0] = neighbors[1] = 169 if (found == null) {
171 neighbors[2] = neighbors[3] = null; 170 continue;
172 171 }
173 /* bit code of neighbors 172
174 0---1 173 XYColumn n0 = (XYColumn)found.p1;
175 | x | 174 XYColumn n1 = (XYColumn)found.p2;
176 2---3 175 XYColumn n2 = (XYColumn)found.p3;
177 */ 176 XYColumn n3 = (XYColumn)found.p4;
178 177
179 int mask = 0; 178 if (n0.prepare(xyDepth)
180 // reversed order is essential here! 179 && n1.prepare(xyDepth)
181 for (int j = potential.size()-1; j >= 0; --j) { 180 && n2.prepare(xyDepth)
182 XYColumn n = (XYColumn)potential.get(j); 181 && n3.prepare(xyDepth)
183 int code = n.x > center.x ? 1 : 0;
184 if (n.y > center.y) code |= 2;
185 neighbors[code] = n;
186 mask |= 1 << code;
187 }
188
189 int numNeighbors = Integer.bitCount(mask);
190
191 // only interpolate if we have all four neighbors
192 // and we do not have any gaps.
193 if (numNeighbors == 4
194 && !neighbors[0].hasIGap(neighbors[1])
195 && !neighbors[1].hasJGap(neighbors[3])
196 && !neighbors[3].hasIGap(neighbors[2])
197 && !neighbors[2].hasJGap(neighbors[0])
198 && neighbors[0].prepare(xyDepth)
199 && neighbors[1].prepare(xyDepth)
200 && neighbors[2].prepare(xyDepth)
201 && neighbors[3].prepare(xyDepth)
202 ) { 182 ) {
203 double y1 = Interpolation2D.interpolate( 183 double y1 = Interpolation2D.interpolate(
204 neighbors[0].x, neighbors[0].y, 184 n0.x, n0.y,
205 neighbors[1].x, neighbors[1].y, 185 n1.x, n1.y,
206 center.x); 186 center.x);
207 double y2 = Interpolation2D.interpolate( 187 double y2 = Interpolation2D.interpolate(
208 neighbors[2].x, neighbors[2].y, 188 n2.x, n2.y,
209 neighbors[3].x, neighbors[3].y, 189 n3.x, n3.y,
210 center.x); 190 center.x);
211 int j = i; 191 int j = i;
212 for (double z = -cellHeight*0.5; 192 for (double z = -cellHeight*0.5;
213 j < raster.length && z >= depth; 193 j < raster.length && z >= depth;
214 z -= cellHeight, j += width) { 194 z -= cellHeight, j += width) {
215 195
216 double v0 = neighbors[0].value(z); 196 double v0 = n0.value(z);
217 double v1 = neighbors[1].value(z); 197 double v1 = n1.value(z);
218 double v2 = neighbors[2].value(z); 198 double v2 = n2.value(z);
219 double v3 = neighbors[3].value(z); 199 double v3 = n3.value(z);
220 200
221 double z1 = Interpolation2D.interpolate( 201 double z1 = Interpolation2D.interpolate(
222 neighbors[0].x, v0, 202 n0.x, v0,
223 neighbors[1].x, v1, 203 n1.x, v1,
224 center.x); 204 center.x);
225 double z2 = Interpolation2D.interpolate( 205 double z2 = Interpolation2D.interpolate(
226 neighbors[2].x, v2, 206 n2.x, v2,
227 neighbors[3].x, v3, 207 n3.x, v3,
228 center.x); 208 center.x);
229 double value = Interpolation2D.interpolate( 209 double value = Interpolation2D.interpolate(
230 y1, z1, 210 y1, z1,
231 y2, z2, 211 y2, z2,
232 center.y); 212 center.y);
http://dive4elements.wald.intevation.org