Mercurial > dive4elements > river
view flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java @ 625:c0c60a611fca
Introduce model to store q values of WST columns efficiently.
flys-artifacts/trunk@1983 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Mon, 23 May 2011 22:33:37 +0000 |
parents | 663aa18bee7f |
children | 2dbbb5be30a1 |
line wrap: on
line source
package de.intevation.flys.artifacts.math; import java.util.ArrayList; import java.util.List; import java.io.Serializable; import org.apache.commons.math.analysis.interpolation.SplineInterpolator; import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; import org.apache.commons.math.ArgumentOutsideDomainException; import org.apache.log4j.Logger; public class BackJumpCorrector implements Serializable { private static Logger log = Logger.getLogger(BackJumpCorrector.class); protected ArrayList<Double> backjumps; protected double [] corrected; public BackJumpCorrector() { backjumps = new ArrayList<Double>(); } public boolean hasBackJumps() { return !backjumps.isEmpty(); } public List<Double> getBackJumps() { return backjumps; } public double [] getCorrected() { return corrected; } private static final boolean isIncreasing(double [] ws) { int inc = 0; int dec = 0; int sweet = (ws.length-1)/2; for (int i = 1; i < ws.length; ++i) { if (ws[i] > ws[i-1]) { if (++inc > sweet) { return true; } } else if (++dec > sweet) { return false; } } return inc > sweet; } public boolean doCorrection(double [] km, double [] ws) { boolean isIncreasing = isIncreasing(ws); if (isIncreasing) { // mirror along x axis to simulate decreasing values ws = (double [])ws.clone(); for (int i = 0; i < ws.length; ++i) { ws[i] = -ws[i]; } } boolean hasBackJumps = doCorrectionClean(km, ws); if (hasBackJumps && isIncreasing) { // mirror back for (int i = 0; i < corrected.length; ++i) { corrected[i] = -corrected[i]; } } return hasBackJumps; } protected boolean doCorrectionClean(double [] km, double [] ws) { int N = km.length; if (N != ws.length) { throw new IllegalArgumentException("km.length != ws.length"); } if (N < 2) { return false; } SplineInterpolator interpolator = null; for (int i = 1; i < N; ++i) { if (ws[i] <= ws[i-1]) { // no back jump continue; } backjumps.add(km[i]); if (corrected == null) { // lazy cloning ws = corrected = (double [])ws.clone(); } int back = i-2; double distance = Math.abs(km[i] - km[i-1]); double rest = 0.0; double ikm = 0.0; while (back > -1) { if (ws[back] < ws[i]) { // under water // continue scanning backwards distance += Math.abs(km[back+1] - km[back]); --back; continue; } if (ws[back] > ws[i]) { // over water // need to calculate intersection log.debug("intersection case"); double m = (km[back+1]-km[back])/(ws[back+1]-ws[back]); double b = km[back]-ws[back]*m; ikm = m*ws[i] + b; distance += Math.abs(ikm - km[back+1]); rest = Math.abs(km[back] - ikm); } else { // equals height at ws[back] log.debug("exact height match"); distance += Math.abs(km[back+1] - km[back]); ikm = km[back]; } break; } if (back < 0) { log.debug("run over left border"); // fill with same as ws[i] for (int j = 0; j < i; ++j) { ws[j] = ws[i]; } continue; } double quarterDistance = 0.25*distance; // Now further seek back for the max height double restDistance = Math.max(0.0, quarterDistance - rest); --back; double mkmw = ws[i]; double mkm = km[0]; while (back > -1) { double d = Math.abs(km[back+1] - km[back]); restDistance -= d; if (restDistance > 0.0) { --back; continue; } if (restDistance < 0.0) { // need to calculate intersection if (km[back] == km[back+1]) { // should not happen mkm = km[back]; mkmw = 0.5*(ws[back] + ws[back+1]); } else { double m = (ws[back+1]-ws[back])/(km[back+1]-km[back]); double b = ws[back] - km[back]*m; mkm = km[back] + restDistance; mkmw = m*mkm + b; } } else { // exact match mkm = km[back]; mkmw = ws[back]; } break; } double factor = back >= 0 && Math.abs(restDistance) < 1e-4 ? 1.0 : 1.0 - Math.min(1, Math.max(0, restDistance/quarterDistance)); double ikmw = factor*0.25*(mkmw-ws[i]) + ws[i]; double end = ikm + quarterDistance*factor; double [] x = { mkm, ikm, end }; double [] y = { mkmw, ikmw, ws[i] }; if (interpolator == null) { interpolator = new SplineInterpolator(); } PolynomialSplineFunction spline = interpolator.interpolate(x, y); try { if (log.isDebugEnabled()) { log.debug("spline points:"); for (int j = 0; j < x.length; ++j) { log.debug(x[j] + " " + y[j] + " " + spline.value(x[j])); } } for (back = Math.max(back, 0); back < i && km[back] < end; ++back ) { // to 3/4 point fill with spline values ws[back] = spline.value(km[back]); } while (back < i) { // fill the rest with ws[i] ws[back++] = ws[i]; } } catch (ArgumentOutsideDomainException aode) { log.error("spline interpolation failed", aode); } } // for all km return !backjumps.isEmpty(); } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :