[Git][debian-gis-team/pykdtree][master] 5 commits: New upstream version 1.3.0
Antonio Valentino
gitlab at salsa.debian.org
Mon Jun 4 20:42:09 BST 2018
Antonio Valentino pushed to branch master at Debian GIS Project / pykdtree
Commits:
9cc94bdf by Antonio Valentino at 2018-06-04T21:18:06+02:00
New upstream version 1.3.0
- - - - -
80e07ece by Antonio Valentino at 2018-06-04T21:18:07+02:00
Update upstream source from tag 'upstream/1.3.0'
Update to upstream version '1.3.0'
with Debian dir 536a068de51f5a19e05b108069271f26433ac174
- - - - -
a9259da8 by Antonio Valentino at 2018-06-04T21:19:06+02:00
New upstream release
- - - - -
c9ac706c by Antonio Valentino at 2018-06-04T19:36:15+00:00
Drop 0002-safe-clean.patch
- - - - -
2bad0fea by Antonio Valentino at 2018-06-04T19:37:19+00:00
Set distribution to unstable
- - - - -
13 changed files:
- PKG-INFO
- − README
- + README
- debian/changelog
- − debian/patches/0002-safe-clean.patch
- debian/patches/series
- pykdtree.egg-info/PKG-INFO
- pykdtree.egg-info/requires.txt
- pykdtree/_kdtree_core.c
- pykdtree/kdtree.c
- pykdtree/test_tree.py
- setup.cfg
- setup.py
Changes:
=====================================
PKG-INFO
=====================================
--- a/PKG-INFO
+++ b/PKG-INFO
@@ -1,6 +1,6 @@
Metadata-Version: 1.1
Name: pykdtree
-Version: 1.2.2
+Version: 1.3.0
Summary: Fast kd-tree implementation with OpenMP-enabled queries
Home-page: UNKNOWN
Author: Esben S. Nielsen
=====================================
README deleted
=====================================
--- a/README
+++ /dev/null
@@ -1 +0,0 @@
-README.rst
\ No newline at end of file
=====================================
README
=====================================
--- /dev/null
+++ b/README
@@ -0,0 +1,146 @@
+.. image:: https://travis-ci.org/storpipfugl/pykdtree.svg?branch=master
+ :target: https://travis-ci.org/storpipfugl/pykdtree
+.. image:: https://ci.appveyor.com/api/projects/status/ubo92368ktt2d25g/branch/master
+ :target: https://ci.appveyor.com/project/storpipfugl/pykdtree
+
+========
+pykdtree
+========
+
+Objective
+---------
+pykdtree is a kd-tree implementation for fast nearest neighbour search in Python.
+The aim is to be the fastest implementation around for common use cases (low dimensions and low number of neighbours) for both tree construction and queries.
+
+The implementation is based on scipy.spatial.cKDTree and libANN by combining the best features from both and focus on implementation efficiency.
+
+The interface is similar to that of scipy.spatial.cKDTree except only Euclidean distance measure is supported.
+
+Queries are optionally multithreaded using OpenMP.
+
+Installation
+------------
+Default build of pykdtree with OpenMP enabled queries using libgomp
+
+.. code-block:: bash
+
+ $ cd <pykdtree_dir>
+ $ python setup.py install
+
+If it fails with undefined compiler flags or you want to use another OpenMP implementation please modify setup.py at the indicated point to match your system.
+
+Building without OpenMP support is controlled by the USE_OMP environment variable
+
+.. code-block:: bash
+
+ $ cd <pykdtree_dir>
+ $ export USE_OMP=0
+ $ python setup.py install
+
+Note evironment variables are by default not exported when using sudo so in this case do
+
+.. code-block:: bash
+
+ $ USE_OMP=0 sudo -E python setup.py install
+
+Usage
+-----
+The usage of pykdtree is similar to scipy.spatial.cKDTree so for now refer to its documentation
+
+ >>> from pykdtree.kdtree import KDTree
+ >>> kd_tree = KDTree(data_pts)
+ >>> dist, idx = kd_tree.query(query_pts, k=8)
+
+The number of threads to be used in OpenMP enabled queries can be controlled with the standard OpenMP environment variable OMP_NUM_THREADS.
+
+The **leafsize** argument (number of data points per leaf) for the tree creation can be used to control the memory overhead of the kd-tree. pykdtree uses a default **leafsize=16**.
+Increasing **leafsize** will reduce the memory overhead and construction time but increase query time.
+
+pykdtree accepts data in double precision (numpy.float64) or single precision (numpy.float32) floating point. If data of another type is used an internal copy in double precision is made resulting in a memory overhead. If the kd-tree is constructed on single precision data the query points must be single precision as well.
+
+Benchmarks
+----------
+Comparison with scipy.spatial.cKDTree and libANN. This benchmark is on geospatial 3D data with 10053632 data points and 4276224 query points. The results are indexed relative to the construction time of scipy.spatial.cKDTree. A leafsize of 10 (scipy.spatial.cKDTree default) is used.
+
+Note: libANN is *not* thread safe. In this benchmark libANN is compiled with "-O3 -funroll-loops -ffast-math -fprefetch-loop-arrays" in order to achieve optimum performance.
+
+================== ===================== ====== ======== ==================
+Operation scipy.spatial.cKDTree libANN pykdtree pykdtree 4 threads
+------------------ --------------------- ------ -------- ------------------
+
+Construction 100 304 96 96
+
+query 1 neighbour 1267 294 223 70
+
+Total 1 neighbour 1367 598 319 166
+
+query 8 neighbours 2193 625 449 143
+
+Total 8 neighbours 2293 929 545 293
+================== ===================== ====== ======== ==================
+
+Looking at the combined construction and query this gives the following performance improvement relative to scipy.spatial.cKDTree
+
+========== ====== ======== ==================
+Neighbours libANN pykdtree pykdtree 4 threads
+---------- ------ -------- ------------------
+1 129% 329% 723%
+
+8 147% 320% 682%
+========== ====== ======== ==================
+
+Note: mileage will vary with the dataset at hand and computer architecture.
+
+Test
+----
+Run the unit tests using nosetest
+
+.. code-block:: bash
+
+ $ cd <pykdtree_dir>
+ $ python setup.py nosetests
+
+Installing on AppVeyor
+----------------------
+
+Pykdtree requires the "stdint.h" header file which is not available on certain
+versions of Windows or certain Windows compilers including those on the
+continuous integration platform AppVeyor. To get around this the header file(s)
+can be downloaded and placed in the correct "include" directory. This can
+be done by adding the `anaconda/missing-headers.ps1` script to your repository
+and running it the install step of `appveyor.yml`:
+
+ # install missing headers that aren't included with MSVC 2008
+ # https://github.com/omnia-md/conda-recipes/pull/524
+ - "powershell ./appveyor/missing-headers.ps1"
+
+In addition to this, AppVeyor does not support OpenMP so this feature must be
+turned off by adding the following to `appveyor.yml` in the
+`environment` section:
+
+ environment:
+ global:
+ # Don't build with openmp because it isn't supported in appveyor's compilers
+ USE_OMP: "0"
+
+Changelog
+---------
+v1.3.0 : Keyword argument "mask" added to "qeury" method. OpenMP compilation now works for MS Visual Studio compiler
+
+v1.2.2 : Build process fixes
+
+v1.2.1 : Fixed OpenMP thread safety issue introduced in v1.2.0
+
+v1.2.0 : 64 and 32 bit MSVC Windows support added
+
+v1.1.1 : Same as v1.1 release due to incorrect pypi release
+
+v1.1 : Build process improvements. Add data attribute to kdtree class for scipy interface compatibility
+
+v1.0 : Switched license from GPLv3 to LGPLv3
+
+v0.3 : Avoid zipping of installed egg
+
+v0.2 : Reduced memory footprint. Can now handle single precision data internally avoiding copy conversion to double precision. Default leafsize changed from 10 to 16 as this reduces the memory footprint and makes it a cache line multiplum (negligible if any query performance observed in benchmarks). Reduced memory allocation for leaf nodes. Applied patch for building on OS X.
+
+v0.1 : Initial version.
=====================================
debian/changelog
=====================================
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,11 +1,14 @@
-pykdtree (1.2.2-5) UNRELEASED; urgency=medium
+pykdtree (1.3.0-1) unstable; urgency=medium
+ * New upstream release.
* Update Vcs-* URLs for Salsa.
* Bump Standards-Version to 4.1.4, no changes.
* Change Homepage & Source URL to GitHub repository.
* Strip trailing whitespace from rules files.
+ * debian/patches
+ - drop 0002-safe-clean.patch, no longer necessary
- -- Antonio Valentino <antonio.valentino at tiscali.it> Sat, 31 Mar 2018 12:50:49 +0200
+ -- Antonio Valentino <antonio.valentino at tiscali.it> Mon, 04 Jun 2018 19:37:00 +0000
pykdtree (1.2.2-4) unstable; urgency=medium
=====================================
debian/patches/0002-safe-clean.patch deleted
=====================================
--- a/debian/patches/0002-safe-clean.patch
+++ /dev/null
@@ -1,17 +0,0 @@
-From: Antonio Valentino <antonio.valentino at tiscali.it>
-Date: Sun, 22 Jun 2014 14:28:53 +0200
-Subject: safe clean
-
----
- pykdtree.egg-info/requires.txt | 3 ++-
- 1 file changed, 2 insertions(+), 1 deletion(-)
-
-diff --git a/pykdtree.egg-info/requires.txt b/pykdtree.egg-info/requires.txt
-index 296d654..b26e81f 100644
---- a/pykdtree.egg-info/requires.txt
-+++ b/pykdtree.egg-info/requires.txt
-@@ -1 +1,2 @@
--numpy
-\ No newline at end of file
-+numpy
-+
=====================================
debian/patches/series
=====================================
--- a/debian/patches/series
+++ b/debian/patches/series
@@ -1,2 +1 @@
0001-fix_egginfo_source.patch
-0002-safe-clean.patch
=====================================
pykdtree.egg-info/PKG-INFO
=====================================
--- a/pykdtree.egg-info/PKG-INFO
+++ b/pykdtree.egg-info/PKG-INFO
@@ -1,6 +1,6 @@
Metadata-Version: 1.1
Name: pykdtree
-Version: 1.2.2
+Version: 1.3.0
Summary: Fast kd-tree implementation with OpenMP-enabled queries
Home-page: UNKNOWN
Author: Esben S. Nielsen
=====================================
pykdtree.egg-info/requires.txt
=====================================
--- a/pykdtree.egg-info/requires.txt
+++ b/pykdtree.egg-info/requires.txt
@@ -1 +1 @@
-numpy
\ No newline at end of file
+numpy
=====================================
pykdtree/_kdtree_core.c
=====================================
--- a/pykdtree/_kdtree_core.c
+++ b/pykdtree/_kdtree_core.c
@@ -1,6 +1,6 @@
/*
pykdtree, Fast kd-tree implementation with OpenMP-enabled queries
-
+
Copyright (C) 2013 - present Esben S. Nielsen
This program is free software: you can redistribute it and/or modify it under
@@ -53,7 +53,7 @@ typedef struct
float *bbox;
int8_t no_dims;
uint32_t *pidx;
- struct Node_float *root;
+ struct Node_float *root;
} Tree_float;
@@ -74,14 +74,14 @@ typedef struct
double *bbox;
int8_t no_dims;
uint32_t *pidx;
- struct Node_double *root;
+ struct Node_double *root;
} Tree_double;
void insert_point_float(uint32_t *closest_idx, float *closest_dist, uint32_t pidx, float cur_dist, uint32_t k);
void get_bounding_box_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, float *bbox);
-int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *bbox, int8_t *cut_dim,
+int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *bbox, int8_t *cut_dim,
float *cut_val, uint32_t *n_lo);
Tree_float* construct_tree_float(float *pa, int8_t no_dims, uint32_t n, uint32_t bsp);
Node_float* construct_subtree_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, float *bbox);
@@ -92,18 +92,20 @@ void print_tree_float(Node_float *root, int level);
float calc_dist_float(float *point1_coord, float *point2_coord, int8_t no_dims);
float get_cube_offset_float(int8_t dim, float *point_coord, float *bbox);
float get_min_dist_float(float *point_coord, int8_t no_dims, float *bbox);
-void search_leaf_float(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
+void search_leaf_float(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
uint32_t k, uint32_t *restrict closest_idx, float *restrict closest_dist);
-void search_splitnode_float(Node_float *root, float *pa, uint32_t *pidx, int8_t no_dims, float *point_coord,
- float min_dist, uint32_t k, float distance_upper_bound, float eps_fac, uint32_t * closest_idx, float *closest_dist);
-void search_tree_float(Tree_float *tree, float *pa, float *point_coords,
- uint32_t num_points, uint32_t k, float distance_upper_bound,
- float eps, uint32_t *closest_idxs, float *closest_dists);
+void search_leaf_float_mask(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
+ uint32_t k, uint8_t *restrict mask, uint32_t *restrict closest_idx, float *restrict closest_dist);
+void search_splitnode_float(Node_float *root, float *pa, uint32_t *pidx, int8_t no_dims, float *point_coord,
+ float min_dist, uint32_t k, float distance_upper_bound, float eps_fac, uint8_t *mask, uint32_t * closest_idx, float *closest_dist);
+void search_tree_float(Tree_float *tree, float *pa, float *point_coords,
+ uint32_t num_points, uint32_t k, float distance_upper_bound,
+ float eps, uint8_t *mask, uint32_t *closest_idxs, float *closest_dists);
void insert_point_double(uint32_t *closest_idx, double *closest_dist, uint32_t pidx, double cur_dist, uint32_t k);
void get_bounding_box_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t n, double *bbox);
-int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *bbox, int8_t *cut_dim,
+int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *bbox, int8_t *cut_dim,
double *cut_val, uint32_t *n_lo);
Tree_double* construct_tree_double(double *pa, int8_t no_dims, uint32_t n, uint32_t bsp);
Node_double* construct_subtree_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, uint32_t bsp, double *bbox);
@@ -114,13 +116,15 @@ void print_tree_double(Node_double *root, int level);
double calc_dist_double(double *point1_coord, double *point2_coord, int8_t no_dims);
double get_cube_offset_double(int8_t dim, double *point_coord, double *bbox);
double get_min_dist_double(double *point_coord, int8_t no_dims, double *bbox);
-void search_leaf_double(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
+void search_leaf_double(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
uint32_t k, uint32_t *restrict closest_idx, double *restrict closest_dist);
-void search_splitnode_double(Node_double *root, double *pa, uint32_t *pidx, int8_t no_dims, double *point_coord,
- double min_dist, uint32_t k, double distance_upper_bound, double eps_fac, uint32_t * closest_idx, double *closest_dist);
-void search_tree_double(Tree_double *tree, double *pa, double *point_coords,
- uint32_t num_points, uint32_t k, double distance_upper_bound,
- double eps, uint32_t *closest_idxs, double *closest_dists);
+void search_leaf_double_mask(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
+ uint32_t k, uint8_t *restrict mask, uint32_t *restrict closest_idx, double *restrict closest_dist);
+void search_splitnode_double(Node_double *root, double *pa, uint32_t *pidx, int8_t no_dims, double *point_coord,
+ double min_dist, uint32_t k, double distance_upper_bound, double eps_fac, uint8_t *mask, uint32_t * closest_idx, double *closest_dist);
+void search_tree_double(Tree_double *tree, double *pa, double *point_coords,
+ uint32_t num_points, uint32_t k, double distance_upper_bound,
+ double eps, uint8_t *mask, uint32_t *closest_idxs, double *closest_dists);
@@ -143,7 +147,7 @@ void insert_point_float(uint32_t *closest_idx, float *closest_dist, uint32_t pid
closest_dist[i] = closest_dist[i - 1];
closest_idx[i] = closest_idx[i - 1];
}
- else
+ else
{
break;
}
@@ -204,7 +208,7 @@ Params:
bbox : bounding box of data points
cut_dim : dimension used for partition (return)
cut_val : value of cutting point (return)
- n_lo : number of point below cutting plane (return)
+ n_lo : number of point below cutting plane (return)
************************************************/
int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *bbox, int8_t *cut_dim, float *cut_val, uint32_t *n_lo)
{
@@ -212,7 +216,7 @@ int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_id
uint32_t p, q, i2;
float size = 0, min_val, max_val, split, side_len, cur_val;
uint32_t end_idx = start_idx + n - 1;
-
+
/* Find largest bounding box side */
for (i = 0; i < no_dims; i++)
{
@@ -223,17 +227,17 @@ int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_id
size = side_len;
}
}
-
+
min_val = bbox[2 * dim];
max_val = bbox[2 * dim + 1];
-
+
/* Check for zero length or inconsistent */
if (min_val >= max_val)
return 1;
-
- /* Use middle for splitting */
+
+ /* Use middle for splitting */
split = (min_val + max_val) / 2;
-
+
/* Partition all data points around middle */
p = start_idx;
q = end_idx;
@@ -246,37 +250,37 @@ int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_id
else if (PA(q, dim) >= split)
{
/* Guard for underflow */
- if (q > 0)
+ if (q > 0)
{
q--;
}
else
{
- break;
- }
+ break;
+ }
}
else
{
PASWAP(p, q);
p++;
- q--;
- }
+ q--;
+ }
}
/* Check for empty splits */
if (p == start_idx)
{
- /* No points less than split.
+ /* No points less than split.
Split at lowest point instead.
Minimum 1 point will be in lower box.
*/
-
+
uint32_t j = start_idx;
split = PA(j, dim);
- for (i2 = start_idx + 1; i2 <= end_idx; i2++)
+ for (i2 = start_idx + 1; i2 <= end_idx; i2++)
{
/* Find lowest point */
- cur_val = PA(i2, dim);
+ cur_val = PA(i2, dim);
if (cur_val < split)
{
j = i2;
@@ -288,11 +292,11 @@ int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_id
}
else if (p == end_idx + 1)
{
- /* No points greater than split.
+ /* No points greater than split.
Split at highest point instead.
Minimum 1 point will be in higher box.
*/
-
+
uint32_t j = end_idx;
split = PA(j, dim);
for (i2 = start_idx; i2 < end_idx; i2++)
@@ -303,12 +307,12 @@ int partition_float(float *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_id
{
j = i2;
split = cur_val;
- }
+ }
}
PASWAP(j, end_idx);
- p = end_idx;
+ p = end_idx;
}
-
+
/* Set return values */
*cut_dim = dim;
*cut_val = split;
@@ -337,9 +341,9 @@ Node_float* construct_subtree_float(float *pa, uint32_t *pidx, int8_t no_dims, u
uint32_t n_lo;
float cut_val, lv, hv;
if (is_leaf)
- {
+ {
/* Make leaf node */
- root->cut_dim = -1;
+ root->cut_dim = -1;
}
else
{
@@ -349,30 +353,30 @@ Node_float* construct_subtree_float(float *pa, uint32_t *pidx, int8_t no_dims, u
if (rval == 1)
{
root->cut_dim = -1;
- return root;
+ return root;
}
root->cut_val = cut_val;
- root->cut_dim = cut_dim;
-
+ root->cut_dim = cut_dim;
+
/* Recurse on both subsets */
lv = bbox[2 * cut_dim];
hv = bbox[2 * cut_dim + 1];
-
+
/* Set bounds for cut dimension */
root->cut_bounds_lv = lv;
root->cut_bounds_hv = hv;
-
+
/* Update bounding box before call to lower subset and restore after */
bbox[2 * cut_dim + 1] = cut_val;
- root->left_child = (struct Node_float *)construct_subtree_float(pa, pidx, no_dims, start_idx, n_lo, bsp, bbox);
+ root->left_child = (struct Node_float *)construct_subtree_float(pa, pidx, no_dims, start_idx, n_lo, bsp, bbox);
bbox[2 * cut_dim + 1] = hv;
-
+
/* Update bounding box before call to higher subset and restore after */
bbox[2 * cut_dim] = cut_val;
- root->right_child = (struct Node_float *)construct_subtree_float(pa, pidx, no_dims, start_idx + n_lo, n - n_lo, bsp, bbox);
- bbox[2 * cut_dim] = lv;
+ root->right_child = (struct Node_float *)construct_subtree_float(pa, pidx, no_dims, start_idx + n_lo, n - n_lo, bsp, bbox);
+ bbox[2 * cut_dim] = lv;
}
- return root;
+ return root;
}
/************************************************
@@ -389,16 +393,16 @@ Tree_float* construct_tree_float(float *pa, int8_t no_dims, uint32_t n, uint32_t
uint32_t i;
uint32_t *pidx;
float *bbox;
-
+
tree->no_dims = no_dims;
-
+
/* Initialize permutation array */
pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
for (i = 0; i < n; i++)
{
pidx[i] = i;
}
-
+
bbox = (float *)malloc(2 * sizeof(float) * no_dims);
get_bounding_box_float(pa, pidx, no_dims, n, bbox);
tree->bbox = bbox;
@@ -414,16 +418,16 @@ Tree_float* construct_tree_float(float *pa, int8_t no_dims, uint32_t n, uint32_t
Create a tree node.
Params:
start_idx : index of first data point to use
- n : number of data points
+ n : number of data points
************************************************/
Node_float* create_node_float(uint32_t start_idx, uint32_t n, int is_leaf)
{
- Node_float *new_node;
+ Node_float *new_node;
if (is_leaf)
{
/*
Allocate only the part of the struct that will be used in a leaf node.
- This relies on the C99 specification of struct layout conservation and padding and
+ This relies on the C99 specification of struct layout conservation and padding and
that dereferencing is never attempted for the node pointers in a leaf.
*/
new_node = (Node_float *)malloc(sizeof(Node_float) - 2 * sizeof(Node_float *));
@@ -490,7 +494,7 @@ Params:
************************************************/
float calc_dist_float(float *point1_coord, float *point2_coord, int8_t no_dims)
{
- /* Calculate squared distance */
+ /* Calculate squared distance */
float dist = 0, dim_dist;
int8_t i;
for (i = 0; i < no_dims; i++)
@@ -511,11 +515,11 @@ Params:
float get_cube_offset_float(int8_t dim, float *point_coord, float *bbox)
{
float dim_coord = point_coord[dim];
-
+
if (dim_coord < bbox[2 * dim])
{
/* Left of cube in dimension */
- return dim_coord - bbox[2 * dim];
+ return dim_coord - bbox[2 * dim];
}
else if (dim_coord > bbox[2 * dim + 1])
{
@@ -526,7 +530,7 @@ float get_cube_offset_float(int8_t dim, float *point_coord, float *bbox)
{
/* Inside cube in dimension */
return 0.;
- }
+ }
}
/************************************************
@@ -544,7 +548,7 @@ float get_min_dist_float(float *point_coord, int8_t no_dims, float *bbox)
for (i = 0; i < no_dims; i++)
{
cube_offset_dim = get_cube_offset_float(i, point_coord, bbox);
- cube_offset += cube_offset_dim * cube_offset_dim;
+ cube_offset += cube_offset_dim * cube_offset_dim;
}
return cube_offset;
@@ -555,21 +559,58 @@ Search a leaf node for closest point
Params:
pa : data points
pidx : permutation index of data points
- no_dims : number of dimensions
+ no_dims : number of dimensions
start_idx : index of first data point to use
size : number of data points
point_coord : query point
closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
+ closest_dist : distance to closest point (return)
************************************************/
-void search_leaf_float(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
+void search_leaf_float(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
uint32_t k, uint32_t *restrict closest_idx, float *restrict closest_dist)
{
float cur_dist;
uint32_t i;
- /* Loop through all points in leaf */
+ /* Loop through all points in leaf */
+ for (i = 0; i < n; i++)
+ {
+ /* Get distance to query point */
+ cur_dist = calc_dist_float(&PA(start_idx + i, 0), point_coord, no_dims);
+ /* Update closest info if new point is closest so far*/
+ if (cur_dist < closest_dist[k - 1])
+ {
+ insert_point_float(closest_idx, closest_dist, pidx[start_idx + i], cur_dist, k);
+ }
+ }
+}
+
+
+/************************************************
+Search a leaf node for closest point with data point mask
+Params:
+ pa : data points
+ pidx : permutation index of data points
+ no_dims : number of dimensions
+ start_idx : index of first data point to use
+ size : number of data points
+ point_coord : query point
+ mask : boolean array of invalid (True) and valid (False) data points
+ closest_idx : index of closest data point found (return)
+ closest_dist : distance to closest point (return)
+************************************************/
+void search_leaf_float_mask(float *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, float *restrict point_coord,
+ uint32_t k, uint8_t *mask, uint32_t *restrict closest_idx, float *restrict closest_dist)
+{
+ float cur_dist;
+ uint32_t i;
+ /* Loop through all points in leaf */
for (i = 0; i < n; i++)
{
+ /* Is this point masked out? */
+ if (mask[start_idx + i])
+ {
+ continue;
+ }
/* Get distance to query point */
cur_dist = calc_dist_float(&PA(start_idx + i, 0), point_coord, no_dims);
/* Update closest info if new point is closest so far*/
@@ -589,17 +630,19 @@ Params:
no_dims : number of dimensions
point_coord : query point
min_dist : minumum distance to nearest neighbour
+ mask : boolean array of invalid (True) and valid (False) data points
closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
+ closest_dist : distance to closest point (return)
************************************************/
void search_splitnode_float(Node_float *root, float *pa, uint32_t *pidx, int8_t no_dims, float *point_coord,
- float min_dist, uint32_t k, float distance_upper_bound, float eps_fac, uint32_t *closest_idx, float *closest_dist)
+ float min_dist, uint32_t k, float distance_upper_bound, float eps_fac, uint8_t *mask,
+ uint32_t *closest_idx, float *closest_dist)
{
int8_t dim;
float dist_left, dist_right;
float new_offset;
float box_diff;
-
+
/* Skip if distance bound exeeded */
if (min_dist > distance_upper_bound)
{
@@ -607,40 +650,47 @@ void search_splitnode_float(Node_float *root, float *pa, uint32_t *pidx, int8_t
}
dim = root->cut_dim;
-
+
/* Handle leaf node */
if (dim == -1)
{
- search_leaf_float(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, closest_idx, closest_dist);
+ if (mask)
+ {
+ search_leaf_float_mask(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, mask, closest_idx, closest_dist);
+ }
+ else
+ {
+ search_leaf_float(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, closest_idx, closest_dist);
+ }
return;
}
-
+
/* Get distance to cutting plane */
new_offset = point_coord[dim] - root->cut_val;
-
+
if (new_offset < 0)
{
/* Left of cutting plane */
dist_left = min_dist;
if (dist_left < closest_dist[k - 1] * eps_fac)
{
- /* Search left subtree if minimum distance is below limit */
- search_splitnode_float((Node_float *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, closest_idx, closest_dist);
+ /* Search left subtree if minimum distance is below limit */
+ search_splitnode_float((Node_float *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
-
- /* Right of cutting plane. Update minimum distance.
+
+ /* Right of cutting plane. Update minimum distance.
See Algorithms for Fast Vector Quantization
Sunil Arya and David M. Mount. */
box_diff = root->cut_bounds_lv - point_coord[dim];
if (box_diff < 0)
- {
+ {
box_diff = 0;
}
dist_right = min_dist - box_diff * box_diff + new_offset * new_offset;
if (dist_right < closest_dist[k - 1] * eps_fac)
{
- /* Search right subtree if minimum distance is below limit*/
- search_splitnode_float((Node_float *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, closest_idx, closest_dist);
+ /* Search right subtree if minimum distance is below limit*/
+ search_splitnode_float((Node_float *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
}
else
@@ -649,23 +699,23 @@ void search_splitnode_float(Node_float *root, float *pa, uint32_t *pidx, int8_t
dist_right = min_dist;
if (dist_right < closest_dist[k - 1] * eps_fac)
{
- /* Search right subtree if minimum distance is below limit*/
- search_splitnode_float((Node_float *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, closest_idx, closest_dist);
+ /* Search right subtree if minimum distance is below limit*/
+ search_splitnode_float((Node_float *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
-
+
/* Left of cutting plane. Update minimum distance.
See Algorithms for Fast Vector Quantization
Sunil Arya and David M. Mount. */
box_diff = point_coord[dim] - root->cut_bounds_hv;
if (box_diff < 0)
- {
+ {
box_diff = 0;
}
dist_left = min_dist - box_diff * box_diff + new_offset * new_offset;
if (dist_left < closest_dist[k - 1] * eps_fac)
{
- /* Search left subtree if minimum distance is below limit*/
- search_splitnode_float((Node_float *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, closest_idx, closest_dist);
+ /* Search left subtree if minimum distance is below limit*/
+ search_splitnode_float((Node_float *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
}
}
@@ -678,29 +728,37 @@ Params:
pidx : permutation index of data points
point_coords : query points
num_points : number of query points
+ mask : boolean array of invalid (True) and valid (False) data points
closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
+ closest_dist : distance to closest point (return)
************************************************/
-void search_tree_float(Tree_float *tree, float *pa, float *point_coords,
+void search_tree_float(Tree_float *tree, float *pa, float *point_coords,
uint32_t num_points, uint32_t k, float distance_upper_bound,
- float eps, uint32_t *closest_idxs, float *closest_dists)
+ float eps, uint8_t *mask, uint32_t *closest_idxs, float *closest_dists)
{
float min_dist;
float eps_fac = 1 / ((1 + eps) * (1 + eps));
int8_t no_dims = tree->no_dims;
float *bbox = tree->bbox;
uint32_t *pidx = tree->pidx;
- uint32_t i, j;
+ uint32_t j = 0;
+#if defined(_MSC_VER) && defined(_OPENMP)
+ int32_t i = 0;
+ int32_t local_num_points = (int32_t) num_points;
+#else
+ uint32_t i;
+ uint32_t local_num_points = num_points;
+#endif
Node_float *root = (Node_float *)tree->root;
/* Queries are OpenMP enabled */
#pragma omp parallel
{
- /* The low chunk size is important to avoid L2 cache trashing
+ /* The low chunk size is important to avoid L2 cache trashing
for spatial coherent query datasets
*/
#pragma omp for private(i, j) schedule(static, 100) nowait
- for (i = 0; i < num_points; i++)
+ for (i = 0; i < local_num_points; i++)
{
for (j = 0; j < k; j++)
{
@@ -709,7 +767,7 @@ void search_tree_float(Tree_float *tree, float *pa, float *point_coords,
}
min_dist = get_min_dist_float(point_coords + no_dims * i, no_dims, bbox);
search_splitnode_float(root, pa, pidx, no_dims, point_coords + no_dims * i, min_dist,
- k, distance_upper_bound, eps_fac, &closest_idxs[i * k], &closest_dists[i * k]);
+ k, distance_upper_bound, eps_fac, mask, &closest_idxs[i * k], &closest_dists[i * k]);
}
}
}
@@ -733,7 +791,7 @@ void insert_point_double(uint32_t *closest_idx, double *closest_dist, uint32_t p
closest_dist[i] = closest_dist[i - 1];
closest_idx[i] = closest_idx[i - 1];
}
- else
+ else
{
break;
}
@@ -794,7 +852,7 @@ Params:
bbox : bounding box of data points
cut_dim : dimension used for partition (return)
cut_val : value of cutting point (return)
- n_lo : number of point below cutting plane (return)
+ n_lo : number of point below cutting plane (return)
************************************************/
int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *bbox, int8_t *cut_dim, double *cut_val, uint32_t *n_lo)
{
@@ -802,7 +860,7 @@ int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_
uint32_t p, q, i2;
double size = 0, min_val, max_val, split, side_len, cur_val;
uint32_t end_idx = start_idx + n - 1;
-
+
/* Find largest bounding box side */
for (i = 0; i < no_dims; i++)
{
@@ -813,17 +871,17 @@ int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_
size = side_len;
}
}
-
+
min_val = bbox[2 * dim];
max_val = bbox[2 * dim + 1];
-
+
/* Check for zero length or inconsistent */
if (min_val >= max_val)
return 1;
-
- /* Use middle for splitting */
+
+ /* Use middle for splitting */
split = (min_val + max_val) / 2;
-
+
/* Partition all data points around middle */
p = start_idx;
q = end_idx;
@@ -836,37 +894,37 @@ int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_
else if (PA(q, dim) >= split)
{
/* Guard for underflow */
- if (q > 0)
+ if (q > 0)
{
q--;
}
else
{
- break;
- }
+ break;
+ }
}
else
{
PASWAP(p, q);
p++;
- q--;
- }
+ q--;
+ }
}
/* Check for empty splits */
if (p == start_idx)
{
- /* No points less than split.
+ /* No points less than split.
Split at lowest point instead.
Minimum 1 point will be in lower box.
*/
-
+
uint32_t j = start_idx;
split = PA(j, dim);
- for (i2 = start_idx + 1; i2 <= end_idx; i2++)
+ for (i2 = start_idx + 1; i2 <= end_idx; i2++)
{
/* Find lowest point */
- cur_val = PA(i2, dim);
+ cur_val = PA(i2, dim);
if (cur_val < split)
{
j = i2;
@@ -878,11 +936,11 @@ int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_
}
else if (p == end_idx + 1)
{
- /* No points greater than split.
+ /* No points greater than split.
Split at highest point instead.
Minimum 1 point will be in higher box.
*/
-
+
uint32_t j = end_idx;
split = PA(j, dim);
for (i2 = start_idx; i2 < end_idx; i2++)
@@ -893,12 +951,12 @@ int partition_double(double *pa, uint32_t *pidx, int8_t no_dims, uint32_t start_
{
j = i2;
split = cur_val;
- }
+ }
}
PASWAP(j, end_idx);
- p = end_idx;
+ p = end_idx;
}
-
+
/* Set return values */
*cut_dim = dim;
*cut_val = split;
@@ -927,9 +985,9 @@ Node_double* construct_subtree_double(double *pa, uint32_t *pidx, int8_t no_dims
uint32_t n_lo;
double cut_val, lv, hv;
if (is_leaf)
- {
+ {
/* Make leaf node */
- root->cut_dim = -1;
+ root->cut_dim = -1;
}
else
{
@@ -939,30 +997,30 @@ Node_double* construct_subtree_double(double *pa, uint32_t *pidx, int8_t no_dims
if (rval == 1)
{
root->cut_dim = -1;
- return root;
+ return root;
}
root->cut_val = cut_val;
- root->cut_dim = cut_dim;
-
+ root->cut_dim = cut_dim;
+
/* Recurse on both subsets */
lv = bbox[2 * cut_dim];
hv = bbox[2 * cut_dim + 1];
-
+
/* Set bounds for cut dimension */
root->cut_bounds_lv = lv;
root->cut_bounds_hv = hv;
-
+
/* Update bounding box before call to lower subset and restore after */
bbox[2 * cut_dim + 1] = cut_val;
- root->left_child = (struct Node_double *)construct_subtree_double(pa, pidx, no_dims, start_idx, n_lo, bsp, bbox);
+ root->left_child = (struct Node_double *)construct_subtree_double(pa, pidx, no_dims, start_idx, n_lo, bsp, bbox);
bbox[2 * cut_dim + 1] = hv;
-
+
/* Update bounding box before call to higher subset and restore after */
bbox[2 * cut_dim] = cut_val;
- root->right_child = (struct Node_double *)construct_subtree_double(pa, pidx, no_dims, start_idx + n_lo, n - n_lo, bsp, bbox);
- bbox[2 * cut_dim] = lv;
+ root->right_child = (struct Node_double *)construct_subtree_double(pa, pidx, no_dims, start_idx + n_lo, n - n_lo, bsp, bbox);
+ bbox[2 * cut_dim] = lv;
}
- return root;
+ return root;
}
/************************************************
@@ -979,16 +1037,16 @@ Tree_double* construct_tree_double(double *pa, int8_t no_dims, uint32_t n, uint3
uint32_t i;
uint32_t *pidx;
double *bbox;
-
+
tree->no_dims = no_dims;
-
+
/* Initialize permutation array */
pidx = (uint32_t *)malloc(sizeof(uint32_t) * n);
for (i = 0; i < n; i++)
{
pidx[i] = i;
}
-
+
bbox = (double *)malloc(2 * sizeof(double) * no_dims);
get_bounding_box_double(pa, pidx, no_dims, n, bbox);
tree->bbox = bbox;
@@ -1004,16 +1062,16 @@ Tree_double* construct_tree_double(double *pa, int8_t no_dims, uint32_t n, uint3
Create a tree node.
Params:
start_idx : index of first data point to use
- n : number of data points
+ n : number of data points
************************************************/
Node_double* create_node_double(uint32_t start_idx, uint32_t n, int is_leaf)
{
- Node_double *new_node;
+ Node_double *new_node;
if (is_leaf)
{
/*
Allocate only the part of the struct that will be used in a leaf node.
- This relies on the C99 specification of struct layout conservation and padding and
+ This relies on the C99 specification of struct layout conservation and padding and
that dereferencing is never attempted for the node pointers in a leaf.
*/
new_node = (Node_double *)malloc(sizeof(Node_double) - 2 * sizeof(Node_double *));
@@ -1080,7 +1138,7 @@ Params:
************************************************/
double calc_dist_double(double *point1_coord, double *point2_coord, int8_t no_dims)
{
- /* Calculate squared distance */
+ /* Calculate squared distance */
double dist = 0, dim_dist;
int8_t i;
for (i = 0; i < no_dims; i++)
@@ -1101,11 +1159,11 @@ Params:
double get_cube_offset_double(int8_t dim, double *point_coord, double *bbox)
{
double dim_coord = point_coord[dim];
-
+
if (dim_coord < bbox[2 * dim])
{
/* Left of cube in dimension */
- return dim_coord - bbox[2 * dim];
+ return dim_coord - bbox[2 * dim];
}
else if (dim_coord > bbox[2 * dim + 1])
{
@@ -1116,7 +1174,7 @@ double get_cube_offset_double(int8_t dim, double *point_coord, double *bbox)
{
/* Inside cube in dimension */
return 0.;
- }
+ }
}
/************************************************
@@ -1134,7 +1192,7 @@ double get_min_dist_double(double *point_coord, int8_t no_dims, double *bbox)
for (i = 0; i < no_dims; i++)
{
cube_offset_dim = get_cube_offset_double(i, point_coord, bbox);
- cube_offset += cube_offset_dim * cube_offset_dim;
+ cube_offset += cube_offset_dim * cube_offset_dim;
}
return cube_offset;
@@ -1145,19 +1203,19 @@ Search a leaf node for closest point
Params:
pa : data points
pidx : permutation index of data points
- no_dims : number of dimensions
+ no_dims : number of dimensions
start_idx : index of first data point to use
size : number of data points
point_coord : query point
closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
+ closest_dist : distance to closest point (return)
************************************************/
-void search_leaf_double(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
+void search_leaf_double(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
uint32_t k, uint32_t *restrict closest_idx, double *restrict closest_dist)
{
double cur_dist;
uint32_t i;
- /* Loop through all points in leaf */
+ /* Loop through all points in leaf */
for (i = 0; i < n; i++)
{
/* Get distance to query point */
@@ -1170,6 +1228,43 @@ void search_leaf_double(double *restrict pa, uint32_t *restrict pidx, int8_t no_
}
}
+
+/************************************************
+Search a leaf node for closest point with data point mask
+Params:
+ pa : data points
+ pidx : permutation index of data points
+ no_dims : number of dimensions
+ start_idx : index of first data point to use
+ size : number of data points
+ point_coord : query point
+ mask : boolean array of invalid (True) and valid (False) data points
+ closest_idx : index of closest data point found (return)
+ closest_dist : distance to closest point (return)
+************************************************/
+void search_leaf_double_mask(double *restrict pa, uint32_t *restrict pidx, int8_t no_dims, uint32_t start_idx, uint32_t n, double *restrict point_coord,
+ uint32_t k, uint8_t *mask, uint32_t *restrict closest_idx, double *restrict closest_dist)
+{
+ double cur_dist;
+ uint32_t i;
+ /* Loop through all points in leaf */
+ for (i = 0; i < n; i++)
+ {
+ /* Is this point masked out? */
+ if (mask[start_idx + i])
+ {
+ continue;
+ }
+ /* Get distance to query point */
+ cur_dist = calc_dist_double(&PA(start_idx + i, 0), point_coord, no_dims);
+ /* Update closest info if new point is closest so far*/
+ if (cur_dist < closest_dist[k - 1])
+ {
+ insert_point_double(closest_idx, closest_dist, pidx[start_idx + i], cur_dist, k);
+ }
+ }
+}
+
/************************************************
Search subtree for nearest to query point
Params:
@@ -1179,17 +1274,19 @@ Params:
no_dims : number of dimensions
point_coord : query point
min_dist : minumum distance to nearest neighbour
+ mask : boolean array of invalid (True) and valid (False) data points
closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
+ closest_dist : distance to closest point (return)
************************************************/
void search_splitnode_double(Node_double *root, double *pa, uint32_t *pidx, int8_t no_dims, double *point_coord,
- double min_dist, uint32_t k, double distance_upper_bound, double eps_fac, uint32_t *closest_idx, double *closest_dist)
+ double min_dist, uint32_t k, double distance_upper_bound, double eps_fac, uint8_t *mask,
+ uint32_t *closest_idx, double *closest_dist)
{
int8_t dim;
double dist_left, dist_right;
double new_offset;
double box_diff;
-
+
/* Skip if distance bound exeeded */
if (min_dist > distance_upper_bound)
{
@@ -1197,40 +1294,47 @@ void search_splitnode_double(Node_double *root, double *pa, uint32_t *pidx, int8
}
dim = root->cut_dim;
-
+
/* Handle leaf node */
if (dim == -1)
{
- search_leaf_double(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, closest_idx, closest_dist);
+ if (mask)
+ {
+ search_leaf_double_mask(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, mask, closest_idx, closest_dist);
+ }
+ else
+ {
+ search_leaf_double(pa, pidx, no_dims, root->start_idx, root->n, point_coord, k, closest_idx, closest_dist);
+ }
return;
}
-
+
/* Get distance to cutting plane */
new_offset = point_coord[dim] - root->cut_val;
-
+
if (new_offset < 0)
{
/* Left of cutting plane */
dist_left = min_dist;
if (dist_left < closest_dist[k - 1] * eps_fac)
{
- /* Search left subtree if minimum distance is below limit */
- search_splitnode_double((Node_double *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, closest_idx, closest_dist);
+ /* Search left subtree if minimum distance is below limit */
+ search_splitnode_double((Node_double *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
-
- /* Right of cutting plane. Update minimum distance.
+
+ /* Right of cutting plane. Update minimum distance.
See Algorithms for Fast Vector Quantization
Sunil Arya and David M. Mount. */
box_diff = root->cut_bounds_lv - point_coord[dim];
if (box_diff < 0)
- {
+ {
box_diff = 0;
}
dist_right = min_dist - box_diff * box_diff + new_offset * new_offset;
if (dist_right < closest_dist[k - 1] * eps_fac)
{
- /* Search right subtree if minimum distance is below limit*/
- search_splitnode_double((Node_double *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, closest_idx, closest_dist);
+ /* Search right subtree if minimum distance is below limit*/
+ search_splitnode_double((Node_double *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
}
else
@@ -1239,23 +1343,23 @@ void search_splitnode_double(Node_double *root, double *pa, uint32_t *pidx, int8
dist_right = min_dist;
if (dist_right < closest_dist[k - 1] * eps_fac)
{
- /* Search right subtree if minimum distance is below limit*/
- search_splitnode_double((Node_double *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, closest_idx, closest_dist);
+ /* Search right subtree if minimum distance is below limit*/
+ search_splitnode_double((Node_double *)root->right_child, pa, pidx, no_dims, point_coord, dist_right, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
-
+
/* Left of cutting plane. Update minimum distance.
See Algorithms for Fast Vector Quantization
Sunil Arya and David M. Mount. */
box_diff = point_coord[dim] - root->cut_bounds_hv;
if (box_diff < 0)
- {
+ {
box_diff = 0;
}
dist_left = min_dist - box_diff * box_diff + new_offset * new_offset;
if (dist_left < closest_dist[k - 1] * eps_fac)
{
- /* Search left subtree if minimum distance is below limit*/
- search_splitnode_double((Node_double *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, closest_idx, closest_dist);
+ /* Search left subtree if minimum distance is below limit*/
+ search_splitnode_double((Node_double *)root->left_child, pa, pidx, no_dims, point_coord, dist_left, k, distance_upper_bound, eps_fac, mask, closest_idx, closest_dist);
}
}
}
@@ -1268,29 +1372,37 @@ Params:
pidx : permutation index of data points
point_coords : query points
num_points : number of query points
+ mask : boolean array of invalid (True) and valid (False) data points
closest_idx : index of closest data point found (return)
- closest_dist : distance to closest point (return)
+ closest_dist : distance to closest point (return)
************************************************/
-void search_tree_double(Tree_double *tree, double *pa, double *point_coords,
+void search_tree_double(Tree_double *tree, double *pa, double *point_coords,
uint32_t num_points, uint32_t k, double distance_upper_bound,
- double eps, uint32_t *closest_idxs, double *closest_dists)
+ double eps, uint8_t *mask, uint32_t *closest_idxs, double *closest_dists)
{
double min_dist;
double eps_fac = 1 / ((1 + eps) * (1 + eps));
int8_t no_dims = tree->no_dims;
double *bbox = tree->bbox;
uint32_t *pidx = tree->pidx;
- uint32_t i, j;
+ uint32_t j = 0;
+#if defined(_MSC_VER) && defined(_OPENMP)
+ int32_t i = 0;
+ int32_t local_num_points = (int32_t) num_points;
+#else
+ uint32_t i;
+ uint32_t local_num_points = num_points;
+#endif
Node_double *root = (Node_double *)tree->root;
/* Queries are OpenMP enabled */
#pragma omp parallel
{
- /* The low chunk size is important to avoid L2 cache trashing
+ /* The low chunk size is important to avoid L2 cache trashing
for spatial coherent query datasets
*/
#pragma omp for private(i, j) schedule(static, 100) nowait
- for (i = 0; i < num_points; i++)
+ for (i = 0; i < local_num_points; i++)
{
for (j = 0; j < k; j++)
{
@@ -1299,7 +1411,7 @@ void search_tree_double(Tree_double *tree, double *pa, double *point_coords,
}
min_dist = get_min_dist_double(point_coords + no_dims * i, no_dims, bbox);
search_splitnode_double(root, pa, pidx, no_dims, point_coords + no_dims * i, min_dist,
- k, distance_upper_bound, eps_fac, &closest_idxs[i * k], &closest_dists[i * k]);
+ k, distance_upper_bound, eps_fac, mask, &closest_idxs[i * k], &closest_dists[i * k]);
}
}
}
=====================================
pykdtree/kdtree.c
=====================================
The diff for this file was not included because it is too large.
=====================================
pykdtree/test_tree.py
=====================================
--- a/pykdtree/test_tree.py
+++ b/pykdtree/test_tree.py
@@ -103,9 +103,9 @@ data_pts_real = np.array([[ 790535.062, -369324.656, 6310963.5 ],
[ 750396.188, -49264.852, 6326458.5 ],
[ 750056.375, -46624.227, 6326519. ],
[ 749718.875, -43993.633, 6326578. ]])
-
+
def test1d():
-
+
data_pts = np.arange(1000)
kdtree = KDTree(data_pts, leafsize=15)
query_pts = np.arange(400, 300, -10)
@@ -113,16 +113,16 @@ def test1d():
assert idx[0] == 400
assert dist[0] == 0
assert idx[1] == 390
-
+
def test3d():
-
- #7, 93, 45
+
+ #7, 93, 45
query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
[751763.125, -59925.969, 6326205.5],
[769957.188, -202418.125, 6321069.5]])
-
-
+
+
kdtree = KDTree(data_pts_real)
dist, idx = kdtree.query(query_pts, sqr_dists=True)
@@ -136,13 +136,13 @@ def test3d():
def test3d_float32():
-
- #7, 93, 45
+
+ #7, 93, 45
query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
[751763.125, -59925.969, 6326205.5],
[769957.188, -202418.125, 6321069.5]], dtype=np.float32)
-
-
+
+
kdtree = KDTree(data_pts_real.astype(np.float32))
dist, idx = kdtree.query(query_pts, sqr_dists=True)
epsilon = 1e-5
@@ -156,50 +156,50 @@ def test3d_float32():
def test3d_float32_mismatch():
-
- #7, 93, 45
+
+ #7, 93, 45
query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
[751763.125, -59925.969, 6326205.5],
[769957.188, -202418.125, 6321069.5]], dtype=np.float32)
-
+
kdtree = KDTree(data_pts_real)
dist, idx = kdtree.query(query_pts, sqr_dists=True)
def test3d_float32_mismatch2():
-
- #7, 93, 45
+
+ #7, 93, 45
query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
[751763.125, -59925.969, 6326205.5],
[769957.188, -202418.125, 6321069.5]])
-
+
kdtree = KDTree(data_pts_real.astype(np.float32))
try:
dist, idx = kdtree.query(query_pts, sqr_dists=True)
assert False
except TypeError:
- assert True
-
+ assert True
+
def test3d_8n():
query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
[751763.125, -59925.969, 6326205.5],
[769957.188, -202418.125, 6321069.5]])
-
+
kdtree = KDTree(data_pts_real)
dist, idx = kdtree.query(query_pts, k=8)
-
+
exp_dist = np.array([[ 0.00000000e+00, 4.05250235e+03, 4.07389794e+03, 8.08201128e+03,
8.17063009e+03, 1.20904577e+04, 1.22902057e+04, 1.60775136e+04],
[ 1.73205081e+00, 2.70216896e+03, 2.71431274e+03, 5.39537066e+03,
5.43793210e+03, 8.07855631e+03, 8.17119970e+03, 1.07513693e+04],
[ 1.41424892e+02, 3.25500021e+03, 3.44284958e+03, 6.58019346e+03,
6.81038455e+03, 9.89140135e+03, 1.01918659e+04, 1.31892516e+04]])
-
+
exp_idx = np.array([[ 7, 8, 6, 9, 5, 10, 4, 11],
[93, 94, 92, 95, 91, 96, 90, 97],
[45, 46, 44, 47, 43, 48, 42, 49]])
-
+
assert np.array_equal(idx, exp_idx)
assert np.allclose(dist, exp_dist)
@@ -207,65 +207,65 @@ def test3d_8n_ub():
query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
[751763.125, -59925.969, 6326205.5],
[769957.188, -202418.125, 6321069.5]])
-
+
kdtree = KDTree(data_pts_real)
dist, idx = kdtree.query(query_pts, k=8, distance_upper_bound=10e3, sqr_dists=False)
-
+
exp_dist = np.array([[ 0.00000000e+00, 4.05250235e+03, 4.07389794e+03, 8.08201128e+03,
8.17063009e+03, np.Inf, np.Inf, np.Inf],
[ 1.73205081e+00, 2.70216896e+03, 2.71431274e+03, 5.39537066e+03,
5.43793210e+03, 8.07855631e+03, 8.17119970e+03, np.Inf],
[ 1.41424892e+02, 3.25500021e+03, 3.44284958e+03, 6.58019346e+03,
6.81038455e+03, 9.89140135e+03, np.Inf, np.Inf]])
- n = 100
+ n = 100
exp_idx = np.array([[ 7, 8, 6, 9, 5, n, n, n],
[93, 94, 92, 95, 91, 96, 90, n],
[45, 46, 44, 47, 43, 48, n, n]])
-
+
assert np.array_equal(idx, exp_idx)
assert np.allclose(dist, exp_dist)
-
+
def test3d_8n_ub_leaf20():
query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
[751763.125, -59925.969, 6326205.5],
[769957.188, -202418.125, 6321069.5]])
-
+
kdtree = KDTree(data_pts_real, leafsize=20)
dist, idx = kdtree.query(query_pts, k=8, distance_upper_bound=10e3, sqr_dists=False)
-
+
exp_dist = np.array([[ 0.00000000e+00, 4.05250235e+03, 4.07389794e+03, 8.08201128e+03,
8.17063009e+03, np.Inf, np.Inf, np.Inf],
[ 1.73205081e+00, 2.70216896e+03, 2.71431274e+03, 5.39537066e+03,
5.43793210e+03, 8.07855631e+03, 8.17119970e+03, np.Inf],
[ 1.41424892e+02, 3.25500021e+03, 3.44284958e+03, 6.58019346e+03,
6.81038455e+03, 9.89140135e+03, np.Inf, np.Inf]])
- n = 100
+ n = 100
exp_idx = np.array([[ 7, 8, 6, 9, 5, n, n, n],
[93, 94, 92, 95, 91, 96, 90, n],
[45, 46, 44, 47, 43, 48, n, n]])
-
+
assert np.array_equal(idx, exp_idx)
assert np.allclose(dist, exp_dist)
-
+
def test3d_8n_ub_eps():
query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
[751763.125, -59925.969, 6326205.5],
[769957.188, -202418.125, 6321069.5]])
-
+
kdtree = KDTree(data_pts_real)
dist, idx = kdtree.query(query_pts, k=8, eps=0.1, distance_upper_bound=10e3, sqr_dists=False)
-
+
exp_dist = np.array([[ 0.00000000e+00, 4.05250235e+03, 4.07389794e+03, 8.08201128e+03,
8.17063009e+03, np.Inf, np.Inf, np.Inf],
[ 1.73205081e+00, 2.70216896e+03, 2.71431274e+03, 5.39537066e+03,
5.43793210e+03, 8.07855631e+03, 8.17119970e+03, np.Inf],
[ 1.41424892e+02, 3.25500021e+03, 3.44284958e+03, 6.58019346e+03,
6.81038455e+03, 9.89140135e+03, np.Inf, np.Inf]])
- n = 100
+ n = 100
exp_idx = np.array([[ 7, 8, 6, 9, 5, n, n, n],
[93, 94, 92, 95, 91, 96, 90, n],
[45, 46, 44, 47, 43, 48, n, n]])
-
+
assert np.array_equal(idx, exp_idx)
assert np.allclose(dist, exp_dist)
@@ -291,12 +291,54 @@ def test3d_large_query():
assert np.all(abs(dist[2*n:] - 20001.) < epsilon * dist[2*n:])
def test_scipy_comp():
-
+
query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
[751763.125, -59925.969, 6326205.5],
[769957.188, -202418.125, 6321069.5]])
-
+
kdtree = KDTree(data_pts_real)
assert id(kdtree.data) == id(kdtree.data_pts)
+def test1d_mask():
+ data_pts = np.arange(1000)
+ kdtree = KDTree(data_pts, leafsize=15)
+ query_pts = np.arange(400, 300, -10)
+ query_mask = np.zeros(data_pts.shape[0]).astype(bool)
+ query_mask[400] = True
+ dist, idx = kdtree.query(query_pts, mask=query_mask)
+ assert idx[0] == 399 # would be 400 if no mask
+ assert dist[0] == 1.
+ assert idx[1] == 390
+
+
+def test1d_all_masked():
+ data_pts = np.arange(1000)
+ kdtree = KDTree(data_pts, leafsize=15)
+ query_pts = np.arange(400, 300, -10)
+ query_mask = np.ones(data_pts.shape[0]).astype(bool)
+ dist, idx = kdtree.query(query_pts, mask=query_mask)
+ # all invalid
+ assert np.all(i >= 1000 for i in idx)
+ assert np.all(d >= 1001 for d in dist)
+
+
+def test3d_mask():
+ #7, 93, 45
+ query_pts = np.array([[ 787014.438, -340616.906, 6313018.],
+ [751763.125, -59925.969, 6326205.5],
+ [769957.188, -202418.125, 6321069.5]])
+
+ kdtree = KDTree(data_pts_real)
+ query_mask = np.zeros(data_pts_real.shape[0])
+ query_mask[6:10] = True
+ dist, idx = kdtree.query(query_pts, sqr_dists=True, mask=query_mask)
+
+ epsilon = 1e-5
+ assert idx[0] == 5 # would be 7 if no mask
+ assert idx[1] == 93
+ assert idx[2] == 45
+ # would be 0 if no mask
+ assert abs(dist[0] - 66759196.1053) < epsilon * dist[0]
+ assert abs(dist[1] - 3.) < epsilon * dist[1]
+ assert abs(dist[2] - 20001.) < epsilon * dist[2]
=====================================
setup.cfg
=====================================
--- a/setup.cfg
+++ b/setup.cfg
@@ -5,5 +5,4 @@ release = 1
[egg_info]
tag_build =
tag_date = 0
-tag_svn_revision = 0
=====================================
setup.py
=====================================
--- a/setup.py
+++ b/setup.py
@@ -1,5 +1,5 @@
#pykdtree, Fast kd-tree implementation with OpenMP-enabled queries
-#
+#
#Copyright (C) 2013 - present Esben S. Nielsen
#
# This program is free software: you can redistribute it and/or modify it under
@@ -19,7 +19,7 @@ import os
from setuptools import setup, Extension
from setuptools.command.build_ext import build_ext
-# Get OpenMP setting from environment
+# Get OpenMP setting from environment
try:
use_omp = int(os.environ['USE_OMP'])
except KeyError:
@@ -36,9 +36,9 @@ def set_builtin(name, value):
# Custom builder to handler compiler flags. Edit if needed.
class build_ext_subclass(build_ext):
def build_extensions(self):
- comp = self.compiler.compiler_type
+ comp = self.compiler.compiler_type
if comp in ('unix', 'cygwin', 'mingw32'):
- # Check if build is with OpenMP
+ # Check if build is with OpenMP
if use_omp:
extra_compile_args = ['-std=c99', '-O3', '-fopenmp']
extra_link_args=['-lgomp']
@@ -69,11 +69,11 @@ class build_ext_subclass(build_ext):
set_builtin('__NUMPY_SETUP__', False)
import numpy
self.include_dirs.append(numpy.get_include())
-
-
+
+
setup(
name='pykdtree',
- version='1.2.2',
+ version='1.3.0',
description='Fast kd-tree implementation with OpenMP-enabled queries',
author='Esben S. Nielsen',
author_email='storpipfugl at gmail.com',
@@ -83,7 +83,7 @@ setup(
tests_require=['nose'],
zip_safe=False,
test_suite = 'nose.collector',
- ext_modules = [Extension('pykdtree.kdtree',
+ ext_modules = [Extension('pykdtree.kdtree',
['pykdtree/kdtree.c', 'pykdtree/_kdtree_core.c'])],
cmdclass = {'build_ext': build_ext_subclass },
classifiers=[
@@ -95,6 +95,6 @@ setup(
'Topic :: Scientific/Engineering'
]
)
-
+
View it on GitLab: https://salsa.debian.org/debian-gis-team/pykdtree/compare/6cbc4c35cc7a16ace87d64d2c605986a4beeee71...2bad0fea1659bfc197ce1e62d3c6a0ddc6c1e8db
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pykdtree/compare/6cbc4c35cc7a16ace87d64d2c605986a4beeee71...2bad0fea1659bfc197ce1e62d3c6a0ddc6c1e8db
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20180604/d6fcd8a0/attachment-0001.html>
More information about the Pkg-grass-devel
mailing list