Approximate all nearest neighbor search for high dimensional entropy estimation for image registration

Signal Processing 92 (2012) 1302–1316

Contents lists available at SciVerse ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Approximate all nearest neighbor search for high dimensional entropy estimation for image registration Jan Kybic n, Ivan Vnucˇko Center for Machine Perception, Department of Cybernetics, Faculty of Electrical Engineering, Czech Technical University in Prague, Czech Republic

a r t i c l e in f o

abstract

Article history: Received 4 August 2010 Received in revised form 15 September 2011 Accepted 25 November 2011 Available online 7 December 2011

Information theoretic criteria such as mutual information are often used as similarity measures for inter-modality image registration. For better performance, it is useful to consider vector-valued pixel features. However, this leads to the task of estimating entropy in medium to high dimensional spaces, for which standard histogram entropy estimator is not usable. We have therefore previously proposed to use a nearest neighbor-based Kozachenko–Leonenko (KL) entropy estimator. Here we address the issue of determining a suitable all nearest neighbor (NN) search algorithm for this relatively specific task. We evaluate several well-known state-of-the-art standard algorithms based on k-d trees (FLANN), balanced box decomposition (BBD) trees (ANN), and locality sensitive hashing (LSH), using publicly available implementations. In addition, we present our own method, which is based on k-d trees with several enhancements and is tailored for this particular application. We conclude that all tree-based methods perform acceptably well, with our method being the fastest and most suitable for the all-NN search task needed by the KL estimator on image data, while the ANN and especially FLANN methods being most often the fastest on other types of data. On the other hand, LSH is found the least suitable, with the brute force search being the slowest. & 2011 Elsevier B.V. All rights reserved.

Keywords: Nearest neighbor search Post-office problem Closest-point queries Approximation algorithms Entropy estimation k-d tree Incremental update

1. Introduction Nearest neighbor (NN) search problem [1] consist of preprocessing a set of points S such that for a given query point, a nearest neighbor from S can be found quickly. This problem has a large number of applications, such as knowledge discovery in databases [2], case-based reasoning systems [3], computer aided medical diagnosis [4], protein or sequence search in molecular biology [5], scattered data interpolation [6], building a NN classifier

n

Corresponding author. E-mail address: [email protected] (J. Kybic).

0165-1684/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2011.11.027

(or a k-NN classifier) in pattern recognition [7,8], or a vector quantization encoder in signal processing [9,10]. It is also called a post-office problem [11]. Some of the many existing algorithms will be mentioned in Section 1.2. The task is most often formulated in a vector space using standard p-norms (‘1 , ‘2 or ‘1 ), although more general metrics such as the Levenstein distance for strings is used for some applications. All nearest neighbor (all-NN) search is a less studied problem. It consists of finding a nearest neighbor from the set S for all query points from S [12,13]. The task can be solved by repeating a standard NN search n times, where n is the number of points in S, but there are computational savings to be made given that the query set is identical to the set of points S and both are known in advance. All-NN

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

search is useful in estimating a local probability density from a set of samples, based on the observation that small NN distances correspond to points with high probability density and vice versa [14]. Applications include manifold learning [15], estimation of local dimension [16], estimation of Re´nyi entropy [17], Shannon entropy [18–20], Kullback– Leibler divergence [21], and mutual information [22]. Our target application is estimating information theory based similarity criteria for image registration. Image registration or image alignment is one of the fundamental image analysis tasks [23], especially in the biomedical domain. It is used for movement compensation, motion detection and measurement, image mosaicking, image fusion, atlas building and atlas based segmentation. Many image registration methods are based on optimizing an image similarity criterion and information theoretic criteria such as mutual information (MI) are frequently used, especially in multimodal image registration [24–27], for their ability to evaluate the statistical dependency between images. Mutual information of two random variables X, Y [28] is defined as IðX; YÞ ¼ HðXÞ þ HðYÞ HðX,YÞ, where HðXÞ and HðYÞ are marginal entropies and HðX,YÞ is the joint entropy. The joint entropy itself and normalized versions of mutual information are other often used similarity criteria. Normally, scalar pixel intensities are used as features, i.e. samples of the random variables X, Y, requiring entropy estimation in dimension d ¼2. Only a few attempts to use higher dimensionality vector features have been made, using features such as gradients, color, output of spatial filters, texture descriptors, or intensities of neighborhood pixels [29–32]. The main obstacle is that the predominantly used standard histogram-based estimator [33,34] or the kernel density estimator [35] perform poorly in higher dimensions. For this reason, the above mentioned approaches are limited to low dimensions (e.g. d¼ 4 in [29]) or need to use the normal approximation (d ¼18 in [31]). The little known Kozachenko–Leonenko (KL) estimator [19,20], based on NN distances, is a qualitative improvement [36,37], permitting evaluation of the entropy for truly high-dimensional features such as color (d¼6), texture descriptors, or large pixel neighborhoods (d ¼ 10100). The computational bottleneck of the KL estimator is the all-NN search, since the number of pixels n in an image and the dimension d are large, and the estimator needs to be evaluated many times during the registration in order to find the optimal transformation. On the other hand, we can take advantage of the fact that (i) approximate solution is acceptable; (ii) the query set is known and identical to the set of data points; and (iii) as the optimizer converges, the algorithm is called repeatedly with very similar data. The aim of this paper is to find a suitable algorithm for the all-NN search on image generated features. It is organized as follows: after defining the problem in Section 1.1 and reviewing some of the existing algorithms in Section 1.2, we describe our target application in Section 1.3. We propose our own all-NN search algorithm in Section 2, which is based on k-d trees but is tailored to our particular problem, analyze its performance theoretically in Section 3, and compare it with available existing

1303

methods in Section 4. Additional details about the problem formulation and the implementation, as well as additional experiments can be found in our research report [38]. 1.1. Problem definition and notation Given a set of n points S from Rd and a query point q 2 Rd , then p 2 S is a nearest neighbor (NN) of q iff

Rq ¼ JpqJ ¼ minfJp0 qJ; p0 2 Sg

ð1Þ

and Rq is the NN distance. Since exact NN search is often too time consuming for many application, especially for high dimensional spaces (high d) and many points (high n), an approximate NN problem is considered [39]. Here we judge the approximation quality mainly according to the entropy estimation error (see Section 1.3) but we consider also the relative distance error and the probability of finding the correct NN (Section 3). All NN problem (all-NN) consists of finding, for all query points q from S, a NN from S~ ¼ S\fqg [12,13]. Our generalization of the formulation is to consider S to be a multiset [40,41], i.e. it can contain some points several times. If the query point q has multiplicity one then we want the NN from the rest of the points, otherwise we return q itself, the distance Rq equal to zero, and we report the multiplicity. The dynamic all-NN problem consists of solving the allðiÞ ðiÞ NN problem for m sets of points Si ¼ fpðiÞ 1 ,p2 , . . . ,pn g, with i ¼ 1, . . . ,m, each of them containing n points, assuming that the displacement Jpjði þ 1Þ pðiÞ j J of each point j is small with respect to the inter-point distances JpðiÞ pðiÞ J. In j k other words, pðiÞ is seen as a trajectory of a point p j in time j i and it is assumed that the movement is slow. This allows the datastructure built by preprocessing the set Si to be used with few changes also for the set Si þ 1 , bringing computational savings. As far as we know, this scenario has not yet been described in the literature. For compatibility with previous work, here we work with the standard Euclidean ‘2 norm. However, most NN search algorithms including ours can adapt to any ‘p norm with minimal changes. For example, we use the ‘1 norm for entropy estimation in [38] for its better compatibility with the rectangular space partitioning of the k-d tree, which can lead to up to 50% faster computation (see Section 4.1). 1.2. Existing NN search algorithms Nearest neighbor search in low dimensions is a wellstudied problem. For d¼1, an asymptotically optimal solution consists of sorting the points and applying a binary search. For d ¼2, an asymptotically optimal solution involves a Voronoi diagram [42,39], however, it is slow in practice. A number of algorithms is based on griding or treating each dimension separately, thus trimming the set of candidate points and exhaustively searching the remaining ones [43–45]. Unfortunately, these strategies only seem to be suitable for low dimensions and distributions close to uniform.

