Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java @ 540:80630520e25a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged gnv-artifacts/0.4
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:13:49 +0200
parents 44415ae01ddb
children 9a828e5a2390
comparison
equal deleted inserted replaced
415:9f4a0b990d27 540:80630520e25a
1 package de.intevation.gnv.math;
2
3 import com.vividsolutions.jts.geom.Coordinate;
4 import com.vividsolutions.jts.geom.Envelope;
5
6 import com.vividsolutions.jts.index.strtree.STRtree;
7
8 import java.util.List;
9
10 import org.apache.log4j.Logger;
11
12 /**
13 * @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
14 */
15 public final class Interpolation2D
16 {
17 private static Logger log = Logger.getLogger(Interpolation2D.class);
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
24 public interface Consumer {
25 void interpolated(Coordinate point, boolean success);
26 } // interface Consumer
27
28 private Interpolation2D() {
29 }
30
31 public static final Envelope relevantArea(
32 List<? extends Coordinate> path,
33 List<? extends Coordinate> points
34 ) {
35 return relevantArea(path, points, CULL_POINT_THRESHOLD);
36 }
37
38 public static final Envelope relevantArea(
39 Envelope pathBBox,
40 List<? extends Coordinate> points
41 ) {
42 return relevantArea(pathBBox, points, CULL_POINT_THRESHOLD);
43 }
44
45 public static final Envelope relevantArea(
46 Envelope pathBBox,
47 List<? extends Coordinate> points,
48 int threshold
49 ) {
50 return points.size() < threshold
51 ? null
52 : relevantArea(
53 pathBBox,
54 pointsBoundingBox(points));
55 }
56
57 public static final Envelope relevantArea(
58 List<? extends Coordinate> path,
59 List<? extends Coordinate> points,
60 int threshold
61 ) {
62 return points.size() < threshold || path.isEmpty()
63 ? null
64 : relevantArea(
65 pointsBoundingBox(path),
66 pointsBoundingBox(points));
67 }
68
69 public static final Envelope relevantArea(
70 Envelope pathBBox,
71 Envelope pointsBBox
72 ) {
73 double pathArea = pathBBox.getWidth()*pathBBox.getHeight();
74 double pointsArea = pointsBBox.getWidth()*pointsBBox.getHeight();
75
76 if (pathArea > 0.8d*pointsArea) { return null; }
77
78 double nArea = 1.44d * pathArea;
79 if (nArea < 0.1d*pointsArea) nArea = 0.1d*pointsArea;
80 double w = pathBBox.getWidth();
81 double h = pathBBox.getHeight();
82 double [] d = solveQuadratic(1d, w+h, pathArea - nArea);
83
84 if (d == null) { return null; }
85
86 double extra = pos(d);
87
88 pathBBox.expandBy(extra);
89
90 return pathBBox;
91 }
92
93 public static final double [] solveQuadratic(
94 double a, double b, double c
95 ) {
96 double d = b*b - 4d*a*c;
97 if (d < 0d) { return null; }
98
99 d = Math.sqrt(d);
100 a = 1d/(2d*a);
101 b = -b;
102
103 return new double [] { a*(b + d), a*(b - d) };
104 }
105
106 public static final double pos(double [] x) {
107 return x[0] >= 0 ? x[0] : x[1];
108 }
109
110
111 public static Envelope pointsBoundingBox(
112 List<? extends Coordinate> path
113 ) {
114 int N = path.size();
115 Envelope area = new Envelope(path.get(N-1));
116
117 for (int i = N-2; i >= 0; --i) {
118 area.expandToInclude(path.get(i));
119 }
120
121 return area;
122 }
123
124 public static void interpolate(
125 List<? extends Coordinate> path,
126 List<? extends Point2d> points,
127 double from,
128 double to,
129 int steps,
130 Metrics metrics,
131 Consumer consumer
132 ) {
133 boolean debug = log.isDebugEnabled();
134
135 int N = path.size();
136 int M = points.size();
137
138 if (debug) {
139 log.debug("Size of path: " + N);
140 log.debug("Size of points: " + M);
141 }
142
143 if (M < 1 || N < 2) { // nothing to do
144 return;
145 }
146
147 List<GridCell> cells = GridCell.pointsToGridCells(
148 points, relevantArea(path, points));
149
150 if (cells.isEmpty()) {
151 log.warn("no cells found");
152 return;
153 }
154
155 // put into spatial index to speed up finding neighbors.
156 STRtree spatialIndex = new STRtree();
157
158 for (GridCell cell: cells) {
159 spatialIndex.insert(cell.getEnvelope(), cell);
160 }
161
162 LinearToMap linearToMap = new LinearToMap(
163 path, from, to, metrics);
164
165 double dP = (to - from)/steps;
166
167 Coordinate center = new Coordinate();
168 Envelope queryBuffer = new Envelope();
169 GridCell.CellFinder finder = new GridCell.CellFinder();
170
171 int missedInterpolations = 0;
172 int interpolations = 0;
173
174 for (double p = from; p <= to; p += dP) {
175 if (!linearToMap.locate(p, center)) {
176 continue;
177 }
178 queryBuffer.init(
179 center.x - EPS, center.x + EPS,
180 center.y - EPS, center.y + EPS);
181
182 finder.prepare(center);
183 spatialIndex.query(queryBuffer, finder);
184
185 GridCell found = finder.found;
186
187 if (found == null) {
188 consumer.interpolated(center, false);
189 ++missedInterpolations;
190 continue;
191 }
192
193 Point2d n0 = found.p1;
194 Point2d n1 = found.p2;
195 Point2d n2 = found.p3;
196 Point2d n3 = found.p4;
197
198 double z1 = interpolate(
199 n0.x, n0.z,
200 n1.x, n1.z,
201 center.x);
202 double z2 = interpolate(
203 n2.x, n2.z,
204 n3.x, n3.z,
205 center.x);
206 double y1 = interpolate(
207 n0.x, n0.y,
208 n1.x, n1.y,
209 center.x);
210 double y2 = interpolate(
211 n2.x, n2.y,
212 n3.x, n3.y,
213 center.x);
214 center.z = interpolate(
215 y1, z1,
216 y2, z2,
217 center.y);
218 consumer.interpolated(center, true);
219 ++interpolations;
220 }
221
222 if (debug) {
223 log.debug("interpolations: " +
224 interpolations + " / " + missedInterpolations);
225 }
226 }
227
228 public static final double interpolate(
229 double x1, double y1,
230 double x2, double y2,
231 double x
232 ) {
233 if (x2 == x1) {
234 return (y1 + y2)*0.5d;
235 }
236 double m = (y2-y1)/(x2-x1);
237 double b = y1 - m*x1;
238 return m*x + b;
239 }
240 }
241 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:
http://dive4elements.wald.intevation.org