Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/geom/Lines.java @ 3468:f37e7e8907cb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged flys-artifacts/2.8.1
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:14:39 +0200
parents ed07dd55f487
children
comparison
equal deleted inserted replaced
3387:5ffad8bde8ad 3468:f37e7e8907cb
1 package de.intevation.flys.artifacts.geom;
2
3 import java.util.ArrayList;
4 import java.util.List;
5 import java.util.Iterator;
6
7 import java.awt.geom.Point2D;
8 import java.awt.geom.Line2D;
9
10 import de.intevation.flys.artifacts.math.Linear;
11
12 import org.apache.log4j.Logger;
13
14 import gnu.trove.TDoubleArrayList;
15
16 /**
17 * Utility to create lines (intersect water with cross-section etc).
18 */
19 public class Lines
20 {
21 private static Logger log = Logger.getLogger(Lines.class);
22
23 public static final double EPSILON = 1e-4;
24
25 public static enum Mode { UNDEF, WET, DRY };
26
27
28 /** Never instantiate Lines, use static functions instead. */
29 protected Lines() {
30 }
31
32
33 /**
34 * Calculate area of polygon with four vertices.
35 * @return area of polygon with four vertices.
36 */
37 public static double area(Point2D p1, Point2D p2, Point2D p3, Point2D p4) {
38 double[] x = new double[] {p1.getX(), p2.getX(), p3.getX(), p4.getX(), p1.getX() };
39 double[] y = new double[] {p1.getY(), p2.getY(), p3.getY(), p4.getY(), p1.getY() };
40 double area = 0d;
41 for (int i=0; i <4; i++) {
42 area += (x[i] * y[i+1]) - (x[i+1] * y[i]);
43 }
44 return Math.abs(area * 0.5d);
45 }
46
47
48 /**
49 * Calculate the 'length' of the given lines.
50 * @param lines lines of which to calculate length.
51 */
52 public static double length(List<Line2D> lines) {
53 double sum = 0d;
54 for (Line2D line: lines) {
55 double xDiff = line.getX1() - line.getX2();
56 double yDiff = line.getY1() - line.getY2();
57 sum += Math.sqrt(xDiff*xDiff + yDiff*yDiff);
58 }
59 return sum;
60 }
61
62
63 /** List of lines and a double-precision area. */
64 private static class ListWithArea {
65 public List<Line2D> lines;
66 public double area;
67 public ListWithArea(List<Line2D> lines, double area) {
68 this.lines = lines;
69 this.area = area;
70 }
71 }
72
73
74 /**
75 * For a cross section given as points and a waterlevel (in meters),
76 * create a set of lines that represent the water surface, assuming it
77 * is distributed horizontally equally.
78 * @param points the points describing the river bed.
79 * @param waterLevel the height of the horizontal water line.
80 * @return A list of Lines representing the water surface and the
81 * calculated area between water surface and river bed.
82 */
83 public static ListWithArea fillWater(List<Point2D> points, double waterLevel) {
84
85 boolean debug = log.isDebugEnabled();
86
87 if (debug) {
88 log.debug("fillWater");
89 log.debug("----------------------------");
90 }
91
92 List<Line2D> result = new ArrayList();
93
94 int N = points.size();
95
96 if (N == 0) {
97 return new ListWithArea(result, 0d);
98 }
99
100 if (N == 1) {
101 Point2D p = points.get(0);
102 // Only generate point if over profile
103 if (waterLevel > p.getY()) {
104 result.add(new Line2D.Double(
105 p.getX(), waterLevel,
106 p.getX(), waterLevel));
107 }
108 // TODO continue calculating area.
109 return new ListWithArea(result, 0d);
110 }
111
112 double minX = Double.MAX_VALUE;
113 double minY = Double.MAX_VALUE;
114 double maxX = -Double.MAX_VALUE;
115 double maxY = -Double.MAX_VALUE;
116
117 // To ensure for sequences of equals x's that
118 // the original index order is preserved.
119 for (Point2D p: points) {
120 double x = p.getX(), y = p.getY();
121 if (x < minX) minX = x;
122 if (x > maxX) maxX = x;
123 if (y < minY) minY = y;
124 if (y > maxY) maxY = y;
125 }
126
127 if (minY > waterLevel) { // profile completely over water level
128 log.debug("complete over water");
129 return new ListWithArea(result, 0d);
130 }
131
132 if (waterLevel > maxY) { // water floods profile
133 log.debug("complete under water");
134 result.add(new Line2D.Double(minX, waterLevel, maxX, waterLevel));
135 return new ListWithArea(result, 0d);
136 }
137
138 // Water is sometimes above, sometimes under profile.
139 Mode mode = Mode.UNDEF;
140
141 double startX = minX;
142
143 double area = 0d;
144 // Walking along the profile.
145 for (int i = 1; i < N; ++i) {
146 Point2D p1 = points.get(i-1);
147 Point2D p2 = points.get(i);
148
149 if (p1.getY() < waterLevel && p2.getY() < waterLevel) {
150 // completely under water
151 if (debug) {
152 log.debug("under water: " + p1 + " " + p2);
153 }
154 if (mode != Mode.WET) {
155 startX = p1.getX();
156 mode = Mode.WET;
157 }
158 area += area(p1, p2,
159 new Point2D.Double(p2.getX(), waterLevel),
160 new Point2D.Double(p1.getX(), waterLevel));
161 continue;
162 }
163
164 // TODO trigger area calculation
165 if (p1.getY() > waterLevel && p2.getY() > waterLevel) {
166 if (debug) {
167 log.debug("over water: " + p1 + " " + p2);
168 }
169 // completely over water
170 if (mode == Mode.WET) {
171 log.debug("over/wet");
172 result.add(new Line2D.Double(
173 startX, waterLevel,
174 p1.getX(), waterLevel));
175 }
176 mode = Mode.DRY;
177 continue;
178 }
179
180 // TODO trigger area calculation
181 if (Math.abs(p1.getX() - p2.getX()) < EPSILON) {
182 // vertical line
183 switch (mode) {
184 case WET:
185 log.debug("vertical/wet");
186 mode = Mode.DRY;
187 result.add(new Line2D.Double(
188 startX, waterLevel,
189 p1.getX(), waterLevel));
190 break;
191 case DRY:
192 log.debug("vertical/dry");
193 mode = Mode.WET;
194 startX = p2.getX();
195 break;
196 default: // UNDEF
197 log.debug("vertical/undef");
198 if (p2.getY() < waterLevel) {
199 mode = Mode.WET;
200 startX = p2.getX();
201 }
202 else {
203 mode = Mode.DRY;
204 }
205 }
206 continue;
207 }
208
209 // check if waterlevel directly hits the vertices;
210
211 boolean p1W = Math.abs(waterLevel - p1.getY()) < EPSILON;
212 boolean p2W = Math.abs(waterLevel - p2.getY()) < EPSILON;
213
214 // TODO trigger area calculation
215 if (p1W || p2W) {
216 if (debug) {
217 log.debug("water hits vertex: " + p1 + " " + p2 + " " + mode);
218 }
219 if (p1W && p2W) { // parallel to water -> dry
220 log.debug("water hits both vertices");
221 if (mode == Mode.WET) {
222 result.add(new Line2D.Double(
223 startX, waterLevel,
224 p1.getX(), waterLevel));
225 }
226 mode = Mode.DRY;
227 }
228 else if (p1W) { // p1 == waterlevel
229 log.debug("water hits first vertex");
230 if (p2.getY() > waterLevel) { // --> dry
231 if (mode == Mode.WET) {
232 result.add(new Line2D.Double(
233 startX, waterLevel,
234 p1.getX(), waterLevel));
235 }
236 mode = Mode.DRY;
237 }
238 else { // --> wet
239 if (mode != Mode.WET) {
240 startX = p1.getX();
241 mode = Mode.WET;
242 }
243 area += area(p1, p2,
244 new Point2D.Double(p2.getX(), waterLevel),
245 new Point2D.Double(p2.getX(), waterLevel));
246 }
247 }
248 else { // p2 == waterlevel
249 log.debug("water hits second vertex");
250 if (p1.getY() > waterLevel) { // --> wet
251 if (mode != Mode.WET) {
252 startX = p2.getX();
253 mode = Mode.WET;
254 }
255 }
256 else { // --> dry
257 if (mode == Mode.WET) {
258 result.add(new Line2D.Double(
259 startX, waterLevel,
260 p2.getX(), waterLevel));
261 }
262 mode = Mode.DRY;
263 area += area(p1, p2,
264 new Point2D.Double(p1.getX(), waterLevel),
265 new Point2D.Double(p1.getX(), waterLevel));
266 }
267 }
268 if (debug) {
269 log.debug("mode is now: " + mode);
270 }
271 continue;
272 }
273
274 // TODO trigger area calculation
275 // intersection case
276 double x = Linear.linear(
277 waterLevel,
278 p1.getY(), p2.getY(),
279 p1.getX(), p2.getX());
280
281 if (debug) {
282 log.debug("intersection p1:" + p1);
283 log.debug("intersection p2:" + p2);
284 log.debug("intersection at x: " + x);
285 }
286
287 // Add area of that part of intersection that is 'wet'.
288 if (p1.getY() > waterLevel) {
289 area += area(new Point2D.Double(x, waterLevel),
290 p2,
291 new Point2D.Double(p2.getX(), waterLevel),
292 new Point2D.Double(x, waterLevel));
293 }
294 else {
295 area += area(new Point2D.Double(x, waterLevel),
296 p1,
297 new Point2D.Double(p1.getX(), waterLevel),
298 new Point2D.Double(x, waterLevel));
299 }
300
301 switch (mode) {
302 case WET:
303 log.debug("intersect/wet");
304 mode = Mode.DRY;
305 result.add(new Line2D.Double(
306 startX, waterLevel,
307 x, waterLevel));
308 break;
309
310 case DRY:
311 log.debug("intersect/dry");
312 mode = Mode.WET;
313 startX = x;
314 break;
315
316 default: // UNDEF
317 log.debug("intersect/undef");
318 if (p2.getY() > waterLevel) {
319 log.debug("intersect/undef/over");
320 mode = Mode.DRY;
321 result.add(new Line2D.Double(
322 p1.getX(), waterLevel,
323 x, waterLevel));
324 }
325 else {
326 mode = Mode.WET;
327 startX = x;
328 }
329 } // switch mode
330 } // for all points p[i] and p[i-1]
331
332 if (mode == Mode.WET) {
333 result.add(new Line2D.Double(
334 startX, waterLevel,
335 maxX, waterLevel));
336 }
337
338 return new ListWithArea(result, area);
339 }
340
341
342 /**
343 * Class holding points that form lines and the calculated length.
344 */
345 public static class LineData {
346 public double [][] points;
347 public double width;
348 public double area;
349 public LineData(double[][] points, double width, double area) {
350 this.points = points;
351 this.width = width;
352 this.area = area;
353 }
354 }
355
356
357 /** Return length of a single line. */
358 public static double lineLength(Line2D line) {
359 double xDiff = line.getX1() - line.getX2();
360 double yDiff = line.getY1() - line.getY2();
361 return Math.sqrt(xDiff*xDiff + yDiff*yDiff);
362 }
363
364
365 /**
366 * @param points the riverbed.
367 */
368 public static LineData createWaterLines(
369 List<Point2D> points,
370 double waterlevel
371 ) {
372 ListWithArea listAndArea = fillWater(points, waterlevel);
373 List<Line2D> lines = listAndArea.lines;
374
375 TDoubleArrayList lxs = new TDoubleArrayList();
376 TDoubleArrayList lys = new TDoubleArrayList();
377 double linesLength = 0.0f;
378
379 for (Iterator<Line2D> iter = lines.iterator(); iter.hasNext();) {
380 Line2D line = iter.next();
381 Point2D p1 = line.getP1();
382 Point2D p2 = line.getP2();
383 lxs.add(p1.getX());
384 lys.add(p1.getY());
385 lxs.add(p2.getX());
386 lys.add(p2.getY());
387
388 // Length calculation.
389 linesLength += lineLength(line);
390
391 if (iter.hasNext()) {
392 lxs.add(Double.NaN);
393 lys.add(Double.NaN);
394 }
395 }
396
397 return new LineData(
398 new double [][] { lxs.toNativeArray(), lys.toNativeArray() },
399 linesLength, listAndArea.area
400 );
401 }
402 }
403 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org