sascha@625: package de.intevation.flys.artifacts.model;
sascha@625: 
sascha@625: import java.io.Serializable;
sascha@625: 
sascha@625: import java.math.BigDecimal;
sascha@625: 
sascha@625: import java.util.List;
sascha@625: 
sascha@625: import de.intevation.flys.model.WstQRange;
sascha@625: import de.intevation.flys.model.Range;
sascha@625: 
sascha@625: public class QRangeTree
sascha@625: implements   Serializable
sascha@625: {
sascha@625:     public static class Node
sascha@625:     implements          Serializable
sascha@625:     {
sascha@625:         Node left;
sascha@625:         Node right;
sascha@625:         Node prev;
sascha@625:         Node next;
sascha@625: 
sascha@625:         double a;
sascha@625:         double b;
sascha@625:         double q;
sascha@625: 
sascha@625:         public Node() {
sascha@625:         }
sascha@625: 
sascha@625:         public Node(double a, double b, double q) {
sascha@625:             this.a = a;
sascha@625:             this.b = b;
sascha@625:             this.q = q;
sascha@625:         }
sascha@625: 
sascha@625:         protected final double interpolatePrev(double pos) {
sascha@625:             /*
sascha@625:             f(prev.b) = prev.q
sascha@625:             f(a)      = q
sascha@625: 
sascha@625:             prev.q = m*prev.b + n
sascha@625:             q      = m*a      + n <=> n = q - m*a
sascha@625: 
sascha@625:             q - prev.q = m*(a - prev.b)
sascha@625: 
sascha@625:             m = (q - prev.q)/(a - prev.b) # a != prev.b
sascha@625:             */
sascha@625: 
sascha@625:             if (a == prev.b) {
sascha@625:                 return 0.5*(q + prev.q);
sascha@625:             }
sascha@625:             double m = (q - prev.q)/(a - prev.b);
sascha@625:             double n = q - m*a;
sascha@625:             return m*pos + n;
sascha@625:         }
sascha@625: 
sascha@625:         protected final double interpolateNext(double pos) {
sascha@625:             /*
sascha@625:             f(next.a) = next.q
sascha@625:             f(b)      = q
sascha@625: 
sascha@625:             next.q = m*next.a + n
sascha@625:             q      = m*b      + n <=> n = q - m*b
sascha@625: 
sascha@625:             q - next.q = m*(b - next.a)
sascha@625:             m = (q - next.q)/(b - next.a) # b != next.a
sascha@625:             */
sascha@625: 
sascha@625:             if (b == next.a) {
sascha@625:                 return 0.5*(q + next.q);
sascha@625:             }
sascha@625:             double m = (q - next.q)/(b - next.a);
sascha@625:             double n = q - m*b;
sascha@625:             return m*pos + n;
sascha@625:         }
sascha@625: 
sascha@625:         public double findQ(double pos) {
sascha@625:             Node current = this;
sascha@625:             for (;;) {
sascha@625:                 if (pos < current.a) {
sascha@625:                     if (current.left != null) {
sascha@625:                         current = current.left;
sascha@625:                         continue;
sascha@625:                     }
sascha@625:                     return current.left != null
sascha@625:                         ? current.interpolatePrev(pos)
sascha@625:                         : Double.NaN;
sascha@625:                 }
sascha@625:                 if (pos > current.b) {
sascha@625:                     if (current.right != null) {
sascha@625:                         current = current.right;
sascha@625:                         continue;
sascha@625:                     }
sascha@625:                     return current.right != null
sascha@625:                         ? current.interpolateNext(pos)
sascha@625:                         : Double.NaN;
sascha@625:                 }
sascha@625:                 return current.q;
sascha@625:             }
sascha@625:         }
sascha@625:     } // class Node
sascha@625: 
sascha@625:     protected Node root;
sascha@625: 
sascha@625:     public QRangeTree() {
sascha@625:     }
sascha@625: 
sascha@625:     /** wstQRanges need to sorted by range.a */
sascha@625:     public QRangeTree(List<WstQRange> wstQRanges) {
sascha@625: 
sascha@625:         if (wstQRanges.isEmpty()) {
sascha@625:             return;
sascha@625:         }
sascha@625: 
sascha@625:         Node [] nodes = new Node[wstQRanges.size()];
sascha@625:         for (int i = 0; i < nodes.length; ++i) {
sascha@625:             WstQRange wstQRange = wstQRanges.get(i);
sascha@625:             Range     range     = wstQRange.getRange();
sascha@625:             BigDecimal a = range.getA();
sascha@625:             BigDecimal b = range.getB();
sascha@625:             nodes[i] = new Node(
sascha@625:                 a != null ? a.doubleValue() : -Double.MAX_VALUE,
sascha@625:                 b != null ? b.doubleValue() :  Double.MAX_VALUE,
sascha@625:                 wstQRange.getQ().doubleValue());
sascha@625:         }
sascha@625: 
sascha@625:         for (int i = 0; i < nodes.length; ++i) {
sascha@625:             Node node = nodes[i];
sascha@625:             if (i > 0             ) node.prev = nodes[i-1];
sascha@625:             if (i < nodes.length-1) node.next = nodes[i+1];
sascha@625:         }
sascha@625: 
sascha@625:         root = buildTree(nodes, 0, nodes.length-1);
sascha@625:     }
sascha@625: 
sascha@625:     protected static Node buildTree(Node [] nodes, int lo, int hi) {
sascha@625:         int mid = (lo + hi) >> 1;
sascha@625:         Node parent = nodes[mid];
sascha@625: 
sascha@625:         if (hi > lo) {
sascha@625:             parent.left  = buildTree(nodes, lo, mid-1);
sascha@625:             parent.right = buildTree(nodes, mid+1, hi);
sascha@625:         }
sascha@625: 
sascha@625:         return parent;
sascha@625:     }
sascha@625: 
sascha@625:     public double findQ(double pos) {
sascha@625:         return root != null ? root.findQ(pos) : Double.NaN;
sascha@625:     }
sascha@625: }
sascha@625: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :