z-order space filling curve encoding/decoding/search

index/README.md

@@ -0,0 +1,4 @@
+# geotrellis-index
+-----------------
+
+Provides space z-order space filling curve for building multi-dimentional indices.

index/src/main/scala/geotrellis/index/zcurve/MergeQueue.scala

@@ -0,0 +1,127 @@
+package geotrellis.index.zcurve
+
+class MergeQueue(initialSize: Int = 1) {
+  private var array = if(initialSize <= 1) { Array.ofDim[(Long, Long)](1) } else { Array.ofDim[(Long, Long)](initialSize) }
+  private var _size = 0
+
+  def size = _size
+
+  private def removeElement(i: Int): Unit = {
+    if(i < _size - 1) {
+      val result = array.clone
+      System.arraycopy(array, i + 1, result, i, _size - i - 1)
+      array = result
+    }
+    _size = _size - 1
+  }
+
+  private def insertElement(range: (Long, Long), i: Int): Unit = {
+    ensureSize(_size + 1)
+    if(i == _size) {
+      array(i) = range
+    } else {
+      val result = array.clone
+      System.arraycopy(array, 0, result, 0, i)
+      System.arraycopy(array, i, result, i + 1, _size - i)
+      result(i) = range
+      array = result
+    }
+    _size += 1
+  }
+
+
+  /** Ensure that the internal array has at least `n` cells. */
+  protected def ensureSize(n: Int) {
+    // Use a Long to prevent overflows
+    val arrayLength: Long = array.length
+    if (n > arrayLength - 1) {
+      var newSize: Long = arrayLength * 2
+      while (n > newSize) {
+        newSize = newSize * 2
+      }
+      // Clamp newSize to Int.MaxValue
+      if (newSize > Int.MaxValue) newSize = Int.MaxValue
+
+      val newArray: Array[(Long, Long)] = new Array(newSize.toInt)
+      scala.compat.Platform.arraycopy(array, 0, newArray, 0, _size)
+      array = newArray
+    }
+  }
+
+  val ordering = implicitly[Ordering[(Long, Long)]]
+
+  /** Inserts a single range into the priority queue.
+   *
+   *  @param  range        the element to insert.
+   */
+  @annotation.tailrec
+  final def +=(range: (Long, Long)): Unit = {
+    val res = if(_size == 0) { -1 } else { java.util.Arrays.binarySearch(array, 0, _size, range, ordering) }
+    if(res < 0) {
+      val i = -(res + 1)
+      var (thisStart, thisEnd) = range
+      var removeLeft = false
+
+      var removeRight = false
+      var rightRemainder: Option[(Long, Long)] = None
+
+      // Look at the left range
+      if(i != 0) {
+        val (prevStart, prevEnd) = array(i - 1)
+        if(prevStart == thisStart) {
+          removeLeft = true
+        }
+        if (prevEnd + 1 >= thisStart) { 
+          removeLeft = true
+          thisStart = prevStart
+          if(prevEnd > thisEnd) {
+            thisEnd = prevEnd
+          }
+        }
+      }
+
+      // Look at the right range
+      if(i < _size  && _size > 0) {
+        val (nextStart, nextEnd) = array(i)
+        if(thisStart == nextStart) {
+          removeRight = true
+          thisEnd = nextEnd
+        } else {
+          if(thisEnd + 1 >= nextStart) {
+            removeRight = true
+            if(nextEnd - 1 >= thisEnd) {
+              thisEnd = nextEnd
+            } else if (nextEnd < thisEnd - 1) {
+              rightRemainder = Some((nextEnd + 1, thisEnd))
+              thisEnd = nextEnd
+            }
+          }
+        }
+      }
+
+      if(removeRight) { 
+        if(!removeLeft) {
+          array(i) = (thisStart, thisEnd)
+        } else {
+          array(i-1) = (thisStart, thisEnd)
+          removeElement(i)
+        }
+      } else if(removeLeft) {
+        array(i-1) = (thisStart, thisEnd)
+      } else {
+        insertElement(range, i)
+      }
+
+      rightRemainder match {
+        case Some(r) => this += r
+        case None =>
+      }
+    }
+  }
+
+  def toSeq: Seq[(Long, Long)] = {
+    val result = Array.ofDim[(Long, Long)](size)
+    System.arraycopy(array, 0, result, 0, size)
+    result
+  }
+}

