Mercurial > dive4elements > gnv-client
comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java @ 657:af3f56758f59
merged gnv-artifacts/0.5
| author | Thomas Arendsen Hein <thomas@intevation.de> |
|---|---|
| date | Fri, 28 Sep 2012 12:13:53 +0200 |
| parents | b248531fa20b |
| children | c4156275c1e1 |
comparison
equal
deleted
inserted
replaced
| 590:5f5f273c8566 | 657:af3f56758f59 |
|---|---|
| 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.awt.Dimension; | |
| 9 | |
| 10 import java.io.Serializable; | |
| 11 | |
| 12 import java.util.Arrays; | |
| 13 import java.util.List; | |
| 14 | |
| 15 import org.apache.log4j.Logger; | |
| 16 | |
| 17 /** | |
| 18 * @author Sascha L. Teichmann (sascha.teichmann@intevation.de) | |
| 19 */ | |
| 20 public class AreaInterpolation | |
| 21 implements Serializable | |
| 22 { | |
| 23 private static Logger log = Logger.getLogger(AreaInterpolation.class); | |
| 24 | |
| 25 protected double [] raster; | |
| 26 | |
| 27 protected int width; | |
| 28 protected int height; | |
| 29 | |
| 30 public AreaInterpolation() { | |
| 31 } | |
| 32 | |
| 33 public int getWidth() { | |
| 34 return width; | |
| 35 } | |
| 36 | |
| 37 public int getHeight() { | |
| 38 return height; | |
| 39 } | |
| 40 | |
| 41 public double [] getRaster() { | |
| 42 return raster; | |
| 43 } | |
| 44 | |
| 45 public static final double EPS = 1e-6d; | |
| 46 | |
| 47 | |
| 48 public boolean interpolate( | |
| 49 List<? extends Point2d> points, | |
| 50 Envelope boundingBox, | |
| 51 Dimension samples, | |
| 52 XYDepth depth, | |
| 53 int extrapolationRounds | |
| 54 ) { | |
| 55 boolean debug = log.isDebugEnabled(); | |
| 56 | |
| 57 if (points == null || points.isEmpty()) { | |
| 58 log.warn("no points to interpolate"); | |
| 59 return false; | |
| 60 } | |
| 61 | |
| 62 List<GridCell> cells = GridCell.pointsToGridCells( | |
| 63 points, | |
| 64 Interpolation2D.relevantArea( | |
| 65 boundingBox, | |
| 66 points), | |
| 67 extrapolationRounds); | |
| 68 | |
| 69 if (cells.isEmpty()) { | |
| 70 log.warn("no cells to interpolate"); | |
| 71 return false; | |
| 72 } | |
| 73 | |
| 74 int W = samples.width; | |
| 75 int H = samples.height; | |
| 76 | |
| 77 double cellWidth = boundingBox.getWidth() / W; | |
| 78 double cellHeight = boundingBox.getHeight() / H; | |
| 79 | |
| 80 STRtree spatialIndex = new STRtree(); | |
| 81 | |
| 82 for (GridCell cell: cells) { | |
| 83 spatialIndex.insert(cell.getEnvelope(), cell); | |
| 84 } | |
| 85 | |
| 86 if (debug) { | |
| 87 log.debug("width: " + boundingBox.getWidth()); | |
| 88 log.debug("height: " + boundingBox.getHeight()); | |
| 89 log.debug("sample width: " + W); | |
| 90 log.debug("sample height: " + H); | |
| 91 log.debug("cell width: " + cellWidth); | |
| 92 log.debug("cell height: " + cellHeight); | |
| 93 } | |
| 94 | |
| 95 Envelope queryBuffer = new Envelope(); | |
| 96 Coordinate center = new Coordinate(); | |
| 97 GridCell.CellFinder finder = new GridCell.CellFinder(); | |
| 98 | |
| 99 double [] raster = new double[W*H]; | |
| 100 Arrays.fill(raster, Double.NaN); | |
| 101 | |
| 102 double minX = boundingBox.getMinX(); | |
| 103 double minY = boundingBox.getMinY(); | |
| 104 | |
| 105 long startTime = System.currentTimeMillis(); | |
| 106 | |
| 107 int pos = 0; | |
| 108 for (int j = 0; j < H; ++j) { | |
| 109 | |
| 110 double y = j*cellHeight + 0.5d*cellHeight + minY; | |
| 111 double x = 0.5d*cellWidth + minX; | |
| 112 | |
| 113 for (int end = pos + W; pos < end; ++pos, x += cellWidth) { | |
| 114 | |
| 115 queryBuffer.init(x - EPS, x + EPS, y - EPS, y + EPS); | |
| 116 center.x = x; center.y = y; | |
| 117 finder.prepare(center); | |
| 118 spatialIndex.query(queryBuffer, finder); | |
| 119 | |
| 120 GridCell found = finder.found; | |
| 121 | |
| 122 if (found == null || depth.depth(center) > 0d) { | |
| 123 continue; | |
| 124 } | |
| 125 | |
| 126 double z1 = Interpolation2D.interpolate( | |
| 127 found.p1.x, found.p1.z, | |
| 128 found.p2.x, found.p2.z, | |
| 129 center.x); | |
| 130 double z2 = Interpolation2D.interpolate( | |
| 131 found.p3.x, found.p3.z, | |
| 132 found.p4.x, found.p4.z, | |
| 133 center.x); | |
| 134 double y1 = Interpolation2D.interpolate( | |
| 135 found.p1.x, found.p1.y, | |
| 136 found.p2.x, found.p2.y, | |
| 137 center.x); | |
| 138 double y2 = Interpolation2D.interpolate( | |
| 139 found.p3.x, found.p3.y, | |
| 140 found.p4.x, found.p4.y, | |
| 141 center.x); | |
| 142 raster[pos] = Interpolation2D.interpolate( | |
| 143 y1, z1, | |
| 144 y2, z2, | |
| 145 center.y); | |
| 146 } | |
| 147 } | |
| 148 | |
| 149 long stopTime = System.currentTimeMillis(); | |
| 150 | |
| 151 if (debug) { | |
| 152 log.debug("interpolation took: " + (stopTime - startTime)/1000f + " secs"); | |
| 153 } | |
| 154 | |
| 155 this.raster = raster; | |
| 156 this.width = W; | |
| 157 this.height = H; | |
| 158 | |
| 159 return true; | |
| 160 } | |
| 161 } | |
| 162 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : |