Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java @ 376:d8f3ef441bf2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged gnv-artifacts/0.3
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:13:47 +0200
parents 086e3af38b96
children 04a242c67fe6
comparison
equal deleted inserted replaced
293:6b0ef2324d02 376:d8f3ef441bf2
1 package de.intevation.gnv.math;
2
3 import java.util.ArrayList;
4 import java.util.List;
5 import java.util.HashMap;
6 import java.util.Collections;
7
8 import com.vividsolutions.jts.geom.Coordinate;
9 import com.vividsolutions.jts.geom.Envelope;
10
11 import com.vividsolutions.jts.index.quadtree.Quadtree;
12
13 import org.apache.log4j.Logger;
14
15 /**
16 * @author Sascha L. Teichmann
17 */
18 public final class Interpolation2D
19 {
20 private static Logger log = Logger.getLogger(Interpolation2D.class);
21
22 public interface Consumer {
23 void interpolated(Coordinate point);
24 } // interface Consumer
25
26 private Interpolation2D() {
27 }
28
29 public static void interpolate(
30 List<? extends Coordinate> path,
31 List<? extends Point2d> points,
32 double from,
33 double to,
34 int steps,
35 Metrics metrics,
36 Consumer consumer
37 ) {
38 int N = path.size();
39 int M = points.size();
40
41 log.debug("Size of path: " + N);
42 log.debug("Size of points: " + M);
43
44 if (M < 1 || N < 2) { // nothing to do
45 return;
46 }
47
48 HashMap<Integer, ArrayList<Point2d>> map = new HashMap<Integer, ArrayList<Point2d>>();
49
50 for (int k = M-1; k >= 0; --k) {
51 Point2d p = points.get(k);
52
53 ArrayList<Point2d> list = map.get(p.j);
54
55 if (list == null) {
56 map.put(p.j, list = new ArrayList<Point2d>());
57 }
58 list.add(p);
59 }
60
61 double dxMax = -Double.MAX_VALUE;
62
63 for (ArrayList<Point2d> v: map.values()) {
64 Collections.sort(v, Point2d.X_COMPARATOR);
65 for (int i = 1, L = v.size(); i < L; ++i) {
66 double dx = Math.abs(v.get(i).x - v.get(i-1).x);
67 if (dx > dxMax) {
68 dxMax = dx;
69 }
70 }
71 }
72
73 dxMax = dxMax + 1e-5d;
74
75 map.clear();
76
77 for (int k = M-1; k >= 0; --k) {
78 Point2d p = points.get(k);
79
80 ArrayList<Point2d> list = map.get(p.i);
81
82 if (list == null) {
83 map.put(p.i, list = new ArrayList<Point2d>());
84 }
85 list.add(p);
86 }
87
88 double dyMax = -Double.MAX_VALUE;
89
90 for (ArrayList<Point2d> v: map.values()) {
91 Collections.sort(v, Point2d.Y_COMPARATOR);
92 for (int i = 1, L = v.size(); i < L; ++i) {
93 double dy = Math.abs(v.get(i).y - v.get(i-1).y);
94 if (dy > dyMax) {
95 dyMax = dy;
96 }
97 }
98 }
99
100 dyMax = dyMax + 1e-5d;
101
102 map = null;
103
104 log.debug("buffer size: " + dxMax + " / " + dyMax);
105
106 // put into spatial index to speed up finding neighbors.
107 Quadtree spatialIndex = new Quadtree();
108
109 for (int i = 0; i < M; ++i) {
110 Point2d p = points.get(i);
111 spatialIndex.insert(p.envelope(), p);
112 }
113
114 LinearToMap linearToMap = new LinearToMap(
115 path, from, to, metrics);
116
117 double dP = (to - from)/steps;
118
119 Coordinate center = new Coordinate();
120
121 Envelope queryBuffer = new Envelope();
122
123 Point2d [] neighbors = new Point2d[4];
124
125 int missedInterpolations = 0;
126 int interpolations = 0;
127
128 for (double p = from; p <= to; p += dP) {
129 if (!linearToMap.locate(p, center)) {
130 continue;
131 }
132 queryBuffer.init(
133 center.x - dxMax, center.x + dxMax,
134 center.y - dyMax, center.y + dyMax);
135
136 List potential = spatialIndex.query(queryBuffer);
137
138 L1Comparator invL1 = new L1Comparator(center);
139 Collections.sort(potential, invL1);
140
141 neighbors[0] = neighbors[1] =
142 neighbors[2] = neighbors[3] = null;
143
144 /* bit code of neighbors
145 0---1
146 | x |
147 2---3
148 */
149
150 int mask = 0;
151 // reversed order is essential here!
152 for (int i = potential.size()-1; i >= 0; --i) {
153 Point2d n = (Point2d)potential.get(i);
154 int code = n.x > center.x ? 1 : 0;
155 if (n.y > center.y) code |= 2;
156 neighbors[code] = n;
157 mask |= 1 << code;
158 }
159
160 int numNeighbors = Integer.bitCount(mask);
161
162 // only interpolate if we have all four neighbors
163 // and we do not have any gaps.
164 if (numNeighbors == 4
165 && !neighbors[0].hasIGap(neighbors[1])
166 && !neighbors[1].hasJGap(neighbors[3])
167 && !neighbors[3].hasIGap(neighbors[2])
168 && !neighbors[2].hasJGap(neighbors[0])
169 ) {
170 double z1 = interpolate(
171 neighbors[0].x, neighbors[0].z,
172 neighbors[1].x, neighbors[1].z,
173 center.x);
174 double z2 = interpolate(
175 neighbors[2].x, neighbors[2].z,
176 neighbors[3].x, neighbors[3].z,
177 center.x);
178 double y1 = interpolate(
179 neighbors[0].x, neighbors[0].y,
180 neighbors[1].x, neighbors[1].y,
181 center.x);
182 double y2 = interpolate(
183 neighbors[2].x, neighbors[2].y,
184 neighbors[3].x, neighbors[3].y,
185 center.x);
186 center.z = interpolate(
187 y1, z1,
188 y2, z2,
189 center.y);
190 consumer.interpolated(center);
191 ++interpolations;
192 }
193 else {
194 ++missedInterpolations;
195 }
196 }
197
198 log.debug("interpolations: " + interpolations + " / " + missedInterpolations);
199 }
200
201 public static final double interpolate(
202 double x1, double y1,
203 double x2, double y2,
204 double x
205 ) {
206 if (x2 == x1) {
207 return (y1 + y2)*0.5d;
208 }
209 double m = (y2-y1)/(x2-x1);
210 double b = y1 - m*x1;
211 return m*x + b;
212 }
213 }
214 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:
http://dive4elements.wald.intevation.org