diff --git a/index/rtree/rtree.go b/index/rtree/rtree.go index 590165a7..53f41c8f 100644 --- a/index/rtree/rtree.go +++ b/index/rtree/rtree.go @@ -1,22 +1,5 @@ -// Package rtree - A 2d Implementation of RTree, a bounding rectangle tree. -// -// This file is derived from the work done by Toni Gutman. R-Trees: A Dynamic Index Structure for -// Spatial Searching, Proc. 1984 ACM SIGMOD International Conference on Management of Data, pp. -// 47-57. The implementation found in SQLite is a refinement of Guttman's original idea, commonly -// called "R*Trees", that was described by Norbert Beckmann, Hans-Peter Kriegel, Ralf Schneider, -// Bernhard Seeger: The R*-Tree: An Efficient and Robust Access Method for Points and Rectangles. -// SIGMOD Conference 1990: 322-331 -// -// The original C code can be found at "http://www.superliminal.com/sources/sources.htm". -// -// And the website carries this message: "Here are a few useful bits of free source code. You're -// completely free to use them for any purpose whatsoever. All I ask is that if you find one to -// be particularly valuable, then consider sending feedback. Please send bugs and suggestions too. -// Enjoy" package rtree -import "math" - // Item is an rtree item type Item interface { Rect() (minX, minY, maxX, maxY float64) @@ -32,626 +15,58 @@ func (item *Rect) Rect() (minX, minY, maxX, maxY float64) { return item.MinX, item.MinY, item.MaxX, item.MaxY } -func min(a, b float64) float64 { - if a < b { - return a - } - return b -} -func max(a, b float64) float64 { - if a > b { - return a - } - return b -} - -const ( - unitSphereVolume1 = 2.000000 - unitSphereVolume2 = 3.141593 - unitSphereVolume3 = 4.188790 - unitSphereVolume4 = 4.934802 -) -const ( - maxNodes = 16 - minNodes = maxNodes / 2 - useSphericalVolume = true - unitSphereVolume = unitSphereVolume2 -) - -/// Minimal bounding rectangle (n-dimensional) -type rectT struct { - min [2]float64 ///< Min dimensions of bounding box - max [2]float64 ///< Max dimensions of bounding box -} - -/// May be data or may be another subtree -/// The parents level determines this. -/// If the parents level is 0, then this is data -type branchT struct { - rect rectT ///< Bounds - child *nodeT ///< Child node - item Item ///< Data ID or Ptr -} - -/// nodeT for each branch level -type nodeT struct { - count int ///< Count - level int ///< Leaf is zero, others positive - branch [maxNodes]branchT ///< branchT -} - -func (node *nodeT) isInternalNode() bool { return node.level > 0 } // Not a leaf, but a internal node - -/// A link list of nodes for reinsertion after a delete operation -type listNodeT struct { - next *listNodeT ///< Next in list - node *nodeT ///< nodeT -} - -/// Variables for finding a split partition -type partitionVarsT struct { - partition [maxNodes + 1]int - total int - minFill int - taken [maxNodes + 1]bool - count [2]int - cover [2]rectT - area [2]float64 - branchBuf [maxNodes + 1]branchT - branchCount int - coverSplit rectT - coverSplitArea float64 -} - // RTree is an implementation of an rtree type RTree struct { - root *nodeT -} - -func itemRect(item Item) (rect rectT) { - minX, minY, maxX, maxY := item.Rect() - return rectT{ - min: [2]float64{minX, minY}, - max: [2]float64{maxX, maxY}, - } + tr *d2RTree } // New creates a new RTree func New() *RTree { - return &RTree{} + return &RTree{ + tr: d2New(), + } } // Insert inserts item into rtree func (tr *RTree) Insert(item Item) { - if tr.root == nil { - tr.root = &nodeT{} - } - insertRect(itemRect(item), item, &tr.root, 0) + minX, minY, maxX, maxY := item.Rect() + tr.tr.Insert([2]float64{minX, minY}, [2]float64{maxX, maxY}, item) } // Remove removes item from rtree func (tr *RTree) Remove(item Item) { - if tr.root == nil { - tr.root = &nodeT{} - } - removeRect(itemRect(item), item, &tr.root) + minX, minY, maxX, maxY := item.Rect() + tr.tr.Remove([2]float64{minX, minY}, [2]float64{maxX, maxY}, item) } // Search finds all items in bounding box. func (tr *RTree) Search(minX, minY, maxX, maxY float64, iterator func(item Item) bool) { - if iterator == nil { - return - } - rect := rectT{ - min: [2]float64{minX, minY}, - max: [2]float64{maxX, maxY}, - } - // NOTE: May want to return search result another way, perhaps returning the number of found elements here. - if tr.root == nil { - tr.root = &nodeT{} - } - search(tr.root, rect, iterator) + tr.tr.Search([2]float64{minX, minY}, [2]float64{maxX, maxY}, func(data interface{}) bool { + return iterator(data.(Item)) + }) } // Count return the number of items in rtree. func (tr *RTree) Count() int { - return countRec(tr.root, 0) + return tr.tr.Count() } // RemoveAll removes all items from rtree. func (tr *RTree) RemoveAll() { - tr.root = nil + tr.tr.RemoveAll() } func (tr *RTree) Bounds() (minX, minY, maxX, maxY float64) { - var rect rectT - if tr.root != nil { - if tr.root.count > 0 { - rect = tr.root.branch[0].rect - for i := 1; i < tr.root.count; i++ { - rect = combineRect(rect, tr.root.branch[i].rect) + var rect d2rectT + if tr.tr.root != nil { + if tr.tr.root.count > 0 { + rect = tr.tr.root.branch[0].rect + for i := 1; i < tr.tr.root.count; i++ { + rect2 := tr.tr.root.branch[i].rect + rect = d2combineRect(&rect, &rect2) } } } minX, minY, maxX, maxY = rect.min[0], rect.min[1], rect.max[0], rect.max[1] return } - -func countRec(node *nodeT, counter int) int { - if node.isInternalNode() { // not a leaf node - for index := 0; index < node.count; index++ { - counter = countRec(node.branch[index].child, counter) - } - } else { // A leaf node - if node.count > 256 { - println(node.count) - } - counter += node.count - } - return counter -} - -// Inserts a new data rectangle into the index structure. -// Recursively descends tree, propagates splits back up. -// Returns 0 if node was not split. Old node updated. -// If node was split, returns 1 and sets the pointer pointed to by -// new_node to point to the new node. Old node updated to become one of two. -// The level argument specifies the number of steps up from the leaf -// level to insert; e.g. a data rectangle goes in at level = 0. -func insertRectRec(rect rectT, item Item, node *nodeT, newNode **nodeT, level int) bool { - var index int - var branch branchT - var otherNode *nodeT - // Still above level for insertion, go down tree recursively - if node == nil { - return false - } - if node.level > level { - index = pickBranch(rect, node) - if !insertRectRec(rect, item, node.branch[index].child, &otherNode, level) { - // Child was not split - node.branch[index].rect = combineRect(rect, node.branch[index].rect) - return false - } // Child was split - node.branch[index].rect = nodeCover(node.branch[index].child) - branch.child = otherNode - if branch.child == nil { - println(">> child assigned is nil") - } - branch.rect = nodeCover(otherNode) - return addBranch(&branch, node, newNode) - } else if node.level == level { // Have reached level for insertion. Add rect, split if necessary - branch.rect = rect - branch.item = item - // Child field of leaves contains id of data record - return addBranch(&branch, node, newNode) - } else { - // Should never occur - return false - } -} - -// Insert a data rectangle into an index structure. -// InsertRect provides for splitting the root; -// returns 1 if root was split, 0 if it was not. -// The level argument specifies the number of steps up from the leaf -// level to insert; e.g. a data rectangle goes in at level = 0. -// InsertRect2 does the recursion. -func insertRect(rect rectT, item Item, root **nodeT, level int) bool { - var newRoot *nodeT - var newNode *nodeT - var branch branchT - if *root == nil { - println(">> root is nil") - } - if insertRectRec(rect, item, *root, &newNode, level) { // Root split - newRoot = &nodeT{} // Grow tree taller and new root - newRoot.level = (*root).level + 1 - branch.rect = nodeCover(*root) - branch.child = *root - if branch.child == nil { - println(">> child assigned is nil 2") - } - addBranch(&branch, newRoot, nil) - branch.rect = nodeCover(newNode) - branch.child = newNode - if branch.child == nil { - println(">> child assigned is nil 3") - } - addBranch(&branch, newRoot, nil) - *root = newRoot - return true - } - return false -} - -// Find the smallest rectangle that includes all rectangles in branches of a node. -func nodeCover(node *nodeT) rectT { - var firstTime = true - var rect rectT - for index := 0; index < node.count; index++ { - if firstTime { - rect = node.branch[index].rect - firstTime = false - } else { - rect = combineRect(rect, node.branch[index].rect) - } - } - return rect -} - -// Add a branch to a node. Split the node if necessary. -// Returns 0 if node not split. Old node updated. -// Returns 1 if node split, sets *new_node to address of new node. -// Old node updated, becomes one of two. -func addBranch(branch *branchT, node *nodeT, newNode **nodeT) bool { - if node.count < maxNodes { // Split won't be necessary - node.branch[node.count] = *branch - node.count++ - return false - } - splitNode(node, branch, newNode) - return true -} - -// Disconnect a dependent node. -// Caller must return (or stop using iteration index) after this as count has changed -func disconnectBranch(node *nodeT, index int) { - // Remove element by swapping with the last element to prevent gaps in array - node.branch[index] = node.branch[node.count-1] - node.count-- -} - -// Pick a branch. Pick the one that will need the smallest increase -// in area to accommodate the new rectangle. This will result in the -// least total area for the covering rectangles in the current node. -// In case of a tie, pick the one which was smaller before, to get -// the best resolution when searching. -func pickBranch(rect rectT, node *nodeT) int { - var firstTime = true - var increase float64 - var bestIncr float64 = -1 - var area float64 - var bestArea float64 - var best int - var tempRect rectT - for index := 0; index < node.count; index++ { - curRect := node.branch[index].rect - area = calcRectVolume(curRect) - tempRect = combineRect(rect, curRect) - increase = calcRectVolume(tempRect) - area - if (increase < bestIncr) || firstTime { - best = index - bestArea = area - bestIncr = increase - firstTime = false - } else if (increase == bestIncr) && (area < bestArea) { - best = index - bestArea = area - bestIncr = increase - } - } - return best -} - -// Combine two rectangles into larger one containing both -func combineRect(rectA, rectB rectT) rectT { - var newRect rectT - for index := 0; index < 2; index++ { - newRect.min[index] = min(rectA.min[index], rectB.min[index]) - newRect.max[index] = max(rectA.max[index], rectB.max[index]) - } - return newRect -} - -// Split a node. -// Divides the nodes branches and the extra one between two nodes. -// Old node is one of the new ones, and one really new one is created. -// Tries more than one method for choosing a partition, uses best result. - -func splitNode(node *nodeT, branch *branchT, newNode **nodeT) { - // Could just use local here, but member or external is faster since it is reused - var localVars partitionVarsT - var parVars = &localVars - var level int - // Load all the branches into a buffer, initialize old node - level = node.level - getBranches(node, branch, parVars) - // Find partition - choosePartition(parVars, minNodes) - // Put branches from buffer into 2 nodes according to chosen partition - *newNode = &nodeT{} - node.level = level - (*newNode).level = node.level - loadNodes(node, *newNode, parVars) -} - -// Calculate the n-dimensional volume of a rectangle -func rectVolume(rect rectT) float64 { - var volume float64 = 1 - for index := 0; index < 2; index++ { - volume *= rect.max[index] - rect.min[index] - } - return volume -} - -// The exact volume of the bounding sphere for the given rectT -func rectSphericalVolume(rect rectT) float64 { - var sumOfSquares float64 - var radius float64 - for index := 0; index < 2; index++ { - var halfExtent = (rect.max[index] - rect.min[index]) * 0.5 - sumOfSquares += halfExtent * halfExtent - } - radius = math.Sqrt(sumOfSquares) - // Pow maybe slow, so test for common dims like 2,3 and just use x*x, x*x*x. - if 2 == 3 { - return radius * radius * radius * unitSphereVolume - } else if 2 == 2 { - return radius * radius * unitSphereVolume - } else { - return math.Pow(radius, 2) * unitSphereVolume - } -} - -// Use one of the methods to calculate retangle volume -func calcRectVolume(rect rectT) float64 { - if useSphericalVolume { - return rectSphericalVolume(rect) // Slower but helps certain merge cases - } - return rectVolume(rect) // Faster but can cause poor merges -} - -// Load branch buffer with branches from full node plus the extra branch. -func getBranches(node *nodeT, branch *branchT, parVars *partitionVarsT) { - // Load the branch buffer - for index := 0; index < maxNodes; index++ { - parVars.branchBuf[index] = node.branch[index] - } - parVars.branchBuf[maxNodes] = *branch - parVars.branchCount = maxNodes + 1 - // Calculate rect containing all in the set - parVars.coverSplit = parVars.branchBuf[0].rect - for index := 1; index < maxNodes+1; index++ { - parVars.coverSplit = combineRect(parVars.coverSplit, parVars.branchBuf[index].rect) - } - parVars.coverSplitArea = calcRectVolume(parVars.coverSplit) - node.count = 0 - node.level = -1 -} - -// Method #0 for choosing a partition: -// As the seeds for the two groups, pick the two rects that would waste the -// most area if covered by a single rectangle, i.e. evidently the worst pair -// to have in the same group. -// Of the remaining, one at a time is chosen to be put in one of the two groups. -// The one chosen is the one with the greatest difference in area expansion -// depending on which group - the rect most strongly attracted to one group -// and repelled from the other. -// If one group gets too full (more would force other group to violate min -// fill requirement) then other group gets the rest. -// These last are the ones that can go in either group most easily. -func choosePartition(parVars *partitionVarsT, minFill int) { - var biggestDiff float64 - var group, chosen, betterGroup int - initParVars(parVars, parVars.branchCount, minFill) - pickSeeds(parVars) - for ((parVars.count[0] + parVars.count[1]) < parVars.total) && - (parVars.count[0] < (parVars.total - parVars.minFill)) && - (parVars.count[1] < (parVars.total - parVars.minFill)) { - biggestDiff = -1 - for index := 0; index < parVars.total; index++ { - if !parVars.taken[index] { - var curRect = parVars.branchBuf[index].rect - rect0 := combineRect(curRect, parVars.cover[0]) - rect1 := combineRect(curRect, parVars.cover[1]) - growth0 := calcRectVolume(rect0) - parVars.area[0] - growth1 := calcRectVolume(rect1) - parVars.area[1] - diff := growth1 - growth0 - if diff >= 0 { - group = 0 - } else { - group = 1 - diff = -diff - } - if diff > biggestDiff { - biggestDiff = diff - chosen = index - betterGroup = group - } else if (diff == biggestDiff) && (parVars.count[group] < parVars.count[betterGroup]) { - chosen = index - betterGroup = group - } - } - } - classify(chosen, betterGroup, parVars) - } - // If one group too full, put remaining rects in the other - if (parVars.count[0] + parVars.count[1]) < parVars.total { - if parVars.count[0] >= parVars.total-parVars.minFill { - group = 1 - } else { - group = 0 - } - for index := 0; index < parVars.total; index++ { - if !parVars.taken[index] { - classify(index, group, parVars) - } - } - } -} - -// Copy branches from the buffer into two nodes according to the partition. -func loadNodes(nodeA *nodeT, nodeB *nodeT, parVars *partitionVarsT) { - for index := 0; index < parVars.total; index++ { - if parVars.partition[index] == 0 { - addBranch(&parVars.branchBuf[index], nodeA, nil) - } else if parVars.partition[index] == 1 { - addBranch(&parVars.branchBuf[index], nodeB, nil) - } - } -} - -// Initialize a partitionVarsT structure. -func initParVars(parVars *partitionVarsT, maxRects int, minFill int) { - parVars.count[1] = 0 - parVars.count[0] = parVars.count[1] - parVars.area[1] = 0 - parVars.area[0] = parVars.area[1] - parVars.total = maxRects - parVars.minFill = minFill - for index := 0; index < maxRects; index++ { - parVars.taken[index] = false - parVars.partition[index] = -1 - } -} - -func pickSeeds(parVars *partitionVarsT) { - var seed0, seed1 int - var worst, waste float64 - var area [maxNodes + 1]float64 - for index := 0; index < parVars.total; index++ { - area[index] = calcRectVolume(parVars.branchBuf[index].rect) - } - worst = -parVars.coverSplitArea - 1 - for indexA := 0; indexA < parVars.total-1; indexA++ { - for indexB := indexA + 1; indexB < parVars.total; indexB++ { - var oneRect = combineRect(parVars.branchBuf[indexA].rect, parVars.branchBuf[indexB].rect) - waste = calcRectVolume(oneRect) - area[indexA] - area[indexB] - if waste > worst { - worst = waste - seed0 = indexA - seed1 = indexB - } - } - } - classify(seed0, 0, parVars) - classify(seed1, 1, parVars) -} - -// Put a branch in one of the groups. -func classify(index int, group int, parVars *partitionVarsT) { - parVars.partition[index] = group - parVars.taken[index] = true - - if parVars.count[group] == 0 { - parVars.cover[group] = parVars.branchBuf[index].rect - } else { - parVars.cover[group] = combineRect(parVars.branchBuf[index].rect, parVars.cover[group]) - } - parVars.area[group] = calcRectVolume(parVars.cover[group]) - parVars.count[group]++ -} - -// Delete a data rectangle from an index structure. -// Pass in a pointer to a rectT, the tid of the record, ptr to ptr to root node. -// Returns 1 if record not found, 0 if success. -// RemoveRect provides for eliminating the root. -func removeRect(rect rectT, item Item, root **nodeT) bool { - var tempNode *nodeT - var reInsertList *listNodeT - if !removeRectRec(rect, item, *root, &reInsertList) { - // Found and deleted a data item - // Reinsert any branches from eliminated nodes - for reInsertList != nil { - tempNode = reInsertList.node - for index := 0; index < tempNode.count; index++ { - insertRect(tempNode.branch[index].rect, - tempNode.branch[index].item, - root, - tempNode.level) - } - reInsertList = reInsertList.next - } - // Check for redundant root (not leaf, 1 child) and eliminate - if (*root).count == 1 && (*root).isInternalNode() { - tempNode = (*root).branch[0].child - *root = tempNode - } - return false - } - return true -} - -// Delete a rectangle from non-root part of an index structure. -// Called by RemoveRect. Descends tree recursively, -// merges branches on the way back up. -// Returns 1 if record not found, 0 if success. -func removeRectRec(rect rectT, item Item, node *nodeT, listNode **listNodeT) bool { - if node == nil { - return true - } - if node.isInternalNode() { // not a leaf node - for index := 0; index < node.count; index++ { - if overlap(rect, node.branch[index].rect) { - if !removeRectRec(rect, item, node.branch[index].child, listNode) { - if node.branch[index].child.count >= minNodes { - // child removed, just resize parent rect - node.branch[index].rect = nodeCover(node.branch[index].child) - } else { - // child removed, not enough entries in node, eliminate node - reInsert(node.branch[index].child, listNode) - disconnectBranch(node, index) // Must return after this call as count has changed - } - return false - } - } - } - return true - } - // A leaf node - for index := 0; index < node.count; index++ { - if node.branch[index].item == item { - disconnectBranch(node, index) // Must return after this call as count has changed - return false - } - } - return true -} - -// Decide whether two rectangles overlap. -func overlap(rectA rectT, rectB rectT) bool { - for index := 0; index < 2; index++ { - if rectA.min[index] > rectB.max[index] || - rectB.min[index] > rectA.max[index] { - return false - } - } - return true -} - -// Add a node to the reinsertion list. All its branches will later -// be reinserted into the index structure. -func reInsert(node *nodeT, listNode **listNodeT) { - *listNode = &listNodeT{ - node: node, - next: *listNode, - } -} - -// Search in an index tree or subtree for all data retangles that overlap the argument rectangle. -func search(node *nodeT, rect rectT, iterator func(item Item) bool) bool { - if node == nil { - return true - } - if node.isInternalNode() { // This is an internal node in the tree - for index := 0; index < node.count; index++ { - nrect := node.branch[index].rect - if overlap(rect, nrect) { - if !search(node.branch[index].child, rect, iterator) { - return false // Don't continue searching - } - } - } - } else { // This is a leaf node - for index := 0; index < node.count; index++ { - if overlap(rect, node.branch[index].rect) { - // NOTE: There are different ways to return results. Here's where to modify - if !iterator(node.branch[index].item) { - return false // Don't continue searching - } - } - } - } - return true // Continue searching -} diff --git a/index/rtree/rtreed.go b/index/rtree/rtreed.go new file mode 100644 index 00000000..d7153ba4 --- /dev/null +++ b/index/rtree/rtreed.go @@ -0,0 +1,658 @@ +package rtree + +import "math" + +func d2fmin(a, b float64) float64 { + if a < b { + return a + } + return b +} +func d2fmax(a, b float64) float64 { + if a > b { + return a + } + return b +} + +const ( + d2numDims = 2 + d2maxNodes = 8 + d2minNodes = d2maxNodes / 2 + d2useSphericalVolume = true // Better split classification, may be slower on some systems +) + +var d2unitSphereVolume = []float64{ + 0.000000, 2.000000, 3.141593, // Dimension 0,1,2 + 4.188790, 4.934802, 5.263789, // Dimension 3,4,5 + 5.167713, 4.724766, 4.058712, // Dimension 6,7,8 + 3.298509, 2.550164, 1.884104, // Dimension 9,10,11 + 1.335263, 0.910629, 0.599265, // Dimension 12,13,14 + 0.381443, 0.235331, 0.140981, // Dimension 15,16,17 + 0.082146, 0.046622, 0.025807, // Dimension 18,19,20 +}[d2numDims] + +type d2RTree struct { + root *d2nodeT ///< Root of tree +} + +/// Minimal bounding rectangle (n-dimensional) +type d2rectT struct { + min [d2numDims]float64 ///< Min dimensions of bounding box + max [d2numDims]float64 ///< Max dimensions of bounding box +} + +/// May be data or may be another subtree +/// The parents level determines this. +/// If the parents level is 0, then this is data +type d2branchT struct { + rect d2rectT ///< Bounds + child *d2nodeT ///< Child node + data interface{} ///< Data Id or Ptr +} + +/// d2nodeT for each branch level +type d2nodeT struct { + count int ///< Count + level int ///< Leaf is zero, others positive + branch [d2maxNodes]d2branchT ///< Branch +} + +func (node *d2nodeT) isInternalNode() bool { + return (node.level > 0) // Not a leaf, but a internal node +} +func (node *d2nodeT) isLeaf() bool { + return (node.level == 0) // A leaf, contains data +} + +/// A link list of nodes for reinsertion after a delete operation +type d2listNodeT struct { + next *d2listNodeT ///< Next in list + node *d2nodeT ///< Node +} + +const d2notTaken = -1 // indicates that position + +/// Variables for finding a split partition +type d2partitionVarsT struct { + partition [d2maxNodes + 1]int + total int + minFill int + count [2]int + cover [2]d2rectT + area [2]float64 + + branchBuf [d2maxNodes + 1]d2branchT + branchCount int + coverSplit d2rectT + coverSplitArea float64 +} + +func d2New() *d2RTree { + // We only support machine word size simple data type eg. integer index or object pointer. + // Since we are storing as union with non data branch + return &d2RTree{ + root: &d2nodeT{}, + } +} + +/// Insert entry +/// \param a_min Min of bounding rect +/// \param a_max Max of bounding rect +/// \param a_dataId Positive Id of data. Maybe zero, but negative numbers not allowed. +func (tr *d2RTree) Insert(min, max [d2numDims]float64, dataId interface{}) { + var branch d2branchT + branch.data = dataId + for axis := 0; axis < d2numDims; axis++ { + branch.rect.min[axis] = min[axis] + branch.rect.max[axis] = max[axis] + } + d2insertRect(&branch, &tr.root, 0) +} + +/// Remove entry +/// \param a_min Min of bounding rect +/// \param a_max Max of bounding rect +/// \param a_dataId Positive Id of data. Maybe zero, but negative numbers not allowed. +func (tr *d2RTree) Remove(min, max [d2numDims]float64, dataId interface{}) { + var rect d2rectT + for axis := 0; axis < d2numDims; axis++ { + rect.min[axis] = min[axis] + rect.max[axis] = max[axis] + } + d2removeRect(&rect, dataId, &tr.root) +} + +/// Find all within d2search rectangle +/// \param a_min Min of d2search bounding rect +/// \param a_max Max of d2search bounding rect +/// \param a_searchResult d2search result array. Caller should set grow size. Function will reset, not append to array. +/// \param a_resultCallback Callback function to return result. Callback should return 'true' to continue searching +/// \param a_context User context to pass as parameter to a_resultCallback +/// \return Returns the number of entries found +func (tr *d2RTree) Search(min, max [d2numDims]float64, resultCallback func(data interface{}) bool) int { + var rect d2rectT + for axis := 0; axis < d2numDims; axis++ { + rect.min[axis] = min[axis] + rect.max[axis] = max[axis] + } + foundCount, _ := d2search(tr.root, rect, 0, resultCallback) + return foundCount +} + +/// Count the data elements in this container. This is slow as no internal counter is maintained. +func (tr *d2RTree) Count() int { + var count int + d2countRec(tr.root, &count) + return count +} + +/// Remove all entries from tree +func (tr *d2RTree) RemoveAll() { + // Delete all existing nodes + tr.root = &d2nodeT{} +} + +func d2countRec(node *d2nodeT, count *int) { + if node.isInternalNode() { // not a leaf node + for index := 0; index < node.count; index++ { + d2countRec(node.branch[index].child, count) + } + } else { // A leaf node + *count += node.count + } +} + +// Inserts a new data rectangle into the index structure. +// Recursively descends tree, propagates splits back up. +// Returns 0 if node was not split. Old node updated. +// If node was split, returns 1 and sets the pointer pointed to by +// new_node to point to the new node. Old node updated to become one of two. +// The level argument specifies the number of steps up from the leaf +// level to insert; e.g. a data rectangle goes in at level = 0. +func d2insertRectRec(branch *d2branchT, node *d2nodeT, newNode **d2nodeT, level int) bool { + // recurse until we reach the correct level for the new record. data records + // will always be called with a_level == 0 (leaf) + if node.level > level { + // Still above level for insertion, go down tree recursively + var otherNode *d2nodeT + //var newBranch d2branchT + + // find the optimal branch for this record + index := d2pickBranch(&branch.rect, node) + + // recursively insert this record into the picked branch + childWasSplit := d2insertRectRec(branch, node.branch[index].child, &otherNode, level) + + if !childWasSplit { + // Child was not split. Merge the bounding box of the new record with the + // existing bounding box + node.branch[index].rect = d2combineRect(&branch.rect, &(node.branch[index].rect)) + return false + } else { + // Child was split. The old branches are now re-partitioned to two nodes + // so we have to re-calculate the bounding boxes of each node + node.branch[index].rect = d2nodeCover(node.branch[index].child) + var newBranch d2branchT + newBranch.child = otherNode + newBranch.rect = d2nodeCover(otherNode) + + // The old node is already a child of a_node. Now add the newly-created + // node to a_node as well. a_node might be split because of that. + return d2addBranch(&newBranch, node, newNode) + } + } else if node.level == level { + // We have reached level for insertion. Add rect, split if necessary + return d2addBranch(branch, node, newNode) + } else { + // Should never occur + return false + } +} + +// Insert a data rectangle into an index structure. +// d2insertRect provides for splitting the root; +// returns 1 if root was split, 0 if it was not. +// The level argument specifies the number of steps up from the leaf +// level to insert; e.g. a data rectangle goes in at level = 0. +// InsertRect2 does the recursion. +// +func d2insertRect(branch *d2branchT, root **d2nodeT, level int) bool { + var newNode *d2nodeT + + if d2insertRectRec(branch, *root, &newNode, level) { // Root split + + // Grow tree taller and new root + newRoot := &d2nodeT{} + newRoot.level = (*root).level + 1 + + var newBranch d2branchT + + // add old root node as a child of the new root + newBranch.rect = d2nodeCover(*root) + newBranch.child = *root + d2addBranch(&newBranch, newRoot, nil) + + // add the split node as a child of the new root + newBranch.rect = d2nodeCover(newNode) + newBranch.child = newNode + d2addBranch(&newBranch, newRoot, nil) + + // set the new root as the root node + *root = newRoot + + return true + } + return false +} + +// Find the smallest rectangle that includes all rectangles in branches of a node. +func d2nodeCover(node *d2nodeT) d2rectT { + rect := node.branch[0].rect + for index := 1; index < node.count; index++ { + rect = d2combineRect(&rect, &(node.branch[index].rect)) + } + return rect +} + +// Add a branch to a node. Split the node if necessary. +// Returns 0 if node not split. Old node updated. +// Returns 1 if node split, sets *new_node to address of new node. +// Old node updated, becomes one of two. +func d2addBranch(branch *d2branchT, node *d2nodeT, newNode **d2nodeT) bool { + if node.count < d2maxNodes { // Split won't be necessary + node.branch[node.count] = *branch + node.count++ + return false + } else { + d2splitNode(node, branch, newNode) + return true + } +} + +// Disconnect a dependent node. +// Caller must return (or stop using iteration index) after this as count has changed +func d2disconnectBranch(node *d2nodeT, index int) { + // Remove element by swapping with the last element to prevent gaps in array + node.branch[index] = node.branch[node.count-1] + node.branch[node.count-1].data = nil + node.branch[node.count-1].child = nil + node.count-- +} + +// Pick a branch. Pick the one that will need the smallest increase +// in area to accomodate the new rectangle. This will result in the +// least total area for the covering rectangles in the current node. +// In case of a tie, pick the one which was smaller before, to get +// the best resolution when searching. +func d2pickBranch(rect *d2rectT, node *d2nodeT) int { + var firstTime bool = true + var increase float64 + var bestIncr float64 = -1 + var area float64 + var bestArea float64 + var best int + var tempRect d2rectT + + for index := 0; index < node.count; index++ { + curRect := &node.branch[index].rect + area = d2calcRectVolume(curRect) + tempRect = d2combineRect(rect, curRect) + increase = d2calcRectVolume(&tempRect) - area + if (increase < bestIncr) || firstTime { + best = index + bestArea = area + bestIncr = increase + firstTime = false + } else if (increase == bestIncr) && (area < bestArea) { + best = index + bestArea = area + bestIncr = increase + } + } + return best +} + +// Combine two rectangles into larger one containing both +func d2combineRect(rectA, rectB *d2rectT) d2rectT { + var newRect d2rectT + + for index := 0; index < d2numDims; index++ { + newRect.min[index] = d2fmin(rectA.min[index], rectB.min[index]) + newRect.max[index] = d2fmax(rectA.max[index], rectB.max[index]) + } + + return newRect +} + +// Split a node. +// Divides the nodes branches and the extra one between two nodes. +// Old node is one of the new ones, and one really new one is created. +// Tries more than one method for choosing a partition, uses best result. +func d2splitNode(node *d2nodeT, branch *d2branchT, newNode **d2nodeT) { + // Could just use local here, but member or external is faster since it is reused + var localVars d2partitionVarsT + parVars := &localVars + + // Load all the branches into a buffer, initialize old node + d2getBranches(node, branch, parVars) + + // Find partition + d2choosePartition(parVars, d2minNodes) + + // Create a new node to hold (about) half of the branches + *newNode = &d2nodeT{} + (*newNode).level = node.level + + // Put branches from buffer into 2 nodes according to the chosen partition + node.count = 0 + d2loadNodes(node, *newNode, parVars) +} + +// Calculate the n-dimensional volume of a rectangle +func d2rectVolume(rect *d2rectT) float64 { + var volume float64 = 1 + for index := 0; index < d2numDims; index++ { + volume *= rect.max[index] - rect.min[index] + } + return volume +} + +// The exact volume of the bounding sphere for the given d2rectT +func d2rectSphericalVolume(rect *d2rectT) float64 { + var sumOfSquares float64 = 0 + var radius float64 + + for index := 0; index < d2numDims; index++ { + halfExtent := (rect.max[index] - rect.min[index]) * 0.5 + sumOfSquares += halfExtent * halfExtent + } + + radius = math.Sqrt(sumOfSquares) + + // Pow maybe slow, so test for common dims just use x*x, x*x*x. + if d2numDims == 5 { + return (radius * radius * radius * radius * radius * d2unitSphereVolume) + } else if d2numDims == 4 { + return (radius * radius * radius * radius * d2unitSphereVolume) + } else if d2numDims == 3 { + return (radius * radius * radius * d2unitSphereVolume) + } else if d2numDims == 2 { + return (radius * radius * d2unitSphereVolume) + } else { + return (math.Pow(radius, d2numDims) * d2unitSphereVolume) + } +} + +// Use one of the methods to calculate retangle volume +func d2calcRectVolume(rect *d2rectT) float64 { + if d2useSphericalVolume { + return d2rectSphericalVolume(rect) // Slower but helps certain merge cases + } else { // RTREE_USE_SPHERICAL_VOLUME + return d2rectVolume(rect) // Faster but can cause poor merges + } // RTREE_USE_SPHERICAL_VOLUME +} + +// Load branch buffer with branches from full node plus the extra branch. +func d2getBranches(node *d2nodeT, branch *d2branchT, parVars *d2partitionVarsT) { + // Load the branch buffer + for index := 0; index < d2maxNodes; index++ { + parVars.branchBuf[index] = node.branch[index] + } + parVars.branchBuf[d2maxNodes] = *branch + parVars.branchCount = d2maxNodes + 1 + + // Calculate rect containing all in the set + parVars.coverSplit = parVars.branchBuf[0].rect + for index := 1; index < d2maxNodes+1; index++ { + parVars.coverSplit = d2combineRect(&parVars.coverSplit, &parVars.branchBuf[index].rect) + } + parVars.coverSplitArea = d2calcRectVolume(&parVars.coverSplit) +} + +// Method #0 for choosing a partition: +// As the seeds for the two groups, pick the two rects that would waste the +// most area if covered by a single rectangle, i.e. evidently the worst pair +// to have in the same group. +// Of the remaining, one at a time is chosen to be put in one of the two groups. +// The one chosen is the one with the greatest difference in area expansion +// depending on which group - the rect most strongly attracted to one group +// and repelled from the other. +// If one group gets too full (more would force other group to violate min +// fill requirement) then other group gets the rest. +// These last are the ones that can go in either group most easily. +func d2choosePartition(parVars *d2partitionVarsT, minFill int) { + var biggestDiff float64 + var group, chosen, betterGroup int + + d2initParVars(parVars, parVars.branchCount, minFill) + d2pickSeeds(parVars) + + for ((parVars.count[0] + parVars.count[1]) < parVars.total) && + (parVars.count[0] < (parVars.total - parVars.minFill)) && + (parVars.count[1] < (parVars.total - parVars.minFill)) { + biggestDiff = -1 + for index := 0; index < parVars.total; index++ { + if d2notTaken == parVars.partition[index] { + curRect := &parVars.branchBuf[index].rect + rect0 := d2combineRect(curRect, &parVars.cover[0]) + rect1 := d2combineRect(curRect, &parVars.cover[1]) + growth0 := d2calcRectVolume(&rect0) - parVars.area[0] + growth1 := d2calcRectVolume(&rect1) - parVars.area[1] + diff := growth1 - growth0 + if diff >= 0 { + group = 0 + } else { + group = 1 + diff = -diff + } + + if diff > biggestDiff { + biggestDiff = diff + chosen = index + betterGroup = group + } else if (diff == biggestDiff) && (parVars.count[group] < parVars.count[betterGroup]) { + chosen = index + betterGroup = group + } + } + } + d2classify(chosen, betterGroup, parVars) + } + + // If one group too full, put remaining rects in the other + if (parVars.count[0] + parVars.count[1]) < parVars.total { + if parVars.count[0] >= parVars.total-parVars.minFill { + group = 1 + } else { + group = 0 + } + for index := 0; index < parVars.total; index++ { + if d2notTaken == parVars.partition[index] { + d2classify(index, group, parVars) + } + } + } +} + +// Copy branches from the buffer into two nodes according to the partition. +func d2loadNodes(nodeA, nodeB *d2nodeT, parVars *d2partitionVarsT) { + for index := 0; index < parVars.total; index++ { + targetNodeIndex := parVars.partition[index] + targetNodes := []*d2nodeT{nodeA, nodeB} + + // It is assured that d2addBranch here will not cause a node split. + d2addBranch(&parVars.branchBuf[index], targetNodes[targetNodeIndex], nil) + } +} + +// Initialize a d2partitionVarsT structure. +func d2initParVars(parVars *d2partitionVarsT, maxRects, minFill int) { + parVars.count[0] = 0 + parVars.count[1] = 0 + parVars.area[0] = 0 + parVars.area[1] = 0 + parVars.total = maxRects + parVars.minFill = minFill + for index := 0; index < maxRects; index++ { + parVars.partition[index] = d2notTaken + } +} + +func d2pickSeeds(parVars *d2partitionVarsT) { + var seed0, seed1 int + var worst, waste float64 + var area [d2maxNodes + 1]float64 + + for index := 0; index < parVars.total; index++ { + area[index] = d2calcRectVolume(&parVars.branchBuf[index].rect) + } + + worst = -parVars.coverSplitArea - 1 + for indexA := 0; indexA < parVars.total-1; indexA++ { + for indexB := indexA + 1; indexB < parVars.total; indexB++ { + oneRect := d2combineRect(&parVars.branchBuf[indexA].rect, &parVars.branchBuf[indexB].rect) + waste = d2calcRectVolume(&oneRect) - area[indexA] - area[indexB] + if waste > worst { + worst = waste + seed0 = indexA + seed1 = indexB + } + } + } + + d2classify(seed0, 0, parVars) + d2classify(seed1, 1, parVars) +} + +// Put a branch in one of the groups. +func d2classify(index, group int, parVars *d2partitionVarsT) { + parVars.partition[index] = group + + // Calculate combined rect + if parVars.count[group] == 0 { + parVars.cover[group] = parVars.branchBuf[index].rect + } else { + parVars.cover[group] = d2combineRect(&parVars.branchBuf[index].rect, &parVars.cover[group]) + } + + // Calculate volume of combined rect + parVars.area[group] = d2calcRectVolume(&parVars.cover[group]) + + parVars.count[group]++ +} + +// Delete a data rectangle from an index structure. +// Pass in a pointer to a d2rectT, the tid of the record, ptr to ptr to root node. +// Returns 1 if record not found, 0 if success. +// d2removeRect provides for eliminating the root. +func d2removeRect(rect *d2rectT, id interface{}, root **d2nodeT) bool { + var reInsertList *d2listNodeT + + if !d2removeRectRec(rect, id, *root, &reInsertList) { + // Found and deleted a data item + // Reinsert any branches from eliminated nodes + for reInsertList != nil { + tempNode := reInsertList.node + + for index := 0; index < tempNode.count; index++ { + // TODO go over this code. should I use (tempNode->m_level - 1)? + d2insertRect(&tempNode.branch[index], root, tempNode.level) + } + reInsertList = reInsertList.next + } + + // Check for redundant root (not leaf, 1 child) and eliminate TODO replace + // if with while? In case there is a whole branch of redundant roots... + if (*root).count == 1 && (*root).isInternalNode() { + tempNode := (*root).branch[0].child + *root = tempNode + } + return false + } else { + return true + } +} + +// Delete a rectangle from non-root part of an index structure. +// Called by d2removeRect. Descends tree recursively, +// merges branches on the way back up. +// Returns 1 if record not found, 0 if success. +func d2removeRectRec(rect *d2rectT, id interface{}, node *d2nodeT, listNode **d2listNodeT) bool { + if node.isInternalNode() { // not a leaf node + for index := 0; index < node.count; index++ { + if d2overlap(*rect, node.branch[index].rect) { + if !d2removeRectRec(rect, id, node.branch[index].child, listNode) { + if node.branch[index].child.count >= d2minNodes { + // child removed, just resize parent rect + node.branch[index].rect = d2nodeCover(node.branch[index].child) + } else { + // child removed, not enough entries in node, eliminate node + d2reInsert(node.branch[index].child, listNode) + d2disconnectBranch(node, index) // Must return after this call as count has changed + } + return false + } + } + } + return true + } else { // A leaf node + for index := 0; index < node.count; index++ { + if node.branch[index].data == id { + d2disconnectBranch(node, index) // Must return after this call as count has changed + return false + } + } + return true + } +} + +// Decide whether two rectangles d2overlap. +func d2overlap(rectA, rectB d2rectT) bool { + for index := 0; index < d2numDims; index++ { + if rectA.min[index] > rectB.max[index] || + rectB.min[index] > rectA.max[index] { + return false + } + } + return true +} + +// Add a node to the reinsertion list. All its branches will later +// be reinserted into the index structure. +func d2reInsert(node *d2nodeT, listNode **d2listNodeT) { + newListNode := &d2listNodeT{} + newListNode.node = node + newListNode.next = *listNode + *listNode = newListNode +} + +// d2search in an index tree or subtree for all data retangles that d2overlap the argument rectangle. +func d2search(node *d2nodeT, rect d2rectT, foundCount int, resultCallback func(data interface{}) bool) (int, bool) { + if node.isInternalNode() { + // This is an internal node in the tree + for index := 0; index < node.count; index++ { + if d2overlap(rect, node.branch[index].rect) { + var ok bool + foundCount, ok = d2search(node.branch[index].child, rect, foundCount, resultCallback) + if !ok { + // The callback indicated to stop searching + return foundCount, false + } + } + } + } else { + // This is a leaf node + for index := 0; index < node.count; index++ { + if d2overlap(rect, node.branch[index].rect) { + id := node.branch[index].data + foundCount++ + if !resultCallback(id) { + return foundCount, false // Don't continue searching + } + + } + } + } + return foundCount, true // Continue searching +}