ingo@427: package de.intevation.gnv.math;
ingo@427: 
ingo@427: import java.util.ArrayList;
sascha@779: import java.util.Collections;
ingo@427: import java.util.List;
sascha@779: 
sascha@779: import org.apache.commons.math.FunctionEvaluationException;
sascha@779: import org.apache.commons.math.MathException;
sascha@779: 
sascha@779: import org.apache.commons.math.analysis.UnivariateRealFunction;
sascha@433: 
sascha@433: import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
sascha@433: import org.apache.commons.math.analysis.interpolation.UnivariateRealInterpolator;
sascha@433: 
sascha@433: import org.apache.log4j.Logger;
ingo@427: 
ingo@427: /**
sascha@807:  * A column of discrete attributed height values located at a point.
sascha@807:  * Values between the discrete height values are interpolated.
sascha@807:  *
sascha@780:  * @author <a href="mailto:ingo.weinzierl@intevation.de">Ingo Weinzierl</a>
sascha@780:  * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
ingo@427:  */
ingo@429: public class XYColumn
sascha@431: extends      Point2d
sascha@433: implements   UnivariateRealFunction
ingo@429: {
sascha@434:     private static Logger log = Logger.getLogger(XYColumn.class);
sascha@433: 
sascha@807:     /**
sascha@807:      * The list of discrete height values.
sascha@807:      */
sascha@432:     protected List<HeightValue> values;
ingo@427: 
sascha@807:     /**
sascha@807:      * The curve used to interpolate the points between the
sascha@807:      * discrete height values.
sascha@807:      */
sascha@433:     protected transient UnivariateRealFunction curve;
sascha@433: 
sascha@807:     /**
sascha@807:      * Default constructor.
sascha@807:      */
sascha@431:     public XYColumn() {
sascha@433:         values = new ArrayList<HeightValue>();
sascha@431:     }
ingo@427: 
sascha@807:     /**
sascha@807:      * Constructor to create an XYColumn with a given (x, y) coordinate
sascha@807:      * and an (i, j) index tuple.
sascha@807:      * @param x The x coordinate.
sascha@807:      * @param y The y coordinate.
sascha@807:      * @param i The i component.
sascha@807:      * @param j The j component.
sascha@807:      */
ingo@427:     public XYColumn(double x, double y, int i, int j) {
sascha@431:         super(x, y, i, j);
sascha@433:         values = new ArrayList<HeightValue>();
ingo@427:     }
ingo@427: 
sascha@807:     /**
sascha@807:      * Adds a given height value to the list of height values.
sascha@807:      * @param value The height value.
sascha@807:      */
ingo@427:     public void add(HeightValue value) {
ingo@427:         values.add(value);
ingo@427:     }
ingo@427: 
sascha@807:     /**
sascha@807:      * Returns the list of height values.
sascha@807:      * @return The list of height values.
sascha@807:      */
sascha@432:     public List<HeightValue> getValues() {
sascha@432:         return values;
ingo@427:     }
sascha@433: 
sascha@433:     public double value(double depth) {
sascha@433:         try {
sascha@433:             if (curve != null) {
sascha@434:                 HeightValue h = values.get(0);
sascha@434:                 // extrapolate beyond boundaries by repeating
sascha@434:                 if (depth > h.z) return h.v;
sascha@434:                 h = values.get(values.size()-1);
sascha@434:                 if (depth < h.z) return h.v;
sascha@433:                 return curve.value(depth);
sascha@433:             }
sascha@433:         }
sascha@433:         catch (FunctionEvaluationException fee) {
sascha@433:             log.error("evaluation failed", fee);
sascha@433:         }
sascha@433: 
sascha@433:         return Double.NaN;
sascha@433:     }
sascha@433: 
sascha@807:     /**
sascha@807:      * Prepares this XYColumn to be queried. A given XYDepth
sascha@807:      * object is used to fill the values to a certain depth.
sascha@807:      * @param xyDepth To figure out the depth a the cordinate.
sascha@807:      * @return true if preparation succeeds else false.
sascha@807:      */
sascha@434:     public boolean prepare(XYDepth xyDepth) {
sascha@445:         if (curve == null) {
sascha@451:             int N = values.size();
sascha@445:             if (N == 0) {
sascha@445:                 log.error("no points for interpolation");
sascha@445:                 return false;
sascha@445:             }
sascha@445: 
sascha@433:             if (N == 1) {
sascha@433:                 // only one value -> constant function
sascha@433:                 curve = new ConstantFunction(values.get(0).v);
sascha@433:             }
sascha@433:             else { // more than on value
sascha@433:                 double depth = xyDepth.depth(this);
sascha@433:                 Collections.sort(values, HeightValue.INV_Z_COMPARATOR);
sascha@433: 
sascha@433:                 // if there is no value at 0 repeat first value
sascha@433:                 HeightValue first = values.get(0);
sascha@433:                 if (first.z < 0d) {
sascha@448:                     values.add(0, new HeightValue(0d, first.v, first.k-1));
sascha@445:                     ++N;
sascha@433:                 }
sascha@433: 
sascha@433:                 // if there is no value at depth repeat last value
sascha@433:                 HeightValue last = values.get(N-1);
sascha@433:                 if (last.z > depth) {
sascha@448:                     values.add(new HeightValue(depth, last.v, last.k+1));
sascha@445:                     ++N;
sascha@433:                 }
sascha@433:                 if (N < 3) { // interpolate linear
sascha@433:                     first = values.get(0);
sascha@433:                     last  = values.get(N-1);
sascha@433:                     curve = new LinearFunction.Univariate(
sascha@433:                         first.z, first.v,
sascha@433:                         last.z,  last.v);
sascha@433:                 }
sascha@433:                 else { // higher degree interpolation
sascha@433:                     double [] z = new double[N];
sascha@433:                     double [] v = new double[N];
sascha@433:                     for (int i = 0; i < N; ++i) {
sascha@445:                         HeightValue h = values.get(N-1-i);
sascha@433:                         z[i] = h.z;
sascha@433:                         v[i] = h.v;
sascha@433:                     }
sascha@433:                     try {
sascha@433:                         curve = getInterpolator().interpolate(z, v);
sascha@433:                     }
sascha@433:                     catch (MathException me) {
sascha@433:                         log.error("interpolation failed", me);
sascha@434:                         return false;
sascha@433:                     }
sascha@433:                 }
sascha@433:             }
sascha@433:         }
sascha@434:         return true;
sascha@433:     }
sascha@433: 
sascha@807:     /**
sascha@807:      * Returns the interpolator used to interpolate the values between
sascha@807:      * the discrete height values. This class returns an instance of
sascha@807:      * {@link org.apache.commons.math.analysis.interpolation.SplineInterpolator}.
sascha@807:      * Override this if you want to use another kind of interpolation.
sascha@807:      * @return The interpolator to be used.
sascha@807:      */
sascha@433:     protected UnivariateRealInterpolator getInterpolator() {
sascha@433:         return new SplineInterpolator();
sascha@433:     }
ingo@427: }
sascha@807: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :