A Simple Algorithm for Nearest Neighbor Search in High Dimensions.
AbstractThe problem of finding the closest point in high-dimensional spaces is common in pattern recognition. Unfortunately, the complexity of most existing search algorithms, such as k-d tree and R-tree, grows exponentially with dimension, making them impractical for dimensionality above 15. In nearly all applications, the closest point is of interest only if it lies within a user-specified distance $\epsilon.$ We present a simple and practical algorithm to efficiently search for the nearest neighbor within Euclidean distance $\epsilon.$ The use of projection search combined with a novel data structure dramatically improves performance in high dimensions. A complexity analysis is presented which helps to automatically determine $\epsilon$ in structured problems. A comprehensive set of benchmarks clearly shows the superiority of the proposed algorithm for a variety of structured and unstructured search problems. Object recognition is demonstrated as an example application. The simplicity of the algorithm makes it possible to construct an inexpensive hardware search engine which can be 100 times faster than its software equivalent. A C++ implementation of our algorithm is available upon request to search@cs.columbia.edu/CAVE/.
Introduction
Searching for nearest neighbors continues to prove itself as an important problem in many
fields of science and engineering. The nearest neighbor problem in multiple dimensions is
stated as follows: given a set of n points and a novel query point Q in a d-dimensional
space, "Find a point in the set such that its distance from Q is lesser than, or equal to, the
distance of Q from any other point in the set" [ 21 ] . A variety of search algorithms have been
advanced since Knuth first stated this (post-office) problem. Why then, do we need a new
algorithm? The answer is that existing techniques perform very poorly in high dimensional
spaces. The complexity of most techniques grows exponentially with the dimensionality,
d. By high dimensional, we mean when, say d ? 25. Such high dimensionality occurs
commonly in applications that use eigenspace based appearance matching, such as real-time
object recognition [ tracking and inspection [ 26 ] , and feature detection
. Moreover, these techniques require that nearest neighbor search be performed using the
Euclidean distance (or L 2 ) norm. This can be a hard problem, especially when dimensionality
is high. High dimensionality is also observed in visual correspondence problems such as
motion estimation in MPEG coding (d ? estimation in binocular stereo
(d=25-81), and optical flow computation in structure from motion (also d=25-81).
In this paper, we propose a simple algorithm to efficiently search for the nearest neighbor
within distance ffl in high dimensions. We shall see that the complexity of the proposed
algorithm, for small ffl, grows very slowly with d. Our algorithm is successful because it
does not tackle the nearest neighbor problem as originally stated; it only finds points within
distance ffl from the novel point. This property is sufficient in most pattern recognition
problems (and for the problems stated above), because a "match" is declared with high
confidence only when a novel point is sufficiently close to a training point. Occasionally, it
is not possible to assume that ffl is known, so we suggest a method to automatically choose
ffl. We now briefly outline the proposed algorithm.
Our algorithm is based on the projection search paradigm first used by Friedman
Friedman's simple technique works as follows. In the preprocessing step, d dimensional
training points are ordered in d different ways by individually sorting each of their coordi-
nates. Each of the d sorted coordinate arrays can be thought of as a 1-D axis with the entire
d dimensional space collapsed (or projected) onto it. Given a novel point Q, the nearest
neighbor is found as follows. A small offset ffl is subtracted from and added to each of Q's
coordinates to obtain two values. Two binary searches are performed on each of the sorted
arrays to locate the positions of both the values. An axis with the minimum number of
points in between the positions is chosen. Finally, points in between the positions on the
chosen axis are exhaustively searched to obtain the closest point. The complexity of this
technique is roughly O(ndffl) and is clearly inefficient in high d.
This simple projection search was improved upon by Yunck utilizes a precomputed
data structure which maintains a mapping from the sorted to the unsorted (original)
coordinate arrays. In addition to this mapping, an indicator array of n elements is used.
Each element of the indicator array, henceforth called an indicator, corresponds to a point.
At the beginning of a search, all indicators are initialized to the number '1'. As before, a
small offset ffl is subtracted from and added to each of the novel point Q's coordinates to
obtain two values. Two binary searches are performed on each of the d sorted arrays to
locate the positions of both the values. The mapping from sorted to unsorted arrays is used
to find the points corresponding to the coordinates in between these values. Indicators corresponding
to these points are (binary) shifted to the left by one bit and the entire process
repeated for each of the d dimensions. At the end, points whose indicators have the value
must lie within an 2ffl hypercube. An exhaustive search can now be performed on the
hypercube points to find the nearest neighbor.
With the above data structure, Yunck was able to find points within the hypercube using
primarily integer operations. However, the total number of machine operations required (in-
teger and floating point) to find points within the hypercube are similar to that of Friedman's
algorithm (roughly O(ndffl)). Due to this and the fact that most modern CPUs do not significantly
penalize floating point operations, the improvement is only slight (benchmarked in a
later section). We propose a data structure that significantly reduces the total number of machine
operations required to locate points within the hypercube to roughly O(nffl
)).
Moreover, this data structure facilitates a very simple hardware implementation which can
result in a further increase in performance by two orders of magnitude.
Previous Work
Search algorithms can be divided into the following broad categories: (a) Exhaustive search,
(b) hashing and indexing, (c) static space partitioning, (d) dynamic space partitioning, and
randomized algorithms. The algorithm described in this paper falls in category (d). The
algorithms can be further categorized into those that work in vector spaces and those that
work in metric spaces. Categories (b)-(d) fall in the former, while category (a) falls in the
later. Metric space search techniques are used when it is possible to somehow compute a
distance measure between sample "points" or pieces of data but the space in which the points
reside lacks an explicit coordinate structure. In this paper, we focus only on vector space
techniques. For a detailed discussion on searching in metric spaces, refer to [
Exhaustive search, as the term implies, involves computing the distance of the novel
point from each and every point in the set and finding the point with the minimum distance.
This approach is clearly inefficient and its complexity is O(nd). Hashing and indexing are
the fastest search techniques and run in constant time. However, the space required to store
an index table increases exponentially with d. Hence, hybrid schemes of hashing from a
high dimensional space to a low (1 or 2) dimensional space and then indexing in this low
dimensional space have been proposed. Such a dimensionality reduction is called geometric
hashing . The problem is that, with increasing dimensionality, it becomes difficult
to construct a hash function that distributes data uniformly across the entire hash table
(index). An added drawback arises from the fact that hashing inherently partitions space
into bins. If two points in adjacent bins are closer to each other than a third point within
the same bin. A search algorithm that uses a hash table, or an index, will not correctly find
the point in the adjacent bin. Hence, hashing and indexing are only really effective when
the novel point is exactly equal to one of the database points.
Space partitioning techniques have led to a few elegant solutions to multi-dimensional
search problems. A method of particular theoretical significance divides the search space into
polygons. A Voronoi polygon is a geometrical construct obtained by intersecting
perpendicular bisectors of adjacent points. In a 2-D search space, Voronoi polygons allow
the nearest neighbor to be found in O(log 2
n) operations, where, n is the number of points in
the database. Unfortunately, the cost of constructing and storing Voronoi diagrams grows
exponentially with the number of dimensions. Details can be found in [ 3
Another algorithm of interest is the 1-D binary search generalized to d dimensions [ 11 ] . This
runs in O(log 2 n) time but requires storage O(n 4 ), which makes it impractical for n ? 100.
Perhaps the most widely used algorithm for searching in multiple dimensions is a static
space partitioning technique based on a k dimensional binary search tree, called the k-d tree
. The k-d tree is a data structure which partitions space using hyperplanes placed
perpendicular to the coordinate axes. The partitions are arranged hierarchically to form a
tree. In its simplest form, a k-d tree is constructed as follows. A point in the database
is chosen to be the root node. Points lying on one side of a hyperplane passing through
the root node are added to the left child and the points on the other side are added to
the right child. This process is applied recursively on the left and right children until a
small number of points remain. The resulting tree of hierarchically arranged hyperplanes
induces a partition of space into hyper-rectangular regions, termed buckets, each containing
a small number of points. The k-d tree can be used to search for the nearest neighbor as
follows. The k coordinates of a novel point are used to descend the tree to find the bucket
which contains it. An exhaustive search is performed to determine the closest point within
that bucket. The size of a "query" hypersphere is set to the distance of this closest point.
Information stored at the parent nodes is used to determine if this hypersphere intersects
with any other buckets. If it does, then that bucket is exhaustively searched and the size
of the hypersphere is revised if necessary. For fixed d, and under certain assumptions about
the underlying data, the k-d tree requires O(nlog 2 n) operations to construct and O(log 2 n)
operations to search [
k-d trees are extremely versatile and efficient to use in low dimensions. However, the
performance degrades exponentially 1 with increasing dimensionality. This is because, in
high dimensions, the query hypersphere tends to intersect many adjacent buckets, leading
to a dramatic increase in the number of points examined. k-d trees are dynamic data
structures which means that data can be added or deleted at a small cost. The impact
of adding or deleting data on the search performance is however quite unpredictable and
is related to the amount of imbalance the new data causes in the tree. High imbalance
1 Although this appears contradictory to the previous statement, the claim of O(log 2 n) complexity is
made assuming fixed d and varying n [ . The exact relationship between d and complexity has not
yet been established, but it has been observed by us and many others that it is roughly exponential.
generally means slower searches. A number of improvements to the basic algorithm have
been suggested. Friedman recommends that the partitioning hyperplane be chosen such that
it passes through the median point and is placed perpendicular to the coordinate axis along
whose direction the spread of the points is maximum [ 15 ] . Sproull suggests using a truncated
distance computation to increase efficiency in high dimensions [ 36 ] . Variants of the k-d tree
have been used to address specific search problems
An R-tree is also a space partitioning structure, but unlike k-d trees, the partitioning
element is not a hyperplane but a hyper-rectangular region [ This hierarchical rectangular
structure is useful in applications such as searching by image content [
needs to locate the closest manifold (or cluster) to a novel manifold (or cluster). An R-tree
also addresses some of the problems involved in implementing k-d trees in large disk based
databases. The R-tree is also a dynamic data structure, but unlike the k-d tree, the search
performance is not affected by addition or deletion of data. A number of variants of R-Trees
improve on the basic technique, such as packed R-trees [ 34
Although R-trees are useful in implementing sophisticated queries and managing large
databases, the performance of nearest neighbor point searches in high dimensions is very
similar to that of k-d trees; complexity grows exponentially with d.
Other static space partitioning techniques have been proposed such as branch and bound
none of which significantly improve
performance for high dimensions. Clarkson describes a randomized algorithm which finds
the closest point in d dimensional space in O(log 2 n) operations using a RPO (randomized
post . However, the time taken to construct the RPO tree is O(n dd=2e(1+ffl) )
and the space required to store it is also O(n dd=2e(1+ffl) ). This makes it impractical when the
number of points n is large or or if d ? 3.
3 The Algorithm
3.1 Searching by Slicing
We illustrate the proposed high dimensional search algorithm using a simple example in 3-D
space, shown in Figure 1. We call the set of points in which we wish to search for the closest
point as the point set. Then, our goal is to find the point in the point set that is closest to
2e
2e
2e
2e
x-e x+e
z-e
z+e
Y
Z

Figure

1: The proposed algorithm efficiently finds points inside a cube of size 2ffl around the novel
query point Q. The closest point is then found by performing an exhaustive search within the cube
using the Euclidean distance metric.
a novel query point Q(x; and within a distance ffl. Our approach is to first find all the
points that lie inside a cube (see Figure 1) of side 2ffl centered at Q. Since ffl is typically small,
the number of points inside the cube is also small. The closest point can then be found by
performing an exhaustive search on these points. If there are no points inside the cube, we
know that there are no points within ffl.
The points within the cube can be found as follows. First, we find the points that are
sandwiched between a pair of parallel planes X 1 and X 2 (see Figure 1) and add them to a
list, which we call the candidate list. The planes are perpendicular to the first axis of the
coordinate frame and are located on either side of point Q at a distance of ffl. Next, we trim
the candidate list by discarding points that are not also sandwiched between the parallel
pair of planes Y 1 and Y 2 , that are perpendicular to X 1 and X 2 , again located on either side
of Q at a distance ffl. This procedure is repeated for planes Z 1 and Z 2 , at the end of which,
the candidate list contains only points within the cube of size 2ffl centered on Q.
Since the number of points in the final trimmed list is typically small, the cost of the
exhaustive search is negligible. The major computational cost in our technique is therefore
in constructing and trimming the candidate list.
3.2 Data Structure
Candidate list construction and trimming can done in a variety of ways. Here, we propose
a method that uses a simple pre-constructed data structure along with 1-D binary searches
to efficiently find points sandwiched between a pair of parallel hyperplanes. The data
structure is constructed from the raw point set and is depicted in Figure 2. It is assumed
that the point set is static and hence, for a given point set, the data structure needs to be
constructed only once. The point set is stored as a collection of d 1-D arrays, where the j th
array contains the j th coordinate of the points. Thus, in the point set, coordinates of a point
lie along the same row. This is illustrated by the dotted lines in Figure 2. Now suppose that
novel point Q has coordinates . Recall that in order to construct the candidate
list, we need to find points in the point set that lie between a pair of parallel hyperplanes
separated by a distance 2ffl, perpendicular to the first coordinate axis, and centered at
that is, we need to locate points whose first coordinate lies between the limits
ffl. This can be done with the help of two binary searches, one for each limit, if the
coordinate array were sorted beforehand.
To this end, we sort each of the d coordinate arrays in the point set independently to
obtain the ordered set. Unfortunately, sorting raw coordinates does not leave us with any
information regarding which points in the arrays of the ordered set correspond to any given
point in the point set, and vice versa. For this purpose, we maintain two maps. The backward
map maps a coordinate in the ordered set to the corresponding coordinate in the point set
and, conversely, the forward map maps a point in the point set to a point in the ordered
set. Notice that the maps are simple integer arrays; if P [d][n] is the point set, O[d][n] is
the ordered set, F [d][n] and B[d][n] are the forward and backward maps, respectively, then
Using the backward map, we find the corresponding points in the point set (shown as
dark shaded areas) and add the appropriate points to the candidate list. With this, the
construction of the candidate list is complete. Next, we trim the candidate list by iterating
POINT SET
ORDERED SET
Dimensions
Points
Forward
Map
Forward
Map
Backward
Map
Input
-e
+e
-e
-e
+e
+e
,.,

Figure

2: Data structures used for constructing and trimming the candidate list. The point set
corresponds to the raw list of data points, while in the ordered set each coordinate is sorted. The
forward and backward maps enable efficient correspondence between the point and ordered sets.
on as follows. In iteration k, we check every point in the candidate list, by
using the forward map, to see if its k th coordinate lies within the limits
Each of these limits are also obtained by binary search. Points with k th coordinates that lie
outside this range (shown in light grey) are discarded from the list.
At the end of the final iteration, points remaining on the candidate list are the ones
which lie inside a hypercube of side 2ffl centered at Q. In our discussion, we proposed
constructing the candidate list using the first dimension, and then performing list trimming
using dimensions 2; 3; d, in that order. We wish to emphasize that these operations can
be done in any order and still yield the desired result. In the next section, we shall see that it
is possible to determine an optimal ordering such that the cost of constructing and trimming
the list is minimized.
It is important to note that the only operations used in trimming the list are integer
comparisons and memory lookups. Moreover, by using the proposed data structure, we have
limited the use of floating point operations to just the binary searches needed to find the
row indices corresponding to the hyperplanes. This feature is critical to the efficiency of the
proposed algorithm, when compared with competing ones. It not only facilitates a simple
software implementation, but also permits the implementation of a hardware search engine.
As previously stated, the algorithm needs to be supplied with an "appropriate" ffl prior
to search. This is possible for a large class of problems (in pattern recognition, for instance)
where a match can be declared only if the novel point Q is sufficiently close to a database
point. It is reasonable to assume that ffl is given a priori, however, the choice of ffl can prove
problematic if this is not the case. One solution is to set ffl large, but this might seriously
impact performance. On the other hand, a small ffl could result in the hypercube being
empty. How do we determine an optimal ffl for a given problem? How exactly does ffl affect
the performance of the algorithm? We seek answers to these questions in the following
section.
4 Complexity
In this section, we attempt to analyze the computational complexity of data structure stor-
age, construction and nearest neighbor search. As we saw in the previous section, constructing
the data structure is essentially sorting d arrays of size n. This can be done in
O(dn log 2 n) time. The only additional storage necessary is to hold the forward and bacward
maps. This requires space O(nd). For nearest neighbor search, the major computational cost
is in the process of candidate list construction and trimming. The number of points initially
added to the candidate list depends not only on ffl, but also on the distribution of data in the
point set and the location of the novel point Q. Hence, to facilitate analysis, we structure
the problem by assuming widely used distributions for the point set. The following notation
is used. Random variables are denoted by uppercase letters, for instance, Q. Vectors are in
bold, such as, q. Suffixes are used to denote individual elements of vectors, for instance, Q k
is the k th element of vector Q. Probability density is written as
and as f Q (q) if Q is continuous.
2e
Axis c

Figure

3: The projection of the point set and the novel point onto one of the dimensions of the
search space. The number of points inside bin B is given by the binomial distribution.

Figure

3 shows the novel point Q and a set of n points in 2-D space drawn from a known
distribution. Recall that the candidate list is initialized with points sandwiched between
a hyperplane pair in the first dimension, or more generally, in the c th dimension. This
corresponds to the points inside bin B in Figure 3, where the entire point set and Q are
projected to the c th coordinate axis. The boundaries of bin B are where the hyperplanes
intersect the axis c, at c be the number of points in bin B. In order
to determine the average number of points added to the candidate list, we must compute
c to be the distance between Q c and any point on the candidate list. The
distribution of Z c may be calculated from the the distribution of the point set. Define P c
to be the probability that any projected point in the point set is within distance ffl from Q c ;
that is,
It is now possible to write an expression for the density of M c in terms of P c . Irrespective of
the distribution of the points, M c is binomially distributed
2 This is equivalent to the elementary probability problem: given that a success (a point is within bin
B) can occur with probability P c
, the number of successes that occur in n independent trials (points) is
binomially distributed.
From the above expression, the average number of points in bin B, E[M c j Q c ], is easily
determined to be
Note that E[M c j Q c ] is itself a random variable that depends on c and the location of Q. If
the distribution of Q is known, the expected number of points in the bin can be computed as
we perform one lookup in the backward map for every point
between a hyperplane pair, and this is the main computational effort, equation (3) directly
estimates the cost of candidate list construction.
Next, we derive an expression for the total number of points remaining on the candidate
list as we trim through the dimensions in the sequence c 1 . Recall that in the
iteration k, we perform a forward map lookup for every point in the candidate list and see if
it lies between the c k
th hyperplane pair. How many points on the candidate list lie between
this hyperplane pair? Once again, equation (3) can be used, this time replacing n with the
number of points on the candidate list rather than the entire point set. We assume that the
point set is independently distributed. Hence, if N k is the total number of points on the
candidate list before the iteration k,
Y
Define N to be the total cost of constructing and trimming the candidate list. For each
trim, we need to perform one forward map lookup and two integer comparisons. Hence, if
we assign one cost unit to each of these operations, an expression for N can be written with
the aid of equation (4) as
Y
which, on the average is
Y
Equation (6) suggests that if the distributions f Q (q) and f Z (z) are known, we can compute
the average cost E[N in terms of ffl. In the next section, we shall examine two
cases of particular interest: (a) Z is uniformly distributed, and (b) Z is normally distributed.
Note that we have left out the cost of exhaustive search on points within the final hypercube.
The reason is that the cost of an exhaustive search is dependant on the distance metric used.
This cost is however very small and can be neglected in most cases when n AE d. If it needs
to be considered, it can be added to equation (6).
We end this section by making an observation. We had mentioned earlier that it is of
advantage to examine the dimensions in a specific order. What is this order? By expanding
the summation and product and by factoring terms, equation (5) can be rewritten as
It is immediate that the value of N is minimum when P c 1
. In other
words, should be chosen such that the numbers of sandwiched points between
hyperplane pairs are in ascending order. This can be easily ensured by simply sorting
the numbers of sandwiched points. Note that there are only d such numbers, which can
be obtained in time O(d) by simply taking the difference of the indices to the ordered
set returned by each pair of binarysearchs. Further, the cost of sorting these numbers
is O(d log 2 d) by heapsort these costs are negligible in any problem of
reasonable dimensionality.
4.1 Uniformly Distributed Point Set
We now look at the specific case of a point set that is uniformly distributed. If X is a point
in the point set, we assume an independent and uniform distribution with extent l on each
of it's coordinates as
1=l if \Gammal=2 - x - l=2
1.5e+062.5e+06Cost
e
d=5
e5000001.5e+062.5e+06Cost
n=50000
n=150000
(a) (b)

Figure

4: The average cost of the algorithm is independent of d and grows only linearly for small
ffl. The point set in both cases is assumed to be uniformly distributed with extent l = 1. (a) The
point set contains 100,000 points in 5-D, 10-D, 15-D, 20-D and 25-D spaces. (b) The point set is
15-D and contains 50000, 75000, 100000, 125000 and 150000 points.
Using equation (8) and the fact that Z , an expression for the density of Z c can
be written as
f Zc jQc
P c can now be written as
\Gammaffl
f Zc jQc (z)dz
\Gammaffll
dz
l
Substituting equation (10) in equation (6) and considering the upper bound (worst case),
we get
l
l
l
l
l
l
l
By neglecting constants, we write
d=5
Cost
e
n=50000
Cost
(a) (b)

Figure

5: The average cost of the algorithm is independent of d and grows only linearly for small
ffl. The point set in both cases is assumed to be normally distributed with variance oe = 1. (a) The
point set contains 100,000 points in 5-D, 10-D, 15-D, 20-D and 25-D spaces (b) The point
set is 15-D and contains 50000, 75000, 100000, 125000 and 150000 points
For small ffl, we observe that ffl d - 0, because of which cost is independent of d:
In

Figure

4, equation (11) is plotted against ffl for different d (Figure 4(a)) and different

Figure

1. Observe that as long as ffl ! :25, the cost varies little with d
and is linearly proportional to n. This also means that keeping ffl small is crucial to the
performance of the algorithm. As we shall see later, ffl can in fact be kept small for many
problems. Hence, even though the cost of our algorithm grows linearly with n, ffl is small
enough that in many real problems, it is better to pay this price of linearity, rather than an
exponential dependence on d.
4.2 Normally Distributed Point Set
Next, we look at the case when the point set is normally distributed. If X is a point in the
point set, we assume an independent and normal distribution with variance oe on each of it's
coordinates:
f Xc (x) =p
2-oe
exp
As before, using Z , an expression for the density of Z c can be obtained to get
f Zc jQc
2-oe
P c can then be written as
\Gammaffl
f Zc jQc (z)dz
oe
oe
p!
This expression can be substituted into equation (6) and evaluated numerically to estimate
cost for a given Q. Figure 5 shows the cost as a function of ffl for As
with uniform distribution, we observe that when ffl ! 1, the cost is nearly independent of d
and grows linearly with n. In a variety of pattern classification problems, data take the form
of individual Gaussian clusters or mixtures of Gaussian clusters. In such cases, the above
results can serve as the basis for complexity analysis.
5 Determining ffl
It is apparent from the analysis in the preceding section that the cost of the proposed
algorithm depends critically on ffl. Setting ffl too high results in a huge increase in cost with
d, while setting ffl too small may result in an empty candidate list. Although the freedom
to choose ffl may be attractive in some applications, it may prove non-intuitive and hard in
others. In such cases, can we automatically determine ffl so that the closest point can be
found with high certainty? If the distribution of the point set is known, we can.
We first review well known facts about L p norms. Figure 6 illustrates these norms for
a few selected values of p. All points on these surfaces are equidistant (in the sense of the
respective norm) from the central point. More formally, the L p distance between two vectors
a and b is defined as
These distance metrics are also known as Minkowski-p metrics. So how are these relevant
to determining ffl? The L 2 norm occurs most frequently in pattern recognition problems.
Unfortunately, candidate list trimming in our algorithm does not find points within L 2 , but

Figure

An illustration of various norms, also known as Minkowski p-metrics. All points on these
surfaces are equidistant from the central point. The L1 metric bounds L p for all p.
within L1 (i.e. the hypercube). Since L1 bounds L 2 , one can naively perform an exhaustive
search inside L1 . However, as seen in figure 7(a), this does not always correctly find the
closest point. Notice that P 2
is closer to Q than P 1
, although an exhaustive search within
the cube will incorrectly identify
to be the closest. There is a simple solution to this
problem. When performing an exhaustive search, impose an additional constraint that only
points within an L 2 radius ffl should be considered (see figure 7(b)). This, however, increases
the possibility that the hypersphere is empty. In the above example, for instance, P 1 will
be discarded and we would not be able to find any point. Clearly then, we need to consider
this fact in our automatic method of determining ffl which we describe next.
We propose two methods to automatically determine ffl. The first computes the radius of
the smallest hypersphere that will contain at least one point with some (specified) probability.
ffl is set to this radius and the algorithm proceeds to find all points within a circumscribing
hypercube of side 2ffl. This method is however not efficient in very high dimensions; the reason
being as follows. As we increase dimensionality, the difference between the hypersphere and
hypercube volumes becomes so great that the hypercube "corners" contain far more points
than the inscribed hypersphere. Consequently, the extra effort necessary to perform L 2
distance computations on these corner points is eventually wasted. So rather than find the
circumscribing hypercube, in our second method, we simply find the length of a side of the
smallest hypercube that will contain at least one point with some (specified) probability. ffl
can then be set to half the length of this side. This leads to the problem we described earlier
that, when searching some points outside a hypercube can be closer in the L 2 sense than
points inside. We shall now describe both the methods in detail and see how we can remedy
e
r
2e
e
(a) (b)

Figure

7: An exhaustive search within a hypercube may yield an incorrect result. (a) P 2 is closer
to
, but just an exhaustive search within the cube will incorrectly identify
as the
closest point. (b) This can be remedied by imposing the constraint that the exhaustive search
should consider only points within an L 2 distance ffl from Q (given that the length of a side of the
hypercube is 2ffl).
this problem.
5.1 Smallest Hypersphere Method
Let us now see how to analytically compute the minimum size of a hypersphere given that
we want to be able guarantee that it is non empty with probability p. Let the radius of such
a hypersphere be ffl hs . Let M be the total number of points within this hypersphere. Let Q
be the novel point and define kZk to be the L 2 distance between Q and any point in the
point set. Once again, M is binomially distributed with the density
Now, the probability p that there is at least one point in the hypersphere is simply
e hs
r
(a) (b)

Figure

8: ffl can be computed using two methods: (a) By finding the radius of the smallest hyper-sphere
that will contain at least one point with high probability. A search is performed by setting
ffl to this radius and constraining the exhaustive search within ffl. (b) By finding the size of the
smallest hypercube that will contain at least one point with high probability. When searching, ffl is
set to half the length of a side. Additional searches have to be performed in the areas marked in
bold.
The above equation suggests that if we know Q, the density f Z jQ (z), and the probability
p, we can solve for ffl hs .
For example, consider the case when the point set is uniformly distributed with density
given by equation (9). The cumulative distribution function of kZk is the uniform
distribution integrated within a hypersphere; which is simply it's volume. Thus,
l d d\Gamma(d=2)
Substituting the above in equation 19 and solving for ffl hs , we get
l d d\Gamma(d=2)
Using equation (21), ffl hs is plotted against probability for two cases. In figure 9(a), d is fixed
to different values between 5 to 25 with n is fixed to 100000, and in figure 9(b), n is fixed
to different values between 50000 to 150000 with d fixed to 5. Both the figures illustrate
an important property, which is that large changes in the probability p result in very small
Probability of Success0.20.61Epsilon
d=5
Probability of Success0.050.15
Epsilon
n=150000
n=50000
(a) (b)

Figure

9: The radius ffl necessary to find a point inside a hypersphere varies very little with
probability. This means that ffl can be set to the knee where probability is close to unity. The point
set in both cases is uniformly distributed with extent l = 1. (a) The point set contains 100000
points in 5, 10, 15, 20 and 25 dimensional space. (b) The point is 5-D and contains 50000, 75000,
100000, 125000 and 150000 points.
changes in ffl hs . This suggests that ffl hs can be set to the right hand "knee" of both the curves
where probability is very close to unity. In other words, it is easy to guarantee that at least
one point is within the hypersphere. A search can now be performed by setting the length
of a side of the circumscribing hypercube to 2ffl hs and by imposing an additional constraint
during exhaustive search that only points within an L 2 distance ffl hs be considered.
5.2 Smallest Hypercube Method
As before, we attempt to analytically compute the size of the smallest hypercube given
that we want to be able guarantee that it is non empty with probability p. Let M be the
number of points within a hypercube of size 2ffl hc . Define Z c to be the distance between the
c th coordinate of a point set point and the novel point Q. Once again, M is binomially
distributed with the density
Y
!k/
d
Y
kC A :(22)
Now, the probability p that there is at least one point in the hypercube is simply
Probability of Success0.10.30.5
Epsilon
d=5
Probability of Success0.050.150.250.35
Epsilon n=150000
n=50000
(a) (b)

Figure

10: The value of ffl necessary to find a point inside a hypercube varies very little with
probability. This means that ffl can be set to the knee where probability is close to unity. The point
set in both cases is uniformly distributed with extent l = 1. (a) The point set contains 100000
points in 5, 10, 15, 20 and 25 dimensional space. (b) The point set is 5-D and contains 50000,
75000, 100000, 125000 and 150000 points.
d
Y
!n
Again, above equation suggests that if we know Q, the density f Zc jQc (z), and the probability
p, we can solve for ffl hc . If for the specific case that the point set is uniformly distributed, an
expression for ffl hc can be obtained in closed form as follows. Let the density of the uniform
distribution be given by equation (9). Using equation (10) we get,
d
Y
l
Substituting the above in equation (23) and solving for ffl hc , we get
li
Using equation (25), ffl hc is plotted against probability for two cases. In figure 10(a), d is fixed
to different values between 5 to 25 with n is fixed to 100000, and in figure 10(b), n is fixed to
different values between 50000 to 150000 with d fixed to 5. These are similar to the graphs
obtained in the case of a hypersphere and again, ffl hc can be set to the right hand "knee" of
both the curves where probability is very close to unity. Notice that the value of ffl hc required
for the hypercube is much smaller than that required for the hypersphere, especially in high
d. This is precisely the reason why we prefer the second (smallest hypercube) method.
Recall that it is not sufficient to simply search for the closest point within a hypercube
because a point outside can be closer than a point inside. To remedy this problem, we
suggest the following technique. First, an exhaustive search is performed to compute the
distance to the closest point within the hypercube. Call this distance r. In figure 8(b),
the closest point
within the hypercube is at a distance of r from Q. Clearly, if a closer
point exists, it can only be within a hypersphere of radius r. Since parts of this hypersphere
lie outside the original hypercube, we also search in the hyper-rectangular regions shown in
bold (by performing additional list trimmings). When performing an exhaustive search in
each of these hyper-rectangles, we impose the constraint that a point is considered only if it
is less than distance r from Q. In figure 8(b), P 2
is present in one such hyper-rectangular
region and happens to be closer to Q than P 1 . Although this method is more complicated,
it gives excellent performance in sparsely populated high dimensional spaces (such as a high
dimensional uniform distribution).
To conclude, we wish to emphasize that both the hypercube and hypersphere methods can
be used interchangeably and both are guaranteed to find the closest point within ffl. However,
the choice of which one of these methods to use should depend on the dimensionality of
the space and the local density of points. In densely populated low dimensional spaces,
the hypersphere method performs quite well and searching the hyper-rectangular regions
is not worth the additional overhead. In sparsely populated high dimensional spaces, the
effort needed to exhaustively search the huge circumscribing hypercube is far more than the
overhead of searching the hyper-rectangular regions. It is, however, difficult to analytically
predict which one of these methods suits a particular class of data. Hence, we encourage the
reader to implement both the methods and use the one which performs the best. Finally,
although the above discussion is relevant only for the L 2 norm, an equivalent analysis can
be easily performed for any other norm.
6 Benchmarks
We have performed an extensive set of benchmarks on the proposed algorithm. We looked
at two representative classes of search problems that may benefit from the algorithm. In
the first class, the data has statistical structure. This is the case, for instance, when points
are uniformly or normally distributed. The second class of problems are statistically un-
structured, for instance, when points lie on a high dimensional multivariate manifold, and it
is difficult to say anything about their distribution. In this section, we will present results
for benchmarks performed on statistically structured data. For benchmarks on statistically
unstructured data, we refer the reader to section 7.
We tested two commonly occurring distributions, normal and uniform. The proposed algorithm
was compared with the k-d tree and exhaustive search algorithms. Other algorithms
were not included in this benchmark because they did not yield comparable performance.
For the first set of benchmarks, two normally distributed point sets containing 30,000 and
100,000 points with variance 1.0 were used. To test the per search execution time, another
set of points, which we shall call the test set, was constructed. The test set contained 10,000
points, also normally distributed with variance 1.0. For each algorithm, the execution time
was calculated by averaging over the total time required to perform a nearest neighbor search
on each of the 10,000 points in the test set. To determine ffl, we used the 'smallest hypercube'
method described in section 5.2. Since the point set is normally distributed, we cannot use a
closed form solution for ffl. However, it can be numerically computed as follows. Substituting
equation (16) into equation (23), we get
d
Y
oe
oe
p!!n
By setting p (the probability that there is at least one point in the hypercube) to .99 and
oe (the variance) to 1.0, we computed ffl for each search point Q using the fast and simple
bisection technqiue [

Figures

11(a) and 11(b) show the average execution time per search when the point set
contains 30,000 and 100,000 points respectively. These execution times include the time
taken for search, computation of ffl using equation (26), and the time taken for the few (1%)
additional 3 searches necessary when a point was not found within the hypercube. Although
3 When a point was not found within the hypercube, we incremented ffl by 0.1 and searched again. This
Time
(secs.)
Dimensions
Proposed algorithm
k-d tree
Exhaustive
Time
(secs.)
Dimensions
Proposed algorithm
k-d tree
Exhaustive search
(a) (b)0.050.150.250.35
Time
(secs.)
Dimensions
Proposed algorithm
k-d tree
Exhaustive
Time
(secs.)
Dimensions
Proposed algorithm
k-d tree
Exhaustive search
(c) (d)

Figure

11: The average execution time of the proposed algorithm is benchmarked for statistically
structured problems. (a) The point set is normally distributed with variance 1.0 and contains
30,000 points. (b) The point set is normally distributed with variance 1.0 and contains 100,000
points. The proposed algorithm is clearly faster in high d. (c) The point set is uniformly distributed
with extent 1.0 and contains 30,000 points. (d) The point set is uniformly distributed with extent
1.0 and contains 100,000 points. The proposed algorithm does not perform as well for uniform
distributions due to the extreme sparseness of the point set in high d.
ffl varies for each Q, values of ffl for a few sample points are as follows. For the
values of ffl at the point corresponding to
respectively. At the point the values of ffl were
0:24; 0:61; 0:86; 1:04; 1:17 corresponding to respectively. For
the values of ffl at the point corresponding
to respectively. At the point the values of ffl were
corresponding to
that the proposed algorithm is faster than the k-d tree algorithm for all d in Figure 11(a).
In

Figure

11(b), the proposed algorithm is faster for d ? 12. Also notice that the k-d
tree algorithm actually runs slower than exhaustive search for d ? 15. The reason for this
observation is as follows. In high dimensions, the space is so sparsely populated that the
radius of the query hypersphere is very large. Consequently, the hypersphere intersects
almost all the buckets and thus a large number of points are examined. This, along with the
additional overhead of traversing the tree structure makes it very inefficient to search the
sparse high dimensional space.
For the second set of benchmarks, we used uniformly distributed point sets containing
30,000 and 100,000 points with extent 1.0. The test set contained 10,000 points, also
uniformly distributed with extent 1.0. The execution time per search was calculated by
averaging over the total time required to perform a closest point search on each of the 10,000
points in the test set. As before, to determine ffl, the 'smallest hypercube' method described
in section 5.2 was used. Recall that for uniformly distributed point sets, ffl can be computed
in the closed form using equation (25). Figures 11(c) and 11(d) show execution times
when the point set contains 30,000 and 100,000 points respectively. For the
values of ffl were corresponding to
tively. For the values of ffl were corresponding to
respectively. For uniform distribution, the proposed algorithm does not
perform as well, although, it does appear to be slightly faster than the k-d tree and exhaustive
search algorithms. The reason is that the high dimensional space is very sparsely populated
and hence requires ffl to be quite large. As a result, the algorithm ends up examining almost
all points, thereby approaching exhaustive search.
process was repeated till a point was found.
7 An Example Application: Appearance Matching
We now demonstrate two applications where a fast and efficient high dimensional search
technique is desirable. The first, real time object recognition, requires the closest point to be
found among 36,000 points in a 35-D space. In the second, the closest point is required to be
found from points lying on a multivariate high dimensional manifold. Both these problems
are examples of statistically unstructured data.
Let us briefly review the object recognition technique of Murase and Nayar [ 24 ] . Object
recognition is performed in two phases: 1) appearance learning phase, and 2) appearance
recognition phase. In the learning phase, images of each of the hundred objects in all poses
are captured. These images are used to compute a high dimensional subspace, called the
eigenspace. The images are projected to eigenspace to obtain discrete high dimensional
points. A smooth curve is then interpolated through points that belong to the same object.
In this way, for each object, we get a curve (or a univariate manifold) parameterized by it's
pose. Once we have the manifolds, the second phase, object recognition, is easy. An image
of an object is projected to eigenspace to obtain a single point. The manifold closest to this
point identifies the object. The closest point on the manifold identifies the pose. Note that
the manifold is continuous, so in order to find the closest point on the manifold, we need to
finely sample it to obtain discrete closely spaced points.
For our benchmark, we used the Columbia Object Image Library [ along with the
SLAM software package [ 28 ] to compute 100 univariate manifolds in a 35-D eigenspace.
These manifolds correspond to appearance models of the 100 objects (20 of the 100 objects
shown in Figure 12(a)). Each of the 100 manifolds were sampled at 360 equally spaced
points to obtain 36,000 discrete points in 35-D space. It was impossible to manually capture
the large number of object images that would be needed for a large test set. Hence, we
automatically generated a test set of 100,000 points by sampling the manifolds at random
locations. This is roughly equivalent to capturing actual images, but, without image sensor
noise, lens blurring, and perspective projection effects. It is important to simulate these
effects because they cause the projected point to shift away from the manifold and hence,
substantially affect the performance of nearest neighbor search algorithms 4 .
4 For instance, in the k-d tree, a large query hypersphere would result in a large increase in the number
of adjacent buckets that may have to be searched.
Algorithm Time (secs.)
Proposed Algorithm .0025
k-d tree .0045
Exhaustive Search .1533
Projection Search .2924
(b)

Figure

12: The proposed algorithm was used to recognize and estimate pose of hundred objects
using the Columbia Object Image Library. (a) Twenty of the hundred objects are shown. The point
set consisted of 36,000 points (360 for each object) in 35-D eigenspace. (b) The average execution
time per search is compared with other algorithms.
Unfortunately, it is very difficult to relate image noise, perspective projection and other
distortion effects to the location of points in eigenspace. Hence, we used a simple model
where we add uniformly distributed noise with extent 5 .01 to each of the coordinates of
points in the test set. We found that this approximates real-world data. We determined
that setting gave us good recognition accuracy. Figure 12(b) shows the time taken
per search by the different algorithms. The search time was calculated by averaging the
total time taken to perform 100,000 closest point searches using points in the test set. It
can be seen that the proposed algorithm outperforms all the other techniques. ffl was set to
a predetermined value such that a point was found within the hypersphere all the time. For
object recognition, it is useful to search for the closest point within ffl because this provides
us with a means to reject points that are "far" from the manifold (most likely from objects
not in the database).
Next, we examine another case when data is statistically unstructured. Here, the closest
point is required to be found from points lying on a single smooth multivariate high dimensional
manifold. Such a manifold appears frequently in appearance matching problems such
as visual tracking [ 26 ] , visual inspection [ 26 ] , and parametric feature detection [ 25 ] . As with
object recognition, the manifold is a representation of visual appearance. Given a novel
appearance (point), matching involves finding a point on the manifold closest to that point.
Given that the manifold is continuous, to pose appearance matching as a nearest neighbor
problem, as before, we sample the manifold densely to obtain discrete closely spaced points.
The trivariate manifold we used in our benchmarks was obtained from a visual tracking
experiment conducted by Nayar et al. [ 26 ] . In the first benchmark, the manifold was sampled
to obtain 31,752 discrete points. In the second benchmark, it was sampled to obtain 107,163
points. In both cases, a test set of 10,000 randomly sampled manifold points was used. As
explained previously, noise (with extent .01) was added to each coordinate in the test set.
The execution time per search was averaged over this test set of 10,000 points. For this point
set, it was determined that gave good recognition accuracy. Figure 13(a) shows the
algorithm to be more than two orders of magnitude faster than the other algorithms. Notice
the exponential behaviour of the R-tree algorithm. Also notice that Yunck's algorithm is
5 The extent of the eigenspace is from -1.0 to +1.0. The maximum noise amplitude is hence about 0.5%
of the extent of eigenspace.
Time
(secs.)
Dimensions
Proposed algorithm
R-tree
Exhaustive search
Projection search
Yunck search
(a)0.0050.0150.0250.035
Time
(secs.)
Dimensions
Proposed algorithm
k-d tree0.0050.0150.0250.035
Time
(secs.)
Dimensions
Proposed algorithm
k-d tree
(b) (c)

Figure

13: The average execution time of the proposed algorithm is benchmarked for an unstructured
problem. The point set is constructed by sampling a high dimensional trivariate manifold.
(a) The manifold is sampled to obtain 31,752 points. The proposed algorithm is more than two
orders of magnitude faster than the other algorithms. (b) The manifold is sampled as before to
obtain 31,752 points. (c) The manifold is sampled to obtain 107,163 points. The k-d tree algorithm
is slightly faster in low dimension but degrades rapidly with increase in dimension.
only slightly faster than Friedman's; the difference is due to use of integer operations. We
could only benchmark Yunck's algorithm till due to use of a 32-bit word in the
indicator array. In Figure 13(b), it can be seen that the proposed algorithm is faster than
the k-d tree for all d, while in Figure 13(c), the proposed algorithm is faster for all d ? 21.
8 Hardware Architecture
A major advantage of our algorithm is its simplicity. Recall that the main computations
performed by the algorithm are simple integer map lookups (backward and forward maps)
and two integer comparisons (to see if a point lies within hyperplane boundaries). Conse-
quently, it is possible to implement the algorithm in hardware using off-the-shelf, inexpensive
components. This is hard to envision in the case of any competitive techniques such as k-d
trees or R-trees, given the difficulties involved in constructing parallel stack machines.
The proposed architecture is shown in Figure 14. A Field Programmable Gate Array
acts as an algorithm state machine controller and performs I/O with the CPU. The
Dynamic RAMs (DRAMs) hold the forward and backward maps which are downloaded from
the CPU during initialization. The CPU initiates a search by performing a binary search
to obtain the hyperplane boundaries. These are then passed on to the search engine and
held in the Static RAMs (SRAMs). The FPGA then independently begins the candidate list
construction and trimming. A candidate is looked up in the backward map and each of the
forward maps. The integer comparator returns a true if the candidate is within range, else it
is discarded. After trimming all the candidate points by going through the dimensions, the
final point list (in the form of point set indices) is returned to the CPU for exhaustive search
and/or further processing. Note that although we have described an architecture with a
single comparator, any number of them can be added and run in parallel with a near linear
performance scaling in the number of comparators. While the search engine is trimming the
candidate list, the CPU is of course free to carry out other tasks in parallel.
We have begun implementation of the proposed architecture. The result is intended to
be a small low-cost SCSI based module that can be plugged in to any standard workstation
or PC. We estimate the module to result in a 100 fold speedup over an optimized software
implementation.
ALGORITHM
CONTROL
UNIT
BACKWARD
MAP
FORWARD
MAP
A
A
CPU
COMPARATOR
WITHIN LIMIT FLAG
LOWER LIMIT
UPPER LIMIT
CONTROL BUS
FPGA/PLD

Figure

14: Architecture for an inexpensive hardware search engine that is based on the proposed
algorithm.
9 Discussion
9.1 k Nearest Neighbor Search
In section 5, we saw that it is possible to determine the minimum value of ffl necessary
to ensure that at least one point is found within a hypercube or hypersphere with high
probability. It is possible to extend this notion to ensure that at least k points are found
with high certainty. Recall that the probability that there exists at least one point in a
hypersphere of radius ffl is given by equation (19). Now define p k to be the probability that
there are at least k points within the hypersphere. We can then write p k as
The above expression can now be substituted in equation (18) and, given p k , numerically
solved for ffl hs . Similarly, it can be substituted in equation (22) to compute the minimum
value of ffl hc for a hypercube.
9.2 Dynamic Point Insertion and Deletion
Currently, the algorithm uses d floating point arrays to store the ordered set, and 2d integer
arrays to store the backward and forward maps. As a result, it is not possible to efficiently
insert or delete points in the search space. This limitation can be easily overcome if the
ordered set is not stored as an array but as a set of d binary search trees (BST) (each
BST corresponds to an array of the ordered set). Similarly, the d forward maps have to be
replaced with a single linked list. The backward maps can be done away with completely
as the indices can be made to reside within a node of the BST. Although BSTs would allow
efficient insertion and deletion, nearest neighbor searches would no longer be as efficient as
with integer arrays. Also, in order to get maximum efficiency, the BSTs would have to be
well balanced (see [ 19 ] for a discussion on balancing techniques).
9.3 Searching with Partial Data
Many times, it is required to search for the nearest neighbor in the absence of complete
data. For instance, consider an application which requires features to be extracted from an
image and then matched against other features in a feature space. Now, if it is not possible
to extract all features, then the matching has to be done partially. It is trivial to adapt
our algorithm to such a situation: while trimming the list, you need to only look at the
dimensions for which you have data. This is hard to envision in the case of k-d trees for
example, because the space has been partitioned by hyperplanes in particular dimensions.
So, when traversing the tree to locate the bucket that contains the query point, it is not
possible to choose a traversal direction at a node if data corresponding to the partitioning
dimension at that node is missing from the query point.

Acknowledgements

We wish to thank Simon Baker and Dinkar Bhat for their detailed comments, criticisms and
suggestions which have helped greatly in improving the paper.
This research was conducted at the Center for Research on Intelligent Systems at the
Department of Computer Science, Columbia University. It was supported in parts by ARPA
Contract DACA-76-92-C-007, DOD/ONR MURI Grant N00014-95-1-0601, and a NSF National
Young Investigator Award.



--R

The Design and Analysis of Computer Algorithms.
Nearest neighbor searching and applications.
diagrams - a survey of a fundamental geometric data struc- ture

Multidimensional binary search trees used for associative searching.
Multidimensional binary search trees in database applications.
Data structures for range searching.
Optimal expected-time algorithms for closest point problems
Multidimensional indexing for recognizing visual shapes.
A randomized algorithm for closest-point queries
Multidimensional searching problems.
Algorithms in Combinatorial Geometry.
Fast nearest-neighbor search in dissimilarity spaces
An algorithm for finding nearest neigh- bors
An algorithm for finding best matches in logarithmic expected time.
A branch and bound algorithm for computing k-nearest neighbors
An effective way to represent quadtrees.
A dynamic index structure for spatial searching.
Fundamentals of Data Structures in.
On the complexity of d-dimensional voronoi diagrams
Sorting and Searching

A new version of the nearest-neighbor approximating and eliminating search algorithm (aesa) with linear preprocessing time and memory requirements
Visual learning and recognition of 3d objects from appear- ance
Parametric feature detection.


Slam: A software library for appearance matching.
Digital Pictures
Similarity searching in large image databases.
Computational Geometry: An Introduction.
Numerical Recipes in C.

Direct spatial search on pictorial databases using packed r-trees

Refinements to nearest-neighbor searching in k-dimensional trees
Reducing the overhead of the aesa metric-space nearest neighbour searching algorithm

Data structures and algorithms for nearest neighbor search in general metric spaces.
A technique to identify nearest neighbors.
--TR

--CTR
James McNames, A Fast Nearest-Neighbor Algorithm Based on a Principal Axis Search Tree, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.23 n.9, p.964-976, September 2001
David W. Patterson , Mykola Galushka , Niall Rooney, Characterisation of a Novel Indexing Technique for Case-Based Reasoning, Artificial Intelligence Review, v.23 n.4, p.359-393, June      2005
Philip Quick , David Capson, Subspace position measurement in the presence of occlusion, Pattern Recognition Letters, v.23 n.14, p.1721-1733, December 2002
T. Freeman , Thouis R. Jones , Egon C Pasztor, Example-Based Super-Resolution, IEEE Computer Graphics and Applications, v.22 n.2, p.56-65, March 2002
Bin Zhang , Sargur N. Srihari, Fast k-Nearest Neighbor Classification Using Cluster-Based Trees, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.26 n.4, p.525-528, April 2004
Jim Z. C. Lai , Yi-Ching Liaw , Julie Liu, Fast k-nearest-neighbor search based on projection and triangular inequality, Pattern Recognition, v.40 n.2, p.351-359, February, 2007
Lin Liang , Ce Liu , Ying-Qing Xu , Baining Guo , Heung-Yeung Shum, Real-time texture synthesis by patch-based sampling, ACM Transactions on Graphics (TOG), v.20 n.3, p.127-150, July 2001
Fast texture synthesis using tree-structured vector quantization, Proceedings of the 27th annual conference on Computer graphics and interactive techniques, p.479-488, July 2000
Yong-Sheng Chen , Yi-Ping Hung , Ting-Fang Yen , Chiou-Shann Fuh, Fast and versatile algorithm for nearest neighbor search based on a lower bound tree, Pattern Recognition, v.40 n.2, p.360-375, February, 2007
Mineichi Kudo , Naoto Masuyama , Jun Toyama , Masaru Shimbo, Simple termination conditions for k-nearest neighbor method, Pattern Recognition Letters, v.24 n.9-10, p.1203-1213, 01 June
John T. Favata, Offline General Handwritten Word Recognition Using an Approximate BEAM Matching Algorithm, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.23 n.9, p.1009-1021, September 2001
Edgar Chvez , Gonzalo Navarro, A compact space decomposition for effective metric indexing, Pattern Recognition Letters, v.26 n.9, p.1363-1376, 1 July 2005
Aaron Hertzmann , Steven M. Seitz, Example-Based Photometric Stereo: Shape Reconstruction with General, Varying BRDFs, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.27 n.8, p.1254-1264, August 2005
O. Stahlhut, Extending natural textures with multi-scale synthesis, Graphical Models, v.67 n.6, p.496-517, November 2005
Gonzalo Navarro, Searching in metric spaces by spatial approximation, The VLDB Journal  The International Journal on Very Large Data Bases, v.11 n.1, p.28-46, August 2002
Edgar Chvez , Jos L. Marroqun , Gonzalo Navarro, Fixed Queries Array: A Fast and Economical Data Structure for Proximity Searching, Multimedia Tools and Applications, v.14 n.2, p.113-135, June 2001
Jin-xiang Chai , Jing Xiao , Jessica Hodgins, Vision-based control of 3D facial animation, Proceedings of the ACM SIGGRAPH/Eurographics symposium on Computer animation, July 26-27, 2003, San Diego, California
Edgar Chvez , Gonzalo Navarro , Ricardo Baeza-Yates , Jos Luis Marroqun, Searching in metric spaces, ACM Computing Surveys (CSUR), v.33 n.3, p.273-321, September 2001
Huzefa Neemuchwala , Alfred Hero , Paul Carson, Image matching using alpha-entropy measures and entropic graphs, Signal Processing, v.85 n.2, p.277-296, February 2005
Gisli R. Hjaltason , Hanan Samet, Index-driven similarity search in metric spaces, ACM Transactions on Database Systems (TODS), v.28 n.4, p.517-580, December
Richard Szeliski, Image alignment and stitching: a tutorial, Foundations and Trends in Computer Graphics and Vision, v.2 n.1, p.1-104, January 2006
object recognition;benchmarks;nearest neighbor;searching by slicing;visual correspondence;pattern classification;hardware architecture