index/src/main/scala/geotrellis/index/zcurve/Z2.scala

@@ -0,0 +1,114 @@
+package geotrellis.index.zcurve
+
+class Z2(val z: Long) extends AnyVal {
+  import Z2._
+
+  def < (other: Z2) = z < other.z
+  def > (other: Z2) = z > other.z
+
+  def + (offset: Long) = new Z2(z + offset)
+  def - (offset: Long) = new Z2(z - offset)
+
+  def == (other: Z2) = other.z == z
+
+  def decode: (Int, Int) = {    
+    ( combine(z), combine(z>>1) )      
+  }
+
+  def dim(i: Int) = Z2.combine(z >> i)
+
+  def mid(p: Z2): Z2 = {
+    var ans: Z2  = new Z2(0)
+    if (p.z < z)
+      ans  = new Z2(p.z + (z - p.z)/2)
+    else
+      ans = new Z2(z + (p.z - z)/2)
+    ans
+  }
+
+  def bitsToString = f"(${z.toBinaryString}%16s)(${dim(0).toBinaryString}%8s,${dim(1).toBinaryString}%8s)"
+  override def toString = f"$z ${decode}"
+}
+
+object Z2 {  
+  final val MAX_BITS = 31
+  final val MAX_MASK = 0x7fffffff // ignore the sign bit, using it breaks < relationship
+  final val MAX_DIM = 2
+
+  /** insert 0 between every bit in value. Only first 31 bits can be considred. */
+  def split(value: Long): Long = {
+    var x: Long = value & MAX_MASK  
+    x = (x ^ (x << 32)) & 0x00000000ffffffffL
+    x = (x ^ (x << 16)) & 0x0000ffff0000ffffL
+    x = (x ^ (x <<  8)) & 0x00ff00ff00ff00ffL // 11111111000000001111111100000000..
+    x = (x ^ (x <<  4)) & 0x0f0f0f0f0f0f0f0fL // 1111000011110000
+    x = (x ^ (x <<  2)) & 0x3333333333333333L // 11001100..
+    x = (x ^ (x <<  1)) & 0x5555555555555555L // 1010...
+    x
+  }
+
+  /** combine every other bit to form a value. Maximum value is 31 bits. */
+  def combine(z: Long): Int = {
+    var x = z & 0x5555555555555555L;
+    x = (x ^ (x >>  1)) & 0x3333333333333333L;
+    x = (x ^ (x >>  2)) & 0x0f0f0f0f0f0f0f0fL;
+    x = (x ^ (x >>  4)) & 0x00ff00ff00ff00ffL;
+    x = (x ^ (x >>  8)) & 0x0000ffff0000ffffL;
+    x = (x ^ (x >> 16)) & 0x00000000ffffffffL;       
+    x.toInt
+  }
+
+  /**
+   * Bits of x and y will be encoded as ....y1x1y0x0
+   */
+  def apply(x: Int, y:  Int): Z2 = 
+    new Z2(split(x) | split(y) << 1)  
+
+  def unapply(z: Z2): Option[(Int, Int)] = 
+    Some(z.decode)
+
+  def zdivide(p: Z2, rmin: Z2, rmax: Z2): (Z2, Z2) = {
+    val (litmax,bigmin) = zdiv(load, MAX_DIM)(p.z, rmin.z, rmax.z)
+    (new Z2(litmax), new Z2(bigmin))
+  }
+
+  /** Loads either 1000... or 0111... into starting at given bit index of a given dimention */
+  def load(target: Long, p: Long, bits: Int, dim: Int): Long = {    
+    val mask = ~(Z2.split(MAX_MASK >> (MAX_BITS-bits)) << dim)
+    val wiped = target & mask
+    wiped | (split(p) << dim)
+  }  
+
+  /** Recurse down the quad-tree and report all z-ranges which are contained in the rectangle defined by the min and max points */
+  def zranges(min: Z2, max: Z2): Seq[(Long, Long)] = {
+    val mq = new MergeQueue
+    val sr = Z2Range(min, max)
+
+    var recCounter = 0
+    var reportCounter = 0
+
+    def _zranges(prefix: Long, offset: Int, quad: Long): Unit = {      
+      recCounter += 1
+
+      val min: Long = prefix | (quad << offset) // QR + 000..
+      val max: Long = min | (1L << offset) - 1  // QR + 111..
+
+      val qr = Z2Range(new Z2(min), new Z2(max))
+      if (sr contains qr){                         // whole range matches, happy day
+        mq += (qr.min.z, qr.max.z)
+        reportCounter +=1
+      } else if (offset > 0 && (sr overlaps qr)) { // some portion of this range are excluded
+        _zranges(min, offset - MAX_DIM, 0)
+        _zranges(min, offset - MAX_DIM, 1)
+        _zranges(min, offset - MAX_DIM, 2)
+        _zranges(min, offset - MAX_DIM, 3)        
+        //let our children punt on each subrange
+      }
+    }
+
+    var prefix: Long = 0
+    var offset = MAX_BITS*MAX_DIM                
+    _zranges(prefix, offset, 0) // the entire space
+    mq.toSeq
+  }
+}

