Index: trunk/psLib/src/math/psSort.h =================================================================== --- trunk/psLib/src/math/psSort.h (revision 18144) +++ trunk/psLib/src/math/psSort.h (revision 18145) @@ -30,5 +30,5 @@ -// Use of the PSSORT() macro requires additional macros defined: +// Use of the PSSORT() and PSSELECT() macros require additional macros defined: // COMPAREEXPR(A,B): expression to compare element A with element B; true if element A is smaller than B. // SWAPFUNC(TYPE,A,B): swap element A with element B; the type is provided so that a temporary variable can @@ -74,3 +74,64 @@ + +// The following algorithm for selection was provided by http://en.wikipedia.org/wiki/Selection_algorithm +// (version as of 21 May 2008, at 17:38), and implemented below as PSSELECT(): +// +// function partition(list, left, right, pivotIndex) +// pivotValue := list[pivotIndex] +// swap list[pivotIndex] and list[right] // Move pivot to end +// storeIndex := left +// for i from left to right-1 +// if list[i] < pivotValue +// swap list[storeIndex] and list[i] +// storeIndex := storeIndex + 1 +// swap list[right] and list[storeIndex] // Move pivot to its final place +// return storeIndex +// +// function select(list, k, left, right) +// loop +// select a pivot value list[pivotIndex] +// pivotNewIndex := partition(list, left, right, pivotIndex) +// if k = pivotNewIndex +// return list[k] +// else if k < pivotNewIndex +// right := pivotNewIndex-1 +// else +// left := pivotNewIndex+1 + + +// Select the RANK-th element +// This macro reorders the array so that the RANK-th element is in the correct position +// SIZE: Size of the array +// RANK: RANK to move into the correct position +// COMPAREEXPR: Macro with expression for comparison +// SWAPFUNC: Macro to swap elements +// SWAPTYPE: Type for swapping, passed to SWAPFUNC +#define PSSELECT(SIZE, RANK, COMPAREEXPR, SWAPFUNC, SWAPTYPE) { \ + bool selectContinue = true; /* Continue swapping? */ \ + long selectMax = SIZE - 1; /* Maximum index */ \ + long selectMin = 0; /* Minimum index */ \ + while (selectContinue) { \ + long selectPivot = (selectMin + selectMax) >> 1; /* Pivot index */ \ + SWAPFUNC(SWAPTYPE, selectPivot, selectMax); /* Move pivot to end */ \ + long selectStore = selectMin; /* Index of interest */ \ + for (long i = selectMin; i < selectMax; i++) { \ + if (COMPAREEXPR(i, selectMax)) { /* Note: comparing with the original pivot */ \ + SWAPFUNC(SWAPTYPE, selectStore, i); \ + selectStore++; \ + } \ + } \ + SWAPFUNC(SWAPTYPE, selectMax, selectStore); /* Move pivot to its final place */ \ + if (selectStore == RANK) { \ + selectContinue = false; /* Done */ \ + } else if (selectStore > RANK) { \ + selectMax = selectStore - 1; \ + } else { \ + selectMin = selectStore + 1; \ + } \ + } \ +} + + + #endif Index: trunk/psLib/src/pslib_strict.h =================================================================== --- trunk/psLib/src/pslib_strict.h (revision 18144) +++ trunk/psLib/src/pslib_strict.h (revision 18145) @@ -9,6 +9,6 @@ * @author Eric Van Alst, MHPCC * -* @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ -* @date $Date: 2008-03-04 22:02:09 $ +* @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ +* @date $Date: 2008-06-16 21:58:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -114,4 +114,5 @@ #include "psSparse.h" #include "psSlurp.h" +#include "psTree.h" #endif // #ifndef PS_LIB_STRICT_H Index: trunk/psLib/src/types/Makefile.am =================================================================== --- trunk/psLib/src/types/Makefile.am (revision 18144) +++ trunk/psLib/src/types/Makefile.am (revision 18145) @@ -15,5 +15,6 @@ psMetadataItemCompare.c \ psPixels.c \ - psArguments.c + psArguments.c \ + psTree.c EXTRA_DIST = types.i @@ -30,5 +31,6 @@ psMetadataItemCompare.h \ psPixels.h \ - psArguments.h + psArguments.h \ + psTree.h CLEANFILES = *~ *.bb *.bbg *.da Index: trunk/psLib/src/types/psTree.c =================================================================== --- trunk/psLib/src/types/psTree.c (revision 18145) +++ trunk/psLib/src/types/psTree.c (revision 18145) @@ -0,0 +1,522 @@ +#include +#include + +#include "psAssert.h" +#include "psAbort.h" +#include "psMemory.h" +#include "psType.h" +#include "psArray.h" +#include "psVector.h" +#include "psImage.h" +#include "psList.h" +#include "psSort.h" + +#include "psTree.h" + + +// XXX Upgrades: +// * Could allow different types of coordinates (saving memory) +// * Could make tree "growable" by including a pointer to a new tree within each node. +// * Rebalancing the tree ("pruning"?) is done by collecting all the points and re-planting a new tree. + +static void treeCoordArrayFree(psTreeCoordArray *array) +{ + psFree(array->F64); + psFree(array->raw); +} + +static void treeNodeFree(psTreeNode *node) +{ + // All contents are views only + return; +} + +static void treeFree(psTree *tree) +{ + for (long i = 0; i < tree->numNodes; i++) { + psFree(tree->nodes[i]); + } + psFree(tree->nodes); + psFree(tree->min); + psFree(tree->max); + psFree(tree->division); + // 'root' is a view only + psFree(tree->data); + psFree(tree->contents); +} + +psTreeCoordArray *psTreeCoordArrayAlloc(long size, int dim) +{ + psAssert(size > 0, "Size (%ld) must be positive", size); + psAssert(dim > 0, "Dimension (%d) must be positive", dim); + + psTreeCoordArray *array = psAlloc(sizeof(psTreeCoordArray)); // Array to return + psMemSetDeallocator(array, (psFreeFunc)treeCoordArrayFree); + + array->raw = psAlloc(dim * size * sizeof(psF64)); + array->F64 = psAlloc(size * sizeof(psF64*)); + + psF64 *ptr = (psF64*)array->raw; // Pointer into raw vector + for (long i = 0; i < size; i++, ptr += dim) { + array->F64[i] = ptr; + } + + return array; +} + +psTreeNode *psTreeNodeAlloc(long index) +{ + psAssert(index >= 0, "Index (%ld) must be non-negative", index); + + psTreeNode *node = psAlloc(sizeof(psTreeNode)); // Node to return + psMemSetDeallocator(node, (psFreeFunc)treeNodeFree); + + node->index = index; + node->divideDim = -1; + node->num = 0; + node->contents = NULL; + node->tree = NULL; + node->up = NULL; + node->left = NULL; + node->right = NULL; + + return node; +} + +psTree *psTreeAlloc(int dim, int maxLeafContents, long numNodes) +{ + psAssert(dim > 0, "Dimensionality (%d) must be positive", dim); + psAssert(maxLeafContents > 0, "Maximum number per leaf (%d) must be positive", maxLeafContents); + psAssert(numNodes > 0, "Number of nodes (%ld) must be positive", numNodes); + + psTree *tree = psAlloc(sizeof(psTree)); // Tree to return + psMemSetDeallocator(tree, (psFreeFunc)treeFree); + + tree->dim = dim; + tree->maxLeafContents = maxLeafContents; + + tree->numNodes = numNodes; + tree->nodes = psAlloc(numNodes * sizeof(psTreeNode*)); + for (long i = 0; i < numNodes; i++) { + tree->nodes[i] = psTreeNodeAlloc(i); + } + tree->root = tree->nodes[0]; + + tree->min = psTreeCoordArrayAlloc(numNodes, dim); + tree->max = psTreeCoordArrayAlloc(numNodes, dim); + + tree->division = psAlloc(numNodes * sizeof(double)); + + tree->numData = 0; + tree->data = NULL; + + return tree; +} + + +psTree *psTreePlant(int dim, int maxLeafContents, ...) +{ + PS_ASSERT_INT_POSITIVE(dim, NULL); + PS_ASSERT_INT_POSITIVE(maxLeafContents, NULL); + + // Parse coordinate list + va_list args; // Variable argument list + va_start(args, maxLeafContents); + psArray *coords = psArrayAlloc(dim); // Array of coordinates + long numData; // Number of data points + for (int i = 0; i < dim; i++) { + psVector *data = va_arg(args, psVector*); + if (i == 0) { + numData = data->n; + } else if (data->n != numData) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input vectors have differing sizes: %ld vs %ld", + data->n, numData); + psFree(coords); + va_end(args); + return NULL; + } + coords->data[i] = psMemIncrRefCounter(data); + } + va_end(args); + + long numNodes; + { + // From equation 21.2.2 in Numerical Recipes vol 3, p1102. + long pow2; // Smallest power of two >= numData + for (pow2 = 1; pow2 < numData; pow2 <<= 1); + numNodes = PS_MIN(pow2 - 1, 2 * numData - pow2 / 2 - 1); + } + + psTree *tree = psTreeAlloc(dim, maxLeafContents, numNodes); + tree->data = psTreeCoordArrayAlloc(numData, dim); + tree->numData = numData; + psF64 **data = tree->data->F64; // Dereference data + + // Copy data into the tree +#define PSTREE_COPY_DATA_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + for (long j = 0; j < numData; j++) { \ + data[j][i] = values->data.TYPE[j]; \ + } \ + break; \ + } + + for (int i = 0; i < dim; i++) { + psVector *values = coords->data[i]; // Vector of values for this dimension + switch (values->type.type) { + PSTREE_COPY_DATA_CASE(U8); + PSTREE_COPY_DATA_CASE(U16); + PSTREE_COPY_DATA_CASE(U32); + PSTREE_COPY_DATA_CASE(U64); + PSTREE_COPY_DATA_CASE(S8); + PSTREE_COPY_DATA_CASE(S16); + PSTREE_COPY_DATA_CASE(S32); + PSTREE_COPY_DATA_CASE(S64); + PSTREE_COPY_DATA_CASE(F32); + PSTREE_COPY_DATA_CASE(F64); + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unsupported vector type for dimension %d: %x", + i, values->type.type); + psFree(tree); + psFree(coords); + } + } + psFree(coords); + + long *contents = tree->contents = psAlloc(numData * sizeof(long)); // Contents of tree + for (long i = 0; i < numData; i++) { + contents[i] = i; + } + + psTreeNode *root = tree->root; // Root of the tree + root->tree = tree; + root->contents = contents; + root->num = numData; + root->divideDim = 0; + for (int i = 0; i < dim; i++) { + tree->min->F64[0][i] = -INFINITY; + tree->max->F64[0][i] = INFINITY; + } + + psArray *work = psArrayAlloc(numNodes); // Work queue + work->data[0] = root; + long workNum = 1; + +// Compare points (in a single dimension) by their index +#define PSTREE_SELECT_COMPARE(A,B) (data[contents[A]][divideDim] < data[contents[B]][divideDim]) +// Swap points by their index +#define PSTREE_SELECT_SWAP(TYPE,A,B) { \ + if (A != B) { \ + TYPE temp = contents[A]; \ + contents[A] = contents[B]; \ + contents[B] = temp; \ + } \ +} + + for (long workIndex = 0, nextIndex = 0; workIndex < workNum; workIndex++) { + psTreeNode *node = work->data[workIndex]; // Node to be worked on + long num = node->num; // Number of points in node + long *contents = node->contents;// Indices of points in node + psF64 **data = tree->data->F64; // Dereference the values + long middle = num / 2; // Index of middle point + int divideDim = node->divideDim;// Dimension in which to divide + long nodeIndex = node->index; // Index of node + + PSSELECT(num, middle, PSTREE_SELECT_COMPARE, PSTREE_SELECT_SWAP, long); + + psF64 *nodeMin = tree->min->F64[nodeIndex]; // Minimum bounds for current node + psF64 *nodeMax = tree->max->F64[nodeIndex]; // Maximum bounds for current node + + // Divide to the left and right + double division = tree->division[nodeIndex] = data[contents[middle]][divideDim]; + + long leftIndex = ++nextIndex, rightIndex = ++nextIndex; // Indices for left and right nodes + psTreeNode *left = node->left = tree->nodes[leftIndex]; // Node to the left + psTreeNode *right = node->right = tree->nodes[rightIndex]; // Node to the right + long numLeft = middle, numRight = num - middle; // Number in each node + + for (int i = 0; i < dim; i++) { + tree->min->F64[leftIndex][i] = nodeMin[i]; + tree->min->F64[rightIndex][i] = (i == divideDim ? division : nodeMin[i]); + tree->max->F64[leftIndex][i] = (i == divideDim ? division : nodeMax[i]); + tree->max->F64[rightIndex][i] = nodeMax[i]; + } + + left->num = numLeft; + right->num = numRight; + + left->contents = contents; + right->contents = contents + middle; + +#if 0 + // Check contents + for (int i = 0; i < numLeft; i++) { + long index = left->contents[i]; // Index of point + for (int j = 0; j < dim; j++) { + if (data[index][j] < tree->min->F64[leftIndex][j] || + data[index][j] > tree->max->F64[leftIndex][j]) { + psWarning("Bad division"); + } + } + } + for (int i = 0; i < numRight; i++) { + long index = right->contents[i]; // Index of point + for (int j = 0; j < dim; j++) { + if (data[index][j] < tree->min->F64[rightIndex][j] || + data[index][j] > tree->max->F64[rightIndex][j]) { + psWarning("Bad division"); + } + } + } +#endif + + divideDim++; + if (divideDim == dim) { + divideDim = 0; + } + left->divideDim = right->divideDim = divideDim; + + left->up = right->up = node; + left->tree = right->tree = tree; + + if (numLeft > maxLeafContents) { + work->data[workNum++] = left; + } + if (numRight > maxLeafContents) { + work->data[workNum++] = right; + } + + work->data[workIndex] = NULL; + } + psFree(work); + + return tree; +} + +static void treePrint(FILE *fptr, const psTreeNode *node, int level) +{ + long index = node->index; // Index of node + psTree *tree = node->tree; // Parent tree + int dim = tree->dim; // Dimensions + + for (int i = 0; i < level; i++) { + fprintf(fptr, " "); + } + fprintf(fptr, "%ld: (", index); + for (int i = 0; i < dim; i++) { + fprintf(fptr, "%lf", tree->min->F64[index][i]); + if (i < dim - 1) { + fprintf(fptr, ","); + } + } + fprintf(fptr, ") --> ("); + for (int i = 0; i < dim; i++) { + fprintf(fptr, "%lf", tree->max->F64[index][i]); + if (i < dim - 1) { + fprintf(fptr, ","); + } + } + fprintf(fptr, ")"); + level++; + + if (node->left && node->right) { + fprintf(fptr, " ==> %lf\n", tree->division[index]); + treePrint(fptr, node->left, level); + treePrint(fptr, node->right, level); + return; + } + fprintf(fptr, "\n"); + for (int i = 0; i < level; i++) { + fprintf(fptr, " "); + } + fprintf(fptr, "("); + for (int i = 0; i < node->num; i++) { + psF64 *coords = tree->data->F64[node->contents[i]]; // Coordinates + for (int j = 0; j < dim; j++) { + fprintf(fptr, "%lf", coords[j]); + if (j < dim - 1) { + fprintf(fptr, ","); + } + } + fprintf(fptr, ")"); + if (i < node->num - 1) { + fprintf(fptr, " ("); + } + } + fprintf(fptr, "\n"); + + return; +} + + +void psTreePrint(FILE *fptr, const psTree *tree) +{ + PS_ASSERT_TREE_NON_NULL(tree, ); + + treePrint(fptr, tree->root, 0); + return; +} + +// Given an arbitrary point, find the appropriate leaf +psTreeNode *psTreeLeaf(const psTree *tree, const psVector *coords) +{ +#if 1 // Might be in a tight loop + PS_ASSERT_TREE_NON_NULL(tree, NULL); + PS_ASSERT_VECTOR_NON_NULL(coords, NULL); + PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_SIZE(coords, (long)tree->dim, NULL); +#endif + + psTreeNode *node = tree->root; // Node of interest + for (; node->left; node = coords->data.F64[node->divideDim] < tree->division[node->index] ? + node->left : node->right); + return node; +} + +// Returns the square of the distance between a point on the tree and an aribtary point +static inline double treeContentDistance(const psTree *tree, long index, const psVector *coords) +{ + int dim = tree->dim; // Dimensionality + switch (dim) { + case 2: + return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) + + PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]); + case 3: + return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) + + PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]) + + PS_SQR(coords->data.F64[2] - tree->data->F64[index][2]); + default: { + double distance2 = 0.0; // Distance of interest + for (int i = 0; i < dim; i++) { + distance2 += PS_SQR(coords->data.F64[i] - tree->data->F64[index][i]); + } + return distance2; + } + } + return NAN; +} + +// Returns the square of the distance between a leaf and an arbitrary point +static inline double treeBranchDistance(const psTree *tree, long index, const psVector *coords) +{ + int dim = tree->dim; // Dimensionality + double distance = 0.0; // Distance to box + for (int i = 0; i < dim; i++) { + double minDiff = tree->min->F64[index][i] - coords->data.F64[i]; + if (minDiff > 0) { + distance += PS_SQR(minDiff); + continue; + } + double maxDiff = coords->data.F64[i] - tree->max->F64[index][i]; + if (maxDiff > 0) { + distance += PS_SQR(maxDiff); + continue; + } + } + return distance; +} + +// Returns whether the designated leaf contains an arbitrary point +static inline bool treeBranchInside(const psTree *tree, long index, const psVector *coords) +{ + int dim = tree->dim; // Dimensionality + for (int i = 0; i < dim; i++) { + if (tree->min->F64[index][i] > coords->data.F64[i]) { + return false; + } + if (tree->max->F64[index][i] < coords->data.F64[i]) { + return false; + } + } + return true; +} + +// Search a leaf node for a closer point +static inline void treeLeafSearch(long *bestIndex, + double *bestDistance, + const psTree *tree, + const psTreeNode *leaf, + const psVector *coords + ) +{ + for (int i = 0; i < leaf->num; i++) { + long index = leaf->contents[i]; // Index of point + double distance = treeContentDistance(tree, index, coords); // Distance to point + if (distance < *bestDistance) { + *bestIndex = index; + *bestDistance = distance; + } + } + return; +} + + + +// Given an arbitrary point, return the index of the nearest neighbour +long psTreeNearest(const psTree *tree, const psVector *coords) +{ +#if 1 // Might be in a tight loop + PS_ASSERT_TREE_NON_NULL(tree, -1); + PS_ASSERT_VECTOR_NON_NULL(coords, -1); + PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F64, -1); + PS_ASSERT_VECTOR_SIZE(coords, (long)tree->dim, -1); +#endif + + // Find the closest point in the leaf that contains the point of interest + psTreeNode *leaf = psTreeLeaf(tree, coords); // Leaf containing the point of interest + long bestIndex = -1; // Index of best point + double bestDistance = INFINITY; // Distance to best point + treeLeafSearch(&bestIndex, &bestDistance, tree, leaf, coords); + + // Work up from the leaf containing the point of interest, and for each in turn, search down the other + // branch. Every time we search down a branch, we first check to see if the current best distance rules + // out that branch. By having the check within the work queue loop, we allow the best distance to evolve + // (smaller), which means our elimination test goes up against the smallest possible distance. + psArray *work = psArrayAlloc(tree->numNodes); // Work queue + while (leaf->up) { + psTreeNode *up = leaf->up; // Parent node + + long workIndex = 0; + work->data[workIndex] = (leaf == up->left ? up->right : up->left); // The other node + while (workIndex >= 0) { + psTreeNode *node = work->data[workIndex]; + if (treeBranchDistance(tree, node->index, coords) > bestDistance) { + // No need to investigate + work->data[workIndex] = NULL; + workIndex--; + continue; + } + if (node->left) { + // Branch node + work->data[workIndex] = node->right; + work->data[++workIndex] = node->left; + continue; + } + // Leaf node + treeLeafSearch(&bestIndex, &bestDistance, tree, node, coords); + work->data[workIndex] = NULL; + workIndex--; + } + + leaf = up; + } + psFree(work); + + return bestIndex; +} + + +psVector *psTreeCoords(psVector *out, const psTree *tree, long index) +{ +#if 1 // Might be in a tight loop + PS_ASSERT_TREE_NON_NULL(tree, NULL); + PS_ASSERT_INT_NONNEGATIVE(index, NULL); + PS_ASSERT_INT_LESS_THAN_OR_EQUAL(index, tree->numData, NULL); +#endif + + int dim = tree->dim; // Dimensions + + out = psVectorRecycle(out, dim, PS_TYPE_F64); + memcpy(out->data.F64, tree->data->F64[index], dim * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + return out; +} Index: trunk/psLib/src/types/psTree.h =================================================================== --- trunk/psLib/src/types/psTree.h (revision 18145) +++ trunk/psLib/src/types/psTree.h (revision 18145) @@ -0,0 +1,102 @@ +#ifndef PS_TREE_H +#define PS_TREE_H + +#include + +/// An array of coordinates for the tree +/// +/// The coordinate for point i, dimension j is psTreeCoordArray->F64[i][j]. We have a "raw" vector which +/// points into the operational vector, as is done for psImage. This keeps all the memory together, reducing +/// fragmentation. Currently supports only F64 type, but could be expanded to support multiple types (to save +/// memory) by changing the F64 member into a union. +typedef struct { + psF64 **F64; ///< Coordinates + psU8 *raw; // Raw vector +} psTreeCoordArray; + +/// A simple kd-tree implementation +/// +/// This parent tree structure contains all the main information for the tree. We choose this instead of +/// including everything within the nodes to avoid memory fragmentation. +typedef struct { + int dim; ///< Dimensionality + int maxLeafContents; ///< Maximum number of points on a leaf + long numNodes; ///< Number of nodes + long numData; ///< Number of data points + struct psTreeNode **nodes; ///< Array of nodes + psTreeCoordArray *min, *max; ///< Bounding box: minimum and maximum coordinates + double *division; ///< Division for each node + struct psTreeNode *root; ///< Convenience pointer to root node of tree + psTreeCoordArray *data; ///< Data: coordinates + long *contents; ///< Indices into the data +} psTree; + +/// A node of the kd-tree +/// +/// Most of the information is stored in the parent tree (to reduce memory fragmentation). All the pointers +/// are views only; the pointers to the nodes are contained within the parent tree's "nodes" element. +typedef struct psTreeNode { + long index; ///< Index of node + int divideDim; ///< Dimension in which to divide + int num; ///< Number of points in node + long *contents; ///< Contents of node; view only + psTree *tree; ///< Parent tree; view only + struct psTreeNode *up; ///< Node higher up; view only + struct psTreeNode *left, *right; ///< Nodes below, to the left and right; views only +} psTreeNode; + +/// Assertion for psTree +#define PS_ASSERT_TREE_NON_NULL(TREE, RETVAL) \ +if (!(TREE) || !(TREE)->nodes || !(TREE)->min || !(TREE)->max || !(TREE)->division || !(TREE)->root || \ + !(TREE)->data || !(TREE)->contents) { \ + psError(PS_ERR_UNEXPECTED_NULL, true, "Tree %s or one of its components is NULL.", #TREE); \ + return RETVAL; \ +} + +/// Allocator for psTreeCoordArray +psTreeCoordArray *psTreeCoordArrayAlloc(long size, ///< Size of array + int dim ///< Dimensionality + ); + +/// Allocator for psTreeNode +psTreeNode *psTreeNodeAlloc(long index ///< Index of node + ); + +/// Allocator for psTree +psTree *psTreeAlloc(int dim, ///< Dimensionality + int maxLeafContents,///< Maximum number of points on a leaf + long numNodes ///< Number of nodes in tree + ); + +/// Plant (create) a tree +/// +/// This is the main startup function. +psTree *psTreePlant(int dim, ///< Dimensionality + int maxLeafContents,///< Maximum number of points on a leaf + ... ///< psVector for each coordinate + ); + +/// Return the leaf node containing given coordinates +/// +/// Returns a view to the pointer. +psTreeNode *psTreeLeaf(const psTree *tree, ///< Tree + const psVector *coords ///< Coordinates of interest + ); + +/// Return the index of the nearest neighbour to given coordinates +long psTreeNearest(const psTree *tree, ///< Tree + const psVector *coords ///< Coordinates of interest + ); + +/// Return the coordinates of a point in the tree, specified by its index +psVector *psTreeCoords(psVector *out, ///< Output vector, or NULL + const psTree *tree, ///< Tree + long index ///< Index of point in tree + ); + +/// Print a representation of the tree +void psTreePrint(FILE *fptr, ///< File to which to print + const psTree *tree ///< Tree + ); + +#endif Index: trunk/psLib/test/types/.cvsignore =================================================================== --- trunk/psLib/test/types/.cvsignore (revision 18144) +++ trunk/psLib/test/types/.cvsignore (revision 18145) @@ -77,2 +77,3 @@ tap_psMetadataUpdate tap_psMetadataOverlay +tap_psTree Index: trunk/psLib/test/types/Makefile.am =================================================================== --- trunk/psLib/test/types/Makefile.am (revision 18144) +++ trunk/psLib/test/types/Makefile.am (revision 18145) @@ -37,5 +37,6 @@ tap_psMetadataConfigFormat \ tap_psMetadataConfig_input \ - tap_psHash_845 + tap_psHash_845 \ + tap_psTree if BUILD_TESTS Index: trunk/psLib/test/types/tap_psTree.c =================================================================== --- trunk/psLib/test/types/tap_psTree.c (revision 18145) +++ trunk/psLib/test/types/tap_psTree.c (revision 18145) @@ -0,0 +1,72 @@ +#include +#include "tap.h" +#include "pstap.h" + +#define NUM 10000 // Number of points + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + plan_tests(6); + + { + psMemId id = psMemGetId(); + + psVector *x = psVectorAlloc(NUM, PS_TYPE_F64); + psVector *y = psVectorAlloc(NUM, PS_TYPE_F64); + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); + for (int i = 0; i < NUM; i++) { + x->data.F64[i] = 2.0 * psRandomUniform(rng) - 1.0; + y->data.F64[i] = 2.0 * psRandomUniform(rng) - 1.0; + } + psFree(rng); + + psTree *tree = psTreePlant(2, 2, x, y); + + ok(tree, "Tree planted"); + skip_start(!tree, 4, "tree died"); + { + // psTreePrint(stderr, tree); + + psVector *coords = psVectorAlloc(2, PS_TYPE_F64); + psVectorInit(coords, 0); + + long closeIndex = psTreeNearest(tree, coords); + psFree(coords); + ok(closeIndex >= 0 && closeIndex < tree->numNodes, "found point: %ld", closeIndex); + + long bestIndex = -1; + double bestDist = INFINITY; + for (int i = 0; i < NUM; i++) { + double dist = PS_SQR(x->data.F64[i]) + PS_SQR(y->data.F64[i]); + if (dist < bestDist) { + bestIndex = i; + bestDist = dist; + } + } + ok(bestIndex == closeIndex, "correct point: %ld vs %ld", closeIndex, bestIndex); + + psVector *closest = psTreeCoords(NULL, tree, closeIndex); + ok(closest, "got coords: %lf,%lf", closest->data.F64[0], closest->data.F64[1]); + ok(closest->data.F64[0] == x->data.F64[bestIndex] && + closest->data.F64[1] == y->data.F64[bestIndex], + "correst coords: %lf,%lf(%lf) vs %lf,%lf(%lf)", + closest->data.F64[0], closest->data.F64[1], + sqrt(PS_SQR(closest->data.F64[0]) + PS_SQR(closest->data.F64[1])), + x->data.F64[bestIndex], y->data.F64[bestIndex], + sqrt(PS_SQR(x->data.F64[bestIndex]) + PS_SQR(y->data.F64[bestIndex]))); + psFree(closest); + } + skip_end(); + + psFree(tree); + psFree(x); + psFree(y); + + ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); + } + + psLibFinalize(); +} +