--- a/gnv-artifacts/ChangeLog	Wed Dec 23 12:20:25 2009 +0000
+++ b/gnv-artifacts/ChangeLog	Wed Dec 23 15:28:40 2009 +0000
@@ -1,3 +1,29 @@
+2009-12-23	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	* src/main/java/de/intevation/gnv/math/Interpolation2D.java: 
+	  Use local variable debug instead of asking log.isDebugEnabled() 
+	  more than once.
+
+	* src/main/java/de/intevation/gnv/math/XYColumn.java: extrapolate
+	  with boundary values.
+
+	* src/main/java/de/intevation/gnv/math/Interpolation3D.java:
+	  Implements a 3D interpolation called 'Profilschnitt' along a 
+	  track similiar to the 'Horizontaler Schnittprofil' which takes 
+	  all k layers into account. 
+
+	  At the interpolated (x, y) points columns of parameter values 
+	  from surface to ground are interpolated. To do so the four 
+	  next neighbor of that columns are figured out. Four
+	  cubic splines are fitted through these parameter values 
+	  of these neighbors. Now its possible to continuous eval
+	  the parameter on each. Every entry in the interpolated column
+	  is interpolated bilinear from the four cubic spline interpolated
+	  neighbor values at the respective depth.
+
+	  The result is stored into a double valued raster. NaN values
+	  indicates interpolation gaps.
+
 2009-12-23	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
 
 	* src/main/java/de/intevation/gnv/math/XYColumn.java: Added a method

--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Wed Dec 23 12:20:25 2009 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Wed Dec 23 15:28:40 2009 +0000
@@ -91,10 +91,12 @@
         Metrics                    metrics,
         Consumer                   consumer
     ) {
+        boolean debug = log.isDebugEnabled();
+
         int N = path.size();
         int M = points.size();
 
-        if (log.isDebugEnabled()) {
+        if (debug) {
             log.debug("Size of path: " + N);
             log.debug("Size of points: " + M);
         }
@@ -107,7 +109,7 @@
         double dxMax = buffer[0];
         double dyMax = buffer[1];
 
-        if (log.isDebugEnabled()) {
+        if (debug) {
             log.debug("buffer size: " + dxMax + " / " + dyMax);
         }
 
@@ -204,7 +206,7 @@
             }
         }
 
-        if (log.isDebugEnabled()) {
+        if (debug) {
             log.debug("interpolations: " + 
                 interpolations + " / " + missedInterpolations);
         }

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java	Wed Dec 23 15:28:40 2009 +0000
@@ -0,0 +1,223 @@
+package de.intevation.gnv.math;
+
+import java.util.List;
+import java.util.Arrays;
+import java.util.Collections;
+
+import com.vividsolutions.jts.geom.Coordinate;
+import com.vividsolutions.jts.geom.Envelope;
+
+import com.vividsolutions.jts.index.quadtree.Quadtree;
+
+import org.apache.log4j.Logger;
+
+public class Interpolation3D
+{
+    private static Logger log = Logger.getLogger(Interpolation3D.class);
+
+    public static final int DEFAULT_WIDTH  = 1024;
+    public static final int DEFAULT_HEIGHT =  768;
+
+    protected int width;
+    protected int height;
+
+    protected double [] raster;
+    protected double [] depths;
+
+    public Interpolation3D() {
+        this(DEFAULT_WIDTH, DEFAULT_HEIGHT);
+    }
+
+    public Interpolation3D(int width, int height) {
+        this.width  = width;
+        this.height = height;
+    }
+
+    public int getHeight() {
+        return height;
+    }
+
+    public int getWidth() {
+        return width;
+    }
+
+    public double [] getRaster() {
+        return raster;
+    }
+
+    public double [] getDepths() {
+        return depths;
+    }
+
+    public boolean interpolate(
+        List<? extends Coordinate> path,
+        List<? extends XYColumn>   points,
+        double                     from,
+        double                     to,
+        Metrics                    metrics,
+        XYDepth                    xyDepth
+    ) {
+        boolean debug = log.isDebugEnabled();
+
+        int N = path.size();
+        int M = points.size();
+
+        if (debug) {
+            log.debug("Size of path: " + N);
+            log.debug("Size of points: " + M);
+        }
+
+        if (M < 1 || N < 2) { // nothing to do
+            return false;
+        }
+
+        double [] buffer = Interpolation2D.calculateBuffer(points);
+        double dxMax = buffer[0];
+        double dyMax = buffer[1];
+
+        if (debug) {
+            log.debug("buffer size: " + dxMax + " / " + dyMax);
+        }
+
+        // put into spatial index to speed up finding neighbors.
+        Quadtree spatialIndex = new Quadtree();
+
+        for (int i = 0; i < M; ++i) {
+            Point2d p = points.get(i);
+            spatialIndex.insert(p.envelope(), p);
+        }
+
+        LinearToMap linearToMap = new LinearToMap(
+            path, from, to, metrics);
+
+        double [] depths = new double[width];
+        Arrays.fill(depths, Double.NaN);
+
+        double cellWidth = (to - from)/width;
+
+        double maxDepth = Double.MAX_VALUE;
+
+        int i = 0;
+        Coordinate center = new Coordinate();
+        for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
+            if (linearToMap.locate(p, center)) {
+                double depth = xyDepth.depth(center);
+                depths[i] = depth;
+                if (depth < maxDepth) {
+                    maxDepth = depth;
+                }
+            }
+        }
+
+        if (maxDepth == Double.MAX_VALUE) {
+            log.warn("no depth found -> no interpolation");
+            return false;
+        }
+
+        if (debug) {
+            log.debug("max depth found: " + maxDepth);
+        }
+
+        double cellHeight = Math.abs(maxDepth)/height;
+
+        if (debug) {
+            log.debug("cell size: " + cellWidth + " x " + cellHeight);
+        }
+
+        double [] raster = new double[width*height];
+        Arrays.fill(raster, Double.NaN);
+
+        Envelope queryBuffer  = new Envelope();
+        XYColumn [] neighbors = new XYColumn[4];
+
+        i = 0;
+        for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
+            double depth = depths[i];
+            if (Double.isNaN(depth)) {
+                continue;
+            }
+            linearToMap.locate(p, center);
+
+            queryBuffer.init(
+                center.x - dxMax, center.x + dxMax, 
+                center.y - dyMax, center.y + dyMax);
+
+            List potential = spatialIndex.query(queryBuffer);
+
+            L1Comparator invL1 = new L1Comparator(center);
+            Collections.sort(potential, invL1);
+
+            neighbors[0] = neighbors[1] = 
+            neighbors[2] = neighbors[3] = null;
+
+            /* bit code of neighbors
+               0---1
+               | x |
+               2---3
+            */
+
+            int mask = 0;
+            // reversed order is essential here!
+            for (int j = potential.size()-1; j >= 0; --j) {
+                XYColumn n = (XYColumn)potential.get(j);
+                int code = n.x > center.x ? 1 : 0;
+                if (n.y > center.y) code |= 2;
+                neighbors[code] = n;
+                mask |= 1 << code;
+            }
+
+            int numNeighbors = Integer.bitCount(mask);
+
+            // only interpolate if we have all four neighbors
+            // and we do not have any gaps.
+            if (numNeighbors == 4
+            && !neighbors[0].hasIGap(neighbors[1])
+            && !neighbors[1].hasJGap(neighbors[3])
+            && !neighbors[3].hasIGap(neighbors[2])
+            && !neighbors[2].hasJGap(neighbors[0])
+            &&  neighbors[0].prepare(xyDepth)
+            &&  neighbors[1].prepare(xyDepth)
+            &&  neighbors[2].prepare(xyDepth)
+            &&  neighbors[3].prepare(xyDepth)
+            ) {
+                double y1 = Interpolation2D.interpolate(
+                    neighbors[0].x, neighbors[0].y,
+                    neighbors[1].x, neighbors[1].y,
+                    center.x);
+                double y2 = Interpolation2D.interpolate(
+                    neighbors[2].x, neighbors[2].y,
+                    neighbors[3].x, neighbors[3].y,
+                    center.x);
+                int j = i;
+                for (double z = -cellHeight*0.5; 
+                    j < raster.length; z -= cellHeight, j += width) {
+
+                    double v0 = neighbors[0].value(z);
+                    double v1 = neighbors[1].value(z);
+                    double v2 = neighbors[2].value(z);
+                    double v3 = neighbors[3].value(z);
+
+                    double z1 = Interpolation2D.interpolate(
+                        neighbors[0].x, v0,
+                        neighbors[1].x, v1,
+                        center.x);
+                    double z2 = Interpolation2D.interpolate(
+                        neighbors[2].x, v2,
+                        neighbors[3].x, v3,
+                        center.x);
+                    double value = Interpolation2D.interpolate(
+                        y1, z1,
+                        y2, z2,
+                        center.y);
+                    raster[j] = value;
+                }
+            } // down the x/y column
+        } // all along the track
+
+        this.depths = depths;
+        this.raster = raster;
+
+        return true;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:

--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/XYColumn.java	Wed Dec 23 12:20:25 2009 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/XYColumn.java	Wed Dec 23 15:28:40 2009 +0000
@@ -24,7 +24,7 @@
 extends      Point2d
 implements   UnivariateRealFunction
 {
-    private static Logger log = Logger.getLogger(Interpolation2D.class);
+    private static Logger log = Logger.getLogger(XYColumn.class);
 
     protected List<HeightValue> values;
 
@@ -50,6 +50,11 @@
     public double value(double depth) {
         try {
             if (curve != null) {
+                HeightValue h = values.get(0);
+                // extrapolate beyond boundaries by repeating
+                if (depth > h.z) return h.v;
+                h = values.get(values.size()-1);
+                if (depth < h.z) return h.v;
                 return curve.value(depth);
             }
         }
@@ -60,7 +65,7 @@
         return Double.NaN;
     }
 
-    public void prepare(XYDepth xyDepth) {
+    public boolean prepare(XYDepth xyDepth) {
         int N = values.size();
         if (curve == null && N > 0) {
             if (N == 1) {
@@ -103,13 +108,16 @@
                     }
                     catch (MathException me) {
                         log.error("interpolation failed", me);
+                        return false;
                     }
                 }
             }
         }
         else {
             log.error("no points for interpolation");
+            return false;
         }
+        return true;
     }
 
     protected UnivariateRealInterpolator getInterpolator() {