LB_Keogh Supports Exact Indexing of Shapes under Rotation Invariance with
Arbitrary Representations and Distance Measures
Eamonn Keogh · Li Wei · Xiaopeng Xi · Michail Vlachos · Sang-Hee Lee · Pavlos Protopapas
Dear Reader: This is an expanded version of our VLDB paper
Abstract Shape matching and indexing is important topic in its
own right, and is a fundamental subroutine in most shape data
mining algorithms. Given the ubiquity of shape, shape matching is
an important problem with applications in domains as diverse as
biometrics, industry, medicine, zoology and anthropology. The
distance/similarity measure for used for shape matching must be
invariant to many distortions, including scale, offset, noise,
articulation, partial occlusion, etc. Most of these distortions are
relatively easy to handle, either in the representation of the data orin the similarity measure used. However rotation invariance isnoted in the literature as being an especially difficult challenge.
Current approaches typically try to achieve rotation invariance inthe representation of the data, at the expense of discrimination
ability, or in the distance measure, at the expense of efficiency. In
this work we show that we can take the slow but accurate
approaches and dramatically speed them up. On real world
problems our technique can take current approaches and makethem four orders of magnitude faster, without false dismissals.
Moreover, our technique can be used with any of the dozens ofexisting shape representations and with all the most popular
distance measures including Euclidean distance, Dynamic Time
Warping and Longest Common Subsequence. We further showthat our indexing technique can be used to index star light curves,
an important type of astronomical data, without modification.

distance/similarity measure for used for shape matching must be
invariant to many distortions, including scale, offset, noise,
articulation, partial occlusion, etc.. Figure 1 gives a visual
intuition of these problems in a familiar domain, butterflies andmoths. Most of these distortions are relatively easy to handle,
particularly if we use the well-known technique of converting theshapes into time series as in Figure 2. However, no matter what
representation is used, rotation invariance seems to be uniquely
difficult to handle. For example [20] notes “rotation is always
something hard to handle compared with translation andscaling”.

Actias maenas

Papilio antimachus

Agrias sardanapalus

Keywords Shape, Indexing, Dynamic Time Warping

1. INTRODUCTION
Shape matching and indexing is important topic in its own right,
and is a fundamental subroutine in most shape data mining
algorithms. Given the ubiquity of shape, shape matching is an
important problem with applications in domains as diverse as
biometrics, industry, medicine, zoology and anthropology. The________________________________________________________
E. Keogh
L. WeiX. Xi
Department of Computer Science & Engineering, UCR.
M. Vlachos
IMB Watson
S-H LeeDepartment of Anthropology, UCR.
P. Protopapas
Harvard-Smithsonian Center for Astrophysics

Papiliorutulus

Sphinx
Ligustri

Figure 1: Examples of the distortions we may be interested in being
invariant to when matching shapes. The left column shows drawings ofinsects dating back to 1734 [32]. The right column shows real insects. The
flexible wingtips of Actias maenas require articulation invariance. One of
the Papilio antimachus must be resized before matching. The Agriassardanapalus need their offsets corrected in order to match. The real
Papilio rutulus has a broken wing which appears as an occlusion to shape
matching algorithms. The real Sphinx Ligustri needs to be rotated to matchthe drawing, achieving this invariance is the focus of this work

Many current approaches try to achieve rotation invariance in the
representation of the data, at the expense of discrimination ability[28], or in the distance measure, at the expense of efficiency[1][2][3][9].

A

B

The rest of this paper is organized as follows. In Section 2 wediscuss background material and related work. In Section 3 weformally introduce the problem and in Section 4 we offer oursolution. Section 5 offers a comprehensive empirical evaluation ofboth the effectiveness and efficiency of our technique. Finally
Section 6 offers some conclusions and directions for future work.

C

2. BACKGROUND AND RELATED WORK
0

200

400

600

800

1000

1200

1400

Figure 2: Shapes can be converted to time series. A) A bitmap of ahuman skull. B) The distance from every point on the profile to the centeris measured and treated as the Y-axis of a time series of length n (C)As an example of the former, the very efficient rotation invariant
technique of [28] cannot differentiate between the shapes of thelowercase letters “d”and “b”. As an example of the latter, thework of Adamek and Connor [1], which is state of the art in termsof accuracy or precision/recall takes an untenable O(n3) for eachshape comparison.
In this work we show that we can take the slow but accurate
approaches and dramatically speed them up. For example we cantake the O(n3) approach of [1] and on real world problems bringthe average complexity down to O(n1.06). This dramatic
improvement in efficiency does not come at the expense of
accuracy; we prove that we will always return the same answer set
as the slower methods.
We achieve speedup over the existing methods in two ways,
dramatically decreasing the CPU requirements, and allowing
indexing. Our technique works by grouping together similarrotations, and defining an admissible lower bound to that group.
Given such a lower bound, we can utilize the many search andindexing techniques known in the database community.
Our technique has the following advantages:

There are dozens of techniques in the literature forconverting shapes to time series [1][3][7][38][39][44],including some that are domain specific [5][31]. Our
approach works for any of these representations.

While there are many distance measures for shapes in theliterature, Euclidean distance, Dynamic Time Warping [2][5][30][31] and Longest Common Subsequence [37] accounts
for the majority of the literature. Our approach works for anyof these distance measures.

Our approach uses the idea of LB_Keogh lower bounding as
its cornerstone. Since the introduction of this idea a few
years ago [16], dozens of researchers world wide haveadopted and extended this framework for applications asdiverse as motion capture indexing [18], P2P searching [13],handwriting retrieval [31], dance indexing, and query byhumming and monitoring streams [40]. This widespread
adoption of LB_Keogh lower bounding has insured that it
has become a mature and widely supported technology, andsuggests that any contributions made here can be rapidly
adopted and expanded.
In some domains it may be useful to express rotation-limited
queries. For example, in order to robustly retrieve examplesof the number “8”, without retrieving infinity symbols “”,
we can issue a query such as: “Find the best match to thisshape allowing a maximum rotation of  15 degrees”. Our
framework supports such rotation-limited queries.

The literature on shape matching is vast; we refer the interested
reader to [7][36] and [44] for excellent surveys. While not all
work on shape matching uses a 1D representation of the 2D
shapes, an increasingly large majority of the literature does. We
therefore only consider such approaches here. Note that we lose
little by this omission. The two most popular measures that
operate directly in the image space, the Chamfer [6] andHausdorff [27] distance measures, require O(n2logn) time1 andrecent experiments (including some in this work) suggest that 1D
representations can achieve comparable or superior accuracy.
In essence there are three major techniques for dealing withrotation invariance, landmarking, rotation invariant features andbrute force rotation alignment. We consider each below.

2.1 Landmarking
The idea of “landmarking”is to find the one “true”rotation andonly use that particular alignment as the input to the distancemeasure. The idea comes in two flavors, domain dependent anddomain independent.
In domain dependent landmarking, we attempt to find a single (orvery few) fixed feature to use as a starting point for conversion of
the shape to a time series. For example, in face profile recognition
the most commonly used landmarks (fiducial points) are the chinor nose [5]. In limited domains this may be useful, but it requires
building special purpose feature extractors. For example, even in adomain as intuitively well understood as human profiles,
accurately locating the nose is a non-trivial problem, even if we
discount the possibility of mustaches and glasses. Probably the
only reason any progress has been made in this area is that mostwork reasonably assumes that faces presented in an image arelikely to be upright. For shape matching in skulls, the canonicallandmark is called the Frankfurt Horizontal [41], which is definedby the right and left porion (the highest point on the margin of theexternal auditory meatus) and the left orbitale (the lowest point on
the orbital margin). However, a skull can be missing the relevantbones to determine this orientation and still have enough global
information to match its shape to similar examples. Indeed thefamous Skhul V skull shown in Figure 14 is such an example.
Other examples of domain dependent landmarking include [39]who use the “sharpest corner”of leafs as landmarks. This idea
appears meaningful in the subset of leaf shapes they considered,but in orbicular (circular) leafs the “sharpest corner”is not welldefined.
In domain independent landmarking, we align all the shapes tosome cardinal orientation, typically the major axis. This approach
may be useful for the limited domains in which there is a welldefined major axis, perhaps the indexing of hand tools. However
1More precisely the time complexity is O(Rplogp), where p is the numberof pixels in the perimeter and R is the number of rotations that need to
be executed. Here p = n, and while R is a user defined parameter, it
should be approximately equal n to guarantee all rotations (up to thelimit of rasterization) are considered.

there is increasing recognition that the “… major axis is sensitive
to noise and unreliable”[44]. For example a recent paper showsthat under some circumstances, a single extra pixel can change therotation by  90 degrees [45].
To show how brittle landmarking can be we performed a simple
clustering experiment where we clustered three primate skulls
using Euclidean distance with both the major axis technique, andthe minimum distance of all possible rotations (as found by brute
force). Figure 3 shows the result. It is clear that the major axes do
not have any biological meaning: the points connecting each axis
for each specimen are not homologous (of shared evolutionary
origin). Therefore, the resulting cluster is meaningless in terms of
biology and morphology [41].
Best Rotation Alignment
Owl Monkey(species unknown)A

Owl Monkey
Northern Gray-Necked

B

OrangutanC

Landmark Alignment

Landmark Alignment

Best Rotation Alignment

A
B

A
B

Figure 3: Top) Three primate skulls, two of them from the same
genus, are clustered using both the landmark rotation beginning at the
major axis, and the best rotation. Bottom) The landmark-basedalignment of A and B explains why the landmark-based clustering isincorrect: a small amount of rotation error results in a large differencein the distance measure

Note that testing all rotations does not guarantee homology.
However by considering all possible alignments, we would expect(under the principle of parsimony) to achieve an approximately
homologous alignment. However, we note that what we aremaking claims here about subjective similarity, notmorphogenesis [41].
The most important lesson we learned from this experiment (anddozens of other similar experiments on diverse domains [14]) isthat rotation (mis)alignment is the most important invariance forshape matching, unless we have the best rotation then nothing elsematters.

2.2 Rotation Invariant Features
A large number of papers achieve fast rotation invariant matchingby extracting only rotation invariant features and indexing themwith a feature vector [7]. This feature vector is often called theshapes “signature”. There are literally dozens of rotation invariant
features including ratio of perimeter to area, fractal measures,
elongatedness, circularity, min/max/mean curvature, entropy,
perimeter of convex hull etc. In addition many researchers haveattempted to frame the shape-matching problem as a more familiar
histogram-matching problem. For example in [28] the authorsbuild a histogram containing the distances between two randomly
chosen points on the perimeter of the shapes in question. The

approach seems to be attractive, for example it can trivially also
handle 3D shapes, however it suffers from extremely poor
precision. For example, it cannot differentiate between the shapesof the lowercase letters “d”and “b”, or “p”and “q”, since thesepairs of shapes have identical histograms. In general, all thesemethods suffer from very poor discrimination ability [7]. In
retrospect this is hardly surprising. In order to achieve rotation
invariance, all information that contains rotation information must
be discarded; inevitably, some useful information may also bediscarded in this process. Our experience with these methods
suggests that they can be useful for making quick coarse
discriminations, for example differentiating between skulls andvertebrae. However we could not get these methods to distinguishbetween the skulls of humans and orangutan, a trivial problem forhuman or the brute force algorithm discussed in the next section.

2.3 Brute Force Rotation Alignment
There are a handful of papers that recognize that the above
attempts at approximating rotation invariance are unsatisfactoryfor most domains, and they achieve true rotation invariance by
exhaustive brute force search over all possible rotations, but onlyat the expense of computational efficiency and indexability[1][2][3][9][23][39]. For example, paper [1] uses DTW to handle
nonrigid shapes in the time series domain, while they note thatmost invariances are trivial to handle in this representation, theystate “rotation invariance can (only) be obtained by checking all
possible circular shifts for the optimal diagonal path.”This step
makes the comparison of two shapes O(n3) and forces them toabandon hope of indexing. Similarly paper [39] notes “In order tofind the best matching result, we have to shift one curve n times,
where n is the number of possible start points.”
In [23] the authors discretize the shapes into chain codes, andintroduce a fast dynamic programming method to test all
rotations. They note “The algorithm runs in O(nnlgn) time, where
n is the length of the compared strings.”. Of course this is anexact bound, but we achieve an empirical O(n1.06) on large
datasets (cf. Section 5). Given that most boundaries have about
1,000 datapoints long, this suggests that we are thousands of times
faster while also able to avoid discretization errors and avoiding
the need to set several parameters.
Dozens of papers have suggested that shape matching can bemade faster by sampling the contours. For example, in [19] theauthors note: “it is first necessary to reduce the number of data
points on the contour to a reasonable number that can beevaluated using shape similarity measurement.”These authors areinterested in classifying fish. The fish shapes are reduced down tomere forty data points because they “… found that a reduced data
set of 40 points was sufficient to retain the important shapefeatures for comparison”[19]. This dramatic data reduction did
make the similarity measure more tractable, but we wondered if
the assumption that it “retain(s) the important shape features”
was true. We compared their results, which after considerable
parameter tuning claimed “the highest recognition accuracy of64%”, with rotation invariant Euclidean distance on the raw data.
Surprisingly this simple, parameter-free method achieves 88.57%
accuracy (cf Section 5), which is much greater than the sampling
approach.
Other techniques introduced mitigate the untenable computational
complexity of testing “all” rotations do so at the expense of
introducing false dismissals. Typically they offer some implicit or
explicit trick to find a one (or a small number of) of startingpoint(s) [2][3][9]. For example paper [2] suggests “In order toavoid evaluation of the dissimilarity measure for every possible
pair of starting contour points … we propose to extract a small setof the most likely starting points for each shape.”Furthermore,both the heuristic used and the number of starting points must “beadjusted to a given application”, and it is not obvious how to best
achieve this.
In forceful experiments on publicly available datasets it has been
demonstrated that brute force rotation alignment produces the best
precision/recall and accuracy in diverse domains [1][2]. In
retrospect this is not too surprising. The rival techniques withrotation invariant features are all using some lossy transformationof the data. In contrast the brute force rotation alignment
techniques are using a (potentially) lossless transformation of thedata. With more high quality information to use, any distancemeasures will have an easer time reflecting the true similarity ofthe original images.
The contribution of this work is to speed up these accurate but
slow methods by many orders of magnitude while producing
identical results.

matching community [38]. However this technique does not help
reduce disk accesses for data which does not fit in main memory,and only allows matching under the Euclidean metric.

3. ROTATION INVARIANT MATCHING
We begin by formally defining the rotation invariant matchingproblem. We begin by assuming the Euclidean distance, andgeneralize to other distance measures later. For clarity ofpresentation we will generally refer to “time series”, which thereader will note can be mapped back to the original shapes.
Suppose we have two time series, Q and C of length n, whichwere extracted from shapes by an arbitrary method.
Q = q1,q2,… ,qi,… ,qn
C = c1,c2,… ,cj,… ,cn
As we are interested in large data collections we denote a databaseof m such time series as Q .
Q = {Q1, Q2, ...Qm}
If we wish to compare two time series, and therefore shapes, we
can use the ubiquitous Euclidean distance:

ED Q , C  

2.4 Indexing Star Light Curves
While this paper is focused on the indexing of shapes, it has cometo our attention that our techniques are ideally suited to theindexing of an important type of astronomical data known as starlight curves. We would be remiss not to make this connection
clear, so we briefly discuss the application and provide someexperimental results below.
Globally there are myriads of telescopes covering the entire sky
and constantly recording massive amounts of valuableastronomical data. Having humans to supervise all observations ispractically impossible; hence the increasing interest in computer
aided astronomy. A star light curve is a time series of brightnessof a celestial object as a function of time. The study of light
curves in astronomy is associated with the study of variability ofsources. That led to the discoveries of pulsars, extra solar planets,
supernovae, the rate of expansion of the universe just to namefew. At the Time Series Center at Harvard University Initiative in
Innovative Computing there are more than a 100 million such
curves (with billions more expected by 2009) however none ofthis data is currently searchable (other than by brute force search).

 qi  ci 
n

2

i 1

When using Euclidean distance as a subroutine in a classification
or indexing algorithm, we may be interested in knowing the exact
distance only when it is eventually going to be less than somethreshold r. For example, this threshold can be the “range” inrange search or the “best-so-far”in nearest neighbor search. If thisis the case, we can potentially speed up the calculation by doing
early abandoning [17].
Definition 1. Early Abandon: During the computation of the
Euclidean distance, if we note that the current sum of thesquared differences between each pair of corresponding data
points exceeds r2, then we can stop the calculation, secure in
the knowledge that the exact Euclidean distance had we
calculated it, would exceed r.
While the idea of early abandoning is fairly obvious and intuitive,
it is so important to our work we illustrate it in Figure 5 andprovide pseudocode in Table 1.

OGLE052401.70-Cc a lc u la tio na b a nd o ne da t th is p o in t

Q

0

20

40

OGLE052357.0260
80

0

100

Figure 4: An examples of two similar star light curves.

There is a need to compare the similarity of light curves for basic
astronomical research, for example in [29] researchers discoverunusual light curves worthy of further examination by finding theexamples with the least similarity to other objects. There are twothings which make this difficult. First is the enormous volume ofdata, the second is the fact that while is it possible to extract asingle period of a light curve, there is no natural starting point. In
order the find the similarity of two light curves it is therefore
necessary to compare every possible circular shift of the data [29],which as we show below corresponds exactly to the rotation
invariance matching problem for shapes in the one-dimensional
representation. The astronomical community [29] has mitigatedsome of the CPU effort for circular-shift matching byrediscovering the convolution “trick” long known to the shape

10

20

30

40

50

60

70

80

90

100

Figure 5: A visual intuition of early abandoning. Once the squared
sum of the accumulated gray hatch lines exceeds r2, we can be surethe full Euclidean distance exceeds r

Note that the “num_steps” value returned by the optimized
Euclidean distance in Table 1 is used only to tell us how useful
the optimization was. If its value is significantly less than n thissuggests dramatic speedup.
Table 1: Euclidean
abandonment

distance

optimized

with

early

algorithm [dist, num_steps] = EA_Euclidean_Dist(Q, C, r )accumulator = 0
for i = 1 to length(Q )// Loop over time series

accumulator += (qi - ci)2// Accumulate error contribution

If accumulator > r 2// Can we abandon?

time complexity for this algorithm is O(mn2). This is simply
untenable for large datasets.

disp(‘doing an early abandon’)num_steps = ireturn [ infinity, num_steps ]end// Terminate and return an

Table 3: An algorithm to find the best rotated match toquery from a database of possible matches// infinite error to signal theend// early abandonment.
return [ sqrt(accumulator), length(Q ) ] // Terminate with true distalgorithm: [best_match_loc, bestSoFar]= Search_Database_for_Rotated_Match(C, Q )While the Euclidean distance is a simple distance measure it
produces surprisingly good results for clustering, classificationand query by content of shapes, if the time series in question
happen to be rotation aligned. For example, in an experiment in[30] we manually performed rotation alignment of the time series
extracted from face profiles by explicitly showing the algorithm
the beginning and endpoint of a face (the nape and Adams apple
respectively).
However if the shapes are not rotation aligned, this method can
produce extremely poor results. Recall the results in Figure 1,
where a few degrees of misalignment give objectively andsubjectively incorrect clusterings. To overcome this problem weneed to hold one shape fixed, rotate the other, and record theminimum distance of all possible rotations.
For reasons that will become apparent later, we achieve this byexpanding one time series into a matrix C of size n by n.

c1 , c2 ,, cn1 , cn C  c2 ,, cn1 , cn , c1 





cn , c1 , c2 , , cn1 
Note that each row of the matrix is simply a time series, shifted(rotated) by one from its neighbors. It will be useful below toaddress the time series in each row individually, so we will denotethe ith row as Ci, which allows us to denote the matrix above in the
more compact form of C = {C1, C2,… , Cn}.
We can now define the Rotation invariant Euclidean Distance(RED) as:RED(Q, C) =min  ED Q , C j  
1 j  n  qi  ci 
n

2

i 1

Table 2 shows the pseudocode to calculate this.
Table 2: An algorithm to find the rotated match between twotime series
algorithm: [bestSoFar] = Test_All_Rotations(Q,C,r)bestSoFar = r
for j = 1 to n
distance = EA_Euclidean_Dist(Q, Cj, bestSoFar) // As in Table 1
if distance < bestSoFar
bestSoFar = distance;
end;
end;
return[bestSoFar]

Note that the algorithm tries to take advantage of early
abandoning by passing EA_Euclidean_Dist the value of r, the bestrotation alignment discovered thus far.
If we are simply measuring the distance between two time series
then the algorithm is invoked with r set to infinity, however, as weshall see below, if the algorithm is being used as a subroutine in a
linear scan of a large dataset Q , the calling routine can set thevalue of r to achieve speedup. In particular the calling function
sets r to the value of the best match (under any rotation)
discovered thus far. Table 3 shows the pseudocode. Note that thebest_match_loc = null
bestSoFar= inf
for i = 1 to number_of_time_series_in_database( Q )distance = Test_All_Rotations( Q ,C, bestSoFar);i// As in Table 2

if distance < bestSoFar
best_match_loc = i
bestSoFar = distance
end;
end;
return[best_match_loc, bestSoFar]

We will review the notation introduced thus far in Table 4.
Table 4: Notation Table

C

C
Ci
Q
Q

A time series
c1,c2,… ,cj,… ,cn
A n by n matrix containing every rotation of C
The ith row of the above
Another time series q1,q2,… ,qi,… ,qn
A database containing many time series Q = {Q1,..,Qm}

Note that our notation seems somewhat space inefficient in that it
expands time series C, of length n, to a matrix of size n by n.
However the rest of the database uses the original (arbitraryrotation) time series, and since the size of the database is assumed
to be large, this overhead is asymptotically irrelevant.
There are two simple and useful generalizations of definitions
thus far.
Mirror Image Invariance: Depending on the application we may
wish to retrieve shapes that are enantiomorphic (mirror images)
to the query. For example, in matching skulls, the best match
may simply be facing the opposite direction. In contrast whenmatching letters we don’t want to match a “d” to a “b”. If
enantiomorphic invariance is required we can trivially achieve
this by augmenting matrix C to contain Ci and reverse(Ci ) for 1
 i  n.
Rotation-Limited Invariance: In some domains it may be useful
to express rotation-limited queries. For example, in order torobustly retrieve examples of the number “6”, without retrieving
examples of the number “9”, we can issue a query such as:“Find the best match to this shape allowing a maximum rotationof  15 degrees”. Our framework trivially supports suchrotation-limited queries, by removing from the matrix C all timeseries that correspond to the unwanted rotations.
Thus far we have shown a brute force search algorithm that can
support rotation invariance, rotation-limited invariance and/or
mirror image invariance. We simply put the appropriate timeseries into matrix C and invoke the algorithm in Table 3. Thisalgorithm, even though speeded up by the early abandoning
optimization, is too slow for large datasets. In the next section we
introduce our novel search mechanism.

4. WEDGE BASED ROTATION MATCHING
We will begin by showing how we can efficiently search for thebest match in main memory. Since large datasets may not fit on
disk we will further show how we can index the data.

4.1 Fast and Exact Main Memory Search
We begin by defining time series wedges. Imagine that we takeseveral time series, C1,..,Ck , from our matrix C. We can use thesesequences to form two new sequences U and L:
Ui = max(C1i,..,Cki )Li = min(C1i,..,Cki )U and L stand for Upper and Lower respectively. We can see why
in Figure 6. They form the smallest possible bounding envelopethat encloses all members of the set C1,..,Ck from above andbelow. More formally:
i Ui  C1i,..,Cki  Li
For notational convenience, we will call the combination of U and
L a wedge, and denote a wedge as W:
W = {U, L}
C1

U

W

L
U

W

Q
L
10

20

30

40

50

60

70

80

90

10 0

Figure 6: Top) Two time series C1 and C2. Middle) A time series wedge
W, created from C1 and C2. Bottom) An illustration of LB_Keogh

We can now define a lower bounding measure between an
arbitrary time series Q and the entire set of candidate sequences
contained in a wedge W:

 (qi  U i ) 2 if qi  U in
LB _ Keogh (Q,W )    (qi  Li ) 2 if qi  Lii 1 
 0 otherwise
Below we show a proof of this lower bounding property. A
similar proof appears in [14] and also in [21], where the authors
use this representation for different problem.
Proposition 1: For any sequence Q of length n and a wedge W
containing a set of time series C1,..,Ck of the same length n, thefollowing inequality holds:

LB _ Keogh(Q,W )  ED (Q, C s ), where s = 1, 2, ..., k.
Proof:
Suppose we know that among the k time series C1,..,Ck , Cs has
the minimal Euclidean distance to query Q. And we wish to prove (q i  U i ) 2 if q i  U i2
 ( q i  Li ) if q i  Li i 1 
 0 otherwisen (q  C )
n

i





Below we will show that every term in the left summation can bematched with some greater or equal term in the right summation.
There are three cases to consider, for the moment we will justconsider the case when qi  U i . We want to show:(qi  U i ) 2  (qi  C si ) 2(q i  U i )  (qi  C si )
U i  C si
C si  U i

q  UiSince i, we can takesquare roots on both sides
q
Subtract i from both sides
U  C si
Add i
to both sides
BydefinitionUi=max(C1i,..,Cki )C si  max(C1i ,..., C ki )This is obviously true.
q  Li
The case iyields to a similar argument. The final case is0  ( qi  C si ) 2(q  C si ) 2
simple to show, since clearlybecause imust be nonnegative.
Thus we have shown that each term on the left side is matchedwith an equal or larger term on the right side. Our inequality
holds. ■C20 ( q i  U i ) 2 if q i  U in
2(q i  C si ) 2
 ( q i  Li ) if q i  Li i 1 
i 1
 0 otherwise
n

2

si

i 1

Since the terms under radicals are positive, we can square bothsides:

Note that the LB_Keogh function has been used before to support
DTW [16][30][31][37], uniform scaling [18], and query filtering[40]. For these tasks the lower bounding distance function is thesame, but the definition of U and L are different.
There are two important observations about LB_Keogh. First, inthe special case where W is created from a single candidate
sequence, it degenerates to the Euclidean distance. Second, notonly does LB_ Keogh lower bound all the candidate sequences
C1,..,Ck, but we can also do early abandon with LB_Keogh.
While the latter fact might be obvious, for clarity we make itexplicit in Table 5.
Table 5: LB_Keogh optimized with early abandonment
algorithm [dist, num_steps] = EA_LB_Keogh(Q, W, r )accumulator = 0
for i = 1 to length(Q )// Loop over time series
if qi > W.Ui// Accumulate error contribution
accumulator += (ci - W.Ui )2
elseif qi < W.Li
accumulator += (ci - W.Li )2
end
if accumulator > r 2// Can we abandon?
return [ infinity, i ]// Terminate and return an infinite error
end// to signal the early abandonment.
end
return [ sqrt(accumulator), length(Q ) ] // Terminate with true distNote once again that the value returned in “num_steps”is merely
a bookkeeping device to allow a post mortem evaluation ofefficiency.
Suppose we have just two time series C1 and C2 of length n, andwe know that in future we will be given a time series query Q and
asked if one (or both) of C1 and C2 are within r of the query. We
naturally wish to minimize the number of steps we must perform(“steps”are measured by “num_steps”). We are now in a position
to outline two possible approaches to this problem.



We can simply compare the two sequences, C1 and C2 (ineither order) to the query using the early abandon algorithm
introduce in Table 1. We will call this algorithm, classic.

We can combine the two candidate sequences into a wedge,
and compare Q to the wedge using LB_Keogh. If the
LB_Keogh function early abandons, we are done. We cansay with absolute certainty that neither of the two candidate
sequences is within r of the query. If we cannot early
abandon on the wedge, we need to individually compare thetwo candidate sequences, C1 and C2 (in either order) to thequery. We will call this algorithm, Merge.
Let us consider the best and worst cases for each approach. For
classic the worst case is if both candidate sequences are within rof the query, which will require 2n steps. In the best case, the firstpoint in the query may be radically different to the first point ineither of the candidates, allowing immediate early abandonmentand giving a total cost of 2 steps.
For Merge, the worst case is also if both candidate sequences arewithin r of the query, because we will waste n steps in the lower
bounding test between the query and the wedge, and then n steps
for each individual candidate, for a total of 3n. However the best
case, also if the first point in the query is radically different,
would allow us to abandon with a total cost of 1 step.
Which of the two approaches is better depends on:The shapes of C1 and C2. If they are similar, this greatly
favors Merge.
The shape of Q. If Q is truly similar to one (or both) of thecandidate sequences, this would greatly favor classic.
The matching distance r. Here the effect is non monotonic
and dependent on the two factors above.
We can generalize the notion of wedges by hierarchically nesting
them. Let us begin by augmenting the notation of a wedge toinclude information about the sequences used to form it. For
example, if a wedge is built from C1 and C2, we will denote it asW(1,2). Note that a single sequence is a special case of a wedge, forexample the sequence C1 can also be denoted as W1. We can
combine W(1,2) and W3 into a single wedge by finding maximum
and minimum values for each ith location, from either wedge.
More concretely:Ui = max(W(1,2)i, W3i )Li = min(W(1,2)i, W3i )W((1,2),3) = {U, L}
In Figure 7 we illustrate this notation. We call W(1,2) and W3children of wedge W((1,2),3). Since individual sequences are special
cases of wedges, we can also call C1 and C2 children of W(1,2).

C1 (or W1 )C2 (or W2 )C3 (or W3 )W(1, 2)W((1, 2), 3)

Figure 7: An illustration of hierarchically nested wedges

Given the generalization to hierarchal wedges, we can now also
generalize the Merge approach. Suppose we have a time series Q
and a wedge W((1,2),3). We can compare the query to the wedge
using LB_Keogh. If the LB_Keogh function early abandons, weare done. We know with certainty that none of the three candidate
sequences is within r of Q. If we cannot early abandon on thewedge, we need to compare the two child wedges, W(1,2) and W3 tothe query. Again, if we cannot early abandon on the wedge W(1,2),
we need to individually compare the two candidate sequences, C1
and C2 (in either order) to the query. We call this algorithm HMerge (Hierarchal Merge).
The utility of a wedge is strongly correlated to its area. We can getsome intuition as to why by visually comparing LB_Keogh(Q,W(1,2)) with LB_Keogh(Q, W((1,2),3)) as shown in Figure 8. Notethat the area of W((1,2),3) is much greater than that of W(1,2), and that
this reduces the value returned by the lower bound function andthus the possibility to early abandon.
W (1,2)Q

W ((1,2),3)

Q

Figure 8: Top) An illustration of LB_Keogh(Q, W(1,2)). Bottom) An
illustration of LB_Keogh(Q, W((1,2),3)). Note that the tightness of thelower bound is proportion to the number and length of vertical lines

For some problems, the H-Merge algorithm can give
exceptionally poor performance. If the wedge W(1,2), created fromC1 and C2 has an exceptional large area (i.e. C1 and C2 are very
dissimilar), it is very unlikely to be able to prune off any steps.
At this point we can see that the efficiency of H-Merge isdependent on the candidate sequences and Q itself. In general,
merging similar sequences into a hierarchal wedge is a good idea,
but merging dissimilar sequences is a bad idea.
The observations above motivate a final generalization of H-Merge.
Recall that to achieve rotation invariance we expanded our timeseries C into a matrix with n time series. Given these n sequences,
we can merge them into K hierarchal wedges, where 1  K  n.
This merging forms a partitioning of the data, with each sequencebelonging to exactly one wedge. We will use W to denote a set ofhierarchal wedges:
W = {Wset(1) , Wset(2) ,.., Wset(K)} ,1K nwhere Wset(i) is a (hierarchically nested) subset of the n candidate
sequences. Note that we haveWset(i)  Wset(j) =  if i  j, andC3 (or W3)

C5 (or W5)| Wset(1) ∪Wset(2) ∪..∪ Wset(K) | = n

C2 (or W2)

We will attempt to merge together only similar sequences. We canthen compare this set of wedges against our query. Table 6
formalizes the algorithm.

C4 (or W4)

Table 6: Algorithm H-Merge
algorithm [dist] = H-Merge(Q, W,K, r )C1 (or W1)

S = {empty }// Initialize a stack.

for i = 1 to K// Place all the wedges into the stack.

Figure 9: A dendrogram of five sequences C1, C2,..., C5, clusteredusing group average linkage

enqueue(Wset(i) ,S )end
while not empty(S )Given that the clustering algorithm produces the tentative wedge
sets, all we need to do is to choose the best one. We could attempt
to do this by eye, for example in Figure 10 it is clear that any
sequence that early abandons on W3, will almost certainly also
early abandon on both W2 and W5; similar remarks apply to W1
and W4. At the other extreme, the wedge at K = 1 is so “fat”that it
is likely have poor pruning power. The set W = {W((2,5),3), W(1,4)}
is probably the best compromise. However because the set of timeseries might be very large, such visual inspection is not scalable.

T = dequeue(S )dist = EA_LB_Keogh(Q,T,r )// Note that is early abandon version.

if isfinite(dist)// We did not early abandon.

if cardinality(T ) = 1// T was an individual sequence.

disp(‘The sequence ’,T, ‘is ’, dist, ‘ units from the query’)return[dist]else// T was a wedge, find its children

enqueue(children(T ) ,S )// and push them onto the stack.

end
end
W3

end

Note

that

this

algorithm

isdesignedtoreplace

W3

W3

the

Test_All_Rotations algorithm that is invoked as a subroutine in theW((2,5),3)W2
W(2,5)

Search_Database_for_Rotated_Match algorithm shown in Table 3.

As we shall see in our empirical evaluations, H-Merge canproduce very impressive speedup if we make judicious choices inthe set of hierarchal wedges that make up W. However, thenumber of possible ways to arrange the hierarchal wedges isgreater than KK, and the vast majority of these arrangements will
be very poor, so specifying a good arrangement of W is critical.
A simple observation alleviates the need to invent a newalgorithm to find a good arrangement of W. Note that hierarchal
clustering algorithms have very similar goals to an ideal wedgeproducing algorithm. In particular, hierarchal clustering
algorithms can be seen as attempting to minimize the distances
between objects in each subtree. A wedge-producing algorithm
should attempt to minimize the area of each wedge. However thearea of a wedge is simply the maximum Euclidean distance
between any sequences contained therein (i.e Newton-Cotes rulefrom elementary calculus). This motivates us to derive wedge setsbased on the result of a hierarchal clustering algorithm. Figure 10
shows wedge sets W, of every size from 1 to 5, derived from thedendrogram shown in Figure 9.

W(2,5)

W(((2,5),3), (1,4))

W5

W1

W1
W(1,4)
W(1,4)

W4

W4

K=5

K=4

K=3

K=2

K=1

Figure 10: Wedge sets W, of size 1 to 5, derived from the dendrogram
shown in Figure 9

The problem is actually even more complex, in that the best value
for K also depends on the current value of r (Recall r is the “bestso-far”in nearest neighbor search.). If r is large then very little
early abandoning is possible and this favors a large value for K. In
contrast, if r is small we can do a lot of early abandoning, and weare better off having many sequences in a single wedge so we canearly abandon all of them with a single calculation. Note howeverthat for nearest neighbor search the value of r will get smaller as
we search through the database.
With this in mind, we dynamically choose the wedge set based ona fast empirical test. We start with the wedge set where K = 2.
Each time the bestSoFar value changes, we test a subset of thepossible values of K and choose the most efficient one (asmeasured by num_steps) as the next K to use. Which subset to test

is decided on-the-fly based on the current K value. They are thevalues which evenly divide the ranges [1, current_K] and[current_K, max_K] into 5 intervals. Note that on average the
bestSoFar value only changes log(m) during a linear search, so thisslight overhead in adjusting the parameter is not too burdensome,however we do include this cost in all experiments in Section 5.

4.2 Lower Bounding in Index Space
True rotation invariance has traditionally been so demanding interms of CPU time that little or no effort was made to index it (orit was indexed with the possibility of false dismissals with regard
to the raw shapes [4]). As we shall see in the experiments in
Section 5.2, the ideas presented in the last section produce suchdramatic reductions in CPU time that it is worth consideringindexing the data.
There are several possible techniques we could consider forindexing. Recent years have seen dozens of papers on indexing
time series envelopes that we could attempt to leverage off[16][21][30][31][37]. The only non-trivial adaptation to be madeis that instead of the query being a single envelope, it would benecessary to search for the best match to K envelopes in thewedge set W.
Note however that we do not necessarily have to use theenveloping idea in the indexing phase. So long as we can lower
bound in the index space we can use an arbitrary technique to get(hopefully a small subset of) the data from disk to main memory[8], where our H-Merge can very efficiently find the distance tothe best rotation.

This technique is adapted from [38], and we refer the reader tothis work for a more complete treatment.

4.3 Generalizing to other Distance Measures
As we shall see in Section 5, the Euclidean distance is typicallyvery effective and intuitive as a distance measure for shapes.
However in some domains it may not produce the best possible
precision/recall or classification accuracy [2][30]. The problem isthat even after best rotation alignment, subjectively similar shapes
may produce time series that are globally similar but contain local“distortions”. These distortions may correspond to local featuresin that are present in both shapes but in different proportions. For
example in Figure 11 we can see that the larger brain case of the
Lowland Gorilla changes the locations in which the brow ridgeand jaw map to in a time series relative to the Mountain Gorilla.

Lowland Gorilla
Gorilla gorilla graueri

Table 7: A Vantage Point Tree for Indexing Shapes
Algorithm [BSF] = NNSearch(C)BSF.ID = null;// BSF is the Best-So-Far variable
BSF.distance = infinity;
W = convert_time_series_to_wedge_set(C);Search( Q ,W, BSF);// Invoke subroutine on the root of index tree
root
Subroutine Search(NODE, W, BSF)
if NODE.isLeaf// we are at a leaf node.
for each compressed time-series cT in node
LB = computeLowerBound(cT, W);
queue.push(cT,LB);// sorted by lower bound.
endwhile (not (queue.empty()) and (queue.top().LB < BSF.distance))
if (BSF.distance > queue.top().LB)
retrieve full time series Q of queue.top() from disk;
distance = H-Merge(Q, W, BSF.distance ) // calculate full distance.
if distance < BSF.distance// update the best-so-far
BSF.distance = distance;// distance and location.
BSF.ID = Q;
end
end
end
else
// we are at a vantage point.
LB = computeLowerBound(VP, W);
queue.push(VP,LB);
if LB < (node.median + BSF.distance)
search(NODE.left, W, BSF);// recursive search left.
else
search(NODE.right, W, BSF); // recursive search right.
end
end

One possible method to achieve this indexable lower bound is touse Fourier methods. Many authors have independently noted thattransforming the signal to the Fourier space and calculating the
Euclidean distance between the magnitude of the coefficients
produces a lower bounds to any rotation [4][38]. We can leverageof this lower bound to use a VP-tree to index our time series asshown in Table 7.

Mountain Gorilla
Gorilla gorilla beringei

Figure 11: The Lowland Gorilla and Mountain Gorilla aremorphologically similar, but have slightly different proportions.
Dynamic Time Warping can be used to align homologous features inthe time series representation space

Even if we assume that the database contains the actual object
used as a query, it is possible that the two time series are distorted
versions of each. Here the distortions may be caused by camera
perspective effect, differences in lighting causing shadows whichappear to be features, parallax etc.
Fortunately there is a well-known technique for compensatingsuch local misalignments, Dynamic Time Warping (DTW)[16][30]. While DTW was invented in the context of 1D speech
signals others have noted its utility for matching shapes, including
face profiles [5], hand gestures [25], leafs [30] and handwriting[31].
To align two sequences using DTW, an n-by-n matrix isconstructed, where the (ith, jth) element of the matrix is the
distance d(qi, cj) between the two points qi and cj (i.e. d(qi, cj) = (qi- cj)2 ). Each matrix element (i, j) corresponds to the alignment
between the points qi and cj, as illustrated in Figure 12.

A warping path P is a contiguous set of matrix elements thatdefines a mapping between Q and C. The tth element of P is
defined as pt = (i, j)t so we have:n ≤ T < 2n-1

P = p1, p2, … , pt, … , pT

The warping path that defines the alignment between the two timeseries is subject to several constraints. For example, the warping
path must start and finish in diagonally opposite corner cells of
the matrix; the steps in the warping path are restricted to adjacentcells (including diagonally adjacent cells); the points in thewarping path must be monotonically spaced in time. In addition tothese constraints, virtually all practitioners using DTW also
constrain the warping path in a global sense by limiting how far it
may stray from the diagonal [16][30][31]. A typical constraint is
the Sakoe-Chiba Band which states that the warping path cannot
deviate more than R cells from diagonal [34].

 ( qi  DTW _ U i ) 2 if qi  DTW _ U i

LB _ KeoghDTW (Q,W )    ( qi  DTW _ Li ) 2 if qi  DTW _ Lii 1 
0otherwisen

We will now prove the claim of the lower bounding.
Proposition 2: For any sequence Q of length n and a wedge W
containing a set of time series C1, … , Ck of the same length n, forany global constraint on the warping path of theform j  R  i  j  R , the following inequality holds:

LB _ KeoghDTW (Q,W )  DTW (Q, C s ), where s = 1, 2, ..., k.
Proof:
Suppose we know that among the k time series C1, … , Ck, Cs has
the minimal DTW distance to query Q. And we wish to prove ( q i  DTW _ U i ) 2 if q i  DTW _ U i2
 ( q i  DTW _ Li ) if q i  DTW _ Li i 1 
0otherwise



n

C

C
Q

Q

p
T

st

t 1

Since the terms under radicals are positive, we can square bothsides:Q ( q i  DTW _ U i ) 2 if q i  DTW _ U i2
 (q i  DTW _ Li ) if q i  DTW _ Lii 1 
0otherwise



n

C

R
Figure 12: Left) Two time series sequences with local differences. Right)
To align the sequences we construct a warping matrix, and search for the
optimal warping path, shown with solid squares. Note that Sakoe-Chiba
Band with width R is used to constrain the warping path

The optimal warping path can be found in O(nR) time by dynamic
programming [16].
Based on an arbitrary wedge W and the allowed warping range R,
we define two new sequences, DTW_U and DTW_L:
DTW_Ui = max(Ui-R : Ui+R )DTW_Li = min(Li-R : Li+R )They form an additional envelope above and below the wedge, as
illustrated in Figure 13.
DTW_U
C1

W
W

C2

A

C

DTW_L
DTW_U

W
W

U

L

st

t 1

Recall that that when we stated the definition of the warping pathabove we had, P = p1, p2, … , pt, … , pTn ≤ T < 2n-1. We
therefore have n ≤ T, so our strategy will be to show that every
term in the left summation can be matched with some greater or
equal term in the right summation.
There are three cases to consider, for the moment we will justconsider the case when qi > DTW_Ui . We want to show:(qi  DTW _ U i ) 2  p st
(q i  DTW _ U i ) 2  (qi  C sj ) 2

By Definition 3(q i  DTW _ U i )  (qi  C sj )Since qi > DTW_Ui , we
can take square roots on
both sides
Subtract qi from bothsides
Add DTW_Ui + Csj toboth sides
By definition DTW_Ui =max(Ui-R : Ui+R ) DTW _ U i  C sj
C sj  DTW _ U i
C sj  max(U i  R : U i  R )Since the query sequence Q and all the candidate sequences C1,… , Ck are of the same length and j-R ≤ i ≤ j+R, we know i-R ≤ j ≤i+R. So we can rewrite the right side and the inequality becomesWW

C sj  max(U i  R , U ( i 1)  R ,..., U j ,..., U i  R )

Q

B

p
T



D

DTW_L

Figure 13: The idea of bounding envelopes introduced in Figure 6 isgeneralized to allow DTW. A) Two time series C1 and C2. B) A timeseries wedge W, created from C1 and C2. C) In order to allow lower
bounding of DTW, an additional envelope is created above and below
the wedge. D) An illustration of LB _ KeoghDTW

We can now define a lower bounding measure for DTW distance
between an arbitrary query Q and the entire set of candidate
sequences contained in a wedge W:

If we remove all terms except Uj from the RHS we are left withC sj  max(U j )which is obviously true since Uj = max(C1j,..,Ckj).
The case qi < DTW_Ui yields to a similar argument. The final case0  (qi  C sj ) 2(q  C sj ) 2
is simple to show, since clearlybecause imust be nonnegative.
Thus we have shown that each term on the left side is matchedwith an equal or larger term on the right side. Our inequality
holds. ■For brevity we omit the very minor modifications required toindex LB_KeoghDTW(Q,W), however [37] contains the necessary
modifications for both DTW and for LCSS which is discussed
below.
To facilitate later efficiency comparisons to Euclidean distanceand other methods, it will be useful to define the time complexityof DTW in terms of “num_steps” as returned by Table 1 andTable 5. The variable “num_steps” is the number of real-value
subtractions that must be performed, and completely dominates
the CPU time, since the square root function is only performed
once (and can be removed, see [17]). If we construct a full n by n
warping matrix, then DTW clearly requires at least n2 steps.
However as we noted above and illustrated in Figure 12, we cantruncate the corners of the matrix to reduce this number toapproximately nR, where R is the width of the Sakoe-Chiba Band.
While nR is the number of steps for a single DTW, we expect the
average number of steps to be less, because some full DTW
calculations will not be needed if the lower bound test fails. Since
the lower bound test requires n steps, the average number of steps
when doing m comparisons should be:m  a ( nR )  m ( n )mWhere a is the fraction of the database that requires the full DTW
calculated. Note that even this is pessimistic, since both DTW2and LB_KeoghDTW are implemented as early abandoning (recall
Table 5). We therefore simply count the “num_steps”required byeach approach and divide it by m to get the average number ofsteps required for one comparison.
In addition to DTW, several researchers have suggested using
Longest Common SubSequence (LCSS) as a distance measure forshapes. The LCSS is very similar to DTW except that while DTW
insists that every point in C maps onto one (or more) point(s) in
Q, LCSS allows some points to go unmatched. The intuition
behind this idea in a time series domain is that subsequences may
contain additions or deletions, for example an extra (or forgotten)
dance move in a motion capture performance, or a missed beat in
ECG data. Rather than forcing DTW to produce an unnaturalalignment between two such sequences, we can use LCSS, whichsimply ignores parts of the time series that are too difficult tomatch. In the image space the missing section of the time series
may correspond to a partial occlusion of an object, or to aphysically missing part of the object, as shown in Figure 14.

A

B

C

This region will
not be matched

LCSS
Alignment

DTW

B

Figure 14: A) The famous Skhul V is generally reproduced with themissing bones extrapolated in epoxy, however the original Skhul V (B)is missing the nose region, which means it will match to a modern
human (C) poorly, even after DTW alignment (inset). In contrast, LCSSalignment will not attempt to match features that are outside a“matching envelope”(heavy gray line) created from the other sequence.

Real world examples of domains that require LCSS abound. For
example anthropologists are interested in exploring large datasetof projectile points (“arrowheads”). At the UCR Lithic
Technology Lab at UCR there are over a million specimens, so
indexing is required for efficient access. While anthropologists
have long been interested in shape, interest in matching suchobjects is further fueled by the availability of computing powerand by a recent movement that notes, “an increasing number ofarchaeologists are showing interest in employing Darwinianevolutionary theory to explain variation in the material record”[26]. Anthropologists have recently used tools from biologicalmorphology to attempt to explain spatial and temporal distributionof projectile points in North America. As we illustrate in Figure
15 many examples are incomplete, missing tip or tangs. LCSS canignore such missing feature to provide more robust matching.

Figure 15: Project points are frequently found with broken tips ortangs. Such objects require LCSS to find meaningful matches tocomplete specimens. From left to right, Edwards, Langtry, andGolondrina projectile points.
2

Note that a recursive implementation of DTW would always require nR
steps, however iterative implementation (as used here) can potentially
early abandon with as few as R steps.

While we considered LCSS for generality, we will not further
explain how to incorporate it into our framework. It has been

shown in [37] that it is trivial to lower bound LCSS using theenvelope-based techniques described above. The minor changes
include reversing some inequality signs since LCSS is a similarity
measure, not a distance measure. Unlike Euclidean distance which
has no parameters, or DTW, which has one intuitive and easy to
set parameter, LCSS requires 2 parameters, and tuning them is
nontrivial. In objective classification experiments we found that
we could sometimes tune LCSS to slightly beat DTW on someproblems, however we did not have large enough datasets to allow
training/test splits that guarded against overfitting to a statisticallysignificant standard. Automatically choosing the correct
parameters for LCSS is a matter for future research.

5. EXPERIMENTAL RESULTS
In this section we empirically evaluate our approach. We begin bystating our experimental philosophy. In a recent paper Veltkamp
and Latecki attempted to reproduce the accuracy claims of several
shape matching papers but discovered to their dismay that they
could not match the claimed accuracy for any approach [36]. One
suggested reason is the observation that many approaches havehighly tuned parameters, a fact which we believe makes Euclidean
distance (zero parameters) and DTW (one parameter) particularly
attractive. Veltkamp and Latecki conclude “It would be good forthe scientific community if the reported test results are madereproducible and verifiable by publishing data sets and software
along with the articles”. We completely concur and have placed
all datasets at the following URL [14].

5.1 Effectiveness of Shape Matching
In general this paper is not making any claims about theeffectiveness of shape matching. Because we are simply speeding
up arbitrary distance calculations on arbitrary 1-dimensional
representations of shapes, we automatically inherit the welldocumented effectiveness of other researchers published work
[1][2][3][9][12][30][38].
Nevertheless, for completeness and in order to justify the extra
computational expense of DTW, we will show the effectiveness ofshape matching on several publicly available datasets.
Table 8 shows the error rate of one-nearest neighbor classification
as measured using leaving-one-out evaluation. Recall that
Euclidean distance has no parameters, DTW has a single
parameter (the warping window width R) which was learned bylooking only at the training data.
Table 8: The Error of Euclidean distance and DTW on
publicly available datasets
Name

Number ofClassesFace

16

Numberof
Instances
2240

Euclidean
Error (%)

DTW
Error (%) {R}

3.839%

3.170% {3}Swedish Leaves

15

1125

13.33%

10.84% {2}Chicken5446

19.96%

19.96% {1}MixedBag

9160

4.375%

4.375% {1}OSU Leaves

6442

33.71%

15.61% {2}Diatoms

37

781

27.53%

27.53% {1}Aircraft

7

210

0.95%

0.0% {3}

Fish

7

350

11.43%

9.71% {1}Light-Curve

3954

14.15%

11.43% {3}

Yoga

2

3300

4.70%

4.85% {1}For the Face and Leaf datasets the (approximate) correct rotation
was known [30]. We removed this information by randomly
rotating the images.

The MixedBag dataset is small enough to run the morecomputationally expensive Chamfer [6] and Hausdorff [27]distance measures. They achieved an error rate of 6.0% and 7.0%
respectively [38], slightly worse than Euclidean distance.
Likewise the Chicken dataset allows us to compare directly to[24], which used identical experiments to test 6 different
algorithms based on discrete sequences extracted from the shapes.
The best of these algorithms had an error rate of 20.5% and took
over a minute for each distance calculation, whereas our approachtakes an average time of 0.0039 seconds for each distancecalculation3. For the Diatom dataset, the results are competitive
with human experts, whose error rates ranged from 57% to 13.5%[12], and only slightly worse than the Morphological Curvature
Scale Spaces (MCSS) approach of [12], which got 26.0%. Notehowever that the Euclidean distance requires zero parameters once
the time series have been extracted, whereas the MCSS has
several parameters to set. On the aircraft dataset the data donors
tried 4 different Hidden Markov Model based approaches andachieve a best error-rate of 0.95% [35]. However, we can achieve
this accuracy with zero-parameter Euclidean distance and get zeroerror with DTW. The data donors of the fish dataset tested many
shape descriptors, such as Fourier descriptors, polygon
approximation and line segments to achieve the best error rate of
36.0% [19], however we get dramatically lower error rates forboth Euclidean distance and DTW.
The yoga dataset is of particular interest. It has previously beenclassified without rotation invariance [15]. In that case humanvolunteers painstakingly located a landmark point in each of
3,300 images. The accuracy for the human annotated dataset was
17.0% and 15.5% for Euclidean distance and DTW respectively.
However the rotation invariant classification here has reduced theerror rate by a factor of three. Since the only thing to differ
between the two experiments is rotation invariance, this strongly
supports the contention made in Section 2.1, that “rotation(mis)alignment is the most important invariance for shape
matching, unless we have the best rotation then nothing elsematters”. Note that here DTW is very slightly worse on thedataset, but the difference is not statistically significant at the 0.01
level.

5.2 Shape Matching Sanity Checks
In general the experiments in the previous section show two
things (which had been noted before), the extra effort of DTW isuseful in some domains, and very simple time series
representations of shapes are completive to other more complex
representations.
We also performed extensive “sanity check” experiments indomains where we could meaningful visualize the results.
Our first example uses a large database of skulls. For all primate
species where we have at least two examples we perform ahierarchical clustering and check to see if both samples of thesame species clustered together. Figure 16 shows a typicalexample.

3

We are aware that one should normally not compare CPU times fromdifferent computers, however here the 4 orders of magnitude offers a
comfortable margin that dwarfs implementation details.

•Orangutan

Alligator
mississippiensis

Crocodylidae•Orangutan(juvenile)•Mantled Howler
Monkey•Red Howler
Monkey•Human Ancestor(Skhul V)•HumanCaimancrocodilus

Crocodylus
cataphractus

Alligatorinae

Tomistoma
schlegelii

Crocodylus
johnstoni•De Brazza monkey(juvenile)•De Brazza monkey
Figure 16: A group average hierarchal clustering of eight primate
skulls based on the lateral view, using Euclidean distance

It is important to recall that Figure 16 shows a phenogram, not aphylogenetic tree. However on larger scale experiments in thisdomain (shown in [14]) we found that large subtrees of thedendrograms did conform to the current consensus on primate
evolution.
While the Euclidean distance works very well on the relatively
simple primate skulls, we found that considering a more
(morphologically) diverse groups of animals, such as all reptiles,
requires DTW as a distance measure. Consider Figure 17 whichshows a hierarchical clustering of a very diverse set of reptiles. As
with the primates, this is not the correct phylogenetic tree forthese animals, once again however, the (uniquely colored)
subtrees do correspond to current consensus on reptiles evolution
based on DNA analysis and/or more complete morphologicalstudies [10][11].
Note that we are not claiming that our shape matching techniques
replace or even complement classic morphometrics in zoology.
The point of these experiments is that if the shape matching
techniques can produce intuitive results in a domain in which weknow the correct relationships by other means, this suggests thatalgorithms may also produce meaningful results in shapeproblems for which there is more uncertainty, including projectile
points (see [26] and Figure 15), petroglyphs, insect bite patternsin leaves [42], mammographic calcifications [43] etc.
It has recently been claimed that shape matching methods that
only look at the contours of shapes (boundary based methods) arebrittle to articulation distortion [33], however we believe thatwhile this may be true for certain boundary based methods (i.e
Hausdorff, Champer etc) the centroid based method we use is very
robust to articulation distortions. To demonstrate this, we
conducted a simple experiment/demonstration. We took three
Lepidoptera, including the very similar and closely related Actias
maenaes and Actias philippinica, and produced a copy of each.
We then took these copies and “bent” the right hindwing. The
clustering of the three originals and three copies under Euclidean
distance, group average linkage is shown in Figure 18.

Cricosaura
typicaXantusia
vigilisAmphisbaenia
Elseyadentata

Chelonia
Glyptemys
muhlenbergii

Phrynosoma
hernandesi

Phrynosoma
douglassii

Iguania

Phrynosoma
taurusPhrynosoma
ditmarsi

Phrynosoma
mcallii

0

200

400

600

800

1000 1200

Figure 17: A group average hierarchal clustering of fourteen reptile
skulls based on the superior view, using DTW distanceChorinea amazon

Actias philippinica

Actias maenaes

returned by Table 1 and Table 5 allows an implementation freemeasure to compare performance.
For Euclidean distance queries we compare to brute force andFourier (FFT) methods, which are the only competitors to also
guarantee no false dismissals. The cost model for the FFT lower
bound is nlogn steps. If the FFT lower bound fails we allow theapproach to avail of our early abandoning techniques discussed in
Section 3.
We tested on two shape datasets, a homogeneous database of
16,000 projectile point images, all of length 251 and aheterogeneous dataset consisting of all the data used in theclassification experiments, plus 1,000 projectile points. In total
the heterogeneous dataset contains 5,844 objects of length 1,024.
To measure the performance we averaged over 50 runs, with thequery object randomly chosen and removed from the dataset.
We measure the average number of steps required by eachapproach for a single comparison of two shapes, divided by the
number of steps require by brute force. For our method, we
include a startup cost of O(n2), which is the time require to build
the wedges. Because the utility of early abandoning depends onthe value of the best-so-far, we expect our method to do better as
we see larger and larger datasets.
Figure 19 shows the results on the projectile points dataset using
Euclidean distance.
Projectile Points
Euclidean

Figure 18: An experiment to demonstrate that the centroid method isreasonably articulation invariant. The gray highlighted areas havebeen randomly “tweaked”in a photo editing program

As we can see, the one dimensional representation has hardly
changed, and the clustering correctly groups the three pairs. We
found that using DTW we can even more radically distort theshapes and achieve similar results.
Given the above, why do boundary based methods have such apoor reputation for domains where articulation is a problem [33]?
We believe the answer is not intrinsic to boundary based methods,
but lies in the measures typically used on them, especially the
Hausdorff distance and its many variants. Consider the following
thought experiment. Imagine we have two identical shapes; solid
automobiles. Assume that they have identical antenna protrudingfrom their roofs. As such, the Hausdorff distance between them iszero, but if we bend the antenna in the spirit of Figure 18, we cantrivially increase the Hausdorff distance to one meter.
In addition to the above, the results in Table 8 already hinted at
the articulation invariance of our chosen representation. Many ofthe datasets considered have significant articulation. For example
the face dataset considers classes with both closed-lip andlaughing/yawning people, the two leaf dataset have significant
amounts of articulation at the stem, and the Yoga dataset features
very flexible people in version poses.
Finally, we note that paper [25] uses the ideas in the conference
version of this work to index hand geometries for biometrics. It isclear that the human hand has a high degree of articulation.

5.3 Main Memory Experiments
There is increasing awareness that comparing two competing
approaches using only CPU time opens the possibility of
implementation bias [17]. As a simple example, while the Haar
wavelet transform is O(n) and DFT is O(nlogn), the DFT is much
faster in the popular language Matlab, simply because it is ahighly optimized subroutine. For this reason many recent papers
compare approaches with some implementation-free metric[16][30][37][38]. As we noted earlier, the variable “num_steps”

1.0

0.8

0.6

0.4

0.2

3264

125 250500 000
1
Num2000 4000
ber of
8000 6000
objec
1
ts in d
ataba
se (m
)

0

Brute
force
F
Early FT
Wedge abandon

Figure 19: The relative performance of four algorithms on the
Projectile Points dataset using the Euclidean distance measure

We can see that for small datasets our approach is slightly worse
than FFT and simple Early abandon because we had to spendsome time building the wedges. However, by the time we haveseen 64 objects we have already broken even, and thereafter
rapidly race towards beating FFT and Early abandon by one
order of magnitude and Brute force by two orders of magnitude.
The results on the projectile points dataset using DTW are shownin Figure 20, and are even more dramatic.

Projectile Points
DTW

0.12

m = 16000
0.1

The Time Series Center at Harvard University Initiative in
Innovative Computing has collected light curves from various
completed surveys (total 100 million examples) and is collecting
data from an on-going survey (TAOS) which is producing few
billion additional examples. While light curves have beencollected for decades it is only recent years in which we havebegun to see serious efforts to index and data mine them. For
consistency we only consider the indexing of the small handlabelled set of examples use in the classification experiments inTable 8. Figure 22 shows the performance of the four rival
methods on the light curves under the Euclidean distance.

1.0

0.08

0.8

0.04

0.02

Brute Force, R =5

Wedge
Early abandon

0.06

0.6

0.4

0

Inset
0.2

3264
125 2500Light Curves- Euclidean

Brute
forc
B
Early rute Force, e
R =5
Wed ge abandon

500 000
1
Numb
2000 4000 00
er of
80
0
objec
1600
ts in d
ataba
se (m
)

1.0

Figure 20: The relative performance of four algorithms on the
Projectile Points dataset using the DTW distance measure. The inset
shows a zoom-in of the 3 best algorithms when m = 16,000

Here the cost of building the wedges is dwarfed by a single brute
force DTW-rotation-invariant comparison, so our approach isfaster even for a database of size 3. By the time we have examinedthe entire database, our approach is more than 5,000 times fasterthan the brute force approach. It is interesting to note that theearly abandoning strategy is by itself quite competitive, yet to ourknowledge no one uses it. We suspect this is because most people
are more familiar with the elegant and terse recursive version of
DTW, which does not allow early abandoning, than the iterative
implementation, which does. Note however that even though ourhighly optimized early abandoning strategy is competitive, ourwedge approach is still an order of magnitude faster once thedataset is larger than 500 objects.
Sometimes indexing methods that work well for highly
homogeneous datasets do not work well for heterogenous
datasets, and vice versa. We consider this possibility by testing on
the heterogenous dataset in Figure 21.
Heterogeneous Dataset

Num
ber

Euclidean

of ob
je

cts in

datab
as e

DTW
(m)

1.0
0.8

Brute force

0.6
0.4
0.2

3264

125 250

3264

500

Numb
er

1000 2000

of obje

EarlyFFT
Wedge abandon

4000 8000

cts indataba

125 250

se (m

500

1000 200 04000 8000

0.8

0.6
0.4

0.2

32
Numb

0

64
er of 125 250
objec
ts in d 500 53
9
ataba
se (m
)

Brute
force
FF
Early T
aband
on
Wedg
e

Figure 22: The relative performance of four algorithms on the light
curve dataset using Euclidean distance

As before, our approach is slightly slower on small datasets due tothe setup overhead. However once the dataset has more than 125
objects the wedge-based approach is slightly faster, and by thetime we see the full dataset it is and order of magnitude better
than the FFT approach.
Recall that in the classification experiments shown in Table 8, theclassification of light curves is significantly faster with the DTW
distance. Figure 23 shows the performance of four rival methods
on the light curves under the Euclidean distance.
Light Curves- DTW

Bru0tefo
EarlyBrute For rce
Wedg abandonce, R =5
e

1.0

)

0.8

Figure 21: The relative performance of four algorithms on theHeterogeneous dataset using Euclidean distance (left) and DTW (right)

In this dataset it takes our wedge approach slightly longer to beat
Early abandon (and FFT for Euclidean search), however by thetime we have seen 8,000 objects our approach is two orders ofmagnitude faster than its Euclidean competitors, and for DTW it
is an order of magnitude faster than Early abandon and 3,976
times faster than brute force.
Recall that our algorithm requires the setting of a single
parameter, the number of intervals to search for a new value for K
every time the bestSoFar variable is updated. In all the experiments
above this value was set to 5. We found that we can change thisvalue to any number in the range 3 to 20 without affecting theperformance of our algorithm by more than 4%, we therefore omit
further discussion of this parameter setting.

0.6
0.4
0.2

32
Numb 64 1250
er of
objec 25 500
3ts in databa 95
se (m
)

0

Brute
fo rc
B
Early rute force, e
a
R =5
bando
Wedg
n
e

Figure 23: The relative performance of four algorithms on the light
curve dataset using DTW distanceAs in the shape dataset, our method is several orders of magnitude
faster.
As a final sanity check we also measured the wall clock time ofour best implementation of all method. The results are essentially
identical to those shown above and are omitted for clarity.

5.4 Disk Access Experiments

Fraction of objects retrieved

The results in the previous section show that we can do true
rotation invariant matching so fast that CPU time is no longer thebottleneck, and we should therefore also attempt to minimize disk
accesses. We will compare to Linear Scan, which is the only othercompetitor that we are aware of that allows exact rotation
invariant indexing under Euclidean distance and DTW with aguarantee of no false dismissals. Recall that the lower bound usedby the VP-tree requires transforming the signal to the Fourier
space and calculating the Euclidean distance between thecoefficient magnitudes [38]. It is well understood that most of theenergy of the signal will be concentrated in a relatively small
number of these coefficients [37] and that using just a few large
valued coefficients is better than using all of them. We therefore
will perform experiments keeping just the first D coefficients,
were D = {4, 8, 16, 32}.
We count the fraction of items that must be retrieved from disk.
Figure 24 illustrates the results for the full projectile points andheterogeneous datasets over a range of dimensionalities.
0.12

Projectile Points

Heterogeneous

0.1
0.08
0.06
0.04
0.02
0

0

4
8
Dime 16 32
nsion
ality

4

8
Wedg
e: DTWWedg
e: Euclidean

16

32

Wedg
e: DTWWedg
e: EuclideanFigure 24: The fraction of items retrieved from disk to answer a 1nearest neighbor query, using dimensionalities D = {4, 8, 16, 32}.

6. CONCLUSIONS AND FUTURE WORK
We have introduced a method to support fast rotation-invariant
search of large shape datasets with arbitrary representations anddistance functions. Our method supports rotation limited queries
and mirror image invariance if desired.
Future work includes both extensions and applications of thecurrent work. We will attempt to extend this approach to theindexing of 3D shapes, and we have begun to use our algorithm asa subroutine in several data mining algorithms which attempt tocluster, classify and discover motifs in a variety ofanthropological datasets, including petroglyph and projectile
point databases.
Reproducible Research Statement: All datasets and images usedin this work are freely available at this URL [14].
Acknowledgements: We would like to acknowledge Chotirat Ann
Ratanamahatana and Longin Jan Latecki for useful suggestions.
We thank Jason Dorff for help with the primate skull images, andDr. Wendy Hodges for assistance with the reptile skulls. In
addition we thank the many donors of datasets.

7. REFERENCES[1] Adamek, T. and O’Connor, N.E. A multiscale representation
method for nonrigid shapes with a single closed contour.
IEEE Circuits and Systems for Video Technology, 14(5):742-753, 2004.
[2] Adamek, T. and O'Connor, N.E. Efficient contour-based
shape representation and matching. Multimedia Information
Retrieval 2003: 138-143.
[3] Attalla, E. and Siy, P. Robust shape similarity retrieval basedon contour segmentation polygonal multiresolution andelastic matching. Pattern Recognition, 38(12): 2229-2241,
2005.
[4] Bartolini, B., Ciaccia, P., and Patella, P. WARP: Accurate
Retrieval of Shapes Using Phase of Fourier Descriptors andTime Warping Distance. IEEE Trans. Pattern Anal. Mach.
Intell. 27(1): 142-147 (2005)[5] Bhanu, B. and Zhou, X. Face recognition from face profile
using dynamic time warping. In Proceedings of International
Conference on Pattern Recognition (ICPN’04), pp. 499-502,
2004.
[6] Borgefors, G. Hierarchical Chamfer Matching: A Parametric
Edge Matching Algorithm. IEEE Transactions on PatternAnalysis and Machine Intelligence, 10(6): 849-865,
November 1988.
[7] Cardone, A., Gupta, S.K., and Karnik, M. A survey of shapesimilarity assessment algorithms for product design andmanufacturing applications. ASME Journal of Computing
and Information Science in Engineering, 3(2): 109-118,
2003.
[8] Ciaccia, P. and Patella, M. Searching in metric spaces withuser-defined and approximate distances. ACM Trans.
Database Syst. 27(4): 398-437[9] Gdalyahu, Y. and Weinshall, D. Flexible syntactic matchingof curves and its application to automatic hierarchicalclassification of silhouettes. IEEE Transactions on PatternAnalysis and Machine Intelligence, 21(12): 1312-1328, Dec.
1999.
[10] Hodges, W. and Zamudio, K. Horned lizard (phrynosoma)
phylogeny inferred from mitochondrial genes andmorphological characters: understanding conflicts using
multiple approaches. Molecular Phylogenetics andEvolution, 31:961–971, 2004.
[11] Iwabe, N. Sister group relationship of turtles to the birdcrocodilian clade revealed by nuclear DNA-coded proteins.
Molecular Biology and Evolution, 22:810–813, 2004.
[12] Jalba, A.C., Wilkinson, M.H.F., Roerdink, J.B.T.M., Bayer,
M.M., and Juggins, S. Automatic Diatom Identification
using Contour Analysis by Morphological Curvature Scale
Spaces. Machine Vision and Applications, 16(4): 217-228,
2005.
[13] Karydis, Y., Nanopoulos, A., Papadopoulos, A.N., andManolopoulos, Y. Evaluation of Similarity Searching
Methods for Music Data in Peer-to-Peer Networks.
International Journal of Business Intelligence and Data
Mining, 1(2): 210-228, 2005.
[14] Keogh, E. www.cs.ucr.edu/~eamonn/shape/shape.htm, 2006.
[15] Keogh, E. http://www.cs.ucr.edu/~eamonn/time_series_data/[16] Keogh, E. Exact indexing of dynamic time warping. In[31] Rath, T. and Manmatha, R. Lower-Bounding of Dynamic

Proceedings of the 28th International Conference on Very
Large Data Bases, Hong Kong. pp 406-417, 2002.
[17] Keogh, E. and Kasetty, S. On the Need for Time Series Data
MiningBenchmarks:
A Survey and
Empirical
Demonstration. In Proceedings of the 8th ACM SIGKDDInternational Conference on Knowledge Discovery and Data
Mining, Edmonton, Canada. pp 102-111, 2002.
[18] Keogh, E., Palpanas, T., Zordan, V., Gunopulos, D., andCardle, M. Indexing Large Human-Motion Databases. In
Proceedings of the 30th International Conference on Very
Large Data Bases, Toronto, Canada, pp 780-791, 2004.
[19] D.J. Lee, R.B. Schoenberger, D.K. Shiozawa, X. Xu, and P.
Zhan, “Contour Matching for a Fish Recognition and
Migration Monitoring System”, SPIE Optics East, Two andThree-Dimensional Vision Systems for Inspection, Control,and Metrology II, vol. 5606-05, p. 37-48, Philadelphia, PA,
USA, October 25-28, 2004.
[20] Li, D. and Simske, S. Shape Retrieval Based on Distance
Ratio Distribution. HP Tech Report. HPL-2002-251, 2002.
[21] Li, Q., Lopez, I., and Moon, B. Skyline Index for Time
Series Data. IEEE Transactions on Knowledge and Data
Engineering. 16(6): pp 669-684, 2004.
[22] Ling, H. and Jacobs, D.W. Using the Inner-Distance forClassification of Articulated Shapes. In Proceedings of IEEE
Conference on Computer Vision and Pattern Recognition(CVPR), Vol. II, pp 719-726, 2005.
[23] Marzal, A. and Palazón V. Dynamic Time Warping of Cyclic
Strings for Shape Matching. ICAPR (2) 2005: 644-652[24] Mollineda, R. A., Vidal, E., and Casacuberta, F. Cyclic
Sequence Alignments: Approximate Versus Optimal
Techniques. International Journal of Pattern Recognition and
Artificial Intelligence (IJPRAI), 16(3): 291-299, 2002.
[25] Niennattrakul, V. Wanichsan, D. and Ratanamahatana, C.
Hand Geometry Verification Using Time Series
Representation. KES (2) 2007: 824-831[26] O’Brien, M. J. & Lyman, R. L. 2003 Resolving phylogeny:
Evolutionary archaeology’s fundamental issue. Essential
Tensions in Archaeological Method and Theory, edited by
VanPool T. L. and VanPool C. S., 115-135, Salt Lake City:
University of Utah Press.
[27] Olson, C. F. and Huttenlocher, D. P. Automatic Target
Recognition by Matching Oriented Edge Pixels. IEEE
Transactions on Image Processing, 6(1): 103-113, January
1997.
[28] Osada, R., Funkhouser, T., Chazelle, B., and Dobkin, D.
Shape Distributions. ACM Transactions on Graphics, 21(4):807-832, October, 2002.
[29] P. Protopapas, J. M. Giammarco, L. Faccioli, M. F. Struble,
R. Dave, C. Alcock (2006) Finding outlier light curves incatalogues of periodic variable stars. Monthly Notices of theRoyal Astronomical Society 369 (2), 677–696[30] Ratanamahatana, C. A. and Keogh, E. Three Myths about
Dynamic Time Warping. In Proceedings of SIAM
International Conference on Data Mining (SDM '05),
Newport Beach, CA, April 21-23, pp 506-510, 2005.

Time Warping Distances for Multivariate Time Series. Tech
Report MM-40, University of Massachusetts Amherst, 2002.
[32] Seba, A, (1734). Locupletissimi rerum naturalium thesauri
accurata descriptio Naaukeurige beschryving van het
schatryke kabinet der voornaamste seldzaamheden der
natuur. - Amsterdam, 1734-1765. 4 vols. 2º. - 394 B 26-29,
vol. 3, plate XXXV.
[33] Sebastian, T and Kimia, B. Curves vs. skeletons in object
recognition. Signal Processing 85(2): 247-263 (2005)[34] Sakoe, H. & Chiba, S. (1978). Dynamic programmingalgorithm optimization for spoken word recognition. IEEE
Trans. Acoustics, Speech, and Signal Proc., Vol. ASSP-26.
pp. 43-49.
[35] Thakoor, N. Gao, J. (2005) Shape classifier based on
generalized probabilistic descent method with hidden
Markov descriptor. Computer Vision, 2005. ICCV 2005.
Tenth IEEE International Conference on 495- 502 Vol. 1[36] Veltkamp, R. C. and Latecki, L. J. Properties andPerformance of Shape Similarity Measures. In Proceedingsof IFCS 2006 Conference: Data Science and Classification.
July, 2006.
[37] Vlachos, M., Hadjieleftheriou, M., Gunopulos, D. and
Keogh. E. Indexing Multi-Dimensional Time-Series with
Support for Multiple Distance Measures. In Proceedings of
the 9th ACM SIGKDD International Conference on
Knowledge Discovery and Data Mining, pp 216-225, August24-27, 2003, Washington, DC, USA.
[38] Vlachos, M., Vagena, Z., Yu, P. S., and Athitsos, V.
Rotation invariant indexing of shapes and line drawings. In
Proceedings of ACM Conference on Information andKnowledge Management (CIKM), pp 131-138, 2005.
[39] Wang, Z., Chi, Z., Feng, D., and Wang, Q. Leaf Image
Retrieval with Shape Features. In Proceedings of the 4thInternational Conference on Advances in Visual Information
Systems, pp 477- 487, 2000.
[40] Wei, L., Keogh, E., Van Herle, H., and Mafra-Neto, A.
Atomic Wedgie: Efficient Query Filtering for Streaming
Time Series. In Proceedings of the 5th IEEE International
Conference on Data Mining (ICDM 2005), pp 490-497,
2005.
[41] White, T. D. Human Osteology. 2nd edition. San Diego:
Academic Press, 2000.
[42] Woodford, B. J., Kasabov, N. K. & Wearing, C. H. Fruit
Image Analysis Using Wavelets, in K. Ko & N. Kasabov,
eds, Proceedings of the ICONIP/ANZIIS/ANNES'99
International Workshop, University of Otago Press, 1999:
pp. 88-91.
[43] Yu, S. Guan, L. A CAD system for the automatic detectionof clustered microcalcifications in digitized mammogram.
IEEE Trans. Med. Imag. 19 (2) (2000) 115–126.
[44] Zhang, D. and Lu, G. Review of shape representation anddescription techniques. Pattern Recognition, 37(1): 1-19,
2004.
[45] Zunic, J., Rosin, p., and Kopanja, L. Shape Orientability.
ACCV (2) 2006: pp 11-20.