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<Line2D> 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<Line2D> lines;
felix@2674:         public double area;
felix@2674:         public ListWithArea(List<Line2D> 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<Point2D> 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<Line2D> 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<Point2D> points,
sascha@1651:         double        waterlevel
sascha@1651:     ) {
felix@2674:         ListWithArea listAndArea = fillWater(points, waterlevel);
felix@2674:         List<Line2D> 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<Line2D> 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 :