1304

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

An often used technique for moderate d (2  100) and large n is hierarchical partitioning of the space or the data points [39,42,46], creating a tree structure. Quad-tree methods [47,48] divide the space to 2d equal hyperrectangles at each level; hyperspheres are also used [49,50]. A prominent example of a data-partitioning approach is the k-d tree [51–55], based on splitting the data points evenly between node children at each level by an axis parallel dividing plane. A balanced box-decomposition (BBD) tree guarantees that both the cardinality and size of the partitions decreases exponentially [39], avoiding degenerate splits. A well-separated pair decomposition [56] constructs a list of point decompositions, speeding up the k-nearest neighbor search. All-NN algorithms with complexity Oðn log nÞ based on quadtree-like decompositions are known [12,57]. For large d and n, solving the exact NN problem is impractical, as it is very likely going to be too slow for many intended applications. In fact, the complexity of most exact NN and all-NN methods is exponential (or worse) in d and there is no known method which is simultaneously approximately linear in storage space and logarithmic in query time [39,58,42,59]. Approximate NN search can guarantee bounds on the error of the reported NN distance [39]. The best-bin-first (BBF) strategy [53] with a limit on the number of visited nodes appears to be faster [60] but without the error guarantee. Muja and Lowe [60] use several k-d trees with randomly chosen splitting dimension, which improves the performance in higher dimensions. They also suggest to use a tree based on recursive k-means clustering. A spill tree [61] method traverses the tree without backtracking while making the boxes overlap, returning the correct NN point with high probability [61]. For very large number of dimensions (d 4 100), the most promising approximative approaches are based on random projections [62,63] or locality sensitive hashing [64–67], especially if the intrinsic dimensionality of the data is large and a large errors are acceptable [68], otherwise spatial tree based approaches are to be preferred [61]. Dynamic NN search algorithms address the case of solving the NN or all-NN problem repeatedly, with a changing point set S. As far as we know, all previously described methods handle only the inserting and deleting of points between queries [69–72] or the case of points moving along a predetermined trajectories [54]. This is different from our problem, where the points move slowly but in an unpredictable way. There are also algorithms for general metric spaces without a vector space structure [73,74] and algorithms for very large datasets which do not fit into memory. These are outside the scope of this paper.

a random variable F in Rd with a probability density distribution f, the Shannon entropy of F Z f ðxÞlog f ðxÞ dx ð2Þ HðFÞ ¼  x2Rd

can be estimated as ! d X ^ HðFÞ ¼ log Rq þ C n q2S

ð3Þ

where n is the number of samples, d is the dimensionality of the space, and Rq is the NN distance for a point q in S~ ¼ S\fqg as defined by (1). The constant C is C ¼ log

ðn1Þpd=2 þ ge Gð1 þ d=2Þ

ð4Þ

if Rq is measured using the ‘2 norm and C ¼ logð2d ðn1ÞÞ þ ge

ð5Þ

