Mercurial > dive4elements > gnv-client
comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java @ 474:ab29e4ff2fda
Added area interpolation needed for "Horizontalschnitt"
gnv-artifacts/trunk@540 c6561f87-3c4e-4783-a992-168aeb5c3f6f
| author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
|---|---|
| date | Thu, 14 Jan 2010 10:34:05 +0000 |
| parents | |
| children | d47b478e662b |
comparison
equal
deleted
inserted
replaced
| 473:a6a33ef35809 | 474:ab29e4ff2fda |
|---|---|
| 1 package de.intevation.gnv.math; | |
| 2 | |
| 3 import java.io.Serializable; | |
| 4 | |
| 5 import java.util.ArrayList; | |
| 6 import java.util.Arrays; | |
| 7 import java.util.List; | |
| 8 import java.util.HashMap; | |
| 9 import java.util.Collections; | |
| 10 | |
| 11 import java.awt.Dimension; | |
| 12 | |
| 13 import com.vividsolutions.jts.geom.Coordinate; | |
| 14 import com.vividsolutions.jts.geom.Envelope; | |
| 15 | |
| 16 import com.vividsolutions.jts.index.quadtree.Quadtree; | |
| 17 | |
| 18 import org.apache.log4j.Logger; | |
| 19 | |
| 20 /** | |
| 21 * @author Sascha L. Teichmann (sascha.teichmann@intevation.de) | |
| 22 */ | |
| 23 public class AreaInterpolation | |
| 24 implements Serializable | |
| 25 { | |
| 26 private static Logger log = Logger.getLogger(AreaInterpolation.class); | |
| 27 | |
| 28 protected double [] raster; | |
| 29 | |
| 30 protected int width; | |
| 31 protected int height; | |
| 32 | |
| 33 public AreaInterpolation() { | |
| 34 } | |
| 35 | |
| 36 public int getWidth() { | |
| 37 return width; | |
| 38 } | |
| 39 | |
| 40 public int getHeight() { | |
| 41 return height; | |
| 42 } | |
| 43 | |
| 44 public double [] getRaster() { | |
| 45 return raster; | |
| 46 } | |
| 47 | |
| 48 public boolean interpolate( | |
| 49 List<? extends Point2d> points, | |
| 50 Envelope boundingBox, | |
| 51 Dimension samples, | |
| 52 XYDepth depth | |
| 53 ) { | |
| 54 if (points == null || points.isEmpty()) { | |
| 55 return false; | |
| 56 } | |
| 57 | |
| 58 double [] buffer = Interpolation2D.calculateBuffer(points); | |
| 59 double dxMax = buffer[0]; | |
| 60 double dyMax = buffer[1]; | |
| 61 | |
| 62 Quadtree spatialIndex = new Quadtree(); | |
| 63 | |
| 64 for (int i = points.size()-1; i >= 0; --i) { | |
| 65 Point2d p = points.get(i); | |
| 66 spatialIndex.insert(p.envelope(), p); | |
| 67 } | |
| 68 | |
| 69 int W = samples.width; | |
| 70 int H = samples.height; | |
| 71 | |
| 72 double cellWidth = boundingBox.getWidth() / W; | |
| 73 double cellHeight = boundingBox.getHeight() / H; | |
| 74 | |
| 75 Envelope queryBuffer = new Envelope(); | |
| 76 Point2d [] neighbors = new Point2d[4]; | |
| 77 Coordinate center = new Coordinate(); | |
| 78 L1Comparator invL1 = new L1Comparator(center); | |
| 79 | |
| 80 double [] raster = new double[W*H]; | |
| 81 Arrays.fill(raster, Double.NaN); | |
| 82 | |
| 83 int row = 0; | |
| 84 for (int j = 0; j < H; ++j, row += W) { | |
| 85 double y = j*cellHeight + 0.5d*cellHeight; | |
| 86 double x = 0.5d*cellWidth; | |
| 87 for (int i = 0; i < W; ++i, x += cellWidth) { | |
| 88 queryBuffer.init( | |
| 89 x - dxMax, x + dxMax, | |
| 90 y - dyMax, y + dyMax); | |
| 91 | |
| 92 List potential = spatialIndex.query(queryBuffer); | |
| 93 | |
| 94 if (potential.isEmpty()) { | |
| 95 continue; | |
| 96 } | |
| 97 | |
| 98 center.x = x; center.y = y; | |
| 99 Collections.sort(potential, invL1); | |
| 100 | |
| 101 neighbors[0] = neighbors[1] = | |
| 102 neighbors[2] = neighbors[3] = null; | |
| 103 | |
| 104 int mask = 0; | |
| 105 // reversed order is essential here! | |
| 106 for (int k = potential.size()-1; k >= 0; --k) { | |
| 107 Point2d n = (Point2d)potential.get(k); | |
| 108 int code = n.x > center.x ? 1 : 0; | |
| 109 if (n.y > center.y) code |= 2; | |
| 110 neighbors[code] = n; | |
| 111 mask |= 1 << code; | |
| 112 } | |
| 113 | |
| 114 int numNeighbors = Integer.bitCount(mask); | |
| 115 | |
| 116 // only interpolate if we have all four neighbors | |
| 117 // we do not have any gaps and we are below surface. | |
| 118 if (numNeighbors == 4 | |
| 119 && !neighbors[0].hasIGap(neighbors[1]) | |
| 120 && !neighbors[1].hasJGap(neighbors[3]) | |
| 121 && !neighbors[3].hasIGap(neighbors[2]) | |
| 122 && !neighbors[2].hasJGap(neighbors[0]) | |
| 123 && depth.depth(center) <= 0d | |
| 124 ) { | |
| 125 double z1 = Interpolation2D.interpolate( | |
| 126 neighbors[0].x, neighbors[0].z, | |
| 127 neighbors[1].x, neighbors[1].z, | |
| 128 center.x); | |
| 129 double z2 = Interpolation2D.interpolate( | |
| 130 neighbors[2].x, neighbors[2].z, | |
| 131 neighbors[3].x, neighbors[3].z, | |
| 132 center.x); | |
| 133 double y1 = Interpolation2D.interpolate( | |
| 134 neighbors[0].x, neighbors[0].y, | |
| 135 neighbors[1].x, neighbors[1].y, | |
| 136 center.x); | |
| 137 double y2 = Interpolation2D.interpolate( | |
| 138 neighbors[2].x, neighbors[2].y, | |
| 139 neighbors[3].x, neighbors[3].y, | |
| 140 center.x); | |
| 141 raster[row+i] = Interpolation2D.interpolate( | |
| 142 y1, z1, | |
| 143 y2, z2, | |
| 144 center.y); | |
| 145 } | |
| 146 } | |
| 147 } | |
| 148 | |
| 149 this.raster = raster; | |
| 150 this.width = W; | |
| 151 this.height = H; | |
| 152 | |
| 153 return true; | |
| 154 } | |
| 155 } | |
| 156 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : |