Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java @ 519:4e347624ee7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Last part to fix gnv/issue153. Now 'Profilschnitte', 'Horizontalschnitte' and 'horizontale Schnittprofile' all use the same x/y interpolation code. gnv-artifacts/trunk@613 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Sat, 23 Jan 2010 21:16:45 +0000
parents d9d933e06875
children 44415ae01ddb
comparison
equal deleted inserted replaced
518:464e03bf786b 519:4e347624ee7c
1 package de.intevation.gnv.math; 1 package de.intevation.gnv.math;
2 2
3 import com.vividsolutions.jts.geom.Coordinate; 3 import com.vividsolutions.jts.geom.Coordinate;
4 import com.vividsolutions.jts.geom.Envelope; 4 import com.vividsolutions.jts.geom.Envelope;
5 5
6 import com.vividsolutions.jts.index.quadtree.Quadtree; 6 import com.vividsolutions.jts.index.strtree.STRtree;
7 7
8 import gnu.trove.TDoubleArrayList;
9
10 import java.util.ArrayList;
11 import java.util.Collections;
12 import java.util.HashMap;
13 import java.util.List; 8 import java.util.List;
14
15 import org.apache.commons.math.stat.descriptive.moment.Mean;
16 import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
17 9
18 import org.apache.log4j.Logger; 10 import org.apache.log4j.Logger;
19 11
20 /** 12 /**
21 * @author Sascha L. Teichmann (sascha.teichmann@intevation.de) 13 * @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
22 */ 14 */
23 public final class Interpolation2D 15 public final class Interpolation2D
24 { 16 {
25 private static Logger log = Logger.getLogger(Interpolation2D.class); 17 private static Logger log = Logger.getLogger(Interpolation2D.class);
26 18
19 public static final int CULL_POINT_THRESHOLD = Integer.getInteger(
20 "gnv.interpolation2d.cull.point.threshold", 1000);
21
22 public static final double EPS = 1e-6d;
23
27 public interface Consumer { 24 public interface Consumer {
28 void interpolated(Coordinate point, boolean success); 25 void interpolated(Coordinate point, boolean success);
29 } // interface Consumer 26 } // interface Consumer
30 27
31 private Interpolation2D() { 28 private Interpolation2D() {
32 } 29 }
33 30
34 private static final double dist(double a, double b) { 31 public static Envelope pathBoundingBox(
35 a -= b; 32 List<? extends Coordinate> path,
36 return Math.sqrt(a*a); 33 double extra
37 } 34 ) {
35 int N = path.size();
36 Envelope area = new Envelope(path.get(N-1));
38 37
39 private static final double OUTLIER_EPS = 0.4d; 38 for (int i = N-2; i >= 0; --i) {
40 39 area.expandToInclude(path.get(i));
41 public static final double [] calculateBufferRemoveOutlier(
42 List <? extends Point2d> points
43 ) {
44 HashMap<Integer, ArrayList<Point2d>> iMap =
45 new HashMap<Integer, ArrayList<Point2d>>();
46
47 HashMap<Integer, ArrayList<Point2d>> jMap =
48 new HashMap<Integer, ArrayList<Point2d>>();
49
50 for (int k = points.size()-1; k >= 0; --k) {
51 Point2d p = points.get(k);
52
53 ArrayList<Point2d> jList = jMap.get(p.j);
54 ArrayList<Point2d> iList = iMap.get(p.i);
55
56 if (jList == null) {
57 jMap.put(p.j, jList = new ArrayList<Point2d>());
58 }
59 jList.add(p);
60
61 if (iList == null) {
62 iMap.put(p.i, iList = new ArrayList<Point2d>());
63 }
64 iList.add(p);
65 } 40 }
66 41
67 TDoubleArrayList deltas = new TDoubleArrayList(); 42 area.expandBy(
43 extra*area.getWidth(),
44 extra*area.getHeight());
68 45
69 Mean mean = new Mean(); 46 return area;
70 StandardDeviation sd = new StandardDeviation();
71
72 for (ArrayList<Point2d> v: jMap.values()) {
73 Collections.sort(v, Point2d.Y_COMPARATOR);
74 for (int i = 1, L = v.size(); i < L; ++i) {
75 // double dy = Math.abs(v.get(i).y - v.get(i-1).y);
76 double dy = v.get(i).distance(v.get(i-1));
77 deltas.add(dy);
78 mean.increment(dy);
79 sd .increment(dy);
80 }
81 }
82
83 deltas.sort();
84
85 double m = mean.getResult();
86 double s = sd .getResult();
87 double eps = OUTLIER_EPS*s;
88
89 int yi = deltas.size() - 1;
90 while (yi > 0 && dist(deltas.get(yi), m) > eps) {
91 --yi;
92 }
93
94 double dyMax = deltas.get(yi) + 1e-5d;
95
96 if (log.isDebugEnabled()) {
97 log.debug("mean y: " + m);
98 log.debug("sigma y: " + s);
99 }
100
101 deltas.reset();
102 mean.clear();
103 sd .clear();
104
105 for (ArrayList<Point2d> v: iMap.values()) {
106 Collections.sort(v, Point2d.X_COMPARATOR);
107 for (int i = 1, L = v.size(); i < L; ++i) {
108 //double dx = Math.abs(v.get(i).x - v.get(i-1).x);
109 double dx = v.get(i).distance(v.get(i-1));
110 deltas.add(dx);
111 mean.increment(dx);
112 sd .increment(dx);
113 }
114 }
115
116 deltas.sort();
117
118 m = mean.getResult();
119 s = sd .getResult();
120 eps = OUTLIER_EPS*s;
121
122 int xi = deltas.size() - 1;
123
124 while (xi > 0 && dist(deltas.get(xi), m) > eps) {
125 --xi;
126 }
127
128 double dxMax = deltas.get(xi) + 1e-5d;
129
130 if (log.isDebugEnabled()) {
131 log.debug("mean y: " + m);
132 log.debug("sigma y: " + s);
133 }
134
135 return new double [] { dxMax, dyMax };
136 }
137
138 public static final double [] calculateBuffer(
139 List <? extends Point2d> points
140 ) {
141 HashMap<Integer, ArrayList<Point2d>> iMap =
142 new HashMap<Integer, ArrayList<Point2d>>();
143
144 HashMap<Integer, ArrayList<Point2d>> jMap =
145 new HashMap<Integer, ArrayList<Point2d>>();
146
147 for (int k = points.size()-1; k >= 0; --k) {
148 Point2d p = points.get(k);
149
150 ArrayList<Point2d> jList = jMap.get(p.j);
151 ArrayList<Point2d> iList = iMap.get(p.i);
152
153 if (jList == null) {
154 jMap.put(p.j, jList = new ArrayList<Point2d>());
155 }
156 jList.add(p);
157
158 if (iList == null) {
159 iMap.put(p.i, iList = new ArrayList<Point2d>());
160 }
161 iList.add(p);
162 }
163
164 double dxMax = -Double.MAX_VALUE;
165 double dyMax = -Double.MAX_VALUE;
166
167 for (ArrayList<Point2d> v: jMap.values()) {
168 Collections.sort(v, Point2d.Y_COMPARATOR);
169 for (int i = 1, L = v.size(); i < L; ++i) {
170 double dy = Math.abs(v.get(i).y - v.get(i-1).y);
171 if (dy > dyMax) {
172 dyMax = dy;
173 }
174 }
175 }
176
177 dyMax += 1e-5d;
178
179 for (ArrayList<Point2d> v: iMap.values()) {
180 Collections.sort(v, Point2d.X_COMPARATOR);
181 for (int i = 1, L = v.size(); i < L; ++i) {
182 double dx = Math.abs(v.get(i).x - v.get(i-1).x);
183 if (dx > dxMax) {
184 dxMax = dx;
185 }
186 }
187 }
188
189 dxMax += 1e-5d;
190
191 return new double [] { dxMax, dyMax };
192 } 47 }
193 48
194 public static void interpolate( 49 public static void interpolate(
195 List<? extends Coordinate> path, 50 List<? extends Coordinate> path,
196 List<? extends Point2d> points, 51 List<? extends Point2d> points,
212 67
213 if (M < 1 || N < 2) { // nothing to do 68 if (M < 1 || N < 2) { // nothing to do
214 return; 69 return;
215 } 70 }
216 71
217 double [] buffer = calculateBuffer(points); 72 Envelope relevantArea = M > CULL_POINT_THRESHOLD
218 double dxMax = buffer[0]; 73 ? pathBoundingBox(path, 0.05d)
219 double dyMax = buffer[1]; 74 : null;
75
76 List<GridCell> cells = GridCell.pointsToGridCells(
77 points, relevantArea);
220 78
221 if (debug) { 79 if (cells.isEmpty()) {
222 log.debug("buffer size: " + dxMax + " / " + dyMax); 80 log.warn("no cells found");
81 return;
223 } 82 }
224 83
225 // put into spatial index to speed up finding neighbors. 84 // put into spatial index to speed up finding neighbors.
226 Quadtree spatialIndex = new Quadtree(); 85 STRtree spatialIndex = new STRtree();
227 86
228 for (int i = 0; i < M; ++i) { 87 for (GridCell cell: cells) {
229 Point2d p = points.get(i); 88 spatialIndex.insert(cell.getEnvelope(), cell);
230 spatialIndex.insert(p.envelope(), p);
231 } 89 }
232 90
233 LinearToMap linearToMap = new LinearToMap( 91 LinearToMap linearToMap = new LinearToMap(
234 path, from, to, metrics); 92 path, from, to, metrics);
235 93
236 double dP = (to - from)/steps; 94 double dP = (to - from)/steps;
237 95
238 Coordinate center = new Coordinate(); 96 Coordinate center = new Coordinate();
239 97 Envelope queryBuffer = new Envelope();
240 Envelope queryBuffer = new Envelope(); 98 GridCell.CellFinder finder = new GridCell.CellFinder();
241
242 Point2d [] neighbors = new Point2d[4];
243 99
244 int missedInterpolations = 0; 100 int missedInterpolations = 0;
245 int interpolations = 0; 101 int interpolations = 0;
246
247 L1Comparator invL1 = new L1Comparator(center);
248 102
249 for (double p = from; p <= to; p += dP) { 103 for (double p = from; p <= to; p += dP) {
250 if (!linearToMap.locate(p, center)) { 104 if (!linearToMap.locate(p, center)) {
251 continue; 105 continue;
252 } 106 }
253 queryBuffer.init( 107 queryBuffer.init(
254 center.x - dxMax, center.x + dxMax, 108 center.x - EPS, center.x + EPS,
255 center.y - dyMax, center.y + dyMax); 109 center.y - EPS, center.y + EPS);
256 110
257 List potential = spatialIndex.query(queryBuffer); 111 finder.prepare(center);
112 spatialIndex.query(queryBuffer, finder);
258 113
259 Collections.sort(potential, invL1); 114 GridCell found = finder.found;
260 115
261 neighbors[0] = neighbors[1] = 116 if (found == null) {
262 neighbors[2] = neighbors[3] = null; 117 consumer.interpolated(center, false);
263 118 ++missedInterpolations;
264 /* bit code of neighbors 119 continue;
265 0---1
266 | x |
267 2---3
268 */
269
270 int mask = 0;
271 // reversed order is essential here!
272 for (int i = potential.size()-1; i >= 0; --i) {
273 Point2d n = (Point2d)potential.get(i);
274 int code = n.x > center.x ? 1 : 0;
275 if (n.y > center.y) code |= 2;
276 neighbors[code] = n;
277 mask |= 1 << code;
278 } 120 }
279 121
280 int numNeighbors = Integer.bitCount(mask); 122 Point2d n0 = found.p1;
123 Point2d n1 = found.p2;
124 Point2d n2 = found.p3;
125 Point2d n3 = found.p4;
281 126
282 // only interpolate if we have all four neighbors 127 double z1 = interpolate(
283 // and we do not have any gaps. 128 n0.x, n0.z,
284 if (numNeighbors == 4 129 n1.x, n1.z,
285 && !neighbors[0].hasIGap(neighbors[1]) 130 center.x);
286 && !neighbors[1].hasJGap(neighbors[3]) 131 double z2 = interpolate(
287 && !neighbors[3].hasIGap(neighbors[2]) 132 n2.x, n2.z,
288 && !neighbors[2].hasJGap(neighbors[0]) 133 n3.x, n3.z,
289 ) { 134 center.x);
290 double z1 = interpolate( 135 double y1 = interpolate(
291 neighbors[0].x, neighbors[0].z, 136 n0.x, n0.y,
292 neighbors[1].x, neighbors[1].z, 137 n1.x, n1.y,
293 center.x); 138 center.x);
294 double z2 = interpolate( 139 double y2 = interpolate(
295 neighbors[2].x, neighbors[2].z, 140 n2.x, n2.y,
296 neighbors[3].x, neighbors[3].z, 141 n3.x, n3.y,
297 center.x); 142 center.x);
298 double y1 = interpolate( 143 center.z = interpolate(
299 neighbors[0].x, neighbors[0].y, 144 y1, z1,
300 neighbors[1].x, neighbors[1].y, 145 y2, z2,
301 center.x); 146 center.y);
302 double y2 = interpolate( 147 consumer.interpolated(center, true);
303 neighbors[2].x, neighbors[2].y, 148 ++interpolations;
304 neighbors[3].x, neighbors[3].y,
305 center.x);
306 center.z = interpolate(
307 y1, z1,
308 y2, z2,
309 center.y);
310 consumer.interpolated(center, true);
311 ++interpolations;
312 }
313 else {
314 consumer.interpolated(center, false);
315 ++missedInterpolations;
316 }
317 } 149 }
318 150
319 if (debug) { 151 if (debug) {
320 log.debug("interpolations: " + 152 log.debug("interpolations: " +
321 interpolations + " / " + missedInterpolations); 153 interpolations + " / " + missedInterpolations);
http://dive4elements.wald.intevation.org