Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 633:d08f77e7f7e8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
WST value table: Qs are now stored in ranges for each column. flys-artifacts/trunk@2006 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 25 May 2011 15:31:25 +0000
parents 07640ab913fd
children 433f67a076aa
comparison
equal deleted inserted replaced
632:07640ab913fd 633:d08f77e7f7e8
23 { 23 {
24 private static Logger log = Logger.getLogger(WstValueTable.class); 24 private static Logger log = Logger.getLogger(WstValueTable.class);
25 25
26 public static final int DEFAULT_Q_STEPS = 500; 26 public static final int DEFAULT_Q_STEPS = 500;
27 27
28 public static class Column 28 public static final class Column
29 implements Serializable 29 implements Serializable
30 { 30 {
31 protected String name; 31 protected String name;
32 32
33 protected QRangeTree qRangeTree; 33 protected QRangeTree qRangeTree;
34 34
52 } 52 }
53 53
54 public void setQRangeTree(QRangeTree qRangeTree) { 54 public void setQRangeTree(QRangeTree qRangeTree) {
55 this.qRangeTree = qRangeTree; 55 this.qRangeTree = qRangeTree;
56 } 56 }
57 } 57 } // class Column
58 // class Column 58
59 59 public static final class QPosition {
60 public static class QPosition { 60
61
62 protected double q;
63 protected int index; 61 protected int index;
64 protected double weight; 62 protected double weight;
65 63
66 public QPosition() { 64 public QPosition() {
67 } 65 }
68 66
69 public QPosition(double q, int index, double weight) { 67 public QPosition(int index, double weight) {
70 this.index = index; 68 this.index = index;
71 this.weight = weight; 69 this.weight = weight;
72 } 70 }
73 71
72 public QPosition set(int index, double weight) {
73 this.index = index;
74 this.weight = weight;
75 return this;
76 }
77
74 } // class Position 78 } // class Position
75 79
76 public static class Row 80 public static final class Row
77 implements Serializable, Comparable<Row> 81 implements Serializable, Comparable<Row>
78 { 82 {
79 double km; 83 double km;
80 double [] ws; 84 double [] ws;
81 double [] qs;
82 boolean qSorted;
83 85
84 public Row() { 86 public Row() {
85 } 87 }
86 88
87 public Row(double km) { 89 public Row(double km) {
88 this.km = km; 90 this.km = km;
89 } 91 }
90 92
91 public Row(double km, double [] ws, double [] qs) { 93 public Row(double km, double [] ws) {
92 this(km); 94 this(km);
93 this.ws = ws; 95 this.ws = ws;
94 this.qs = qs;
95 }
96
97 public static final double linear(
98 double x,
99 double x1, double x2,
100 double y1, double y2
101 ) {
102 // y1 = m*x1 + b
103 // y2 = m*x2 + b
104 // y2 - y1 = m*(x2 - x1)
105 // m = (y2 - y1)/(x2 - x1) # x2 != x1
106 // b = y1 - m*x1
107
108 if (x1 == x2) {
109 return 0.5*(y1 + y2);
110 }
111 double m = (y2 - y1)/(x2 - x1);
112 double b = y1 - m*x1;
113 return x*m + b;
114 }
115
116 public static final double factor(double x, double p1, double p2) {
117 // 0 = m*p1 + b <=> b = -m*p1
118 // 1 = m*p2 + b
119 // 1 = m*(p2 - p1)
120 // m = 1/(p2 - p1) # p1 != p2
121 // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1)
122
123 return p1 == p2 ? 0.0 : (x-p1)/(p2-p1);
124 }
125
126 public static final double weight(double factor, double a, double b) {
127 return (1.0-factor)*a + factor*b;
128 }
129
130 public void interpolateW(
131 Row other,
132 double km,
133 double [] input,
134 double [] output
135 ) {
136 double factor = factor(km, this.km, other.km);
137
138 for (int i = 0; i < input.length; ++i) {
139 double q = input[i];
140 int idx1 = getQIndex(q);
141 int idx2 = other.getQIndex(q);
142
143 double w1 = idx1 >= 0
144 ? ws[idx1]
145 : interpolateW(-idx1-1, q);
146
147 double w2 = idx2 >= 0
148 ? other.ws[idx2]
149 : other.interpolateW(-idx2-1, q);
150
151 output[i] = weight(factor, w1, w2);
152 }
153 }
154
155 public double interpolateW(int idx, double q) {
156 return idx < 1 || idx >= qs.length
157 ? Double.NaN // do not extrapolate
158 : linear(q, qs[idx-1], qs[idx], ws[idx-1], ws[idx]);
159 }
160
161 public double interpolateW(double q) {
162 if (Double.isNaN(q)) return Double.NaN;
163 int index = getQIndex(q);
164 return index >= 0 ? ws[index] : interpolateW(-index -1, q);
165 }
166
167 public int getQIndex(double q) {
168 return qSorted ? binaryQIndex(q) : linearQIndex(q);
169 }
170
171 public int binaryQIndex(double q) {
172 return Arrays.binarySearch(qs, q);
173 }
174
175 public int linearQIndex(double q) {
176 switch (qs.length) {
177 case 0: break;
178 case 1:
179 if (qs[0] == q) return 0;
180 if (qs[0] < q) return -(1+1);
181 return -(0+1);
182 default:
183 for (int i = 1; i < qs.length; ++i) {
184 double qa = qs[i-1];
185 double qb = qs[i];
186 if (qa == q) return i-1;
187 if (qb == q) return i;
188 if (qa > qb) { double t = qa; qa = qb; qb = t; }
189 if (q > qa && q < qb) return -(i+1);
190 }
191 return -qs.length;
192 }
193
194 return -1;
195 } 96 }
196 97
197 public int compareTo(Row other) { 98 public int compareTo(Row other) {
198 double d = km - other.km; 99 double d = km - other.km;
199 if (d < -0.0001) return -1; 100 if (d < -0.0001) return -1;
200 if (d > 0.0001) return +1; 101 if (d > 0.0001) return +1;
201 return 0; 102 return 0;
202 } 103 }
203 104
204 public double [][] interpolateWQ(Row other, double km, int steps) { 105 public void interpolateW(
205 106 Row other,
107 double km,
108 double [] iqs,
109 double [] ows,
110 WstValueTable table
111 ) {
112 double kmWeight = factor(km, this.km, other.km);
113
114 QPosition qPosition = new QPosition();
115
116 for (int i = 0; i < iqs.length; ++i) {
117 if (table.getQPosition(km, iqs[i], qPosition) != null) {
118 double wt = getW(qPosition);
119 double wo = other.getW(qPosition);
120 ows[i] = weight(kmWeight, wt, wo);
121 }
122 else {
123 ows[i] = Double.NaN;
124 }
125 }
126 }
127
128 public double [][] interpolateWQ(
129 Row other,
130 double km,
131 int steps,
132 WstValueTable table
133 ) {
206 int W = Math.min(ws.length, other.ws.length); 134 int W = Math.min(ws.length, other.ws.length);
207 135
208 if (W < 1) { 136 if (W < 1) {
209 return new double[2][0]; 137 return new double[2][0];
210 } 138 }
217 double minQ = Double.MAX_VALUE; 145 double minQ = Double.MAX_VALUE;
218 double maxQ = -Double.MAX_VALUE; 146 double maxQ = -Double.MAX_VALUE;
219 147
220 for (int i = 0; i < W; ++i) { 148 for (int i = 0; i < W; ++i) {
221 double wws = weight(factor, ws[i], other.ws[i]); 149 double wws = weight(factor, ws[i], other.ws[i]);
222 double wqs = weight(factor, qs[i], other.qs[i]); 150 double wqs = weight(
151 factor,
152 table.getQIndex(i, km),
153 table.getQIndex(i, other.km));
154
223 splineW[i] = wws; 155 splineW[i] = wws;
224 splineQ[i] = wqs; 156 splineQ[i] = wqs;
225 if (wqs < minQ) minQ = wqs; 157 if (wqs < minQ) minQ = wqs;
226 if (wqs > maxQ) maxQ = wqs; 158 if (wqs > maxQ) maxQ = wqs;
227 } 159 }
246 } 178 }
247 179
248 return new double [][] { outWs, outQs }; 180 return new double [][] { outWs, outQs };
249 } 181 }
250 182
251 public double [][] interpolateWQ(int steps) { 183 public double [][] interpolateWQ(int steps, WstValueTable table) {
252 184
253 int W = ws.length; 185 int W = ws.length;
254 186
255 if (W < 1) { 187 if (W < 1) {
256 return new double[2][0]; 188 return new double[2][0];
257 } 189 }
258 190
259 double [] splineW = new double[W];
260 double [] splineQ = new double[W]; 191 double [] splineQ = new double[W];
261 192
262 double minQ = Double.MAX_VALUE; 193 double minQ = Double.MAX_VALUE;
263 double maxQ = -Double.MAX_VALUE; 194 double maxQ = -Double.MAX_VALUE;
264 195
196 QPosition qPosition = new QPosition();
197
265 for (int i = 0; i < W; ++i) { 198 for (int i = 0; i < W; ++i) {
266 splineW[i] = ws[i]; 199 splineQ[i] = table.getQIndex(i, km);
267 splineQ[i] = qs[i]; 200 if (splineQ[i] < minQ) minQ = splineQ[i];
268 if (qs[i] < minQ) minQ = qs[i]; 201 if (splineQ[i] > maxQ) maxQ = splineQ[i];
269 if (qs[i] > maxQ) maxQ = qs[i];
270 } 202 }
271 203
272 double stepWidth = (maxQ - minQ)/steps; 204 double stepWidth = (maxQ - minQ)/steps;
273 205
274 SplineInterpolator interpolator = new SplineInterpolator(); 206 SplineInterpolator interpolator = new SplineInterpolator();
275 207
276 PolynomialSplineFunction spline = 208 PolynomialSplineFunction spline =
277 interpolator.interpolate(splineQ, splineW); 209 interpolator.interpolate(splineQ, ws);
278 210
279 double [] outWs = new double[steps]; 211 double [] outWs = new double[steps];
280 double [] outQs = new double[steps]; 212 double [] outQs = new double[steps];
281 213
282 try { 214 try {
290 } 222 }
291 223
292 return new double [][] { outWs, outQs }; 224 return new double [][] { outWs, outQs };
293 } 225 }
294 226
295 public QPosition getQPosition(double q) { 227
296 int qIndex = getQIndex(q); 228 public double getW(QPosition qPosition) {
297 229 int index = qPosition.index;
298 if (qIndex >= 0) { 230 double weight = qPosition.weight;
299 // direct column match 231
300 return new QPosition(q, qIndex, 1.0); 232 return weight == 1.0
301 } 233 ? ws[index]
302 234 : weight(weight, ws[index-1], ws[index]);
303 qIndex = -qIndex -1; 235 }
304 236
305 if (qIndex < 0 || qIndex >= qs.length-1) { 237 public double getW(
306 // do not extrapolate
307 return null;
308 }
309
310 return new QPosition(
311 q, qIndex, factor(q, qs[qIndex-1], qs[qIndex]));
312 }
313
314 public QPosition getQPosition(Row other, double km, double q) {
315
316 QPosition qpt = getQPosition(q);
317 QPosition qpo = other.getQPosition(q);
318
319 if (qpt == null || qpo == null) {
320 return null;
321 }
322
323 double kmWeight = factor(km, this.km, other.km);
324
325 // XXX: Index interpolation is a bit sticky here!
326
327 int index = (int)Math.round(
328 weight(kmWeight, qpt.index, qpo.index));
329
330 double weight = weight(kmWeight, qpt.weight, qpo.weight);
331
332 return new QPosition(q, index, weight);
333 }
334
335 public void storeWQ(
336 QPosition qPosition,
337 double [] ows,
338 double [] oqs,
339 int index
340 ) {
341 int qIdx = qPosition.index;
342 double qWeight = qPosition.weight;
343
344 if (qWeight == 1.0) {
345 oqs[index] = qs[qIdx];
346 ows[index] = ws[qIdx];
347 }
348 else {
349 oqs[index] = weight(qWeight, qs[qIdx-1], qs[qIdx]);
350 ows[index] = weight(qWeight, ws[qIdx-1], ws[qIdx]);
351 }
352 }
353
354 public void storeWQ(
355 Row other, 238 Row other,
356 double km, 239 double km,
357 QPosition qPosition, 240 QPosition qPosition
358 double [] ows,
359 double [] oqs,
360 int index
361 ) { 241 ) {
362 double kmWeight = factor(km, this.km, other.km); 242 double kmWeight = factor(km, this.km, other.km);
363 243
364 double tq, tw; 244 int index = qPosition.index;
365 double oq, ow; 245 double weight = qPosition.weight;
366 246
367 int qIdx = qPosition.index; 247 double tw, ow;
368 double qWeight = qPosition.weight; 248
369 249 if (weight == 1.0) {
370 if (qWeight == 1.0) { 250 tw = ws[index];
371 tq = qs[qIdx]; 251 ow = other.ws[index];
372 tw = ws[qIdx];
373 oq = other.qs[qIdx];
374 ow = other.ws[qIdx];
375 } 252 }
376 else { 253 else {
377 tq = weight(qWeight, qs[qIdx-1], qs[qIdx]); 254 tw = weight(weight, ws[index-1], ws[index]);
378 tw = weight(qWeight, ws[qIdx-1], ws[qIdx]); 255 ow = weight(weight, other.ws[index-1], other.ws[index]);
379 oq = weight(qWeight, other.qs[qIdx-1], other.qs[qIdx]); 256 }
380 ow = weight(qWeight, other.ws[qIdx-1], other.ws[qIdx]); 257
381 } 258 return weight(kmWeight, tw, ow);
382 259 }
383 oqs[index] = weight(kmWeight, tq, oq); 260 } // class Row
384 ows[index] = weight(kmWeight, tw, ow);
385 }
386
387 public void storeWQ(
388 QPosition qPosition,
389 LinearRemap remap,
390 double [] ows,
391 double [] oqs,
392 int index
393 ) {
394 int qIdx = qPosition.index;
395 double qWeight = qPosition.weight;
396
397 oqs[index] = remap.eval(km, qWeight == 1.0
398 ? qs[qIdx]
399 : weight(qWeight, qs[qIdx-1], qs[qIdx]));
400
401 ows[index] = interpolateW(oqs[index]);
402 }
403
404 public void storeWQ(
405 Row other,
406 double km,
407 QPosition qPosition,
408 LinearRemap remap,
409 double [] ows,
410 double [] oqs,
411 int index
412 ) {
413 double kmWeight = factor(km, this.km, other.km);
414
415 double tq, tw;
416 double oq, ow;
417
418 int qIdx = qPosition.index;
419 double qWeight = qPosition.weight;
420
421 if (qWeight == 1.0) {
422 tq = remap.eval(this.km, qs[qIdx]);
423 oq = remap.eval(other.km, other.qs[qIdx]);
424 }
425 else {
426 tq = remap.eval(
427 this.km,
428 weight(qWeight, qs[qIdx-1], qs[qIdx]));
429 oq = remap.eval(
430 other.km,
431 weight(qWeight, other.qs[qIdx-1], other.qs[qIdx]));
432 }
433
434 tw = interpolateW(tq);
435 ow = other.interpolateW(oq);
436
437 oqs[index] = weight(kmWeight, tq, oq);
438 ows[index] = weight(kmWeight, tw, ow);
439 }
440 }
441 // class Row
442 261
443 protected List<Row> rows; 262 protected List<Row> rows;
444 263
445 protected Column [] columns; 264 protected Column [] columns;
446 265
451 public WstValueTable(Column [] columns) { 270 public WstValueTable(Column [] columns) {
452 this(); 271 this();
453 this.columns = columns; 272 this.columns = columns;
454 } 273 }
455 274
275 public WstValueTable(Column [] columns, List<Row> rows) {
276 this.columns = columns;
277 this.rows = rows;
278 }
279
456 public void sortRows() { 280 public void sortRows() {
457 Collections.sort(rows); 281 Collections.sort(rows);
458 } 282 }
459 283
460
461 public double [] interpolateW(double km, double [] qs) {
462 return interpolateW(km, qs, new double[qs.length]);
463 }
464
465 public double [] interpolateW(double km, double [] qs, double [] ws) { 284 public double [] interpolateW(double km, double [] qs, double [] ws) {
466 285
467 int rowIndex = Collections.binarySearch(rows, new Row(km)); 286 int rowIndex = Collections.binarySearch(rows, new Row(km));
287
288 QPosition qPosition = new QPosition();
468 289
469 if (rowIndex >= 0) { // direct row match 290 if (rowIndex >= 0) { // direct row match
470 Row row = rows.get(rowIndex); 291 Row row = rows.get(rowIndex);
471 for (int i = 0; i < qs.length; ++i) { 292 for (int i = 0; i < qs.length; ++i) {
472 ws[i] = row.interpolateW(qs[i]); 293 ws[i] = getQPosition(km, qs[i], qPosition) != null
294 ? row.getW(qPosition)
295 : Double.NaN;
473 } 296 }
474 } 297 }
475 else { // needs bilinear interpolation 298 else { // needs bilinear interpolation
476 rowIndex = -rowIndex -1; 299 rowIndex = -rowIndex -1;
477 300
478 if (rowIndex < 1 || rowIndex >= rows.size()) { 301 if (rowIndex < 1 || rowIndex >= rows.size()) {
479 // do not extrapolate 302 // do not extrapolate
480 Arrays.fill(ws, Double.NaN); 303 Arrays.fill(ws, Double.NaN);
481 } 304 }
482 else { 305 else {
483 rows.get(rowIndex-1).interpolateW( 306 Row r1 = rows.get(rowIndex-1);
484 rows.get(rowIndex), 307 Row r2 = rows.get(rowIndex);
485 km, qs, ws); 308 r1.interpolateW(r2, km, qs, ws, this);
486 } 309 }
487 } 310 }
488 311
489 return ws; 312 return ws;
490 } 313 }
491
492 314
493 public double [][] interpolateWQ(double km) { 315 public double [][] interpolateWQ(double km) {
494 return interpolateWQ(km, DEFAULT_Q_STEPS); 316 return interpolateWQ(km, DEFAULT_Q_STEPS);
495 } 317 }
496 318
498 320
499 int rowIndex = Collections.binarySearch(rows, new Row(km)); 321 int rowIndex = Collections.binarySearch(rows, new Row(km));
500 322
501 if (rowIndex >= 0) { // direct row match 323 if (rowIndex >= 0) { // direct row match
502 Row row = rows.get(rowIndex); 324 Row row = rows.get(rowIndex);
503 return row.interpolateWQ(steps); 325 return row.interpolateWQ(steps, this);
504 } 326 }
505 327
506 rowIndex = -rowIndex -1; 328 rowIndex = -rowIndex -1;
507 329
508 if (rowIndex < 1 || rowIndex >= rows.size()) { 330 if (rowIndex < 1 || rowIndex >= rows.size()) {
511 } 333 }
512 334
513 Row r1 = rows.get(rowIndex-1); 335 Row r1 = rows.get(rowIndex-1);
514 Row r2 = rows.get(rowIndex); 336 Row r2 = rows.get(rowIndex);
515 337
516 return r1.interpolateWQ(r2, km, steps); 338 return r1.interpolateWQ(r2, km, steps, this);
517 } 339 }
518 340
519 public void interpolate( 341 public void interpolate(
520 double [] kms, 342 double [] kms,
521 double [] ws, 343 double [] ws,
524 LinearRemap remap 346 LinearRemap remap
525 ) { 347 ) {
526 int R1 = rows.size()-1; 348 int R1 = rows.size()-1;
527 349
528 Row kmKey = new Row(); 350 Row kmKey = new Row();
351
352 QPosition nPosition = new QPosition();
353
529 for (int i = 0; i < kms.length; ++i) { 354 for (int i = 0; i < kms.length; ++i) {
530 kmKey.km = kms[i]; 355 kmKey.km = kms[i];
356
357 qs[i] = remap.eval(kms[i], getQ(qPosition, kms[i]));
358
359 if (getQPosition(kms[i], qs[i], nPosition) == null) {
360 ws[i] = Double.NaN;
361 continue;
362 }
363
531 int rowIndex = Collections.binarySearch(rows, kmKey); 364 int rowIndex = Collections.binarySearch(rows, kmKey);
532 365
533 if (rowIndex >= 0) { 366 if (rowIndex >= 0) {
534 // direct row match 367 // direct row match
535 rows.get(rowIndex).storeWQ(qPosition, remap, ws, qs, i); 368 ws[i] = rows.get(rowIndex).getW(nPosition);
536 continue; 369 continue;
537 } 370 }
538 371
539 rowIndex = -rowIndex -1; 372 rowIndex = -rowIndex -1;
540 373
541 if (rowIndex < 1 || rowIndex >= R1) { 374 if (rowIndex < 1 || rowIndex >= R1) {
542 // do not extrapolate 375 // do not extrapolate
543 ws[i] = qs[i] = Double.NaN; 376 ws[i] = Double.NaN;
544 continue; 377 continue;
545 } 378 }
546 Row r1 = rows.get(rowIndex-1); 379 Row r1 = rows.get(rowIndex-1);
547 Row r2 = rows.get(rowIndex); 380 Row r2 = rows.get(rowIndex);
548 r1.storeWQ(r2, kms[i], qPosition, remap, ws, qs, i); 381 ws[i] = r1.getW(r2, kms[i], nPosition);
549 } 382 }
550 } 383 }
551 384
552 public QPosition interpolate( 385 public QPosition interpolate(
553 double q, 386 double q,
560 393
561 int rowIndex = Collections.binarySearch(rows, kmKey); 394 int rowIndex = Collections.binarySearch(rows, kmKey);
562 395
563 int R1 = rows.size()-1; 396 int R1 = rows.size()-1;
564 397
565 QPosition qPosition = null; 398 QPosition qPosition = getQPosition(kms[referenceIndex], q);
566 399
567 if (rowIndex >= 0) { // direct row match 400 if (qPosition == null) {
568 qPosition = rows.get(rowIndex).getQPosition(q); 401 // we cannot locate q at km
569 } 402 return null;
570 else { 403 }
404
405 for (int i = 0; i < kms.length; ++i) {
406 kmKey.km = kms[i];
407
408 rowIndex = Collections.binarySearch(rows, kmKey);
409
410 qs[i] = getQ(qPosition, kms[i]);
411
412 if (rowIndex >= 0) {
413 // direct row match
414 ws[i] = rows.get(rowIndex).getW(qPosition);
415 continue;
416 }
417
571 rowIndex = -rowIndex -1; 418 rowIndex = -rowIndex -1;
572 419
573 if (rowIndex < 1 || rowIndex >= R1) { 420 if (rowIndex < 1 || rowIndex >= R1) {
574 // do not extrapolate 421 // do not extrapolate
575 return null; 422 ws[i] = Double.NaN;
576 } 423 continue;
577 424 }
578 // interpolation case
579 Row r1 = rows.get(rowIndex-1); 425 Row r1 = rows.get(rowIndex-1);
580 Row r2 = rows.get(rowIndex); 426 Row r2 = rows.get(rowIndex);
581 427
582 qPosition = r1.getQPosition(r2, kms[referenceIndex], q); 428 ws[i] = r1.getW(r2, kms[i], qPosition);
583 } 429 }
584 430
585 if (qPosition == null) { 431 return qPosition;
586 // we cannot locate q in reference row 432 }
433
434 public QPosition getQPosition(double km, double q) {
435 return getQPosition(km, q, new QPosition());
436 }
437
438 public QPosition getQPosition(double km, double q, QPosition qPosition) {
439
440 if (columns.length == 0) {
587 return null; 441 return null;
588 } 442 }
589 443
590 for (int i = 0; i < kms.length; ++i) { 444 double qLast = columns[0].getQRangeTree().findQ(km);
591 kmKey.km = kms[i]; 445
592 446 if (Math.abs(qLast - q) < 0.00001) {
593 rowIndex = Collections.binarySearch(rows, kmKey); 447 return qPosition.set(0, 1d);
594 448 }
595 if (rowIndex >= 0) { 449
596 // direct row match 450 for (int i = 1; i < columns.length; ++i) {
597 rows.get(rowIndex).storeWQ(qPosition, ws, qs, i); 451 double qCurrent = columns[i].getQRangeTree().findQ(km);
598 continue; 452 if (Math.abs(qCurrent - q) < 0.00001) {
599 } 453 return qPosition.set(i, 1d);
600 454 }
601 rowIndex = -rowIndex -1; 455
602 456 double qMin, qMax;
603 if (rowIndex < 1 || rowIndex >= R1) { 457 if (qLast < qCurrent) { qMin = qLast; qMax = qCurrent; }
604 // do not extrapolate 458 else { qMin = qCurrent; qMax = qLast; }
605 ws[i] = qs[i] = Double.NaN; 459
606 continue; 460 if (q > qMin && q < qMax) {
607 } 461 double weight = factor(q, qLast, qCurrent);
608 Row r1 = rows.get(rowIndex-1); 462 return qPosition.set(i, weight);
609 Row r2 = rows.get(rowIndex); 463 }
610 r1.storeWQ(r2, kms[i], qPosition, ws, qs, i); 464 qLast = qCurrent;
611 } 465 }
612 466
613 return qPosition; 467 return null;
468 }
469
470 public double getQIndex(int index, double km) {
471 return columns[index].getQRangeTree().findQ(km);
472 }
473
474 public double getQ(QPosition qPosition, double km) {
475 int index = qPosition.index;
476 double weight = qPosition.weight;
477
478 if (weight == 1d) {
479 return columns[index].getQRangeTree().findQ(km);
480 }
481 double q1 = columns[index-1].getQRangeTree().findQ(km);
482 double q2 = columns[index ].getQRangeTree().findQ(km);
483 return weight(weight, q1, q2);
484 }
485
486 public static final double linear(
487 double x,
488 double x1, double x2,
489 double y1, double y2
490 ) {
491 // y1 = m*x1 + b
492 // y2 = m*x2 + b
493 // y2 - y1 = m*(x2 - x1)
494 // m = (y2 - y1)/(x2 - x1) # x2 != x1
495 // b = y1 - m*x1
496
497 if (x1 == x2) {
498 return 0.5*(y1 + y2);
499 }
500 double m = (y2 - y1)/(x2 - x1);
501 double b = y1 - m*x1;
502 return x*m + b;
503 }
504
505 public static final double factor(double x, double p1, double p2) {
506 // 0 = m*p1 + b <=> b = -m*p1
507 // 1 = m*p2 + b
508 // 1 = m*(p2 - p1)
509 // m = 1/(p2 - p1) # p1 != p2
510 // f(x) = x/(p2-p1) - p1/(p2-p1) <=> (x-p1)/(p2-p1)
511
512 return p1 == p2 ? 0.0 : (x-p1)/(p2-p1);
513 }
514
515 public static final double weight(double factor, double a, double b) {
516 //return (1.0-factor)*a + factor*b;
517 return a + factor*(b-a);
614 } 518 }
615 } 519 }
616 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : 520 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org