Mercurial > dive4elements > gnv-client

comparison gnv-artifacts/src/main/java/de/intevation/gnv/math/XYColumn.java @ 657:af3f56758f59

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged gnv-artifacts/0.5
author Thomas Arendsen Hein <thomas@intevation.de>
date Fri, 28 Sep 2012 12:13:53 +0200
parents bc5901bb4525
children b1f5f2a8840f
comparison
equal deleted inserted replaced
590:5f5f273c8566 657:af3f56758f59
1 package de.intevation.gnv.math;
2
3 import java.util.ArrayList;
4 import java.util.List;
5 import java.util.Collections;
6
7 import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
8 import org.apache.commons.math.analysis.interpolation.UnivariateRealInterpolator;
9
10 import org.apache.commons.math.analysis.UnivariateRealFunction;
11
12 import org.apache.commons.math.MathException;
13 import org.apache.commons.math.FunctionEvaluationException;
14
15 import org.apache.log4j.Logger;
16
17 /**
18 * @author Ingo Weinzierl <ingo.weinzierl@intevation.de>
19 * @author Sascha L. Teichmann <sascha.teichmann@intevation.de>
20 */
21 public class XYColumn
22 extends Point2d
23 implements UnivariateRealFunction
24 {
25 private static Logger log = Logger.getLogger(XYColumn.class);
26
27 protected List<HeightValue> values;
28
29 protected transient UnivariateRealFunction curve;
30
31 public XYColumn() {
32 values = new ArrayList<HeightValue>();
33 }
34
35 public XYColumn(double x, double y, int i, int j) {
36 super(x, y, i, j);
37 values = new ArrayList<HeightValue>();
38 }
39
40 public void add(HeightValue value) {
41 values.add(value);
42 }
43
44 public List<HeightValue> getValues() {
45 return values;
46 }
47
48 public double value(double depth) {
49 try {
50 if (curve != null) {
51 HeightValue h = values.get(0);
52 // extrapolate beyond boundaries by repeating
53 if (depth > h.z) return h.v;
54 h = values.get(values.size()-1);
55 if (depth < h.z) return h.v;
56 return curve.value(depth);
57 }
58 }
59 catch (FunctionEvaluationException fee) {
60 log.error("evaluation failed", fee);
61 }
62
63 return Double.NaN;
64 }
65
66 public boolean prepare(XYDepth xyDepth) {
67 if (curve == null) {
68 int N = values.size();
69 if (N == 0) {
70 log.error("no points for interpolation");
71 return false;
72 }
73
74 if (N == 1) {
75 // only one value -> constant function
76 curve = new ConstantFunction(values.get(0).v);
77 }
78 else { // more than on value
79 double depth = xyDepth.depth(this);
80 Collections.sort(values, HeightValue.INV_Z_COMPARATOR);
81
82 // if there is no value at 0 repeat first value
83 HeightValue first = values.get(0);
84 if (first.z < 0d) {
85 values.add(0, new HeightValue(0d, first.v, first.k-1));
86 ++N;
87 }
88
89 // if there is no value at depth repeat last value
90 HeightValue last = values.get(N-1);
91 if (last.z > depth) {
92 values.add(new HeightValue(depth, last.v, last.k+1));
93 ++N;
94 }
95 if (N < 3) { // interpolate linear
96 first = values.get(0);
97 last = values.get(N-1);
98 curve = new LinearFunction.Univariate(
99 first.z, first.v,
100 last.z, last.v);
101 }
102 else { // higher degree interpolation
103 double [] z = new double[N];
104 double [] v = new double[N];
105 for (int i = 0; i < N; ++i) {
106 HeightValue h = values.get(N-1-i);
107 z[i] = h.z;
108 v[i] = h.v;
109 }
110 try {
111 curve = getInterpolator().interpolate(z, v);
112 }
113 catch (MathException me) {
114 log.error("interpolation failed", me);
115 return false;
116 }
117 }
118 }
119 }
120 return true;
121 }
122
123 protected UnivariateRealInterpolator getInterpolator() {
124 return new SplineInterpolator();
125 }
126 }
127 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:
http://dive4elements.wald.intevation.org