Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java @ 514:d9d933e06875

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed gnv/issue153 gnv-artifacts/trunk@608 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 22 Jan 2010 18:22:11 +0000
parents 70adafe2b9d5
children 96a1e92e7ed2
comparison
equal deleted inserted replaced
513:ca5048e4e515 514:d9d933e06875
1 package de.intevation.gnv.math; 1 package de.intevation.gnv.math;
2 2
3 import com.vividsolutions.jts.geom.Coordinate; 3 import com.vividsolutions.jts.geom.Coordinate;
4 import com.vividsolutions.jts.geom.Envelope; 4 import com.vividsolutions.jts.geom.Envelope;
5 5
6 import com.vividsolutions.jts.index.quadtree.Quadtree; 6 import com.vividsolutions.jts.index.strtree.STRtree;
7 7
8 import java.awt.Dimension; 8 import java.awt.Dimension;
9 9
10 import java.io.Serializable; 10 import java.io.Serializable;
11 11
12 import java.util.ArrayList;
13 import java.util.Arrays; 12 import java.util.Arrays;
14 import java.util.Collections;
15 import java.util.List; 13 import java.util.List;
16 14
17 import org.apache.log4j.Logger; 15 import org.apache.log4j.Logger;
18 16
19 /** 17 /**
42 40
43 public double [] getRaster() { 41 public double [] getRaster() {
44 return raster; 42 return raster;
45 } 43 }
46 44
45 public static final double EPS = 1e-6d;
46
47
47 public boolean interpolate( 48 public boolean interpolate(
48 List<? extends Point2d> points, 49 List<? extends Point2d> points,
49 Envelope boundingBox, 50 Envelope boundingBox,
50 Dimension samples, 51 Dimension samples,
51 XYDepth depth 52 XYDepth depth
52 ) { 53 ) {
53 boolean debug = log.isDebugEnabled(); 54 boolean debug = log.isDebugEnabled();
54 55
55 if (points == null || points.isEmpty()) { 56 if (points == null || points.isEmpty()) {
56 if (debug) { 57 log.warn("no points to interpolate");
57 log.debug("no points to interpolate"); 58 return false;
58 } 59 }
60
61 List<GridCell> cells = GridCell.pointsToGridCells(points);
62
63 if (cells.isEmpty()) {
64 log.warn("no cells to interpolate");
59 return false; 65 return false;
60 } 66 }
61 67
62 int W = samples.width; 68 int W = samples.width;
63 int H = samples.height; 69 int H = samples.height;
64 70
65 double cellWidth = boundingBox.getWidth() / W; 71 double cellWidth = boundingBox.getWidth() / W;
66 double cellHeight = boundingBox.getHeight() / H; 72 double cellHeight = boundingBox.getHeight() / H;
67 Envelope gridEnvelope = new Envelope(boundingBox);
68 73
69 gridEnvelope.expandBy(2d*cellWidth, 2d*cellWidth); 74 STRtree spatialIndex = new STRtree();
70 75
71 ArrayList<Point2d>relevantPoints = 76 for (GridCell cell: cells) {
72 new ArrayList<Point2d>(); 77 spatialIndex.insert(cell.getEnvelope(), cell);
73
74 for (int i = points.size()-1; i >= 0; --i) {
75 if (gridEnvelope.contains(points.get(i))) {
76 relevantPoints.add(points.get(i));
77 }
78 }
79
80 double [] buffer = Interpolation2D
81 .calculateBufferRemoveOutlier(relevantPoints);
82 double dxMax = buffer[0];
83 double dyMax = buffer[1];
84
85 Quadtree spatialIndex = new Quadtree();
86
87 for (int i = relevantPoints.size()-1; i >= 0; --i) {
88 Point2d p = relevantPoints.get(i);
89 spatialIndex.insert(p.envelope(), p);
90 } 78 }
91 79
92 if (debug) { 80 if (debug) {
93 log.debug("points: " + relevantPoints.size());
94 log.debug("width: " + boundingBox.getWidth()); 81 log.debug("width: " + boundingBox.getWidth());
95 log.debug("height: " + boundingBox.getHeight()); 82 log.debug("height: " + boundingBox.getHeight());
96 log.debug("sample width: " + W); 83 log.debug("sample width: " + W);
97 log.debug("sample height: " + H); 84 log.debug("sample height: " + H);
98 log.debug("cell width: " + cellWidth); 85 log.debug("cell width: " + cellWidth);
99 log.debug("cell height: " + cellHeight); 86 log.debug("cell height: " + cellHeight);
100 log.debug("buffer x: " + dxMax);
101 log.debug("buffer y: " + dyMax);
102 } 87 }
103 88
104 Envelope queryBuffer = new Envelope(); 89 Envelope queryBuffer = new Envelope();
105 Point2d [] neighbors = new Point2d[4]; 90 Coordinate center = new Coordinate();
106 Coordinate center = new Coordinate(); 91 GridCell.CellFinder finder = new GridCell.CellFinder();
107 L1Comparator invL1 = new L1Comparator(center);
108 92
109 double [] raster = new double[W*H]; 93 double [] raster = new double[W*H];
110 Arrays.fill(raster, Double.NaN); 94 Arrays.fill(raster, Double.NaN);
111 95
112 double minX = boundingBox.getMinX(); 96 double minX = boundingBox.getMinX();
117 int row = 0; 101 int row = 0;
118 for (int j = 0; j < H; ++j, row += W) { 102 for (int j = 0; j < H; ++j, row += W) {
119 double y = j*cellHeight + 0.5d*cellHeight + minY; 103 double y = j*cellHeight + 0.5d*cellHeight + minY;
120 double x = 0.5d*cellWidth + minX; 104 double x = 0.5d*cellWidth + minX;
121 for (int i = 0; i < W; ++i, x += cellWidth) { 105 for (int i = 0; i < W; ++i, x += cellWidth) {
122 queryBuffer.init(
123 x - dxMax, x + dxMax,
124 y - dyMax, y + dyMax);
125 106
126 List potential = spatialIndex.query(queryBuffer); 107 queryBuffer.init(x - EPS, x + EPS, y - EPS, y + EPS);
108 center.x = x; center.y = y;
109 finder.prepare(center);
110 spatialIndex.query(queryBuffer, finder);
127 111
128 if (potential.size() < 4) { 112 GridCell found = finder.found;
113
114 if (found == null) {
129 continue; 115 continue;
130 } 116 }
131 117
132 center.x = x; center.y = y; 118 double z1 = Interpolation2D.interpolate(
133 Collections.sort(potential, invL1); 119 found.p1.x, found.p1.z,
134 120 found.p2.x, found.p2.z,
135 neighbors[0] = neighbors[1] = 121 center.x);
136 neighbors[2] = neighbors[3] = null; 122 double z2 = Interpolation2D.interpolate(
137 123 found.p3.x, found.p3.z,
138 int mask = 0; 124 found.p4.x, found.p4.z,
139 // reversed order is essential here! 125 center.x);
140 for (int k = potential.size()-1; k >= 0; --k) { 126 double y1 = Interpolation2D.interpolate(
141 Point2d n = (Point2d)potential.get(k); 127 found.p1.x, found.p1.y,
142 int code = n.x > center.x ? 1 : 0; 128 found.p2.x, found.p2.y,
143 if (n.y > center.y) code |= 2; 129 center.x);
144 neighbors[code] = n; 130 double y2 = Interpolation2D.interpolate(
145 mask |= 1 << code; 131 found.p3.x, found.p3.y,
146 } 132 found.p4.x, found.p4.y,
147 133 center.x);
148 int numNeighbors = Integer.bitCount(mask); 134 raster[row+i] = Interpolation2D.interpolate(
149 135 y1, z1,
150 // only interpolate if we have all four neighbors 136 y2, z2,
151 // we do not have any gaps and we are below surface. 137 center.y);
152 if (numNeighbors == 4
153 && !neighbors[0].hasIGap(neighbors[1])
154 && !neighbors[1].hasJGap(neighbors[3])
155 && !neighbors[3].hasIGap(neighbors[2])
156 && !neighbors[2].hasJGap(neighbors[0])
157 && depth.depth(center) <= 0d
158 ) {
159 double z1 = Interpolation2D.interpolate(
160 neighbors[0].x, neighbors[0].z,
161 neighbors[1].x, neighbors[1].z,
162 center.x);
163 double z2 = Interpolation2D.interpolate(
164 neighbors[2].x, neighbors[2].z,
165 neighbors[3].x, neighbors[3].z,
166 center.x);
167 double y1 = Interpolation2D.interpolate(
168 neighbors[0].x, neighbors[0].y,
169 neighbors[1].x, neighbors[1].y,
170 center.x);
171 double y2 = Interpolation2D.interpolate(
172 neighbors[2].x, neighbors[2].y,
173 neighbors[3].x, neighbors[3].y,
174 center.x);
175 raster[row+i] = Interpolation2D.interpolate(
176 y1, z1,
177 y2, z2,
178 center.y);
179 }
180 } 138 }
181 } 139 }
182 140
183 long stopTime = System.currentTimeMillis(); 141 long stopTime = System.currentTimeMillis();
184 142
http://dive4elements.wald.intevation.org