--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/LinearInterpolated.java	Sun Jun 02 17:52:53 2013 +0200
@@ -0,0 +1,208 @@
+/* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
+ * Software engineering by Intevation GmbH
+ *
+ * This file is Free Software under the GNU AGPL (>=v3)
+ * and comes with ABSOLUTELY NO WARRANTY! Check out the
+ * documentation coming with Dive4Elements River for details.
+ */
+
+package org.dive4elements.river.artifacts.model;
+
+import gnu.trove.TDoubleArrayList;
+
+import java.io.Serializable;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Set;
+import java.util.TreeSet;
+
+import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
+
+import org.dive4elements.river.artifacts.math.Linear;
+
+import org.dive4elements.river.utils.EpsilonComparator;
+
+public class LinearInterpolated
+implements   Serializable
+{
+    public static final double EPSILON = 1e-5;
+    public static final double MULTIPLE_STD_DEV = 4d;
+
+    public static final EpsilonComparator CMP = new EpsilonComparator(EPSILON);
+
+    private TDoubleArrayList xs;
+    private TDoubleArrayList ys;
+
+    private List<Range> gaps;
+
+    public interface Operator {
+        double evaluate(double y1, double y2);
+    } // interface Operator
+
+    public static final Operator SUB = new Operator() {
+        @Override
+        public double evaluate(double y1, double y2) {
+            return y1 - y2;
+        }
+    };
+
+    public static final Operator MAX = new Operator() {
+        @Override
+        public double evaluate(double y1, double y2) {
+            return Math.max(y1, y2);
+        }
+    };
+
+    public LinearInterpolated() {
+        xs = new TDoubleArrayList();
+        ys = new TDoubleArrayList();
+    }
+
+    public LinearInterpolated(int capacity) {
+        xs = new TDoubleArrayList(capacity);
+        ys = new TDoubleArrayList(capacity);
+    }
+
+    public void add(double x, double y) {
+        xs.add(x);
+        ys.add(y);
+    }
+
+    public int size() {
+        return xs.size();
+    }
+
+    public void pointsInRange(double from, double to, Set<Double> points) {
+        if (from > to) {
+            double t = from;
+            from = to;
+            to = t;
+        }
+        for (int i = 0, S = size(); i < S; ++i) {
+            double x = xs.getQuick(i);
+            if (x >= from && x <= to) {
+                points.add(x);
+            }
+        }
+    }
+
+    public boolean inGap(double x) {
+        if (gaps != null) {
+            for (Range gap: gaps) {
+                if (gap.inside(x)) {
+                    return true;
+                }
+            }
+        }
+        return false;
+    }
+
+    public void detectGaps(double threshold) {
+        List<Range> gabs = new ArrayList<Range>();
+        for (int i = 1, S = size(); i < S; ++i) {
+            double x0 = xs.getQuick(i-1);
+            double x1 = xs.getQuick(i);
+            double minX, maxX;
+            if (x0 < x1) { minX = x0; maxX = x1; }
+            else         { minX = x1; maxX = x0; }
+            double diff = maxX - minX - 2d*EPSILON;
+            if (diff > threshold) {
+                gaps.add(new Range(minX+EPSILON, maxX-EPSILON));
+            }
+        }
+        this.gaps = gaps.isEmpty() ? null : gabs;
+    }
+
+    public void resetGaps() {
+        gaps = null;
+    }
+
+    public double guessGapThreshold() {
+        return guessGapThreshold(MULTIPLE_STD_DEV);
+    }
+
+    public double guessGapThreshold(double scale) {
+        int S = size();
+        if (S < 5) { // Too less points.
+            return Double.MAX_VALUE;
+        }
+
+        StandardDeviation s = new StandardDeviation();
+
+        for (int i = 1; i < S; ++i) {
+            double diff = Math.abs(xs.getQuick(i-1) - xs.getQuick(i));
+            s.increment(diff);
+        }
+
+        return scale*s.getResult();
+    }
+
+    public double value(double x) {
+        for (int i = 0, S = size(); i < S; ++i) {
+            double x1 = xs.getQuick(i);
+            if (Math.abs(x1 - x) < EPSILON) {
+                return ys.getQuick(i);
+            }
+            if (i > 0) {
+                double x0 = xs.getQuick(i-1);
+                double minX, maxX;
+                if (x0 < x1) { minX = x0; maxX = x1; }
+                else         { minX = x1; maxX = x0; }
+                if (x > minX && x < maxX) {
+                    return Linear.linear(
+                        x,
+                        x0, x1,
+                        ys.getQuick(i-1), ys.getQuick(i));
+                }
+            }
+        }
+        return Double.NaN;
+    }
+
+    public LinearInterpolated sub(
+        LinearInterpolated other, 
+        double             from,
+        double             to
+    ) {
+        return apply(SUB, other, from, to);
+    }
+
+    public LinearInterpolated max(
+        LinearInterpolated other, 
+        double             from,
+        double             to
+    ) {
+        return apply(MAX, other, from, to);
+    }
+
+    public LinearInterpolated apply(
+        Operator           operator,
+        LinearInterpolated other,
+        double             from,
+        double             to
+    ) {
+        Set<Double> points = new TreeSet<Double>(CMP);
+        points.add(from);
+        points.add(to);
+
+        this .pointsInRange(from, to, points);
+        other.pointsInRange(from, to, points);
+
+        LinearInterpolated result = new LinearInterpolated();
+
+        for (double x: points) {
+            if (!inGap(x) && !other.inGap(x)) {
+                double y1 = this .value(x);
+                double y2 = other.value(x);
+                double y  = operator.evaluate(y1, y2);
+                if (!Double.isNaN(y)) {
+                    result.add(x, y);
+                }
+            }
+        }
+
+        return result;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :