Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java @ 1674:b5209452f6bb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
flys/issue31: Simplified the back jump correction code a lot. flys-artifacts/trunk@2888 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 04 Oct 2011 16:36:19 +0000
parents 0b6dac664bbb
children 6e840e213fdf
comparison
equal deleted inserted replaced
1673:68260e38029a 1674:b5209452f6bb
106 if (corrected == null) { 106 if (corrected == null) {
107 // lazy cloning 107 // lazy cloning
108 ws = corrected = (double [])ws.clone(); 108 ws = corrected = (double [])ws.clone();
109 } 109 }
110 110
111 int back = i-2; 111 double above = aboveWaterKM(km, ws, i);
112 double distance = Math.abs(km[i] - km[i-1]); 112
113 double rest = 0.0; 113 if (Double.isNaN(above)) { // run over start km
114 double ikm = 0.0; 114 // fill all previous
115
116 while (back > -1) {
117 if (ws[back] < ws[i]) { // under water
118 // continue scanning backwards
119 distance += Math.abs(km[back+1] - km[back]);
120 --back;
121 continue;
122 }
123 if (ws[back] > ws[i]) { // over water
124 // need to calculate intersection
125 log.debug("intersection case");
126 double m = (km[back+1]-km[back])/(ws[back+1]-ws[back]);
127 double b = km[back]-ws[back]*m;
128 ikm = m*ws[i] + b;
129 distance += Math.abs(ikm - km[back+1]);
130 rest = Math.abs(km[back] - ikm);
131 }
132 else {
133 // equals height at ws[back]
134 log.debug("exact height match");
135 distance += Math.abs(km[back+1] - km[back]);
136 ikm = km[back];
137 }
138 break;
139 }
140
141 if (back <= 0) {
142 //log.debug("run over left border");
143 // fill with same as ws[i]
144 for (int j = 0; j < i; ++j) { 115 for (int j = 0; j < i; ++j) {
145 ws[j] = ws[i]; 116 ws[j] = ws[i];
146 } 117 }
147 continue; 118 continue;
148 } 119 }
149 120
121 double distance = Math.abs(km[i] - above);
122
150 double quarterDistance = 0.25*distance; 123 double quarterDistance = 0.25*distance;
151 124
152 // Now further seek back for the max height 125 double start = above - quarterDistance;
153 126
154 double restDistance = Math.max(0.0, quarterDistance - rest); 127 double startHeight = DoubleUtil.interpolateSorted(km, ws, start);
155 --back; 128
156 129 if (Double.isNaN(startHeight)) {
157 double mkmw = ws[i]; 130 // run over start km
158 double mkm = km[0]; 131 startHeight = ws[0];
159 132 }
160 while (back > -1) { 133
161 double d = Math.abs(km[back+1] - km[back]); 134 double between = above + quarterDistance;
162 restDistance -= d; 135
163 if (restDistance > 0.0) { 136 double betweenHeight = ws[i] + 0.25*(startHeight - ws[i]);
164 --back; 137
165 continue; 138 double [] x = { start, between, km[i] };
166 } 139 double [] y = { startHeight, betweenHeight, ws[i] };
167 if (restDistance < 0.0) {
168 // need to calculate intersection
169 if (km[back] == km[back+1]) { // should not happen
170 mkm = km[back];
171 mkmw = 0.5*(ws[back] + ws[back+1]);
172 }
173 else {
174 double m = (ws[back+1]-ws[back])/(km[back+1]-km[back]);
175 double b = ws[back] - km[back]*m;
176 mkm = km[back] + restDistance;
177 mkmw = m*mkm + b;
178 }
179 }
180 else {
181 // exact match
182 mkm = km[back];
183 mkmw = ws[back];
184 }
185 break;
186 }
187
188 double factor = back >= 0 && Math.abs(restDistance) < 1e-4
189 ? 1.0
190 : 1.0 - Math.min(1, Math.max(0, restDistance/quarterDistance));
191
192 double ikmw = factor*0.25*(mkmw-ws[i]) + ws[i];
193
194 double end = ikm + quarterDistance*factor;
195
196 double [] x = { mkm, ikm, end };
197 double [] y = { mkmw, ikmw, ws[i] };
198
199 if (interpolator == null) {
200 interpolator = new SplineInterpolator();
201 }
202 140
203 if (log.isDebugEnabled()) { 141 if (log.isDebugEnabled()) {
204 for (int j = 0; j < x.length; ++j) { 142 for (int j = 0; j < x.length; ++j) {
205 log.debug(" " + x[j] + " -> " + y[j]); 143 log.debug(" " + x[j] + " -> " + y[j]);
206 } 144 }
145 }
146
147 if (interpolator == null) {
148 interpolator = new SplineInterpolator();
207 } 149 }
208 150
209 PolynomialSplineFunction spline; 151 PolynomialSplineFunction spline;
210 152
211 try { 153 try {
224 for (int j = 0; j < x.length; ++j) { 166 for (int j = 0; j < x.length; ++j) {
225 log.debug(x[j] + " " + y[j] + " " + spline.value(x[j])); 167 log.debug(x[j] + " " + y[j] + " " + spline.value(x[j]));
226 } 168 }
227 } 169 }
228 170
229 for (back = Math.max(back, 0); 171 for (int j = i-1; j > 0 && ws[j] < startHeight; --j) {
230 back < i && km[back] < end; 172 ws[j] = spline.value(km[j]);
231 ++back
232 ) {
233 // to 3/4 point fill with spline values
234 ws[back] = spline.value(km[back]);
235 }
236 while (back < i) {
237 // fill the rest with ws[i]
238 ws[back++] = ws[i];
239 } 173 }
240 } 174 }
241 catch (ArgumentOutsideDomainException aode) { 175 catch (ArgumentOutsideDomainException aode) {
242 // TODO: I18N 176 // TODO: I18N
243 errors.addProblem("spline interpolation failed."); 177 errors.addProblem("spline interpolation failed.");
245 } 179 }
246 } // for all km 180 } // for all km
247 181
248 return !backjumps.isEmpty(); 182 return !backjumps.isEmpty();
249 } 183 }
184
185
186 protected static double aboveWaterKM(
187 double [] km,
188 double [] ws,
189 int wIndex
190 ) {
191 double w = ws[wIndex];
192
193 while (--wIndex >= 0) {
194 // still under water
195 if (ws[wIndex] < w) continue;
196
197 if (ws[wIndex] > w) {
198 // f(ws[wIndex]) = km[wIndex]
199 // f(ws[wIndex+1]) = km[wIndex+1]
200 return Linear.linear(
201 w,
202 ws[wIndex], ws[wIndex+1],
203 km[wIndex], km[wIndex+1]);
204 }
205 else {
206 return km[wIndex];
207 }
208 }
209
210 return Double.NaN;
211 }
250 } 212 }
251 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : 213 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org