Lecture 11: Orthogonal Range Searching

Reading: Chapter 5 in the 4M's.

Range Queries: 

We shift our focus from algorithm problems to data structures for the
next few lectures. In general, we consider the question, given a
collection of objects, preprocess them (storing the results in a data
structure of some variety) so that queries of a particular form can be
answered efficiently. Generally we measure data structures in terms of
two quantities, the time needed to answer a query and the amount of
space needed by the data structure.

Often there is a tradeoff between these two quantities, but most of
the structures that we will be interested in will have either linear
or near linear space. Preprocessing time is an issue of secondary
importance, but most of the algorithms we will consider will have
either linear or O(n log n) preprocessing time.

In a range queries we are given a set P of points and region R in
space (e.g., a rectangle, polygon, halfspace, or disk) and are asked
list (or count or compute some accumulation function of) the subset of
P lying within the region. To get the best possible performance, the
design of the data structure is tailored to the particular type of
region, but there are some data structures that can be used for a wide
variety of regions.

An important concept behind all geometric range searching is that the
subsets that can be formed by simple geometric ranges is much smaller
than the set of possible subsets (called the power set) of P.  

  (Cool fact that you will not be tested on, and is only marginally
   relevant: There is something called the vapnik-chernovenkis
   dimension of a "shape" which asks: what is the maximum number N
   such that you can place N points in the plane and seperate any
   possible subset of those points using (scaled or rotated versions
   of ) that shape?  ie, if you put 3 points on the plane, you can
   choose any subset of those three points using an appropriately
   placed half-plane.  But you can't do the same thing for four

We can define any range search problem abstractly as follows. Given a
particular class of ranges, a range space is a pair (P, R) consisting
of the points P and the collection R of all subsets of P that be
formed by ranges of this class. For example, the following figure
shows the range space assuming rectangular ranges for a set of points
in the plane. In particular, note that the sets {1, 4} and {1, 2, 4}
cannot be formed by rectangular ranges.

Figure 42: Rectangular range space. Today we consider orthogonal rectangular range queries, that is, ranges defined by rectangles whose sides are aligned with the coordinate axes. One of the nice things about rectangular ranges is that they can be decomposed into a collection 1 dimensional searches. One dimensional range queries: Before consider how to solve general range queries, let us consider how to answer 1 dimensional range queries, or interval queries. We are given a set of points P = {p(1), p(2), ... p(n)} on the line, and given an interval [x(lo), x(hi)], report all the points lying within the interval. Our goal is to find a data structure that can answer these queries in O(log(n) + k) time, where k is the number of points reported (an output sensitive result). [Why is it important that this be an output sensitive result?] Naive Solutions: What if we don't pre-process our input at all? Then for an input [x(lo), x(hi)] we can take O(n) time to check to see if each point p(i) is in [x(lo), x(hi)], and report it or count it if it is. Ok ok, so sort the points first. apply binary search to find the first point of P that is greater than or equal to x(lo) , and less than or equal to x(hi), and then list all the points between. So, range counting queries can be answered in O(log n) time, with minor modifications, with pre-processing consisting of sorting the sort the points. However, this will not generalize to higher dimensions. (ie, what is "sorting the points" when the points are not one-dimensional?). A basic approach to solving almost all range queries is to represent P as a collection of canonical subsets {S(1), S(2), ... S(k)}, each S(i) is a subset of the set S (where k is generally a function of n and the type of ranges), such that the answer to any query can be expressed as a disjoint union of a small number of canonical subsets. Note that these subsets may generally overlap each other. The trick to solving a range searching problem is to identify a collection of canonical subsets having these properties, and then finding the proper subsets for a given range. Suppose that the points of P are sorted in increasing order and then stored in the leaves of a balanced binary tree. Each internal node of the tree is labeled with the largest key appearing in its left child. We can associate each node of this tree (implicitly or explicitly) with the subset of points stored in the leaves that are descendents of this node. This gives rise to the O(n) canonical subsets. For now, these canonical subsets will not be stored explicitly as part of the data structure, but this will change later when we talk about range trees. We claim that the canonical subsets corresponding to any range can be identified in O(log n) time from this structure. Given any interval [x(lo), x(hi)], we search the tree to find the leftmost leaf u whose key is greater than or equal to x(lo) and the rightmost leaf v whose key is greater than or equal to x(hi). Clearly all the leaves between u and v, together possibly with u and v, constitute the points that lie within the range. If key(u) = x(lo) then we include u's canonical (single point) subset and if key(v) = x(hi) then we do the same for v. To form the remaining canonical subsets, we take the subsets of all the maximal subtrees lying between u and v.
Figure 43: Canonical sets for interval queries. Here is how to compute these subtrees. The search paths to u and v may generally share some common subpath, starting at the root of the tree. Once the paths diverge, as we follow the left path to u, whenever the path goes to the left child of some node, we add the canonical subset associated with its right child. Similarly, as we follow the right path to v, whenever the path goes to the right child, we add the canonical subset associated with its left child. To answer a range reporting query we simply traverse these canonical subtrees, reporting the points of their leaves. Each tree can be traversed in time proportional to the number of leaves in each subtree. To answer a range counting query we store the total number of points in each subtree (as part of the preprocessing) and then sum all of these over all the canonical subtrees. Since the search paths are of length O(log n), it follows that there are O(log n) total canonical subsets. Thus the range counting query can be answered in O(log n) time. For reporting queries, since the leaves of each subtree can be listed in time that is proportional to the number of leaves in the tree (a basic fact about binary trees), it follows that the total time in the search is O(log n + k), where k is the number of points reported. Thus, 1 dimensional range queries can be answered in O(logn) time, using O(n) storage. This concept of finding maximal subtrees that are contained within the range is fundamental to all range search data structures. The only question is how to organize the tree and how to locate the desired sets. Let see next how can we extend this to higher dimensional range queries. Kd trees: The natural question is how to extend 1 dimensional range searching to higher dimen sions. First we will consider kd trees. This data structure is easy to implement and quite practical and useful for many different types of searching problems (nearest neighbor search ing for example). However it is not the asymptotically most efficient solution for the orthogonal range searching, as we will see later. Our terminology is a bit nonstandard. The data structure was designed by Jon Bentley. In his notation, these were called k d trees, short for ``k dimensional trees''. The value k was the dimension, and thus there are 2 d trees, 3 d trees, and so on. However, over time, the term kd-tree has become so standard that we simply drop the pretense of using k as a variable. The idea behind a kd tree is to extend the notion of a one dimension tree, but alternate in using the x or y coordinates to split on. In general dimension, the kd-tree cycles among the various possible splitting dimensions, but you might imagine many variations where the choice of the splitting dimension depends on the structure of the data set. Each internal node of the kd-tree is associated with two quantities, a splitting dimension (either x or y), and a splitting value s. It has two subtrees. Since the splitting dimension alternates between x and y, some implementations do not store this explicitly, but keep track of it while traversing the tree. If the splitting dimension is x, then all points whose x coordinates are less than or equal to s are stored in the left subtree and points whose x coordinates are greater than or equal to s are stored in the right subtree. (If a point's coordinate is equal to s, then we reserve the right to store it on either side. This is done to allow us to balance the number of points in the left and right subtrees.) When a single point (or more generally a small constant number of points) remains, we store it in a leaf node. Note that in our definition the data points are stored only in the leaves of the tree. The internal nodes only contain splitting information. How is the splitting value chosen? To guarantee that the tree is balanced, the most common method is to let the splitting value be the median splitting coordinate. If there are an even number of points in the subtree, we may take either the upper or lower median, or we may simply take the midpoint between these two points. As a result, the tree will have O(log n) height. How is the splitting dimension chosen? As mentioned earlier one simple strategy is to cycle through the dimensions one by one. Another is to look at the data points, and select the dimension along which they have the greatest spread, defined to be the difference between the largest and smallest coordinates. Bentley refered to the tree resulting from this splitting rule to be the optimized kd-tree. In the example below, we simply alternate splitting dimensions. If there are a odd number of points, the median value goes with the left (or lower) subtree.
Figure 44: A kd tree and the associated spatial subdivision. A kd-tree naturally defines a subdivision of space, where each splitting node can be thought of as introducing a vertical or horizontal splitting line. Points lying to one side of the line are stored in the left subtree and the other points are stored in the right subtree. This subdivides space into regions, called cells, which are (possibly unbounded) rectangles. In fact this is a special case of a more general class of data structures, called binary space partition or BSP trees, in which the splitting lines may be oriented in any direction. In the case of BSP trees, the cells are convex polygons. Of course, both data structures can be generalized to higher dimensions, where splitting is done using d-1 dimensional hyperplanes. It is possible to build a kd tree in O(n log n) time by a simple recursive procedure. The most costly step of the process is determining the median coordinate. One trick for speeding up processing time is to presort the points into two cross referenced lists, one sorted by x and the other sorted by y. Using these lists, it is an easy matter to find the median at each step in constant time. The two lists can then be split in O(n) time each, where n is the number of remaining points. This leads to a recurrence of the form T (n) = 2 T(n/2) + n, which solves to O(n log n). KD-tree construction applet Searching the kd tree: To answer an orthogonal range query we proceed as follows. Let R denote the query rectangle. Assume that the current splitting line is vertical. (The horizontal case is similar.) Let v denote the current node in the tree. Each node of the tree is associated with a rectangular region of space (its cell). The search proceeds as follows. If v is a leaf, then we check the point(s) stored in v as to whether it lies in R. If so we report/count them. If v is an internal node, then we first consider the left subtree. If its region lies entirely within R then we call a different routine to enumerate all the points in this subtree (or for a counting query we return a precomputed count of the number of points in this subtree). If the left child's region partially overlaps R then we search it recursively. (If the left child's region does not overlap R at all, then we ignore it.) We do the same for the right child of v. The figure below illustrates the nodes that would be visited for the given query range. The leaf nodes with heavy outlines are the points that are finally reported.
Figure 45: A kd tree and the associated spatial subdivision. How many nodes does this method visit altogether? We claim that the total number is O( sqrt(n) + k), where k is the number of points reported. Theorem: The times needed to answer an orthogonal range quep ry using kd tree with n points is O(sqrt(n) + k) where k is the number of reported points. Proof: First observe that for each subtree with m points, we can enumerate the points of this subtree in O(m) time, by a simple traversal. Thus, it suffices to bound the number of nodes visited in the search (not counting subtree enumeration). To count the number of nodes visited, first observe that if we visit a node, then its associated region must intersect one of the four sides of the rectangle range. (Otherwise it lies entirely inside or outside, and so is not visited.) We will do this separately for each of the four sides, and multiply the result by four. Rather than consider the line segment defining a side of the range rectangle, we consider the number of nodes visited if we had applied an orthogonal range query to the entire line on which this side lies. This will only overestimate the number of nodes visited. Because the kd-tree processes even and odd levels differently, it will be important to analyze two levels at a time (or generally d levels at a time for dimension d). Let Q(n) denote the query time for visiting a node v that has n points descended from it. Consider any orthogonal line that intersects the region associated with v. The key observation is that any horizontal or vertical line can intersect at most two of the four subregions associated with the grandchildren of v. (The other two will either be entirely contained within the halfspace or lie entirely outside the halfspace. In either case the algorithm will not make a recursive call on them.) Since each child has half as many points (because we split at the median), the number of points in each grandchild is roughly n/4. Thus we make at most two recursive calls on subtrees of size n/4. This gives the following recurrence: Q(n) = 1 if n = 1, 2Q(n/4) + 1 otherwise. It is an easy process to expand this recursion, and derive the fact that it solves to O(n * (log base 4 of 2 )) = O(sqrt(n)). (This also follows from the Master Theorem in CLR.) Including the factor of four, we still have an O(sqrt(n)) bound on the total number of nodes visited, and adding the enumeration time the total time is O(sqrt(n) + k) as desired. Other spatial data structures links. A Quad-Tree http://www.cs.umd.edu/~brabec/quadtree/index.html A Binary Space Partition Tree http://symbolcraft.com/graphics/bsp/index.html