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 : |