sascha@938: package de.intevation.flys.geom; sascha@933: sascha@1651: import java.util.ArrayList; sascha@933: import java.util.List; sascha@1651: import java.util.Iterator; sascha@933: sascha@933: import java.awt.geom.Point2D; sascha@933: import java.awt.geom.Line2D; sascha@933: sascha@933: import de.intevation.flys.artifacts.math.Linear; sascha@933: sascha@938: import org.apache.log4j.Logger; sascha@938: sascha@1651: import gnu.trove.TDoubleArrayList; sascha@1651: felix@2652: /** felix@2652: * Utility to create lines (intersect water with cross-section etc). felix@2652: */ sascha@933: public class Lines sascha@933: { sascha@938: private static Logger log = Logger.getLogger(Lines.class); sascha@938: sascha@933: public static final double EPSILON = 1e-4; sascha@933: sascha@933: public static enum Mode { UNDEF, WET, DRY }; sascha@933: felix@2681: felix@2681: /** Never instantiate Lines, use static functions instead. */ sascha@933: protected Lines() { sascha@933: } sascha@933: felix@2673: felix@2673: /** felix@2673: * Calculate area of polygon with four vertices. felix@2673: * @return area of polygon with four vertices. felix@2673: */ felix@2673: public static double area(Point2D p1, Point2D p2, Point2D p3, Point2D p4) { felix@2688: double[] x = new double[] {p1.getX(), p2.getX(), p3.getX(), p4.getX(), p1.getX() }; felix@2688: double[] y = new double[] {p1.getY(), p2.getY(), p3.getY(), p4.getY(), p1.getY() }; felix@2673: double area = 0d; felix@2673: for (int i=0; i <4; i++) { felix@2688: area += (x[i] * y[i+1]) - (x[i+1] * y[i]); felix@2673: } felix@2688: return Math.abs(area * 0.5d); felix@2673: } felix@2673: felix@2673: felix@2652: /** felix@2652: * Calculate the 'length' of the given lines. felix@2681: * @param lines lines of which to calculate length. felix@2652: */ sascha@2650: public static double length(List lines) { sascha@2650: double sum = 0d; sascha@2650: for (Line2D line: lines) { sascha@2650: double xDiff = line.getX1() - line.getX2(); sascha@2650: double yDiff = line.getY1() - line.getY2(); sascha@2650: sum += Math.sqrt(xDiff*xDiff + yDiff*yDiff); sascha@2650: } sascha@2650: return sum; sascha@2650: } sascha@2650: felix@2674: felix@2674: /** List of lines and a double-precision area. */ felix@2674: private static class ListWithArea { felix@2674: public List lines; felix@2674: public double area; felix@2674: public ListWithArea(List lines, double area) { felix@2674: this.lines = lines; felix@2674: this.area = area; felix@2674: } felix@2674: } felix@2674: felix@2674: felix@2679: /** felix@2679: * For a cross section given as points and a waterlevel (in meters), felix@2679: * create a set of lines that represent the water surface, assuming it felix@2679: * is distributed horizontally equally. felix@2679: * @param points the points describing the river bed. felix@2679: * @param waterLevel the height of the horizontal water line. felix@2679: * @return A list of Lines representing the water surface and the felix@2679: * calculated area between water surface and river bed. felix@2679: */ felix@2674: public static ListWithArea fillWater(List points, double waterLevel) { sascha@933: sascha@938: boolean debug = log.isDebugEnabled(); sascha@938: sascha@938: if (debug) { sascha@938: log.debug("fillWater"); sascha@938: log.debug("----------------------------"); sascha@933: } sascha@933: sascha@933: List result = new ArrayList(); sascha@933: sascha@933: int N = points.size(); sascha@933: sascha@933: if (N == 0) { felix@2674: return new ListWithArea(result, 0d); sascha@933: } sascha@933: sascha@933: if (N == 1) { sascha@933: Point2D p = points.get(0); felix@2673: // Only generate point if over profile sascha@933: if (waterLevel > p.getY()) { sascha@933: result.add(new Line2D.Double( sascha@934: p.getX(), waterLevel, sascha@934: p.getX(), waterLevel)); sascha@933: } felix@2679: // TODO continue calculating area. felix@2674: return new ListWithArea(result, 0d); sascha@933: } sascha@933: sascha@933: double minX = Double.MAX_VALUE; sascha@933: double minY = Double.MAX_VALUE; sascha@933: double maxX = -Double.MAX_VALUE; sascha@933: double maxY = -Double.MAX_VALUE; sascha@933: sascha@933: // To ensure for sequences of equals x's that sascha@933: // the original index order is preserved. sascha@938: for (Point2D p: points) { sascha@933: double x = p.getX(), y = p.getY(); sascha@933: if (x < minX) minX = x; sascha@933: if (x > maxX) maxX = x; sascha@933: if (y < minY) minY = y; sascha@933: if (y > maxY) maxY = y; sascha@933: } sascha@933: sascha@933: if (minY > waterLevel) { // profile completely over water level sascha@938: log.debug("complete over water"); felix@2674: return new ListWithArea(result, 0d); sascha@933: } sascha@933: sascha@933: if (waterLevel > maxY) { // water floods profile sascha@938: log.debug("complete under water"); sascha@933: result.add(new Line2D.Double(minX, waterLevel, maxX, waterLevel)); felix@2674: return new ListWithArea(result, 0d); sascha@933: } sascha@933: felix@2673: // Water is sometimes above, sometimes under profile. sascha@933: Mode mode = Mode.UNDEF; sascha@933: sascha@938: double startX = minX; sascha@933: felix@2673: double area = 0d; felix@2673: // Walking along the profile. sascha@933: for (int i = 1; i < N; ++i) { sascha@938: Point2D p1 = points.get(i-1); sascha@938: Point2D p2 = points.get(i); sascha@933: sascha@938: if (p1.getY() < waterLevel && p2.getY() < waterLevel) { sascha@933: // completely under water sascha@938: if (debug) { sascha@938: log.debug("under water: " + p1 + " " + p2); sascha@938: } sascha@938: if (mode != Mode.WET) { sascha@938: startX = p1.getX(); sascha@938: mode = Mode.WET; sascha@938: } felix@2673: area += area(p1, p2, felix@2673: new Point2D.Double(p2.getX(), waterLevel), felix@2673: new Point2D.Double(p1.getX(), waterLevel)); sascha@933: continue; sascha@933: } sascha@933: felix@2673: // TODO trigger area calculation sascha@938: if (p1.getY() > waterLevel && p2.getY() > waterLevel) { sascha@938: if (debug) { sascha@938: log.debug("over water: " + p1 + " " + p2); sascha@938: } sascha@933: // completely over water sascha@933: if (mode == Mode.WET) { sascha@938: log.debug("over/wet"); sascha@933: result.add(new Line2D.Double( sascha@933: startX, waterLevel, sascha@933: p1.getX(), waterLevel)); sascha@933: } sascha@933: mode = Mode.DRY; sascha@933: continue; sascha@933: } sascha@933: felix@2673: // TODO trigger area calculation sascha@933: if (Math.abs(p1.getX() - p2.getX()) < EPSILON) { sascha@933: // vertical line sascha@933: switch (mode) { sascha@933: case WET: sascha@938: log.debug("vertical/wet"); sascha@933: mode = Mode.DRY; sascha@933: result.add(new Line2D.Double( sascha@933: startX, waterLevel, sascha@933: p1.getX(), waterLevel)); sascha@933: break; sascha@933: case DRY: sascha@938: log.debug("vertical/dry"); sascha@933: mode = Mode.WET; sascha@933: startX = p2.getX(); sascha@933: break; sascha@933: default: // UNDEF sascha@938: log.debug("vertical/undef"); sascha@933: if (p2.getY() < waterLevel) { sascha@933: mode = Mode.WET; sascha@933: startX = p2.getX(); sascha@933: } sascha@933: else { sascha@933: mode = Mode.DRY; sascha@933: } sascha@933: } sascha@933: continue; sascha@933: } sascha@933: sascha@938: // check if waterlevel directly hits the vertices; sascha@938: sascha@938: boolean p1W = Math.abs(waterLevel - p1.getY()) < EPSILON; sascha@938: boolean p2W = Math.abs(waterLevel - p2.getY()) < EPSILON; sascha@938: felix@2673: // TODO trigger area calculation sascha@938: if (p1W || p2W) { sascha@938: if (debug) { sascha@938: log.debug("water hits vertex: " + p1 + " " + p2 + " " + mode); sascha@938: } sascha@938: if (p1W && p2W) { // parallel to water -> dry sascha@938: log.debug("water hits both vertices"); sascha@938: if (mode == Mode.WET) { sascha@938: result.add(new Line2D.Double( sascha@938: startX, waterLevel, sascha@938: p1.getX(), waterLevel)); sascha@938: } sascha@938: mode = Mode.DRY; sascha@938: } sascha@938: else if (p1W) { // p1 == waterlevel sascha@938: log.debug("water hits first vertex"); sascha@938: if (p2.getY() > waterLevel) { // --> dry sascha@938: if (mode == Mode.WET) { sascha@938: result.add(new Line2D.Double( sascha@938: startX, waterLevel, sascha@938: p1.getX(), waterLevel)); sascha@938: } sascha@938: mode = Mode.DRY; sascha@938: } sascha@938: else { // --> wet sascha@938: if (mode != Mode.WET) { sascha@938: startX = p1.getX(); sascha@938: mode = Mode.WET; sascha@938: } felix@2687: area += area(p1, p2, felix@2687: new Point2D.Double(p2.getX(), waterLevel), felix@2687: new Point2D.Double(p2.getX(), waterLevel)); sascha@938: } sascha@938: } sascha@938: else { // p2 == waterlevel sascha@938: log.debug("water hits second vertex"); sascha@938: if (p1.getY() > waterLevel) { // --> wet sascha@938: if (mode != Mode.WET) { sascha@938: startX = p2.getX(); sascha@938: mode = Mode.WET; sascha@938: } sascha@938: } sascha@938: else { // --> dry sascha@938: if (mode == Mode.WET) { sascha@938: result.add(new Line2D.Double( sascha@938: startX, waterLevel, sascha@938: p2.getX(), waterLevel)); sascha@938: } sascha@938: mode = Mode.DRY; felix@2687: area += area(p1, p2, felix@2687: new Point2D.Double(p1.getX(), waterLevel), felix@2687: new Point2D.Double(p1.getX(), waterLevel)); sascha@938: } sascha@938: } sascha@938: if (debug) { sascha@938: log.debug("mode is now: " + mode); sascha@938: } sascha@938: continue; sascha@938: } sascha@938: felix@2673: // TODO trigger area calculation sascha@933: // intersection case sascha@933: double x = Linear.linear( sascha@933: waterLevel, sascha@938: p1.getY(), p2.getY(), sascha@938: p1.getX(), p2.getX()); sascha@938: sascha@938: if (debug) { sascha@938: log.debug("intersection p1:" + p1); sascha@938: log.debug("intersection p2:" + p2); sascha@938: log.debug("intersection at x: " + x); sascha@938: } sascha@933: felix@2687: // Add area of that part of intersection that is 'wet'. felix@2687: if (p1.getY() > waterLevel) { felix@2687: area += area(new Point2D.Double(x, waterLevel), felix@2687: p2, felix@2687: new Point2D.Double(p2.getX(), waterLevel), felix@2687: new Point2D.Double(x, waterLevel)); felix@2687: } felix@2687: else { felix@2687: area += area(new Point2D.Double(x, waterLevel), felix@2687: p1, felix@2687: new Point2D.Double(p1.getX(), waterLevel), felix@2687: new Point2D.Double(x, waterLevel)); felix@2687: } felix@2687: sascha@933: switch (mode) { sascha@933: case WET: sascha@938: log.debug("intersect/wet"); sascha@933: mode = Mode.DRY; sascha@933: result.add(new Line2D.Double( sascha@933: startX, waterLevel, sascha@933: x, waterLevel)); sascha@933: break; sascha@933: sascha@933: case DRY: sascha@938: log.debug("intersect/dry"); sascha@933: mode = Mode.WET; sascha@933: startX = x; sascha@933: break; sascha@933: sascha@933: default: // UNDEF sascha@938: log.debug("intersect/undef"); sascha@933: if (p2.getY() > waterLevel) { sascha@938: log.debug("intersect/undef/over"); sascha@933: mode = Mode.DRY; sascha@933: result.add(new Line2D.Double( sascha@933: p1.getX(), waterLevel, sascha@933: x, waterLevel)); sascha@933: } sascha@933: else { sascha@933: mode = Mode.WET; sascha@933: startX = x; sascha@933: } sascha@933: } // switch mode sascha@934: } // for all points p[i] and p[i-1] sascha@933: sascha@933: if (mode == Mode.WET) { sascha@933: result.add(new Line2D.Double( sascha@933: startX, waterLevel, sascha@933: maxX, waterLevel)); sascha@933: } sascha@933: felix@2674: return new ListWithArea(result, area); sascha@933: } sascha@1651: felix@2652: felix@2652: /** felix@2652: * Class holding points that form lines and the calculated length. felix@2652: */ felix@2652: public static class LineData { felix@2652: public double [][] points; felix@2652: public double width; felix@2673: public double area; felix@2673: public LineData(double[][] points, double width, double area) { felix@2652: this.points = points; felix@2652: this.width = width; felix@2673: this.area = area; felix@2652: } felix@2652: } felix@2652: felix@2652: felix@2681: /** Return length of a single line. */ felix@2681: public static double lineLength(Line2D line) { felix@2681: double xDiff = line.getX1() - line.getX2(); felix@2681: double yDiff = line.getY1() - line.getY2(); felix@2681: return Math.sqrt(xDiff*xDiff + yDiff*yDiff); felix@2681: } felix@2681: felix@2681: felix@2681: /** felix@2681: * @param points the riverbed. felix@2681: */ felix@2652: public static LineData createWaterLines( sascha@1651: List points, sascha@1651: double waterlevel sascha@1651: ) { felix@2674: ListWithArea listAndArea = fillWater(points, waterlevel); felix@2674: List lines = listAndArea.lines; sascha@1651: sascha@1651: TDoubleArrayList lxs = new TDoubleArrayList(); sascha@1651: TDoubleArrayList lys = new TDoubleArrayList(); felix@2652: double linesLength = 0.0f; sascha@1651: sascha@1651: for (Iterator iter = lines.iterator(); iter.hasNext();) { felix@2652: Line2D line = iter.next(); felix@2652: Point2D p1 = line.getP1(); felix@2652: Point2D p2 = line.getP2(); sascha@1651: lxs.add(p1.getX()); sascha@1651: lys.add(p1.getY()); sascha@1651: lxs.add(p2.getX()); sascha@1651: lys.add(p2.getY()); felix@2652: felix@2652: // Length calculation. felix@2681: linesLength += lineLength(line); felix@2652: sascha@1651: if (iter.hasNext()) { sascha@1651: lxs.add(Double.NaN); sascha@1651: lys.add(Double.NaN); sascha@1651: } sascha@1651: } sascha@1651: felix@2652: return new LineData( felix@2652: new double [][] { lxs.toNativeArray(), lys.toNativeArray() }, felix@2674: linesLength, listAndArea.area felix@2652: ); sascha@1651: } sascha@933: } sascha@933: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :