Mercurial > dive4elements > gnv-client
comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java @ 434:0eed5749fd63
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: |