Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java @ 875:5e9efdda6894

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged gnv-artifacts/1.0
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:13:56 +0200
parents 05bf8534a35a
children f953c9a559d8
comparison
equal deleted inserted replaced
722:bb3ffe7d719e 875:5e9efdda6894
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 * Does an area interpolation for "Horizontalschnitte".
19 *
20 * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
21 */
22 public class AreaInterpolation
23 implements Serializable
24 {
25 private static Logger log = Logger.getLogger(AreaInterpolation.class);
26
27 /**
28 * The generated raster.
29 */
30 protected double [] raster;
31
32 /**
33 * The width of the raster.
34 */
35 protected int width;
36 /**
37 * The height of the raster.
38 */
39 protected int height;
40
41 /**
42 * Default constructor.
43 */
44 public AreaInterpolation() {
45 }
46
47 /**
48 * Returns the width of the generated raster.
49 * @return the width of the raster.
50 */
51 public int getWidth() {
52 return width;
53 }
54
55 /**
56 * Returns the height of the raster.
57 * @return The height of the raster.
58 */
59 public int getHeight() {
60 return height;
61 }
62
63 /**
64 * The generated raster.
65 * @return the raster.
66 */
67 public double [] getRaster() {
68 return raster;
69 }
70
71 /**
72 * Epsilon for numerical stability.
73 */
74 public static final double EPS = 1e-6d;
75
76
77 /**
78 * Fills a raster by interpolating the input values.
79 * @param points The attributed input values.
80 * @param boundingBox The world area where to interpolate.
81 * @param samples The width and height of the raster.
82 * @param depth The callback to query the depth at a given point.
83 * @param extrapolationRounds Number of extrapolation point layers
84 * to generate before doing the interpolation.
85 * @return true if the interpolation succeeds else false.
86 */
87 public boolean interpolate(
88 List<? extends Point2d> points,
89 Envelope boundingBox,
90 Dimension samples,
91 XYDepth depth,
92 int extrapolationRounds
93 ) {
94 boolean debug = log.isDebugEnabled();
95
96 if (points == null || points.isEmpty()) {
97 log.warn("no points to interpolate");
98 return false;
99 }
100
101 List<GridCell> cells = GridCell.pointsToGridCells(
102 points,
103 Interpolation2D.relevantArea(
104 boundingBox,
105 points),
106 extrapolationRounds);
107
108 if (cells.isEmpty()) {
109 log.warn("no cells to interpolate");
110 return false;
111 }
112
113 int W = samples.width;
114 int H = samples.height;
115
116 double cellWidth = boundingBox.getWidth() / W;
117 double cellHeight = boundingBox.getHeight() / H;
118
119 STRtree spatialIndex = new STRtree();
120
121 for (GridCell cell: cells) {
122 spatialIndex.insert(cell.getEnvelope(), cell);
123 }
124
125 if (debug) {
126 log.debug("width: " + boundingBox.getWidth());
127 log.debug("height: " + boundingBox.getHeight());
128 log.debug("sample width: " + W);
129 log.debug("sample height: " + H);
130 log.debug("cell width: " + cellWidth);
131 log.debug("cell height: " + cellHeight);
132 }
133
134 Envelope queryBuffer = new Envelope();
135 Coordinate center = new Coordinate();
136 GridCell.CellFinder finder = new GridCell.CellFinder();
137
138 double [] raster = new double[W*H];
139 Arrays.fill(raster, Double.NaN);
140
141 double minX = boundingBox.getMinX();
142 double minY = boundingBox.getMinY();
143
144 long startTime = System.currentTimeMillis();
145
146 int pos = 0;
147 for (int j = 0; j < H; ++j) {
148
149 double y = j*cellHeight + 0.5d*cellHeight + minY;
150 double x = 0.5d*cellWidth + minX;
151
152 for (int end = pos + W; pos < end; ++pos, x += cellWidth) {
153
154 queryBuffer.init(x - EPS, x + EPS, y - EPS, y + EPS);
155 center.x = x; center.y = y;
156 finder.prepare(center);
157 spatialIndex.query(queryBuffer, finder);
158
159 GridCell found = finder.found;
160
161 if (found == null || depth.depth(center) > 0d) {
162 continue;
163 }
164
165 double z1 = Interpolation2D.interpolate(
166 found.p1.x, found.p1.z,
167 found.p2.x, found.p2.z,
168 center.x);
169 double z2 = Interpolation2D.interpolate(
170 found.p3.x, found.p3.z,
171 found.p4.x, found.p4.z,
172 center.x);
173 double y1 = Interpolation2D.interpolate(
174 found.p1.x, found.p1.y,
175 found.p2.x, found.p2.y,
176 center.x);
177 double y2 = Interpolation2D.interpolate(
178 found.p3.x, found.p3.y,
179 found.p4.x, found.p4.y,
180 center.x);
181 raster[pos] = Interpolation2D.interpolate(
182 y1, z1,
183 y2, z2,
184 center.y);
185 }
186 }
187
188 long stopTime = System.currentTimeMillis();
189
190 if (debug) {
191 log.debug("interpolation took: " +
192 (stopTime - startTime)/1000f + " secs");
193 }
194
195 this.raster = raster;
196 this.width = W;
197 this.height = H;
198
199 return true;
200 }
201 }
202 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org