Mercurial > dive4elements > river
comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java @ 343:f165c7d5d6db
Implementation of the "Ruecksprungkorrektur" to be done in "W fuer angepassten Abflusslaengschnitt".
flys-artifacts/trunk@1743 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Fri, 22 Apr 2011 09:57:27 +0000 |
parents | |
children | 663aa18bee7f |
comparison
equal
deleted
inserted
replaced
342:f72c63713099 | 343:f165c7d5d6db |
---|---|
1 package de.intevation.flys.artifacts.math; | |
2 | |
3 import java.util.ArrayList; | |
4 import java.util.List; | |
5 | |
6 import java.io.Serializable; | |
7 | |
8 import org.apache.commons.math.analysis.interpolation.SplineInterpolator; | |
9 | |
10 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; | |
11 | |
12 import org.apache.commons.math.ArgumentOutsideDomainException; | |
13 | |
14 import org.apache.log4j.Logger; | |
15 | |
16 public class BackJumpCorrector | |
17 implements Serializable | |
18 { | |
19 private static Logger log = Logger.getLogger(BackJumpCorrector.class); | |
20 | |
21 protected ArrayList<Double> backjumps; | |
22 | |
23 protected double [] corrected; | |
24 | |
25 public BackJumpCorrector() { | |
26 backjumps = new ArrayList<Double>(); | |
27 } | |
28 | |
29 public boolean hasBackJumps() { | |
30 return !backjumps.isEmpty(); | |
31 } | |
32 | |
33 public List<Double> getBackJumps() { | |
34 return backjumps; | |
35 } | |
36 | |
37 public double [] getCorrected() { | |
38 return corrected; | |
39 } | |
40 | |
41 public boolean doCorrection(double [] km, double [] ws) { | |
42 int N = km.length; | |
43 | |
44 if (N != ws.length) { | |
45 throw new IllegalArgumentException("km.length != ws.length"); | |
46 } | |
47 | |
48 if (N < 2) { | |
49 return false; | |
50 } | |
51 | |
52 SplineInterpolator interpolator = null; | |
53 | |
54 for (int i = 1; i < N; ++i) { | |
55 if (ws[i] <= ws[i-1]) { | |
56 // no back jump | |
57 continue; | |
58 } | |
59 backjumps.add(km[i]); | |
60 | |
61 if (corrected == null) { | |
62 // lazy cloning | |
63 ws = corrected = (double [])ws.clone(); | |
64 } | |
65 | |
66 int back = i-2; | |
67 double distance = Math.abs(km[i] - km[i-1]); | |
68 double rest = 0.0; | |
69 double ikm = 0.0; | |
70 | |
71 while (back > -1) { | |
72 if (ws[back] < ws[i]) { // under water | |
73 // continue scanning backwards | |
74 distance += Math.abs(km[back+1] - km[back]); | |
75 --back; | |
76 continue; | |
77 } | |
78 if (ws[back] > ws[i]) { // over water | |
79 // need to calculate intersection | |
80 log.debug("intersection case"); | |
81 double m = (km[back+1]-km[back])/(ws[back+1]-ws[back]); | |
82 double b = km[back]-ws[back]*m; | |
83 ikm = m*ws[i] + b; | |
84 distance += Math.abs(ikm - km[back+1]); | |
85 rest = Math.abs(km[back] - ikm); | |
86 } | |
87 else { | |
88 // equals height at ws[back] | |
89 log.debug("exact height match"); | |
90 distance += Math.abs(km[back+1] - km[back]); | |
91 ikm = km[back]; | |
92 } | |
93 break; | |
94 } | |
95 | |
96 if (back < 0) { | |
97 log.debug("run over left border"); | |
98 // fill with same as ws[i] | |
99 for (int j = 0; j < i; ++j) { | |
100 ws[j] = ws[i]; | |
101 } | |
102 continue; | |
103 } | |
104 | |
105 double quarterDistance = 0.25*distance; | |
106 | |
107 // Now further seek back for the max height | |
108 | |
109 double restDistance = Math.max(0.0, quarterDistance - rest); | |
110 --back; | |
111 | |
112 double mkmw = ws[i]; | |
113 double mkm = km[0]; | |
114 | |
115 while (back > -1) { | |
116 double d = Math.abs(km[back+1] - km[back]); | |
117 restDistance -= d; | |
118 if (restDistance > 0.0) { | |
119 --back; | |
120 continue; | |
121 } | |
122 if (restDistance < 0.0) { | |
123 // need to calculate intersection | |
124 if (km[back] == km[back+1]) { // should not happen | |
125 mkm = km[back]; | |
126 mkmw = 0.5*(ws[back] + ws[back+1]); | |
127 } | |
128 else { | |
129 double m = (ws[back+1]-ws[back])/(km[back+1]-km[back]); | |
130 double b = ws[back] - km[back]*m; | |
131 mkm = km[back] + restDistance; | |
132 mkmw = m*mkm + b; | |
133 } | |
134 } | |
135 else { | |
136 // exact match | |
137 mkm = km[back]; | |
138 mkmw = ws[back]; | |
139 } | |
140 break; | |
141 } | |
142 | |
143 double factor = back >= 0 && Math.abs(restDistance) < 1e-4 | |
144 ? 1.0 | |
145 : 1.0 - Math.min(1, Math.max(0, restDistance/quarterDistance)); | |
146 | |
147 double ikmw = factor*0.25*(mkmw-ws[i]) + ws[i]; | |
148 | |
149 double end = ikm + quarterDistance*factor; | |
150 | |
151 double [] x = { mkm, ikm, end }; | |
152 double [] y = { mkmw, ikmw, ws[i] }; | |
153 | |
154 if (interpolator == null) { | |
155 interpolator = new SplineInterpolator(); | |
156 } | |
157 | |
158 PolynomialSplineFunction spline = interpolator.interpolate(x, y); | |
159 | |
160 try { | |
161 if (log.isDebugEnabled()) { | |
162 log.debug("spline points:"); | |
163 for (int j = 0; j < x.length; ++j) { | |
164 log.debug(x[j] + " " + y[j] + " " + spline.value(x[j])); | |
165 } | |
166 } | |
167 | |
168 for (back = Math.max(back, 0); | |
169 back < i && km[back] < end; | |
170 ++back | |
171 ) { | |
172 // to 3/4 point fill with spline values | |
173 ws[back] = spline.value(km[back]); | |
174 } | |
175 while (back < i) { | |
176 // fill the rest with ws[i] | |
177 ws[back++] = ws[i]; | |
178 } | |
179 } | |
180 catch (ArgumentOutsideDomainException aode) { | |
181 log.error("spline interpolation failed", aode); | |
182 } | |
183 } // for all km | |
184 | |
185 return !backjumps.isEmpty(); | |
186 } | |
187 } | |
188 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : |