Algorithm Alley
by Ron Gutman


Example 1:

(a) 
SELECT * FROM homes WHERE
  price > 140000 AND price < 150000

(b)
SELECT * FROM homes WHERE
  longitude > 122.3 AND longitude < 122.4
  AND latitude > 37.5 AND latitude < 37.6


Listing One
/* space.h - Header definitions for computing spatial keys using */
/* space-filling curves.                                         */

typedef struct {
    unsigned long begin;
    unsigned long end;
    } Range;
typedef struct {
    long x;
    long y;
    } Point;
typedef struct {
     Point lowerLeft;  /* lower left corner of the box */
     Point upperRight; /* upper right corner of the box */
     } Box;
#define DEPTH 8        /* values up to 15 supported */
Box world = {{0,0},{2000000,2000000}};
int queryRanges (Box *pQuery, int depth, Range *pRanges, int maxRanges);


Listing Two
/* space.c - C code for computing spatial keys using space-filling curves. */

#include "drdobbs.h"
typedef struct {
    Box   box;
    char  depth;
    long  keyPrefix;
#ifdef HILBERT
    char  rotation;
#endif
    } Cell;
typedef struct {
    Range *pRanges;
    int   nextRange;  /* index of next range to be stored */
    int   maxRanges;
    } RangeList;
/* Returns 1 if the two boxes overlap                     */
/* Returns 0 otherwise.                                   */
char overlap (Box *pBox1, Box *pBox2) {
    if (pBox1->lowerLeft.x > pBox2->upperRight.x)
        return 0; /* box1 is entirely to the right of box2 */
    if (pBox1->lowerLeft.y > pBox2->upperRight.y)
        return 0; /* box1 is entriely above box2 */    
    if (pBox1->upperRight.x < pBox2->lowerLeft.x)
        return 0; /* box1 is entirely to the left of box2 */
    if (pBox1->upperRight.y < pBox2->lowerLeft.y)
        return 0; /* box1 is entriely below box2 */
    /* If none of the above are true then box1 must  overlap with box2. */
    return 1;
    }
/* Add a range to a list of ranges */
int addRange (RangeList *pRangeList, Range *pRange) {
    int nr;  /* index of next range */
    nr = pRangeList->nextRange;
    if (nr >= pRangeList->maxRanges)
       return -1;
    /* if this range and the last range are adjacent, combine them */
    if ( (nr>0) && (pRange->begin==pRangeList->pRanges[nr-1].end+1) ) 
       {
       pRangeList->pRanges[nr-1].end = pRange->end;
       return 0;
       }
    pRangeList->pRanges[nr] = *pRange;
    pRangeList->nextRange++;
    }
#ifdef HILBERT
/* Gives position of child cell in its parent for Hilbert */
/* order given the orientation of parent and position of child */
/* its parent for Morton order. */
static int hilbertIndexTable [4][4] = {
               {0,3,1,2}, /* right side up U */
               {2,1,3,0}, /* upside down U */
               {0,1,3,2}, /* upside left U */
               {2,3,1,0}  /* upside right U */
               };
/* Gives the orientation of a child cell given the orientation of   */
/* the parent cell and the position of the child within the parent. */
static int hilbertRotationTable [4][4] = {
               {2,0,0,3}, /* right side up U */
               {3,1,1,2}, /* upside down U */
               {0,2,2,1}, /* upside left U */
               {1,3,3,0}  /* upside right U */
               };
static int hilbertIndex(Cell *pCell, int zIndex, Cell *pParent) {
   int hIndex;
   hIndex = hilbertIndexTable[pParent->rotation][zIndex];
   pCell->rotation = hilbertRotationTable[pParent->rotation][hIndex];
   return hIndex;
   }
#endif
/* Recursive function that explores the cell hierarchy beneath   */
/* the given cell to find cells overlapping the query region.    */
/* Returns -1 if the buffer for ranges is exhausted.             */
int subtreeRanges (Box *pQuery, int depth, RangeList *pRangeList,
                   Cell *pCell) {
    Cell  childCell[4];
    Point mid;
    char  i,j;
    Range range;
#ifdef HILBERT
    char hilbertOrder[4];
#endif
    if ( !overlap(pQuery, &(pCell->box)) )
        return 0;
    if (pCell->depth >= depth) {
        /* emit range for this cell */
        range.begin = pCell->keyPrefix << ((DEPTH-pCell->depth)*2);
        range.end = range.begin + (1 << ((DEPTH-pCell->depth)*2)) - 1;
        return addRange (pRangeList, &range);          
        }
    /* create child cells of this cell */
    mid.x = (pCell->box.lowerLeft.x >> 1) + (pCell->box.upperRight.x >> 1);
    mid.y = (pCell->box.lowerLeft.y >> 1) + (pCell->box.upperRight.y >> 1);
    childCell[0].box.lowerLeft.x = pCell->box.lowerLeft.x;
    childCell[0].box.lowerLeft.y = mid.y;
    childCell[0].box.upperRight.x = mid.x;
    childCell[0].box.upperRight.y = pCell->box.upperRight.y;
    childCell[1].box.lowerLeft = mid;
    childCell[1].box.upperRight = pCell->box.upperRight;
    childCell[2].box.lowerLeft = pCell->box.lowerLeft;
    childCell[2].box.upperRight = mid;
    childCell[3].box.lowerLeft.x = mid.x;
    childCell[3].box.lowerLeft.y = pCell->box.lowerLeft.y;
    childCell[3].box.upperRight.x = pCell->box.upperRight.x;
    childCell[3].box.upperRight.y = mid.y;
    for (i = 0; i < 4; i++)
        childCell[i].depth = pCell->depth + 1;
#ifdef HILBERT
    for (i = 0; i < 4; i++) 
        hilbertOrder[hilbertIndex(childCell+i,i,pCell)] = i;
#endif
    /* visit children in space-filling curve order */
    for (i = 0; i < 4; i++) {
#ifdef HILBERT
       j = hilbertOrder[i];
#else
       j = i;
#endif
       childCell[j].keyPrefix = (pCell->keyPrefix<<2) + i;
       if ( subtreeRanges(pQuery,depth,pRangeList,childCell+j) < 0 )
          return -1;
       }
    return 0;
    }                                                       
/* This is the main function.  The real work is done by              */
/* subtreeRanges, but this function sets up data structures to set   */
/* the recursion in motion.  Returns the number of ranges stored at  */
/* pRanges. The depth of recursion is controlled by depth which must */
/* be <= DEPTH.  The number of ranges stored at pRanges is limited   */
/* to maxRanges.                                                     */
int queryRanges (Box *pQuery, int depth, Range *pRanges, int maxRanges) {
    Cell rootCell;
    RangeList rangeList;
    rootCell.box = world;
    rootCell.depth = 0;
    rootCell.keyPrefix = 0;
#ifdef HILBERT
    rootCell.rotation = 0;
#endif
    rangeList.pRanges = pRanges;
    rangeList.nextRange = 0;
    rangeList.maxRanges =  maxRanges;
    if (subtreeRanges (pQuery, depth, &rangeList, &rootCell) < 0)
       return -1;
    return rangeList.nextRange;
    }


4