index/src/main/scala/geotrellis/index/zcurve/Z2Range.scala

@@ -0,0 +1,54 @@
+package geotrellis.index.zcurve
+
+import scala.annotation._ 
+
+/**
+ * Represents a rectangle in defined by min and max as two opposing points
+ * 
+ * @param min: lower-left point 
+ * @param max: upper-right point
+ */
+case class Z2Range(min: Z2, max: Z2){  
+  require(min.z < max.z || min.z == max.z, s"NOT: ${min.z} < ${max.z}")  
+  def mid = min mid max
+
+  def length = max.z - min.z
+
+  def contains(z: Z2) = {
+    val (x, y) = z.decode
+    x >= min.dim(0) &&
+    x <= max.dim(0) &&
+    y >= min.dim(1) &&
+    y <= max.dim(1)  
+  }
+
+  def contains(r: Z2Range): Boolean = 
+    contains(r.min) && contains(r.max)
+
+  def overlaps(r: Z2Range): Boolean = {
+    def _overlaps(a1:Int, a2:Int, b1:Int, b2:Int) = math.max(a1,b1) <= math.min(a2,b2)
+    _overlaps(min.dim(0), max.dim(0), r.min.dim(0), r.max.dim(0)) &&
+    _overlaps(min.dim(1), max.dim(1), r.min.dim(1), r.max.dim(1))
+  }
+
+  /** 
+   * Cuts Z-Range in two, can be used to perform augmented binary search
+   * @parm xd: division point
+   * @param inRange: is xd in query range
+   */
+  def cut(xd: Z2, inRange: Boolean): List[Z2Range] = {    
+    if (min.z == max.z)
+      Nil
+    else if (inRange) {
+      if (xd.z == min.z)      // degenerate case, two nodes min has already been counted
+        Z2Range(max, max) :: Nil
+      else if (xd.z == max.z) // degenerate case, two nodes max has already been counted
+        Z2Range(min, min) :: Nil
+      else
+        Z2Range(min, xd-1) :: Z2Range(xd+1, max) :: Nil
+    } else {      
+      val (litmax, bigmin) = Z2.zdivide(xd, min, max)
+      Z2Range(min, litmax) :: Z2Range(bigmin, max) :: Nil
+    }
+  }
+}