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: |