Mercurial > dive4elements > gnv-client
comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java @ 361:aec85d00d82c
Added code to do 2D interpolations along a digitied track on the map.
Not connected to the transition model, yet.
gnv-artifacts/trunk@435 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 15 Dec 2009 22:25:53 +0000 |
parents | |
children | f66088a43ecc |
comparison
equal
deleted
inserted
replaced
360:ee760729f6b8 | 361:aec85d00d82c |
---|---|
1 package de.intevation.gnv.math; | |
2 | |
3 import java.util.List; | |
4 import java.util.Collections; | |
5 | |
6 import com.vividsolutions.jts.geom.Coordinate; | |
7 import com.vividsolutions.jts.geom.Envelope; | |
8 | |
9 import com.vividsolutions.jts.index.quadtree.Quadtree; | |
10 | |
11 /** | |
12 * @author Sascha L. Teichmann | |
13 */ | |
14 public final class Interpolation2D | |
15 { | |
16 public interface Consumer { | |
17 void interpolated(Coordinate point); | |
18 } // interface Consumer | |
19 | |
20 private Interpolation2D() { | |
21 } | |
22 | |
23 public static void interpolate( | |
24 List<? extends Coordinate> path, | |
25 List<? extends Point2d> points, | |
26 double from, | |
27 double to, | |
28 int steps, | |
29 Metrics metrics, | |
30 Consumer consumer | |
31 ) { | |
32 int N = path.size(); | |
33 int M = points.size(); | |
34 | |
35 if (M < 1 || N < 2) { // nothing to do | |
36 return; | |
37 } | |
38 // figure out max delta(p[i].x, p[i-1].x) | |
39 Collections.sort(points, Point2d.X_COMPARATOR); | |
40 double dxMax = -Double.MAX_VALUE; | |
41 for (int i = 1; i < M; ++i) { | |
42 double dx = Math.abs(path.get(i).x - path.get(i-1).x); | |
43 if (dx > dxMax) { | |
44 dxMax = dx; | |
45 } | |
46 } | |
47 | |
48 dxMax = dxMax*0.5d + 1e-5d; | |
49 | |
50 // figure out max delta(p[i].y, p[i-1].y) | |
51 Collections.sort(path, Point2d.X_COMPARATOR); | |
52 double dyMax = -Double.MAX_VALUE; | |
53 for (int i = 1; i < M; ++i) { | |
54 double dy = Math.abs(path.get(i).y - path.get(i-1).y); | |
55 if (dy > dyMax) { | |
56 dyMax = dy; | |
57 } | |
58 } | |
59 | |
60 dyMax = dyMax*0.5d + 1e-5d; | |
61 | |
62 // put into spatial index to speed up finding neighbors. | |
63 Quadtree spatialIndex = new Quadtree(); | |
64 | |
65 for (int i = 0; i < M; ++i) { | |
66 Point2d p = points.get(i); | |
67 spatialIndex.insert(p.envelope(), p); | |
68 } | |
69 | |
70 LinearToMap linearToMap = new LinearToMap( | |
71 path, from, to, metrics); | |
72 | |
73 double dP = (to - from)/steps; | |
74 | |
75 Coordinate center = new Coordinate(); | |
76 | |
77 Envelope queryBuffer = new Envelope(); | |
78 | |
79 Point2d [] neighbors = new Point2d[4]; | |
80 | |
81 for (double p = to; p <= from; p += dP) { | |
82 if (!linearToMap.locate(p, center)) { | |
83 continue; | |
84 } | |
85 queryBuffer.init( | |
86 center.x - dxMax, center.x + dxMax, | |
87 center.y - dyMax, center.y + dyMax); | |
88 | |
89 List potential = spatialIndex.query(queryBuffer); | |
90 | |
91 L1Comparator invL1 = new L1Comparator(center); | |
92 Collections.sort(potential, invL1); | |
93 | |
94 neighbors[0] = neighbors[1] = | |
95 neighbors[2] = neighbors[3] = null; | |
96 | |
97 /* bit code of neighbors | |
98 0---1 | |
99 | x | | |
100 2---3 | |
101 */ | |
102 | |
103 int mask = 0; | |
104 // reversed order is essential here! | |
105 for (int i = potential.size()-1; i >= 0; --i) { | |
106 Point2d n = (Point2d)potential.get(i); | |
107 int code = n.x > center.x ? 1 : 0; | |
108 if (n.y > center.y) code |= 2; | |
109 neighbors[code] = n; | |
110 mask |= 1 << code; | |
111 } | |
112 | |
113 int numNeighbors = Integer.bitCount(mask); | |
114 | |
115 // only interpolate if we have all four neighbors | |
116 // and we do not have any gaps. | |
117 if (numNeighbors == 4 | |
118 && !neighbors[0].hasIGap(neighbors[0]) | |
119 && !neighbors[1].hasJGap(neighbors[3]) | |
120 && !neighbors[3].hasIGap(neighbors[2]) | |
121 && !neighbors[2].hasJGap(neighbors[0]) | |
122 ) { | |
123 double z1 = interpolate( | |
124 neighbors[0].x, neighbors[0].z, | |
125 neighbors[1].x, neighbors[1].z, | |
126 center.x); | |
127 double z2 = interpolate( | |
128 neighbors[2].x, neighbors[2].z, | |
129 neighbors[3].x, neighbors[3].z, | |
130 center.x); | |
131 double y1 = interpolate( | |
132 neighbors[0].x, neighbors[0].y, | |
133 neighbors[1].x, neighbors[1].y, | |
134 center.x); | |
135 double y2 = interpolate( | |
136 neighbors[2].x, neighbors[2].y, | |
137 neighbors[3].x, neighbors[3].y, | |
138 center.x); | |
139 center.z = interpolate( | |
140 y1, z1, | |
141 y2, z2, | |
142 center.y); | |
143 consumer.interpolated(center); | |
144 } | |
145 } | |
146 } | |
147 | |
148 public static final double interpolate( | |
149 double x1, double y1, | |
150 double x2, double y2, | |
151 double x | |
152 ) { | |
153 if (x2 == x1) { | |
154 return (y1 + y2)*0.5d; | |
155 } | |
156 double m = (y2-y1)/(x2-x1); | |
157 double b = y1 - m*x1; | |
158 return m*x + b; | |
159 } | |
160 } | |
161 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: |