for the ‘1 norm, with ge  0:577 being the Euler– Mascheroni constant. As this estimator does not use binning, it is more suitable for high dimensionality data and/or small sample sizes than the standard histogram approach, mainly because of its small bias [20,36,75]. The estimator (3) is asymptotically unbiased and consistent under very general conditions. Here we note only that in practice the estimator works as long as f is smooth and bounded (no Dirac impulses), for details see [19,76]. In image registration, the intensity values are often quantized and conflicts (i.e. repeated values leading to Rq ¼ 0) happen frequently. Some authors advocate adding a low-amplitude perturbation to the data [22]. We prefer to robustify the estimator by replacing the log Rq function [37] in (3) by ( for Rq Z e d log Rq 0 log Rq ¼ ð6Þ logðed =wS ðqÞÞ for Rq o e where wS ðqÞ is the multiplicity of point q. This can be viewed as switching to a kernel estimator for Rq o e. Instead of a hard threshold, a smooth transition can be used [75]. 2. Method Our method [37] uses a k-d tree structure [51] with the best-bin-first (BBF) search strategy [53] and four main enhancements which make it suitable for our application: (i) tight bounding boxes (Section 2.1); (ii) all-NN search instead of repeated NN search (Section 2.2); (iii) incremental update (Section 2.3); (iv) efficient handling of multiple points (Section 2.1). 2.1. Building extended k-d tree

1.3. Nearest neighbor entropy estimation Our target application is estimating entropy using a little known Kozachenko–Leonenko (KL) NN-based entropy estimator [19,20]. Given a set of n samples S of

A k-d tree [51,77,78] stores data points in its leaves. Each node represents a hyperinterval (an axis aligned hyperrectagle), which we call a loose bounding box (LBB) of the node. The LBBs at one level form a complete partitioning of the space Rd . Each leaf node contains a multiset of

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

points from S that belong to its LBB. The tree building Algorithm 1 is based on recursive splitting. We start with the entire input multiset S as the root node and the whole space Rd as its LBB. Then, recursively from the root, for each node Q with at least L points we choose the splitting dimension m as the longest edge of the TBB and the splitting value x as the median of the coordinate m values [52]. Algorithm 1. BuildTree, for building the extended k-d tree. Input: Multiset of points S from Rd . Output: k-d tree T for S. function BuildTree(S): create root node R with points S and LBBR ¼ Rd SplitNode(R) return tree T with root R function SplitNode(tree node Q with points P and LBBQ ): compute TBBQ as follows: Q

Q

li ¼ minfxi ; x 2 Pg, hi ¼ maxfxi ; x 2 Pg if number of points JPJ 4 L then 6 6 m’arg max1 r i r D ðxH xL Þ i i 6 6 6 x’Medianfxm ; x 2 Pg 6 6 PL ’fxm o x; x 2 Pg 6 6 6 PR ’fxm Z x; x 2 Pg 6 6 LBBL ’fxm o x; x 2 LBBQ g 6 6 6 LBBR ’fxm Z x; x 2 LBBQ g 6 6 if neither P L nor P R is empty then 66 66 6 6 create left child L of Q with points P L and LBBL 6 66 create right child R of Q with points P R and LBBR 66 66 66 SplitNodeðLÞ 44 SplitNodeðRÞ

Duplicate points are guaranteed to fall into the same leaf. If all points in a node are identical, the node is marked as degenerated and is not split further. This reduces the depth of the tree and accelerates the search in case of many identical points. As a novelty, we also maintain and store a tight bounding box (TBB) for each node. The TBB of a node Q is the smallest hyperinterval (in the sense of inclusion) containing all points in the subtree of Q. In Fig. 1a, TBBs are drawn in dashed lines. For the purpose of the dynamic tree update (Section 2.3) we introduce a parameter d 2 /0,0:5S, limiting each child subtree to contain at most ð12 þ dÞnp points, where np is the number of points of its parent subtree [54]. The parameter d controls the trade-off between the update speed and search speed (Section 2.3).

2.2. Exact all nearest neighbor search To find all NNs we loop through all leaves and through all points in each leaf. Each point q 2 S acts as a query ^ The method is points and we find its nearest neighbor q. described in Algorithm 2 (ignore the parameter V for the moment). The novelty is that instead of starting from the global root, the search starts from the leaf Q containing the query point q (see Fig. 2). This way we immediately

1305

have a good upper bound on the NN distance and often the correct NN itself. ^ wS Algorithm 2. allNNSearch finds the mapping q/ðq, ðqÞÞ for each point q in the tree T. If V ¼ 1, then q^ ¼ NN0S ðqÞ is the exact NN, otherwise q^ 2 Se is an approximate NN.

Input: Tree T with points S, maximum number of visited points V ^ and ^ the distance Rq ¼ JqqJ Output: Exact or approximate NN q, multiplicity wS ðqÞ for each query point q 2 Se . function allNNSearch(tree T): foreach leaf Q of the tree T do foreach point q from leaf Q do

$

bcall NNSearchðq,Q Þ and store results function NNSearch(querypoint q from leaf Q): pmin ’arg minx2Q \fqg JqxJ Rmin ’Jqpmin J if Rmin ¼ 0 then b return ðpmin , Rmin , wS ðqÞÞ v’VJQ J PQ ’empty priority queue of node pairs ðZ,Z 0 Þ sorted by ZZ PushIfBetter(Parent(Q), Q) while PQ not empty and v 40 do 6 6 ðZ,Z 0 Þ’pop the top element from PQ ==elementisremoved 6 6 6 if ZZ Z Rmin then 6 6 9 return ðpmin , R ,1Þ min 6 6 6 if Z is a leaf then 6 R ’min JqxJ 6 x2Z 6  6  v’vJZJ 6 6 6  if R o Rmin then 6$ 6 Rmin ’R 6 6 6 pmin ’arg minx2Z JqxJ 6 6 6 else 66 6 6 ==Examine yet unseen nodes accessible from Z 66 66 6 for X in fParent,LeftChild,RightChildgðZÞ\fZ 0 g do 44 b PushIfBetterðX,ZÞ return ðpmin , Rmin ,1Þ function PushIfBetternew node X, previous node Z: if X ¼ ParentðZÞ then 9 ZX ¼ dðq, RD \LBBZ Þ else b ZX ¼ dðq,TBBX Þ if ZX o Rmin then b pushðX,Z, ZX Þ to PQ

While traversing the tree, the points in so far unseen leaf nodes are searched and the currently best NN candidate q^ is updated. The search order is determined according to the best-bin-first strategy [53] using a lower bound ZX of the distance from the query q to yet unexplored points reachable from a node X. The lower bound ZX is calculated as follows: If X is not an ancestor of Q (in the sense of the original edge orientations), then q lies outside of X and ZX is the distance of q to its TBB, ZX ¼ dðq,TBBX Þ. Otherwise, if X is an ancestor of Q, then q lies inside X and we set ZX as the distance of q to the yet unseen part of the X subtree, ZX ¼ dðq, Rd \LBBZ Þ, where Z is an already explored child of X and is either equal to Q or its ancestor. Both cases are illustrated in Fig. 2a as ZH and ZA , respectively. A priority queue of nodes to be explored is maintained, sorted by Z. Nodes with Z higher than the currently best

1306

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

Fig. 1. A k-d tree example in a 2D space. (a) Solid blue lines marked by uppercase letters represent the splitting hyperplanes hierarchically subdividing the space into loose bounding boxes (LBB) corresponding to tree nodes. Dashed lines show the tight bounding boxes (TBB). (b) The k-d tree itself, with round non-leaf nodes marked by the corresponding dividing hyperplanes and rectangular leaf nodes each containing a set of data points denoted by numbers. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

NN distance are pruned, since they may not lead to an improvement. If such a node is found at the top of the queue, the search is terminated. 2.2.1. Approximate search The algorithm described in the previous section guarantees to find the exact NN but can be too slow for some applications, despite the pruning. This is not surprising since so far no algorithm is known achieving simultaneously logarithmic exact NN search time in roughly linear memory [39]. We will therefore consider approximate search. A simple change of the presented algorithm leads to an approximative, but potentially much faster version. It is controlled by a parameter V that bounds the

number of visited points [53]. Then the search is stopped and the best result so far is reported. Algorithm 2 already incorporates this change.

2.3. Dynamic tree update If the changes between point positions in subsequent invocations of the all-NN search algorithm are sufficiently small, updating the tree is faster than rebuilding it from scratch. The Algorithm 3 takes a multiset S and a corresponding tree T and updates the tree to correspond to a new set of points S0 ; we assume a one-to-one mapping between the points in S and S0 .

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

1307

Fig. 2. The same k-d tree as in Fig. 1 with an example of a NN search. In the top image (a) we have marked the query point q (number 4) by a red dot. We have also marked the distance ZH from q to TBBH (in green). The distance from q to the LBBC (in cyan) is denoted ZA , as it is a distance to a yet unexplored part of the tree accessible through A (see the text for more details). On the bottom graph (b) we can observe how a query for the NN of point q (point number 4) proceeds, starting from the leaf I, visiting nodes F, B, A, C, D, G, and K, in that order (marked by red arrows), always going to the node with the currently smallest distance bound Z, until the nearest neighbor, point number 3 in leaf K is finally found. All other paths are pruned. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Algorithm 3. UpdateTree, updating the tree after the points have been moved. 0

Input: Tree T with root R and modified points S Output: Updated tree T. function UpdateTree(tree root R): UpdateNode(R) RebuildNode(R) function UpdateNode(tree node Q): if Q is not a leaf then UpdateNodeðLeftChildðQ ÞÞ

$

UpdateNodeðRightChildðQ ÞÞ move points from Q not inside LBBQ to ParentðQ Þ function RebuildNode(tree node Q): if X is a leaf then 9 compute and store TBBX

else 6 6 foreach point x in Q do 66 6 6 ==m and x are the splitting dimension and 66 66 66 value for Q 66 6 6 if xm o x then 66 66 6 6 9 move x to LeftChildðQ Þ 66 6 6 else 64 6 6 b move x to RightChildðQ Þ 6 6 6 ==J  J denotesthenumberofpointsinasubtree 6 6 if maxðJLeftChildðQ ÞJ,JRightChildðQ ÞJÞ 4ð1 þ dÞJQ J then 6 2 6 6 9 X’SplitNodeðQ Þ ==Algorithm1 6 6 else 6 66 RebuildNodeðLeftChildðQ ÞÞ 66 66 66 RebuildNodeðRightChildðQ ÞÞ 44 TBBQ ’combine TBBs of children of X

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

1308

We first detect points from the tree T which have moved out of the LBBs of their original leaves and attribute them to the appropriate new leaves. The TBBs of affected nodes are updated. Finally, parts of the tree violating the balance condition (Section 2.1) are rebuilt.

In the best case for tree update, the points are checked but no other work is done, with a complexity O(dn). In the worst case, the whole tree needs to be rebuilt and the complexity is the complexity of the build operation, Oðdn logðn=LÞÞ.

3. Theoretical performance analysis

3.1. All nearest neighbors search performance

The asymptotically most time consuming step in building the tree (Section 2.1) is calculating the TBBs which takes O(Ld) time per node, where L is the number of points in a node; median calculation takes O(L) time per node. The total asymptotic time complexity of building the tree is therefore Oðdn log ðn=LÞÞ, where log2 ðn=LÞ is the depth of the tree. However, note that for our application the tree building time is usually negligible compared to the all-NN search time. It is easy to establish lower and upper bounds on the complexity of a single exact NN search, given a query point and a cell Q it belongs to. In the worst case, the point needs to be compared with every other point with complexity O(nd), the same as for the brute force algorithm. In the best case, the NN is always found inside the same leaf as the query, so only OðLdÞ  OðdÞ operations are needed, since L is a small constant and L5 n. For the approximate search (Section 2.2.1), the worst-case complexity is reduced to O(Vd), since at most V points are searched. Using TBBs instead of LBBs will not change any of those bounds.

Iterating over all points for the all-NN search (Section 2.2) takes time O(n), because each tree node is visited once. This leads to the total best-case/worst-case complexity of the exact all-NN search to be O(nLd) and Oðn2 dÞ, respectively. In the approximative case, the complexity is O(nVd). This is better than the standard approach by repeated NN search, where iterating over all points takes time Oðn log n þ LdÞ  Oðn log nÞ, leading to a total complexity Oðn log n þ nLdÞ in the best case for the exact search and Oðn log n þnVdÞ for the approximate search (Section 2.2.1). The worst-case bound is unchanged. From our experience, about 50% speedup can be expected (see Table 3). 3.2. Basic average case analysis Average case asymptotic complexity of a particular simple k-d tree-based exact NN search algorithm was analyzed [79,59] and was found to be Oðnr Þ with r  0:06. The multiplicative factors and the possibly exponential dependence on the dimension d were not analyzed.

Table 1 Comparison of the tree build, search, and total times for an all-NN search for the BBF and ANN methods depending on the point distribution, for n ¼105. The shortest total time for each configuration (dimension d and a distribution type) is typeset in bold. Type

Normal Normal Normal Normal Uniform Uniform Uniform Uniform Corrnormal Corrnormal Corrnormal Corrnormal Corruniform Corruniform Corruniform Corruniform Clustered Clustered Clustered Clustered Subspace Subspace Subspace Subspace Imgneighb Imgneighb Imgneighb Imgneighb

d

2 5 10 20 2 5 10 20 2 5 10 20 2 5 10 20 2 5 10 20 2 5 10 20 2 5 10 20

BBF

ANN

FLANN

Bld

Search

Tot

Bld

Search

Tot

Bld

Search

Tot

0.1 0.1 0.2 0.3 0.1 0.1 0.2 0.3 0.1 0.1 0.2 0.3 0.1 0.1 0.2 0.3 0.1 0.1 0.2 0.3 0.1 0.1 0.2 0.3 0.1 0.1 0.2 0.2

0.1 2.1 116.6 2979.8 0.1 1.5 41.8 2994.4 0.1 2.0 77.1 1115.9 0.1 2.0 75.6 1177.6 0.1 1.7 29.8 290.5 0.1 2.0 71.9 1406.7 0.1 0.7 3.8 19.8

0.2 2.2 116.8 2980.1 0.2 1.7 42.0 2994.7 0.2 2.0 77.3 1116.1 0.2 2.1 75.8 1177.9 0.2 1.9 30.0 290.9 0.2 2.1 72.1 1407.0 0.2 0.8 4.0 20.0

0.1 0.2 0.4 0.6 0.1 0.2 0.4 0.8 0.1 0.2 0.4 0.5 0.1 0.2 0.5 0.5 0.1 0.2 0.6 1.0 0.1 0.3 0.4 0.8 0.1 0.3 0.7 1.7

0.3 2.5 88.6 1597.3 0.3 1.9 41.1 1280.5 0.3 2.9 49.5 402.2 0.5 2.3 49.4 464.3 0.4 2.5 39.9 214.4 0.4 3.7 76.6 954.8 0.2 1.2 11.4 71.3

0.4 2.7 89.0 1597.9 0.4 2.1 41.6 1281.2 0.4 3.1 49.9 402.7 0.5 2.5 49.8 464.8 0.5 2.7 40.4 215.4 0.5 4.0 77.0 955.6 0.3 1.5 12.2 73.0

0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.2

0.5 3.7 76.8 2783.0 0.6 3.3 38.9 2350.9 0.5 3.6 43.4 670.1 0.6 3.7 46.2 700.1 0.6 3.8 36.6 444.8 0.6 4.1 60.5 1414.0 0.5 1.9 13.3 94.1

0.6 3.8 77.0 2783.1 0.6 3.3 39.0 2351.1 0.6 3.7 43.6 670.3 0.6 3.8 46.4 700.3 0.7 3.9 36.7 445.0 0.7 4.2 60.6 1414.2 0.5 2.0 13.4 94.3

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

1309

Table 2 We show the build (Tb) and all-NN search (Ts) times for the perturbed dataset as well as the time (Tu) needed to update a tree created for the first dataset and the search time (T s0 ) using the updated tree. All times are in seconds. We also show the sums T b þ s ¼ T b þ T s and T u þ s0 ¼ T u þ T s0 . The best of the two times for the pairs T b þ s ,T u þ s0 and T b ,T u is typeset in bold. The parameter d controls the allowed unbalance of the tree.

s

d

Tb

Ts

Tu

T s0

Tb þ s

T u þ s0

0.001 0.001 0.001 0.001 0.001 0.010 0.010 0.010 0.010 0.010 0.100 0.100 0.100 0.100 0.100

0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4

0.14 0.14 0.14 0.14 0.15 0.14 0.13 0.14 0.14 0.14 0.14 0.14 0.13 0.14 0.14

1.48 1.50 1.49 1.50 1.50 1.49 1.48 1.50 1.47 1.48 1.64 1.63 1.62 1.62 1.64

0.16 0.05 0.05 0.04 0.05 0.19 0.05 0.06 0.05 0.05 0.24 0.09 0.09 0.09 0.08

1.57 1.49 1.49 1.50 1.49 1.59 1.51 1.49 1.47 1.50 1.80 1.66 1.66 1.66 1.64

1.62 1.64 1.63 1.63 1.64 1.63 1.61 1.64 1.61 1.62 1.78 1.78 1.75 1.75 1.78

1.73 1.55 1.54 1.54 1.53 1.79 1.56 1.55 1.52 1.55 2.04 1.76 1.74 1.75 1.72

Table 3 Means and standard deviations of total times (build and all-NN search) in seconds for evaluating joint entropy for the example in Fig. 5 for neighborhood features with neighborhood size h for different all-NN search methods. ‘bbf-tbb’ stands for BBF without TBBs, ‘bbf-ann’ stands for BBF with repeated NN search. Best times (at significance level 5%) are shown in bold.

bbf bbf-tbb bbf-ann flann ann lsh brute

bbf bbf-tbb bbf-ann flann ann lsh brute

h ¼ 1, d ¼2

h ¼ 2, d ¼8

h ¼ 3, d ¼18

0:08 7 0:003 0:12 7 0:015 0:26 7 0:019 0:17 7 0:009 0:08 7 0:003 83:67 7 1:316 57:20 7 0:015

0:35 7 0:025 0:42 7 0:028 0:61 7 0:030 1:11 7 0:072 0:38 7 0:037 18:63 7 2:069 119:09 7 0:491

0:64 7 0:032 0:61 7 0:037 1:03 7 0:047 1:85 7 0:074 2:02 7 0:138 32:37 7 1:268 246:89 7 4:143

h ¼ 4, d ¼32

h ¼ 5, d ¼50

h ¼ 6, d ¼72

1:82 7 0:089 2:08 7 0:113 3:017 0:083 2:41 7 0:066 4:53 7 0:362 115:34 7 6:019 476:26 7 22:841

2:45 7 0:055 2:69 7 0:130 3:76 7 0:100 2:77 7 0:057 9:86 7 1:343 296:62 7 16:898 639:20 7 24:961

3:03 7 0:111 4:03 7 0:213 5:26 7 0:182 3:32 7 0:079 18:69 7 2:335 603:607 263:985 894:30 7 87:120

Furthermore, due to the nature of the equations which do not seem to admit closed-form solutions [80,79] we do not expect the analysis to be extended to the approximate NN search algorithm studied here. Indeed, the BBF approach [53] we are using has been proposed 14 years ago and yet according to our best knowledge it has not been theoretically analyzed. We have therefore chosen to make a number of simplifying assumptions, following and extending the idealized but tractable analysis of Friedman et al. [52], which we recall here briefly for completeness. Let as consider the smallest ball B centered around one of the data points q 2 S which contains exactly k other points. It can be shown [52] that its probability contents R uk ¼ B f ðxÞ dx, where f ðxÞ is the probability density, is a random variable independent of f and following a beta distribution with expected value k=n, where n is the total number of points including the query. Assuming that the density f is locally constant, from uk  vk f we see that the 1 expected volume vk of B is kf =n.

By the same argument, the expected volume of the 1 LBBs of the tree leaves is Lf =n, where L is the number of points per leaf. Since the longest edges tend to be split when constructing the trees, we conclude that the LBBs are approximately hypercubical with edge length 1=d

eL ¼ vL ¼ ðLf =nÞ1=d . We are now ready to estimate the number of leaf nodes that need to be considered in the exact NN search with query point q. It is bound from above by the number of leaves intersecting with B for k¼1, which is in turn bound by the expected number nL of leaves intersecting with a hyperrectangle containing B, calculated as nL ¼

1

 d eB þ1 ¼ eL



GðdÞ L

1=d

!d þ1

¼

2Gðd=2 þ 1Þ1=d þ1 pffiffiffiffi 1=d pL

!d

ð7Þ 1=d

with k ¼1 is the hyperrectangle where eB ¼ ðvk GðdÞÞ edge length and GðdÞ ¼ 2d Gðd=2 þ1Þpd=2 is the ratio of a

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

1310

hyperrectangle and an enclosed hypersphere volumes. Round-off errors are neglected. The number of points to examine is therefore nc ¼ minfnL L,ng1 (not counting the query point). We set L¼1 since it minimizes nc, but note that in practice larger values are usually optimal due to overhead associated with the tree traversal. Then the number of points to examine is nc ¼ minfgd ,ng1

with

g ¼ 2p1=2 Gðd=2þ 1Þ1=d þ1  1 þ

rffiffiffiffiffiffi pffiffiffi 2d  1 þ0:4839 d pe

ð8Þ

where the approximation comes from the Stirling’s formula. The constant g can be interpreted as the number of leaf bounding boxes that the search region intersects in each dimension; g ¼ 2 for d ¼1 and g  6 for d¼100. Then Oðnc Þ and Oðnnc Þ are the total complexities of a single exact NN search and an all-NN exact search, respectively. For the approximate search, the complexities are reduced to OðminfV,nc gÞ and OðnminfV,nc gÞ, respectively, since the search is aborted after having examined V points. Several remarks can be made: first, recall that this is only an upper bound on the expected complexity, we have for example not taken into account that the leaves are examined in the best-bin-first order and the pruning it allows. From practical experiments, it seems that nc indeed behaves like gd but with rather smaller g  1  1:5, depending on the distribution (see e.g. Fig. 4 and Table 3). Second, there is an even simpler and more intuitive argument leading to a formula with similar asymptotic behavior (in d): assuming that a leaf LBB is a hypercube with edge length eL, a 1=2 maximum ‘2 distance inside is eL d , hence we need to examine neighboring leaves intersecting a hypersphere with 1=2 radius eL d ; the volume of the hypersphere and therefore d=2 the number of such leaves is proportional to d , so d=2 nc  d (compare with (8)). Third, the bounds remain valid even if the points are located on a manifold with a local 0 0 dimension d od, with d replacing d; while the edge size will vary depending on the orientation of the tangent space with respect to coordinate axes. Fourth, complexity is reduced by using the ‘1 norm instead of the l2 norm. In this case the hypercube diameter is equal to its edge length eL and in a regular grid such a hypercube can intersect with at most 2d leaves, so nc  2d , i.e. g ¼ 2. Using tight bounding boxes (TBBs) is beneficial if the number of points in a leaf (governed by L) is small or if the point density is inhomogeneous, making the TBBs significantly smaller than the enclosing LBBs. This way the number of nodes to be searched is reduced, since it is less likely for the ball B to intersect the smaller TBB than the larger enclosing LBB. Based on the analysis above, we can estimate the speed improvement as follows: let us assume that the ratio between the size of TBBs and LBBs is b 2 ½0; 1 in each dimension. A search region will on the average overlap g LBBs in each dimension (8). The search region will also overlap the TBBs associated with these LBBs, except the extremal ones, which we assume will be only overlapped with probability b. Hence g~ ¼ g þ b1 TBBs will be overlapped in each dimension, with the total number of points being searched nc ¼ minfg~ d ,ng1. We therefore expect the search time to be reduced by a

constant multiplicative factor, which is consistent with our experimental results (see Table 3). 3.3. Approximation error analysis Based on the derived estimate of nc, the approximation error caused by limiting the number of examined points to V will be evaluated in terms of two quantities: the fraction of correctly identified nearest neighbors x and the relative error e. The fraction of correct NNs, i.e. the probability of finding the correct NN, can be estimated as follows: in order to be sure to find the NN, we should examine a set U of nc points but instead at most V points from a set V D U are examined. Assuming V r nc , the probability of finding the correct NN is

x¼

V nc

ð9Þ

Note that this is a conservative estimate, assuming that the points are examined in a random order, instead of taking the most promising leaves first, as it is done by the BBF strategy. An often used performance metric for NN search algorithms is the mean relative error over all points

e ¼ Eq2S ½eq  with eq ¼

R^ q 1 Rq

ð10Þ

where R^ q is the estimated and Rq the true NN distance for point q. We first express the probability that a ball or radius r does not contain any of N random sample points Q ðr; NÞ ¼ ð1ar d ÞN 1r r r max U r max ¼ a1=d ,

with a ¼ fcd ,

and

cd ¼

pd=2 Gðd=2 þ 1Þ

ð11Þ

where 1  U is the Iverson bracket, f is the locally constant density and cd is the unit hypersphere volume. Note that for large N we have Q  expðar d NÞ, which is another formula used in the literature [20], which would however lead to more difficult integration later on. The corresponding probability density of a NN being found at distance r is qðr; NÞ ¼ Nð1ar d ÞN1 a dr

d1

1r r r max U

ð12Þ

Let us define the NN distance in the complete set U as

rU ¼ minfJpqJ; p 2 Ug and similarly for V and its complement W ¼ U\V. This decomposition has the advantage that RV and RW are independent. Since RU ¼ minfRV , RW g, we can express the expected relative error as Z rmax e ¼ EW ½EV ½e9RW  ¼ qðRW ; nc VÞEV ½e9RW  dRW ð13Þ 0

with EV ½e9RW  ¼

Z

rmax 0

RV minfRV , RW g

qðRV ; VÞ dRV 1

ð14Þ

Substituting (12) into (14) and (13) allows this double integral to be calculated numerically. To obtain analytical formulas, we decompose (14) Z RW Z rmax RV EV ½e9RW  ¼ qðRV ; VÞ dRV þ qðRV ; VÞ dRV 1 0

RW

RW

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

integrate the first part and substitute (12) Z V ad rmax ð1aRdV ÞV1 RdV dRV ¼ Q ðRV ; VÞ þ

RW

RW

Further substitution of s ¼ aRdV leads to Z V 1=d 1 ¼ Q ðRV ; VÞ þ a ð1sÞV1 s1=d ds

RW

ð15Þ

aRdW

where we can recognize the complete and incomplete beta function B      V 1=d 1 1 þ1,V B aRdW ; þ 1,V ¼ Q ðRV ; VÞ þ a B d d RW ð16Þ This way we can calculate the inner integral (14) analytically. The outer integral (13) still has to be calculated numerically. This can be avoided by a polynomial approximation M X

s1=d  gðsÞ ¼

am ð1sÞm

ð17Þ

m¼0

for some moderate M (e.g. 3–10) where the coefficients am are chosen to minimize the least squares error M 1 X

1=d 2

ðgðsi Þsi

Þ

for si ¼

i¼0

i M1

ð18Þ

which can be done by solving a small linear system of equations. It allows us to approximate the incomplete beta function as   M X 1 am ð1ð1rÞm þ V Þ B r; þ1,V  ð19Þ d m þV m¼0 Substituting into (15) and (13) yields, after lengthy calculations, an analytical formula for the mean relative error       V 1 1 e ¼ þ ðnc VÞV B 1 ,nc V B 1þ ,V I 1 nc d d with I ¼

M X

am m þV m¼0

     1 1 B 1 ,nc V B 1 ,nc þ m d d ð20Þ

The behavior of both error measures e and x can be observed in Fig. 3 where the predictions can be compared with measured results on simulated uniformly distributed data with n ¼104. We can see that the global shape of the curves is correct even though the absolute numbers are not, our analysis is too pessimistic. This is probably because we were not able to take into account the order in which the k-d tree nodes are being searched during the BBF procedure. This appears to be an open research problem. 4. Experiments We have implemented our algorithm in Cþþ and performed several experiments on a standard PC with an Intel i7 processor. We are comparing our algorithm (denoted BBF) against the brute force algorithm (referred to as ‘brute’) and three state-of-the-art approximate NN

1311

search implementations: (i) the ANN library1 by Arya and Mount [39] (referred to as ANN), which uses a balanced box decomposition (BBD) tree; (ii) the FLANN library2 by Muja and Lowe [60], which uses randomized k-d or hierarchical k-means trees; and (iii) a locality sensitive hashing (LSH) as implemented in the Caltech Large Scale Image Search.3 All libraries were compiled according to provided instructions. For each method we have identified the main parameters and tried to set them to optimal values for given size, dimensionality, and class of data, by first exhaustively trying a number of parameter combinations from operator-suggested ranges and then optimizing locally, similarly to [60]. However, for larger datasets we could only afford to test a small number of parameter combinations. The ANN library was the easiest to handle in this aspect, as we only needed to specify the approximation factor e, which was very intuitive. Interestingly, sometimes we could set e as high as 10 and still obtain useful entropy estimates with much lower error. All other parameters were left to default or suggested automatically. In our BBF algorithm, we varied the leaf size L and the maximum number of visited points V. Increasing V monotonously decreased the error and increased computation time. For L, there was a shallow data-dependent optimum between L¼15 and L¼150. The FLANN library is capable of automatic parameter selection. It is too time consuming to be used at every run or on the larger datasets, but we have applied it to a number of representative cases to bootstrap our parameter tuning procedure. We have mostly used a single k-d tree with the ‘checks’ parameter determining the number of leaves to check and thus playing the role of V in BBF. Only at high dimensions and if high accuracy was required, we used 4–10 k-d trees. In the LSH method we have used the ‘2 hash function, ‘2 distance and bin width h ¼0.25, which proved optimal. We have then varied the number of hash functions and the number of tables. The number of hash functions influences the number of points per bin and therefore the trade-off between computation time and recall. Increasing the number of tables increases accuracy at the expense of computational time. Typically, we would use 50–500 hash functions and 1–20 tables. 4.1. Exact search The first set of experiments (Fig. 4) demonstrates the dependence of the total time to find all exact NNs on the number of points n for two different dimensions d for normally and isotropically distributed data. LSH was not included in this comparison because it cannot find the exact solution. For d¼5, we can see that the brute force method is quadratic as expected, while the complexity of the other methods is only slightly worse than linear. 1

http://www.cs.umd.edu/  mount/ANN/ http://www.cs.ubc.ca/ mariusm/index.php/FLANN/FLANN http://www.vision.caltech.edu/malaa/software/research/ image-search 2 3

1312

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

Fig. 3. Comparison of predicted (left) and measured (right) dependencies of the ratio of correct NNs found x and relative error e on the dimension d and the number of nodes searched V. We show the means of 10 repetitions for n ¼104 points uniformly distributed inside a unit hypercube. The horizontal scale of the right and left graphs is the same, except for the bottom left graph.

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

1313

Fig. 4. The total time in seconds to find all exact NNs for a normally and isotropically distributed points as a function of the number of points n for dimensions d¼5 and d¼ 15. Thick lines show the use of an ‘2 norm, thin dashed line the ‘1 norm. (Results for the highest d and n are not available for all methods.)

Unfortunately, for the higher dimension d ¼15 and the ‘2 norm, the performance of all methods is close to the brute force search. It demonstrates that in higher dimensions, it is indeed not practical to solve the all-NN problem exactly. The total time is dominated by the search time, except for the smallest n¼103, where the search and tree building times are comparable. Observe that using the ‘1 norm is several times faster compared to using the ‘2 norm for the BBF and ANN methods; FLANN does not support ‘1 norm when using kd trees. 4.2. Effect of different data distributions In Table 1, we show the all-NN build and search times of the BBF, ANN, and FLANN methods for different distributions, mostly taken from [39]: normal, uniform, correlated normal, correlated uniform, clustered, subspace, and image neighborhood (see [38] for more details). There are important differences between the performance of the evaluated algorithms on different types of data. Observe that our method performs especially well for the ‘image neighborhood’ distribution, for which it was designed. This distribution is generated from the 8-bit green channel of the 512  512 Lena image. For each data point to be generated, we randomly select a pixel location and the d features are the intensities on a discrete spiral centered at the chosen pixel. We have also added a small uniformly distributed perturbation (with amplitude 0.01) for the benefit of ANN which would otherwise fail because it cannot handle multiple points. 4.3. Incremental update To evaluate the effectiveness of the tree update operation (Section 2.3), we have generated n ¼106 uniformly distributed points with d ¼5 and perturbed them by adding another set of uniformly distributed points from

½s, sd . The we compare (a) the time needed to build the tree directly from the perturbed dataset and perform an all-NN search with (b) the time needed to update the tree from the unperturbed dataset to the perturbed dataset using Algorithm 3 and then doing the all-NN search. We can see (Table 2) that it is indeed advantageous to update the tree instead of rebuilding it anew, as long as we allow some unbalance d 4 0.

4.4. Sequential entropy estimation The final experiment attempts to simulate using the all-NN search methods for image similarity criteria evaluation. Features are pixel values in a square neighborhood of size h  h. Here we evaluate for different rotation angles the similarity of an image with the image of its smoothed gradient (Fig. 5a,b) but this technique is applicable to any local degradation or transform of an image. As a motivation example, observe (Fig. 5c,d) that the standard scalar joint entropy (h ¼1, d ¼2) is much less informative than the higher order neighborhood entropy (d) for h¼3, d ¼18. Table 3 compares the mean times needed for one joint entropy estimation by all methods, for different neighborhood sizes h and dimensions 2 d ¼ 2h . We have evaluated joint entropy for angles between  51 and 51 with step 11. The parameters of all fast methods were set so that the relative entropy estimation error is at most 1% and the time is as short as possible, using the brute force method as a reference. Besides the full BBF method, we also show the timings for method BBF-TBB, which is BBF without TBBs, and method BBF-ANN, which implements repeated NN search instead of the all-NN search. From their comparison we can see that both TBBs and the all-NN search strategy are helpful. On the other hand, there was no significant effect of using the tree update strategy (not shown). The BBF method was one of the best performers for all tested neighborhood sizes h.

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

1314

0

250

−2 50

200

100

150

150

100

50

−4 −6

100

−8 150

−10 −12

200

50

250

200

−14 −16

250 50

100

150

200

250

50

100

150

200

250

Fig. 5. We evaluate the similarity of image (a) with different rotations of its smoothed gradient image (b). We show the standard scalar joint entropy (c) for h ¼1 (d ¼2) and the higher order neighborhood entropy (b) for h ¼3 (d ¼18) as a function of the rotation angle between the two images.

On the other hand, the LSH method performed poorly, its use being an improvement only with respect to the brute force method. The problem is that to achieve an acceptable accuracy, many hash functions (hundreds) and many tables (4–10) are needed. 5. Conclusions We have compared four (BBF, ANN, FLANN, brute force) approximate all nearest neighbor search algorithms for use with a Kozachenko–Leonenko NN entropy estimator for applications in image registration. An all-NN search formulation was extended to multisets and a new dynamic all-NN search problem was defined. We have described an algorithm (BBF) to find exact and approximate NNs based on k-d trees, with several modifications including TBBs, dedicated all-NN strategy, incremental update and correct handling of multiplicities. We have confirmed that the brute force search as well as any exact NN search are indeed prohibitively slow for our

intended problem size. We can also rule our the evaluated LSH method, which does not work efficiently in low-error mode, when too many random projections are needed to achieve the desired accuracy. The latest version of the freely available FLANN and ANN libraries provide readily available tools for nearest neighbor searches with excellent performance. The FLANN library is slightly faster than ANN on the average, especially in higher dimensions. However, note that ANN and FLANN cannot be directly applied to our image data as they do not explicitly and efficiently handle the case of many identical datapoints, moreover FLANN does not handle ‘1 norm with k-d trees. Our method (BBF) could outperform ANN and FLANN on image data which it was designed to handle. All proposed enhancements of the k-d tree search are shown to be beneficial, even though the gain is sometimes small. In the future, it should be relatively easy to incorporate these enhancement into e.g. the FLANN algorithm, to combine their advantages.

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

The BBF algorithm development was motivated by evaluating entropy using a Kozachenko–Leonenko (KL) NN entropy estimator [19]. The KL estimator has very beneficial properties and we hope that a sufficiently fast all-NN search algorithm such as ours could trigger its more widespread use.

Acknowledgements This work was supported by the Czech Science Foundation project P202/11/0111. References [1] P. Indyk, Nearest Neighbors in High-Dimensional Spaces, 2nd ed., CRC Press, 2004. [2] U. Fayyad, G. Piatetsky-Shapiro, P. Smyth, From data mining to knowledge discovery in databases, AI Magazine 17 (1996) 37–54. [3] S. Wess, K.-D. Althoff, G. Derwand, Using k-d trees to improve the retrieval step in case-based reasoning, in: European Winter Conference on Brain Research, 1993, pp. 167–181. [4] F. Korn, N. Sidiropoulos, C. Faloutsos, E. Siegel, Z. Protopapas, Fast nearest neighbor search in medical image databases, The VLDB Journal (1996) 215–226. [5] M.S. Waterman, Introduction to Computational Biology, Chapman and Hall, 2000. [6] S.W. Park, L. Linsen, O. Kreylos, J.D. Owens, B. Hamann, Discrete Sibson interpolation, IEEE Transactions on Visualization and Computer Graphics 12 (2) (2006) 243–253. [7] R.O. Duda, P.E. Hart, D.G. Stork, Pattern Classification, 2nd ed., Wiley Interscience Publication, John Wiley, New York, 2001. [8] T.M. Cover, P.E. Hart, Nearest neighbor pattern classification, IEEE Transactions on Information Theory 13 (1967) 21–27. [9] A. Gersho, R.M. Gray, Vector Quantization and Signal Compression, Kluwer Academic Publishers, Norwell, MA, USA, 1991. [10] R.M. Gray, D.L. Neuhoff, Quantization, IEEE Transactions on Information Theory 44 (1993) 2325–2383. [11] D.E. Knuth, The Art of Computer Programming, Vol 3: Sorting and Searching, 2nd ed., Addison-Wesley, Reading, MA, 1998. [12] P.M. Vaidya, An Oðn log nÞ algorithm for the all-nearest-neighbors problem, Discrete & Computational Geometry 4 (2) (1989) 101–115. [13] K.L. Clarkson, Fast algorithms for the all nearest neighbors problem, in: Proceedings of 24th Annual IEEE Symposium on Foundations of Computer Science, 1983, pp. 226–232. [14] L.J. Buturovic, Improving k-nearest neighbor density and error estimates, Pattern Recognition 26 (4) (1993) 611–616. [15] J. Costa, A. Hero, Manifold learning using Euclidean k-nearest neighbor graphs, Proceedings of IEEE International Conference on Acoustic Speech and Signal Processing, vol. 3, 2004, pp. 988–991. [16] K.M. Carter, R. Raich, A.O. Hero, On local dimension estimation and its applications, IEEE Transactions on Signal Processing 58 (2) (2010) 650–663, doi:10.1109/TSP.2009.2031722. [17] A.O. Hero, B. Ma, O. Michel, J. Gorman, Applications of entropic spanning graphs, IEEE Signal Processing Magazine 19 (5) (2002) 85–95. ¨ [18] J. Beirlant, E.J. Dudewicz, L. Gyorfi, E.C. van der Meulen, Nonparametric entropy estimation: an overview, International Journal of Mathematics & Statistics (6) (1997) 17–39. [19] L.F. Kozachenko, N.N. Leonenko, On statistical estimation of entropy of random vector, Problems of Information Transmission 23 (9) (1987) in Russian. [20] J.D. Victor, Binless strategies for estimation of information from neural data, Physical Review E 66 (5) (2002) 051903 (15). [21] Q. Wang, S. Kulkarni, S. Verdu, A nearest-neighbor approach to estimating divergence between continuous random vectors, in: IEEE International Symposium on Information Theory, 2006, pp. 242–246. ¨ [22] A. Kraskov, H. Stogbauer, P. Grassberger, Estimating mutual information, Physical Review E 69 (2004) 066138. [23] B. Zitova´, J. Flusser, Image registration methods: a survey, Image and Vision Computing (21) (2003) 977–1000.

1315

[24] P. Viola, W.M.I. Wells, Alignment by maximization of mutual information, International Journal of Computer Vision (2) (1997) 137–154. [25] J. Pluim, J.B.A. Maintz, M.A. Viergever, Mutual-information-based registration of medical images: a survey, IEEE Transactions on Medical Imaging 22 (8) (2003) 986–1004. [26] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, P. Suetens, Multimodality image registration by maximization of mutual information, IEEE Transactions on Medical Imaging 16 (2) (1997) 187–198. [27] P. The´venaz, M. Unser, Optimization of mutual information for multiresolution image registration, IEEE Transactions on Image Processing 9 (12) (2000) 2083–2099. [28] D.J.C. MacKay, Information Theory, Inference, and Learning Algorithms, Cambridge University Press, 2003. [29] D. Rueckert, M.J. Clarkson, D.L.G. Hill, D.J. Hawkes, Non-rigid registration using higher-order mutual information, in: Proceedings of SPIE Medical Imaging 2000: Image Processing, 2000, pp. 438–447. [30] J.P.W. Pluim, J.B.A. Maintz, M.A. Viergever, Image registration by maximization of combined mutual information and gradient information, IEEE Transactions on Medical Imaging 19 (8) (2000). [31] D.B. Russakoff, C. Tomasi, T. Rohlfing, C.R.J. Maurer, Image similarity using mutual information of regions, in: T. Pajdla, J. Matas (Eds.), Proceedings of the Eighth European Conference on Computer Vision (ECCV), Lecture Notes in Computer Science, vol. 3023, Springer, 2004, pp. 596–607. [32] M.R. Sabuncu, P.J. Ramadge, Spatial information in entropy-based image registration, in: J.C. Gee, J.B.A. Maintz, M.W. Vannier (Eds.), WBIR, Lecture Notes in Computer Science, vol. 2717, Springer, 2003, pp. 132–141. [33] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, P. Suetens, Multi-modality image registration maximization of mutual information, in: MMBIA ’96: Proceedings of the 1996 Workshop on Mathematical Methods in Biomedical Image Analysis (MMBIA ’96), IEEE Computer Society, Washington, DC, USA, 1996, p. 14. [34] D.B. Russakoff, C. Tomasi, T. Rohlfing, C.R.M. Jr., Image similarity using mutual information of regions, in: ECCV, vol. 3, 2004, pp. 596–607. [35] P. Viola, W.M. Wells III, Alignment by maximization of mutual information, in: ICCV ’95: Proceedings of the Fifth International Conference on Computer Vision, IEEE Computer Society, Washington, DC, USA, 1995, p. 16. [36] J. Kybic, High-dimensional mutual information estimation for image registration, in: ICIP’04: Proceedings of the 2004 IEEE International Conference on Image Processing, IEEE Computer Society, 445 Hoes Lane, Piscataway, NJ, USA, 2004. [37] J. Kybic, Incremental updating of nearest neighbor-based highdimensional entropy estimation, in: P. Duhamel, L. Vandendorpe (Eds.), ICASSP2006, IEEE, Toulouse, France, 2006, p. III-804 (DVD proceedings). [38] J. Kybic, I. Vnucˇko, Approximate Best Bin First k-d Tree All Nearest Neighbor Search with Incremental Updates, Research Report K33340/10, CTU-CMP-2010-10, Department of Cybernetics, Faculty of Electrical Engineering Czech Technical University, Prague, Czech Republic, July 2010. [39] S. Arya, D.M. Mount, N.S. Netanyahu, R. Silverman, A.Y. Wu, An optimal algorithm for approximate nearest neighbor searching in fixed dimensions, Journal of the ACM 45 (6) (1998) 891–923. [40] W.D. Blizzard, Multiset theory, Notre Dame Journal of Formal Logic (1) (1988) 36–66. [41] D.E. Knuth, The Art of Computer programming. Vol. 2: Seminumerical Algorithms, Addison-Wesley, 1998. [42] P. Tsaparas, Nearest neighbor search in multidimensional spaces, Technical Report 31902, Department of Computer Science, University of Toronto, Canada, Qualifying Depth Oral Report, 1999. [43] J.H. Friedman, J.L. Baskett, L.J. Shustek, An algorithm for finding nearest neighbors, IEEE Transactions on Computing 24 (1975) 1000–1006. [44] T.P. Yunck, A technique to identify nearest neighbors, IEEE Transactions on Systems Man and Cybernetics 6 (10) (1976) 678–683. [45] S.A. Nene, S.K. Nayar, A simple algorithm for nearest neighbor search in high dimensions, IEEE Transactions on Pattern Analysis and Machine Intelligence 19 (9) (1997) 989–1003. [46] M. Smid, Closest-Point Problems in Computational Geometry, North Holland, 2000. [47] H. Samet, The quadtree and related hierarchical data structures, ACM Computing Surveys 16 (2) (1984) 187–260.

1316

J. Kybic, I. Vnucˇko / Signal Processing 92 (2012) 1302–1316

[48] D. Eppstein, M.T. Goodrich, J.Z. Sun, The skip quadtree: a simple dynamic data structure for multidimensional data, in: SCG ’05: Proceedings of the Twenty-First Annual Symposium on Computational Geometry, ACM Press, New York, NY, USA, 2005, pp. 296–305. [49] D.A. White, R. Jain, Similarity indexing with the SS-tree, in: ICDE ’96: Proceedings of the Twelfth International Conference on Data Engineering, IEEE Computer Society, Washington, DC, USA, 1996, pp. 516–523. [50] N. Katayama, S. Satoh, SR-tree: an index structure for nearest neighbor searching of high-dimensional point data, Transactions of the Institute of Electronics, Information and Communication Engineers J80-D-I (8) (1997) 703–717. [51] J.L. Bentley, Multidimensional binary search trees used for associative searching, Communications of the ACM 18 (9) (1975) 509–517. [52] J. Friedman, J.L. Bentley, R.A. Finkel, An algorithm for finding best matches in logarithmic expected time, ACM Transactions on Mathematical Software 3 (3) (1977) 209–226. [53] J.S. Beis, D.G. Lowe, Shape indexing using approximate nearestneighbour search in high-dimensional spaces, in: Proceedings of Conference on Computer Vision and Pattern Recognition, 1997, pp. 1000–1006. [54] P.K. Agarwal, J. Gao, L.J. Guibas, Kinetic medians and k-d trees, in: Proceedings of 10th European Symposium on Algorithms, 2002, pp. 5–16. [55] N. Sample, M. Haines, M. Arnold, T. Purcell, Optimizing search strategies in k-d trees, in: Fifth WSES/IEEE World Multiconference on Circuits, Systems, Communications & Computers (CSCC 2001), 2001. [56] P.B. Callahan, S.R. Kosaraju, A decomposition of multi-dimensional point-sets with applications to k-nearest-neighbors and n-body potential fields, in: Proceedings 24th Annual AMC Symposium on the Theory of Computing, 1992, pp. 546–556. [57] K.L. Clarkson, Nearest neighbor queries in metric spaces, in: STOC ’97: Proceedings of the Twenty-Ninth Annual ACM Symposium on Theory of Computing, ACM Press, New York, NY, USA, 1997, pp. 609–617. /http://dx.doi.org/http://doi.acm.org/10.1145/ 258533.258655S. [58] P. Indyk, Nearest neighbors in high-dimensional spaces, in: J.E. Goodman, J. O’Rourke (Eds.), Handbook of Discrete and Computational Geometry, 2nd ed., CRC Press, 2004 (Chapter 39). [59] P. Chanzy, Range Search and Nearest Neighbor Search in k-d Trees, Ph.D. Thesis, School of Computer Science, McGill University, 1993. [60] M. Muja, D.G. Lowe, Fast approximate nearest neighbors with automatic algorithm configuration, in: International Conference on Computer Vision Theory and Application VISSAPP’09, INSTICC Press, 2009, pp. 331–340. [61] T. Liu, A.W. Moore, A. Gray, K. Yang, An investigation of practical approximate nearest neighbor algorithms, in: L.K. Saul, Y. Weiss, L. Bottou (Eds.), Advances in Neural Information Processing Systems, vol. 17, MIT Press, Cambridge, MA, 2005, pp. 825–832. [62] P.N. Yianilos, Locally lifting the curse of dimensionality for nearest neighbor search (extended abstract), in: 11TH ACM-SIAM Symposium on Discrete Algorithms (SODA’00), 1999, pp. 361–370. [63] E. Kushilevitz, R. Ostrovsky, Y. Rabani, Efficient search for approximate nearest neighbor in high dimensional spaces, SIAM Journal on Computing 30 (2) (2000) 457–474.

[64] P. Indyk, R. Motwani, Approximate nearest neighbors: towards removing the curse of dimensionality, in: Proceedings of 30th STOC, 1998, pp. 604–613. [65] R. Panigrahy, Entropy based nearest neighbor search in high dimensions, in: SODA ’06: Proceedings of the Seventeenth Annual ACM-SIAM Symposium on Discrete Algorithm, ACM Press, New York, NY, USA, 2006, pp. 1186–1195. [66] A. Andoni, P. Indyk, Near-optimal hashing algorithms for approximate nearest neighbor in high dimensions, Communications of the ACM 51 (1) (2008) 117–122. [67] M. Aly, M. Munich, P. Perona, Indexing in large scale image collections: scaling properties and benchmark, in: IEEE Workshop on Applications of Computer Vision (WACV), 2011. [68] M. Datar, N. Immorlica, P. Indyk, V.S. Mirrokni, Locality-sensitive hashing scheme based on p-stable distributions, in: SCG ’04: Proceedings of the Twentieth Annual Symposium on Computational Geometry, ACM Press, New York, NY, USA, 2004, pp. 253–262. [69] P.K. Agarwal, D. Eppstein, J. Matouˇsek, Dynamic half-space reporting, geometric optimization, and minimum spanning trees, in: IEEE Symposium on Foundations of Computer Science, 1992, pp. 80–89. [70] Y. Chiang, R. Tamassia, Dynamic algorithms in computational geometry, Proceedings of the IEEE 80 (9) (1992) 362–381 (special Issue on Computational Geometry). [71] P. Bozanis, A. Nanopoulos, Y. Manopoulos, LR-tree: a logarithmic decomposable spatial index method, The Computer Journal 46 (3) (2003) 319–331. [72] T. Kanungo, D.M. Mount, N.S. Netanyahu, C.D. Piatko, R. Silverman, A.Y. Wu, Computing nearest neighbors for moving points and applications to clustering, in: Symposium on Discrete Algorithms, 1999, pp. 931–932. [73] J.K. Uhlmann, Satisfying general proximity/similarity queries with metric trees, Information Processing Letters 40 (1991) 175–179. [74] P.N. Yianilos, Data structures and algorithms for nearest neighbor search in general metric spaces, in: Proceedings of the Fifth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), 1993. [75] J. Kybic, High-dimensional entropy estimation for finite accuracy data: R-NN entropy estimator, in: N. Karssemeijer, B. Lelieveldt (Eds.), IPMI2007: Information Processing in Medical Imaging, 20th International Conference, Springer, Berlin, Heidelberg, Germany, 2007, pp. 569–580. [76] V. Mergel, On some properties of Kozachenko–Leonenko estimates and maximum entropy principle in goodness of fit test constructions, AIP Conference Proceedings on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, vol. 617, 2002, pp. 174–191, doi:10.1063/1.1477047. [77] F.P. Preparata, M.I. Shamos, Computational Geometry: An Introduction, Texts and Monographs in Computer Science, Springer-Verlag, 1985. [78] R. Sedgewick, Algorithms, Addison-Wesley, Reading, MA, 1988. [79] P. Chanzy, L. Devroye, C. Zamora-cura, Analysis of range search for random k-d trees, Acta Informatica 37 (2001) 355–383. [80] P. Flajolet, C. Puech, Partial match retrieval of multidimensional data, Journal of ACM 33 (1986) 371–407.