1 /***********************************************************************
2  * Copyright (c) 2013-2024 Commonwealth Computer Research, Inc.
3  * All rights reserved. This program and the accompanying materials
4  * are made available under the terms of the Apache License, Version 2.0
5  * which accompanies this distribution and is available at
6  * http://www.opensource.org/licenses/apache2.0.php.
7  ***********************************************************************/
8 
9 package org.locationtech.geomesa.curve
10 
11 import org.locationtech.geomesa.curve.TimePeriod.TimePeriod
12 import org.locationtech.geomesa.curve.XZ3SFC.{QueryWindow, XElement}
13 import org.locationtech.geomesa.zorder.sfcurve.IndexRange
14 
15 import scala.collection.mutable.ArrayBuffer
16 
17 /**
18   * Extended Z-order curve implementation used for efficiently storing polygons.
19   *
20   * Based on 'XZ-Ordering: A Space-Filling Curve for Objects with Spatial Extension'
21   * by Christian Böhm, Gerald Klump  and Hans-Peter Kriegel, expanded to 3 dimensions
22   *
23   * @param g resolution level of the curve - i.e. how many times the space will be recursively split into eighths
24 
25   */
26 class XZ3SFC(val g: Short, val xBounds: (Double, Double), val yBounds: (Double, Double), val zBounds: (Double, Double)) {
27 
28   // TODO see what the max value of g can be where we can use Ints instead of Longs and possibly refactor to use Ints
29 
30   private val xLo = xBounds._1
31   private val xHi = xBounds._2
32   private val yLo = yBounds._1
33   private val yHi = yBounds._2
34   private val zLo = zBounds._1
35   private val zHi = zBounds._2
36 
37   private val xSize = xHi - xLo
38   private val ySize = yHi - yLo
39   private val zSize = zHi - zLo
40 
41   /**
42     * Index a polygon by it's bounding box
43     *
44     * @param xmin min x value in xBounds
45     * @param ymin min y value in yBounds
46     * @param zmin min z value in zBounds
47     * @param xmax max x value in xBounds, must be >= xmin
48     * @param ymax max y value in yBounds, must be >= ymin
49     * @param zmax max z value in zBounds, must be >= tmin
50     * @param lenient standardize boundaries to valid values, or raise an exception
51     * @return z value for the bounding box
52     */
53   def index(xmin: Double, ymin: Double, zmin: Double, xmax: Double, ymax: Double, zmax: Double, lenient: Boolean = false): Long = {
54     // normalize inputs to [0,1]
55     val (nxmin, nymin, nzmin, nxmax, nymax, nzmax) = normalize(xmin, ymin, zmin, xmax, ymax, zmax, lenient)
56 
57     // calculate the length of the sequence code (section 4.1 of XZ-Ordering paper)
58 
59     val maxDim = math.max(math.max(nxmax - nxmin, nymax - nymin), nzmax - nzmin)
60 
61     // l1 (el-one) is a bit confusing to read, but corresponds with the paper's definitions
62     val l1 = math.floor(math.log(maxDim) / XZSFC.LogPointFive).toInt
63 
64     // the length will either be (l1) or (l1 + 1)
65     val length = if (l1 >= g) { g } else {
66       val w2 = math.pow(0.5, l1 + 1) // width of an element at resolution l2 (l1 + 1)
67 
68       // predicate for checking how many axis the polygon intersects
69       // math.floor(min / w2) * w2 == start of cell containing min
70       def predicate(min: Double, max: Double): Boolean = max <= (math.floor(min / w2) * w2) + (2 * w2)
71 
72       if (predicate(nxmin, nxmax) && predicate(nymin, nymax) && predicate(nzmin, nzmax)) l1 + 1 else l1
73     }
74 
75     sequenceCode(nxmin, nymin, nzmin, length)
76   }
77 
78   /**
79     * Determine XZ-curve ranges that will cover a given query window
80     *
81     * @param query a window to cover in the form (xmin, ymin, zmin, xmax, ymax, zmax) where all values are in user space
82     * @return
83     */
84   def ranges(query: (Double, Double, Double, Double, Double, Double)): Seq[IndexRange] = ranges(Seq(query))
85 
86   /**
87     * Determine XZ-curve ranges that will cover a given query window
88     *
89     * @param query a window to cover in the form (xmin, ymin, zmin, xmax, ymax, zmax) where all values are in user space
90     * @param maxRanges a rough upper limit on the number of ranges to generate
91     * @return
92     */
93   def ranges(query: (Double, Double, Double, Double, Double, Double), maxRanges: Option[Int]): Seq[IndexRange] =
94     ranges(Seq(query), maxRanges)
95 
96   /**
97     * Determine XZ-curve ranges that will cover a given query window
98     *
99     * @param xmin min x value in user space
100     * @param ymin min y value in user space
101     * @param zmin min z value in user space
102     * @param xmax max x value in user space, must be >= xmin
103     * @param ymax max y value in user space, must be >= ymin
104     * @param zmax max z value in user space, must be >= zmin
105     * @return
106     */
107   def ranges(xmin: Double, ymin: Double, zmin: Double, xmax: Double, ymax: Double, zmax: Double): Seq[IndexRange] =
108     ranges(Seq((xmin, ymin, zmin, xmax, ymax, zmax)))
109 
110   /**
111     * Determine XZ-curve ranges that will cover a given query window
112     *
113     * @param xmin min x value in user space
114     * @param ymin min y value in user space
115     * @param zmin min z value in user space
116     * @param xmax max x value in user space, must be >= xmin
117     * @param ymax max y value in user space, must be >= ymin
118     * @param zmax max z value in user space, must be >= zmin
119     * @param maxRanges a rough upper limit on the number of ranges to generate
120     * @return
121     */
122   def ranges(xmin: Double,
123              ymin: Double,
124              zmin: Double,
125              xmax: Double,
126              ymax: Double,
127              zmax: Double,
128              maxRanges: Option[Int]): Seq[IndexRange] =
129     ranges(Seq((xmin, ymin, zmin, xmax, ymax, zmax)), maxRanges)
130 
131   /**
132     * Determine XZ-curve ranges that will cover a given query window
133     *
134     * @param queries a sequence of OR'd windows to cover. Each window is in the form
135     *                (xmin, ymin, zmin, xmax, ymax, zmax) where all values are in user space
136     * @param maxRanges a rough upper limit on the number of ranges to generate
137     * @return
138     */
139   def ranges(queries: Seq[(Double, Double, Double, Double, Double, Double)],
140              maxRanges: Option[Int] = None): Seq[IndexRange] = {
141     // normalize inputs to [0,1]
142     val windows = queries.map { case (xmin, ymin, zmin, xmax, ymax, zmax) =>
143       val (nxmin, nymin, nzmin, nxmax, nymax, nzmax) = normalize(xmin, ymin, zmin, xmax, ymax, zmax, lenient = false)
144       QueryWindow(nxmin, nymin, nzmin, nxmax, nymax, nzmax)
145     }
146     ranges(windows.toArray, maxRanges.getOrElse(Int.MaxValue))
147   }
148 
149   /**
150     * Determine XZ-curve ranges that will cover a given query window
151     *
152     * @param query a sequence of OR'd windows to cover, normalized to [0,1]
153     * @param rangeStop a rough max value for the number of ranges to return
154     * @return
155     */
156   private def ranges(query: Array[QueryWindow], rangeStop: Int): Seq[IndexRange] = {
157 
158     import XZ3SFC.{LevelOneElements, LevelTerminator}
159 
160     // stores our results - initial size of 100 in general saves us some re-allocation
161     val ranges = new java.util.ArrayList[IndexRange](100)
162 
163     // values remaining to process - initial size of 100 in general saves us some re-allocation
164     val remaining = new java.util.ArrayDeque[XElement](100)
165 
166     // checks if a quad is contained in the search space
167     def isContained(oct: XElement): Boolean = {
168       var i = 0
169       while (i < query.length) {
170         if (oct.isContained(query(i))) {
171           return true
172         }
173         i += 1
174       }
175       false
176     }
177 
178     // checks if a quad overlaps the search space
179     def isOverlapped(oct: XElement): Boolean = {
180       var i = 0
181       while (i < query.length) {
182         if (oct.overlaps(query(i))) {
183           return true
184         }
185         i += 1
186       }
187       false
188     }
189 
190     // checks a single value and either:
191     //   eliminates it as out of bounds
192     //   adds it to our results as fully matching, or
193     //   adds it to our results as partial matching and queues up it's children for further processing
194     def checkValue(oct: XElement, level: Short): Unit = {
195       if (isContained(oct)) {
196         // whole range matches, happy day
197         val (min, max) = sequenceInterval(oct.xmin, oct.ymin, oct.zmin, level, partial = false)
198         ranges.add(IndexRange(min, max, contained = true))
199       } else if (isOverlapped(oct)) {
200         // some portion of this range is excluded
201         // add the partial match and queue up each sub-range for processing
202         val (min, max) = sequenceInterval(oct.xmin, oct.ymin, oct.zmin, level, partial = true)
203         ranges.add(IndexRange(min, max, contained = false))
204         oct.children.foreach(remaining.add)
205       }
206     }
207 
208     // initial level
209     LevelOneElements.foreach(remaining.add)
210     remaining.add(LevelTerminator)
211 
212     // level of recursion
213     var level: Short = 1
214 
215     while (level < g && !remaining.isEmpty && ranges.size < rangeStop) {
216       val next = remaining.poll
217       if (next.eq(LevelTerminator)) {
218         // we've fully processed a level, increment our state
219         if (!remaining.isEmpty) {
220           level = (level + 1).toShort
221           remaining.add(LevelTerminator)
222         }
223       } else {
224         checkValue(next, level)
225       }
226     }
227 
228     // bottom out and get all the ranges that partially overlapped but we didn't fully process
229     while (!remaining.isEmpty) {
230       val oct = remaining.poll
231       if (oct.eq(LevelTerminator)) {
232         level = (level + 1).toShort
233       } else {
234         val (min, max) = sequenceInterval(oct.xmin, oct.ymin, oct.zmin, level, partial = false)
235         ranges.add(IndexRange(min, max, contained = false))
236       }
237     }
238 
239     // we've got all our ranges - now reduce them down by merging overlapping values
240     // note: we don't bother reducing the ranges as in the XZ paper, as accumulo handles lots of ranges fairly well
241     ranges.sort(IndexRange.IndexRangeIsOrdered)
242 
243     var current = ranges.get(0) // note: should always be at least one range
244     val result = ArrayBuffer.empty[IndexRange]
245     var i = 1
246     while (i < ranges.size()) {
247       val range = ranges.get(i)
248       if (range.lower <= current.upper + 1) {
249         // merge the two ranges
250         current = IndexRange(current.lower, math.max(current.upper, range.upper), current.contained && range.contained)
251       } else {
252         // append the last range and set the current range for future merging
253         result.append(current)
254         current = range
255       }
256       i += 1
257     }
258     // append the last range - there will always be one left that wasn't added
259     result.append(current)
260 
261     result.toSeq
262   }
263 
264   /**
265     * Computes the sequence code for a given point - for polygons this is the lower-left corner.
266     *
267     * Based on Definition 2 from the XZ-Ordering paper
268     *
269     * @param x normalized x value [0,1]
270     * @param y normalized y value [0,1]
271     * @param z normalized z value [0,1]
272     * @param length length of the sequence code that will be generated
273     * @return
274     */
275   private def sequenceCode(x: Double, y: Double, z: Double, length: Int): Long = {
276     var xmin = 0.0
277     var ymin = 0.0
278     var zmin = 0.0
279     var xmax = 1.0
280     var ymax = 1.0
281     var zmax = 1.0
282 
283     var cs = 0L
284 
285     var i = 0
286     while (i < length) {
287       val xCenter = (xmin + xmax) / 2.0
288       val yCenter = (ymin + ymax) / 2.0
289       val zCenter = (zmin + zmax) / 2.0
290       (x < xCenter, y < yCenter, z < zCenter) match {
291         case (true,  true, true)   => cs += 1L                                             ; xmax = xCenter; ymax = yCenter; zmax = zCenter
292         case (false, true, true)   => cs += 1L + 1L * (math.pow(8, g - i).toLong - 1L) / 7L; xmin = xCenter; ymax = yCenter; zmax = zCenter
293         case (true,  false, true)  => cs += 1L + 2L * (math.pow(8, g - i).toLong - 1L) / 7L; xmax = xCenter; ymin = yCenter; zmax = zCenter
294         case (false, false, true)  => cs += 1L + 3L * (math.pow(8, g - i).toLong - 1L) / 7L; xmin = xCenter; ymin = yCenter; zmax = zCenter
295         case (true,  true, false)  => cs += 1L + 4L * (math.pow(8, g - i).toLong - 1L) / 7L; xmax = xCenter; ymax = yCenter; zmin = zCenter
296         case (false, true, false)  => cs += 1L + 5L * (math.pow(8, g - i).toLong - 1L) / 7L; xmin = xCenter; ymax = yCenter; zmin = zCenter
297         case (true,  false, false) => cs += 1L + 6L * (math.pow(8, g - i).toLong - 1L) / 7L; xmax = xCenter; ymin = yCenter; zmin = zCenter
298         case (false, false, false) => cs += 1L + 7L * (math.pow(8, g - i).toLong - 1L) / 7L; xmin = xCenter; ymin = yCenter; zmin = zCenter
299       }
300       i += 1
301     }
302 
303     cs
304   }
305 
306   /**
307     * Computes an interval of sequence codes for a given point - for polygons this is the lower-left corner.
308     *
309     * @param x normalized x value [0,1]
310     * @param y normalized y value [0,1]
311     * @param length length of the sequence code that will used as the basis for this interval
312     * @param partial true if the element partially intersects the query window, false if it is fully contained
313     * @return
314     */
315   private def sequenceInterval(x: Double, y: Double, z: Double, length: Short, partial: Boolean): (Long, Long) = {
316     val min = sequenceCode(x, y, z, length)
317     // if a partial match, we just use the single sequence code as an interval
318     // if a full match, we have to match all sequence codes starting with the single sequence code
319     val max = if (partial) { min } else {
320       // from lemma 3 in the XZ-Ordering paper
321       min + (math.pow(8, g - length + 1).toLong - 1L) / 7L
322     }
323     (min, max)
324   }
325 
326   /**
327     * Normalize user space values to [0,1]
328     *
329     * @param xmin min x value in user space
330     * @param ymin min y value in user space
331     * @param zmin min z value in user space
332     * @param xmax max x value in user space, must be >= xmin
333     * @param ymax max y value in user space, must be >= ymin
334     * @param zmax max z value in user space, must be >= zmin
335     * @param lenient standardize boundaries to valid values, or raise an exception
336     * @return
337     */
338   protected def normalize(
339       xmin: Double,
340       ymin: Double,
341       zmin: Double,
342       xmax: Double,
343       ymax: Double,
344       zmax: Double,
345       lenient: Boolean): (Double, Double, Double, Double, Double, Double) = {
346     require(xmin <= xmax && ymin <= ymax && zmin <= zmax,
347       s"Bounds must be ordered: [$xmin $xmax] [$ymin $ymax] [$zmin $zmax]")
348 
349     try {
350       require(xmin >= xLo && xmax <= xHi && ymin >= yLo && ymax <= yHi && zmin >= zLo && zmax <= zHi,
351         s"Values out of bounds ([$xLo $xHi] [$yLo $yHi] [$zLo $zHi]): [$xmin $xmax] [$ymin $ymax] [$zmin $zmax]")
352 
353       val nxmin = (xmin - xLo) / xSize
354       val nymin = (ymin - yLo) / ySize
355       val nzmin = (zmin - zLo) / zSize
356       val nxmax = (xmax - xLo) / xSize
357       val nymax = (ymax - yLo) / ySize
358       val nzmax = (zmax - zLo) / zSize
359 
360       (nxmin, nymin, nzmin, nxmax, nymax, nzmax)
361     } catch {
362       case _: IllegalArgumentException if lenient =>
363 
364         val bxmin = if (xmin < xLo) { xLo } else if (xmin > xHi) { xHi } else { xmin }
365         val bymin = if (ymin < yLo) { yLo } else if (ymin > yHi) { yHi } else { ymin }
366         val bzmin = if (zmin < zLo) { zLo } else if (zmin > zHi) { zHi } else { zmin }
367         val bxmax = if (xmax < xLo) { xLo } else if (xmax > xHi) { xHi } else { xmax }
368         val bymax = if (ymax < yLo) { yLo } else if (ymax > yHi) { yHi } else { ymax }
369         val bzmax = if (zmax < zLo) { zLo } else if (zmax > zHi) { zHi } else { zmax }
370 
371         val nxmin = (bxmin - xLo) / xSize
372         val nymin = (bymin - yLo) / ySize
373         val nzmin = (bzmin - zLo) / zSize
374         val nxmax = (bxmax - xLo) / xSize
375         val nymax = (bymax - yLo) / ySize
376         val nzmax = (bzmax - zLo) / zSize
377 
378         (nxmin, nymin, nzmin, nxmax, nymax, nzmax)
379     }
380   }
381 }
382 
383 object XZ3SFC {
384 
385   // the initial level of octs
386   private val LevelOneElements = XElement(0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0).children
387 
388   // indicator that we have searched a full level of the oct tree
389   private val LevelTerminator = XElement(-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.0)
390 
391   private val cache = new java.util.concurrent.ConcurrentHashMap[(Short, TimePeriod), XZ3SFC]()
392 
393   def apply(g: Short, period: TimePeriod): XZ3SFC = {
394     var sfc = cache.get((g, period))
395     if (sfc == null) {
396       sfc = new XZ3SFC(g, (-180.0, 180.0), (-90.0, 90.0), (0.0, BinnedTime.maxOffset(period).toDouble))
397       cache.put((g, period), sfc)
398     }
399     sfc
400   }
401 
402   /**
403     * Region being queried. Bounds are normalized to [0-1].
404     *
405     * @param xmin x lower bound in [0-1]
406     * @param ymin y lower bound in [0-1]
407     * @param zmin z lower bound in [0-1]
408     * @param xmax x upper bound in [0-1], must be >= xmin
409     * @param ymax y upper bound in [0-1], must be >= ymin
410     * @param zmax z upper bound in [0-1], must be >= zmin
411     */
412   private case class QueryWindow(xmin: Double, ymin: Double, zmin: Double, xmax: Double, ymax: Double, zmax: Double)
413 
414   /**
415     * An extended Z curve element. Bounds refer to the non-extended z element for simplicity of calculation.
416     *
417     * An extended Z element refers to a normal Z curve element that has it's upper bounds expanded by double it's
418     * width/length/height. By convention, an element is always a cube.
419     *
420     * @param xmin x lower bound in [0-1]
421     * @param ymin y lower bound in [0-1]
422     * @param zmin z lower bound in [0-1]
423     * @param xmax x upper bound in [0-1], must be >= xmin
424     * @param ymax y upper bound in [0-1], must be >= ymin
425     * @param zmax z upper bound in [0-1], must be >= zmin
426     * @param length length of the non-extended side (note: by convention width should be equal to height and depth)
427     */
428   private case class XElement(xmin: Double,
429                               ymin: Double,
430                               zmin: Double,
431                               xmax: Double,
432                               ymax: Double,
433                               zmax: Double,
434                               length: Double) {
435 
436     // extended x and y bounds
437     lazy val xext = xmax + length
438     lazy val yext = ymax + length
439     lazy val zext = zmax + length
440 
441     def isContained(window: QueryWindow): Boolean =
442       window.xmin <= xmin && window.ymin <= ymin && window.zmin <= zmin &&
443           window.xmax >= xext && window.ymax >= yext && window.zmax >= zext
444 
445     def overlaps(window: QueryWindow): Boolean =
446       window.xmax >= xmin && window.ymax >= ymin && window.zmax >= zmin &&
447           window.xmin <= xext && window.ymin <= yext && window.zmin <= zext
448 
449     def children: Seq[XElement] = {
450       val xCenter = (xmin + xmax) / 2.0
451       val yCenter = (ymin + ymax) / 2.0
452       val zCenter = (zmin + zmax) / 2.0
453       val len = length / 2.0
454       val c0 = copy(xmax = xCenter, ymax = yCenter, zmax = zCenter, length = len)
455       val c1 = copy(xmin = xCenter, ymax = yCenter, zmax = zCenter, length = len)
456       val c2 = copy(xmax = xCenter, ymin = yCenter, zmax = zCenter, length = len)
457       val c3 = copy(xmin = xCenter, ymin = yCenter, zmax = zCenter, length = len)
458       val c4 = copy(xmax = xCenter, ymax = yCenter, zmin = zCenter, length = len)
459       val c5 = copy(xmin = xCenter, ymax = yCenter, zmin = zCenter, length = len)
460       val c6 = copy(xmax = xCenter, ymin = yCenter, zmin = zCenter, length = len)
461       val c7 = copy(xmin = xCenter, ymin = yCenter, zmin = zCenter, length = len)
462       Seq(c0, c1, c2, c3, c4, c5, c6, c7)
463     }
464   }
465 }