Mercurial > dive4elements > gnv-client
comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java @ 515:234d9892e497
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); |