# A Linear Approximate Algorithm for Earth Mover's Distance with Thresholded Ground Distance

Authors: Longjie Li; Min Ma; Peng Lei; Xiaoping Wang; Xiaoyun Chen
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Effective and efficient image comparison plays a vital role in content-based image retrieval (CBIR). The earth mover’s distance (EMD) is an enticing measure for image comparison, offering intuitive geometric interpretation and modelling the human perceptions of similarity. Unfortunately, computing EMD, using the simplex method, has cubic complexity. FastEMD, based on min-cost flow, reduces the complexity to (O(N
2
log
⁡
N)). Although both methods can obtain the optimal result, the high complexity prevents the application of EMD on large-scale image datasets. Thresholding the ground distance can make EMD faster and more robust, since it can decrease the impact of noise and reduce the range of transportation. In this paper, we present a new image distance metric, EM
D
+, which applies a threshold to the ground distance. To compute  EM
D
+, the FastEMD approach can be employed. We also propose a novel linear approximation algorithm. Our algorithm achieves O
N complexity with the benefit of qualified bins. Experimental results show that (1) our method is 2 to 3 orders of magnitude faster than EMD (computed by FastEMD) and 2 orders of magnitude faster than FastEMD and (2) the precision of our approximation algorithm is no less than the precision of FastEMD.


---

## Body

## 1. Introduction

The development of advanced multimedia technology increases the importance of image retrieval, since large-scale image datasets are proliferating. Effective and efficient methods for computing the similarity between images is vitally important to content-based image retrieval (CBIR) systems.The earth mover’s distance (EMD) is a famous measure for image retrieval [1]. It not only takes the corresponding bins into account, but also considers the correlations between noncorresponding bins and thus is robust to histogram shift and rotation. The effectiveness of EMD and/or its variants is witnessed in many application domains, such as image retrieval [1–4], common pattern discovery [5, 6], mathematical symbol retrieval [7], contour matching [8], texture retrieval [1, 9], shape matching [10], graph matching [11], visual event recognition [12], face verification [13, 14], visual tracker [15], and uncertain or probability data [16, 17].The EMD defines the comparison of two histograms as thetransportation problem. (In [1], EMD is used to compute the distance between two signatures. For the convenience of statement, we use histogram in this paper. In fact, our algorithm can be used to both histogram and signature.) The minimum cost to transport one histogram to another one is taken as the distance of the two histograms. Rubner et al. [1] employ the transportationsimplex algorithm [18] to compute the EMD. Thus the computation of EMD is very expensive, with cubic time complexity (O
(
N
3
log
⁡
N
)). Pele and Werman [19] present the FastEMD algorithm which is faster than the simplex algorithm. FastEMD is based on the min-cost flow algorithm, thus it still suffers from quadratic complexity (O
(
N
2
log
⁡
N
)). Ling and Okada [10] propose EMD-L
1 which employs L
1-norm as the ground distance. EMD-L
1 can reduce the number of variables in the transportation simplex algorithm from O
(
N
2
) to O
(
N
), so the time complexity can be reduced to O
(
N
2
). But only using L
1-norm as ground distance limits the application of EMD-L
1. Rahimi and Kiram [20] propose a greedy algorithm with complexity of O
(
N
2
). Shirdhonkar and Jacobs [21] present a linear time approximate EMD algorithm,Wavelet EMD, based on a weighted wavelet transform. To approximate EMD, wavelet EMD uses the weighted wavelet coefficients of the difference histogram. Another linear time approximate EMD algorithm proposed by Jang et al. [22] is MHC which is used for dominant color descriptor (DCD). This algorithm employs Hilbert curves to fill the multidimensional color space and computes cost according to the Hilbert order of dominant color. To improve accuracy, MHC uses (
2
k
·
k
!
)
/
2 Hilbert curves to fill the k-dimensional space, then computes (
2
k
·
k
!
)
/
2 cost values, and takes the minimum cost as the distance between two images.Furthermore, EMD is not a metric if the two histograms have different total weight. To solve this drawback, Ljosa et al. [4] introduce a special bin, called thebank, for each image. This extra bin allows the two histograms to exchange mass. The bank has the same ground distance to all other bins. Pele and Werman [3] propose a variant of EMD, named EMD
^, which takes the difference of two weights as an extra mass. The moving distance of the extra mass is θ times maximum ground distance. If the ground distance is a metric and θ
≥
0.5, EMD
^ is a metric for all histograms.The EMD with thresholded ground distance has proven to be more robust, since it can reduce the interference of noise and correspond to the way humans perceive distance [2, 3, 19]. With the help of thresholded ground distance, FastEMD can reduce the number of edges of flow network from O
(
N
2
) to O
(
N
) [19] and thus makes an improvement in both accuracy and speed.In this paper, we present a new image distance metric,EM
D
+, based on thresholded ground distance. We prove that EM
D
+ is a metric by deducing that EM
D
+ equals EMD
^ with the same ground distance threshold. Therefore, we can employ FastEMD [19] to compute the value of EM
D
+. However, its quadratic complexity makes FastEMD a nonideal choice. According to the definition of EM
D
+, we derive a novel linear approximation algorithm which can compute EM
D
+ very fast. In experiments we find that our algorithm is 2 to 3 orders of magnitude faster than EMD (computed by FastEMD) and 2 orders of magnitude faster than FastEMD. Moreover, the results of our approximation algorithm are no worse than the exact results.In the remainder of the paper, we first review the earth mover’s distance andEMD
^ in Section 2. We then formalize the EM
D
+ in Section 3. Section 4 proposes the details of our novel linear approximation algorithm. In Section 5, we experimentally evaluate the practical performance of our algorithm. We conclude this paper in Section 6.


## 2. Earth Mover’s Distance

The earth mover’s distance is a perceptual and intuitive measure between two histograms. EMD which defines the distance between two histograms as the minimum cost that must be paid to transform one histogram to another one is based on the famoustransportation problem. We consider one histogram as asupplier with a given amount of products distributed at some places and another histogram as ademander that has some warehouses with a given limited capacity. To transport products from thesupplier to thedemander, the optimal solution is required to minimize the transportation cost. The following is the formal definition of EMD.Definition 1 (earth mover’s distance).
Given two histogramsP
=
{
(
p
1
,
w
p
1
)
,
…
,
(
p
m
,
w
p
m
)
} and Q
=
{
(
q
1
,
w
q
1
)
,
…
,
(
q
n
,
w
q
n
)
} with sizes m and n respectively, where w
p
i is the weight of bin p
i and w
q
j is the weight of bin q
j. Let d
i
j be the distance between p
i and q
j and let f
i
j represent the amount of weight transferred from p
i to q
j. The total cost of matching P and Q is defined as
(1)
WORK
=
∑
i
=
1
m
‍
∑
j
=
1
n
‍
f
i
j
d
i
j
,
with the following constraints:
(2)
f
i
j
≥
0
,
∑
j
=
1
n
‍
f
i
j
≤
w
p
i
,
∑
i
=
1
m
‍
f
i
j
≤
w
q
j
,
∑
i
=
1
m
‍
∑
j
=
1
n
‍
f
i
j
=
min
⁡
(
∑
i
=
1
m
‍
w
p
i
,
∑
j
=
1
n
‍
w
q
j
)
.
The earth mover’s distance is defined as the minimum total cost for solving this problem; that is,(3)
EMD
(
P
,
Q
)
=
min
⁡
{
f
i
j
}
∑
i
=
1
m
∑
j
=
1
n
f
i
j
d
i
j
∑
i
=
1
m
∑
j
=
1
n
f
i
j
.The EMD is a many-to-many matching method which can handle variable-size structures. If the ground distance is a metric and the total weights of two histograms are equal, that is,∑
i
=
1
m
‍
w
p
i
=
∑
j
=
1
n
‍
w
q
j, EMD is a metric [1]. However, in many applications, the total weights of two histograms cannot be equal. To solve this drawback, Pele and Werman [3] present a new EMD variant, EMD
^, which is defined as
(4)
EMD
^
θ
(
P
,
Q
)
=
min
⁡
{
f
i
j
}
∑
i
=
1
m
‍
∑
j
=
1
n
‍
f
i
j
d
i
j
+
|
∑
i
=
1
m
‍
w
p
i
-
∑
j
=
1
n
‍
w
q
j
|
×
θ
max
⁡
i
,
j
{
d
i
j
}
,
subject to the same constraints as the EMD.If the two histograms are probability distributions (i.e., total weight is one) EMD andEMD
^ are equivalent. Pele and Werman proved that if θ
≥
0.5 and the ground distance is a metric, EMD
^ is a metric. The proof can be found in [3].EMD
^ is more appropriate for two cases. One case is when the total weight of histograms is important. And the other one is when difference of total weight between histograms is distinguished [3].


## 3. New Distance Metric: EMD+

Thresholded distances are distances that are limited by a threshold. Letd
i
j be the distance between p
i and q
j and let μ
>
0 be a threshold; the thresholded distance between p
i and q
j is defined as d
i
j
μ
=
min
⁡
(
d
i
j
,
μ
).As mentioned in [2, 3, 19], thresholding the ground distance can improve the performance of EMD family. Motivated by this thought, we design a new variant of EMD with thresholded ground distance. We call this new measure EM
D
+. If the ground distance is a metric, the EM
D
+ is a metric. Firstly, we formally define EM
D
+ in Definition 2.Definition 2 (distance measure:E
M
D
+).
Given two histogramsP
=
{
(
p
1
,
w
p
1
)
,
…
,
(
p
m
,
w
p
m
)
} and Q
=
{
(
q
1
,
w
q
1
)
,
…
,
(
q
n
,
w
q
n
)
} with sizes m and n, respectively, where w
p
i is the weight of bin p
i and w
q
j is the weight of bin q
j. Let d
i
j be the ground distance between p
i and q
j, let f
i
j be the amount of weight transferred from p
i to q
j, and let μ
>
0 be a threshold of ground distance. If d
i
j
<
μ, EM
D
+ moves as much mass as possible from P to Q; otherwise, EM
D
+ does not move any mass from P to Q. Support ω is the maximum amount that can be transported from P to Q (ω
≤
min
⁡
{
∑
i
=
1
m
‍
w
p
i
,
∑
j
=
1
n
‍
w
q
j
}). The distance between P and Q is defined as
(5)
EM
D
+
(
P
,
Q
)
=
min
⁡
{
f
i
j
}
∑
i
=
1
m
‍
∑
j
=
1
,
d
i
j
<
μ
n
‍
f
i
j
d
i
j
+
(
max
⁡
{
∑
i
=
1
m
‍
w
p
i
,
∑
j
=
1
n
‍
w
q
j
}
-
ω
)
×
μ
,
s.t
(6)
f
i
j
≥
0
,
∑
i
=
1
m
‍
∑
j
=
1
,
d
i
j
<
t
n
‍
f
i
j
=
ω
,
∑
j
=
1
n
‍
f
i
j
≤
w
p
i
,
∑
i
=
1
m
‍
f
i
j
≤
w
q
j
.We denote the total weight ofP and Q as w
(
P
) and w
(
Q
) respectively. Without loss of generality, we assume w
(
P
)
≥
w
(
Q
). The condition d
i
j
<
μ in Definition 2 limits the range of bins that can be matched by p
i. Only the bins whose distances to p
i are smaller than the threshold have the chance to receive mass from p
i. The optimal solution is the minimum total cost of moving ω amount of mass from P to Q. After that, we ignore the limitation of capacity of Q and transport all the remainder weight of P to Q with cost of μ.Figure1 shows two examples of EM
D
+. From Figure 1, we obtain w
(
P
)
=
20, w
(
Q
)
=
18, and w
(
R
)
=
9. Set L
2-norm to be the ground distance and μ
=
2 to be the threshold. We get
(7)
EM
D
+
(
P
,
Q
)
=
(
10
+
6
×
2
)
+
(
20
-
16
)
×
2
≈
26.5
,
EM
D
+
(
P
,
R
)
=
(
6
+
3
×
2
)
+
(
20
-
9
)
×
2
≈
32.2
.
If we calculate the EMD values, we get
(8)
EMD
(
P
,
Q
)
=
10
+
6
×
2
+
2
×
2
18
≈
1.25
,
EMD
(
P
,
R
)
=
6
+
3
×
2
9
≈
1.14
.
That is to say, EMD
(
P
,
Q
)
>
EMD
(
P
,
R
), while EM
D
+
(
P
,
Q
)
<
EM
D
+
(
P
,
R
). Actually, Q is more similar to P than R. Thus, EM
D
+ has more discernibility.Three histogramsP, Q, and R. The light cyan circles indicate histogram P, the light yellow squares indicate histogram Q, and the light yellow hexagons indicate histogram R. In both figures, P is the same histogram. One circle, square, and hexagon denote one bin of P, Q, and R, respectively. The number inside the circle or square or hexagon is the weight of that bin. The blue arrows describe the flow of mass; the number beside the arrow denotes the weight of flow.
(a)
Pand Q
(b)
P and RIn fact, theEM
D
+ is equivalent to the EMD
^ with a thresholded ground distance and θ
=
1 [19]. Because the EMD
^ is a metric, we can deduce that EM
D
+ is a metric.Theorem 3.
If the ground distanced is a metric, then EMD+ is a metric.Proof.
A thresholded ground distance was introduced in [19]. d
i
j is the ground distance between p
i and q
j, and μ
>
0 is a threshold. The thresholded ground distance is defined as d
i
j
μ
=
min
⁡
(
d
i
j
,
μ
). If d is a metric, d
μ is a metric; a proof can be found in [19]. Without loss of generality, we assume w
(
P
)
≥
w
(
Q
). Set θ
=
1; we rewrite EMD
^ as(9)
EMD
^
1
(
P
,
Q
)
=
min
⁡
{
f
i
j
}
∑
i
=
1
m
‍
∑
j
=
1
n
‍
f
i
j
d
i
j
μ
+
(
w
(
P
)
-
w
(
Q
)
)
×
μ
=
min
⁡
{
f
i
j
}
∑
i
=
1
m
‍
∑
j
=
1
,
d
i
j
μ
<
μ
n
‍
f
i
j
d
i
j
μ
+
∑
i
=
1
m
‍
∑
j
=
1
,
d
i
j
μ
=
μ
n
‍
f
i
j
×
μ
+
(
w
(
P
)
-
w
(
Q
)
)
×
μ
=
min
⁡
{
f
i
j
}
∑
i
=
1
m
‍
∑
j
=
1
,
d
i
j
<
μ
n
‍
f
i
j
d
i
j
+
(
w
(
Q
)
-
∑
i
=
1
m
‍
∑
j
=
1
,
d
i
j
<
μ
n
‍
f
i
j
)
×
μ
+
(
w
(
P
)
-
w
(
Q
)
)
×
μ
=
min
⁡
{
f
i
j
}
∑
i
=
1
m
‍
∑
j
=
1
,
d
i
j
<
μ
n
‍
f
i
j
d
i
j
+
(
w
(
P
)
-
∑
i
=
1
m
‍
∑
j
=
1
,
d
i
j
<
μ
n
‍
f
i
j
)
×
μ
=
min
⁡
{
f
i
j
}
∑
i
=
1
m
‍
∑
j
=
1
,
d
i
j
<
μ
n
‍
f
i
j
d
i
j
+
(
w
(
P
)
-
ω
)
×
μ
=
EM
D
+
(
P
,
Q
)
.
Due to the fact thatEMD
^
1 is a metric with d
μ
(
·
,
·
), thus EM
D
+ is a metric.


## 4. Fast Computation of EMD+

Intrinsically,EM
D
+ has the same optimization problem as EMD and EMD
^. Hence, EM
D
+ can be computed by either transportation simplex algorithm [18] (e.g., EMD [1]) or min-cost flow algorithm (e.g., FastEMD [19]). However, the high complexity of simplex and min-cost flow algorithms drives us to develope an optimal method for efficient image retrieval. In this section, we proceed to develop a more efficient algorithm to obtain the value of EM
D
+.In the definition ofEM
D
+ (Definition 2), there is a very important condition, d
i
j
<
μ, which restricts that each bin in one histogram can only match a small part of bins in another histogram (Figure 1). Taking advantage of this characteristic, we derive a linear approximate algorithm to compute EM
D
+. We will explain our basic idea by an example (Figure 2).Figure 2
Red crosses indicate the distribution of histogramP and blue circles indicate the distribution of histogram Q. The radii of the green solid and cyan dash circles equal the threshold μ.In Figure2, it is clear that only the bins that fall into the solid or dash circle have the qualification of receiving mass from p
(
5,3
) or p
(
7,4
), respectively. We call these bins asqualified bins. If b
k is one qualified bin of p
(
5,3
), it is easy to get the qualified bin of p
(
7,4
) at the corresponding position by paralleling b
k. Inspired by this heuristic, we present to precompute all the possible qualified bins of a reference bin (e.g., (
0,0
)) and then utilize the possible qualified bins to improve the computation of EM
D
+.To precompute the qualified bins, we propose a novel method which is composed of functionQ
u
a
l
i
f
i
e
d
B
i
n
s (Algorithm 1) and g
e
t
B
i
n
s (Algorithm 2) (in Algorithm 2, function getBins() uses L
p-norm as the ground distance, in other cases it is distance dependent). Q
u
a
l
i
f
i
e
d
B
i
n
s first calls the function g
e
t
B
i
n
s to enumerate all possible qualified bins of the reference bin (e.g., (
0,0
)) (line 3) and then sorts the bins in ascending order according the distances to the reference bin (line 4). Let b
k be the kth qualified bin to reference bin. For each bin p
i, we can get its kth possible qualified bin by paralleling b
k. Notice that the reference bin itself does not appear in qualified bin set; the reason will be shown later.Algorithm 1:Get the qualified bin set of reference bin.
Input: μ: threshold; d: dimensionality of histogram
Output: Qualified bin set B
(1)B
←
∅; /
/  q
ualified bin set
(2)b
←
(
0
)
d; /
/  reference bin, if d
=
2, b
=
(
0,0
)
(3)B
← getBins(b
,
B
,
μ
,
d
,
1); /
/  A
lgorithm  2
(4) Sort(B);Algorithm 2:Recursive function getBins, Enumerate all possible qualified bins.
Input: b: one qualified bin; B: qualified bin set;
μ: threshold; d: dimensionality;
curr: current dimension.
Output: Qualified bin set B
(1)for  i
=
curr to d  do
(2)b
′
←
b;
(3)if  b
′
[
i
]
=
0  then
(4)
b
′
[
i
]
←
b
′
[
i
]
-
1;
(5)
if  ∥
b
′
∥
p
<
μ  then
(6)B
.
a
d
d
(
b
′
);
(7)        getBins(b
′
,
B
,
μ
,
d
,
i);
(8)
b
′
←
b,  b
′
[
i
]
←
b
′
[
i
]
+
1;
(9)
if  ∥
b
′
∥
p
<
μ  then
(10) B
.
a
d
d
(
b
′
);
(11)  getBins(b
′
,
B
,
μ
,
d
,
i);
(12)else
(13)if  b
′
[
i
]
<
0  then
(14)  b
′
[
i
]
←
b
′
[
i
]
-
1;
(15)else
(16)  b
′
[
i
]
←
b
′
[
i
]
+
1;
(17)if  ∥
b
′
∥
p
<
μ  then
(18)  B
.
a
d
d
(
b
′
);
(19)  getBins(b
′
,
B
,
μ
,
d
,
i);LetK be the size of qualified bin set B, K
≪
2
d
·
μ
2 (d is the dimensionality of histogram). Thus thegetBins() (Algorithm 2) can be done in O
(
K
) time andSort() (line 4) can be done in O
(
K
log
⁡
K
) time. In fact, K is very small corresponding to the size of histogram and can be treated as a constant. Thus the time of Algorithm 1 can be ignored (in our experiments (Section 5), the time of precomputing qualified bins is no more than 0.001 second).Next thing is to compute the value ofEM
D
+, that is, the core of our idea. The details are depicted in the following.(1)
The first step is to saturate all the zero-cost flows, since matching corresponding bins does not cause cost. After this step, many bins will be empty (Figure2).
(2)
The second step focuses on cross-bin matching. We plan to transport as much mass as possible when the cost is as low as possible. We employ the qualified bin setB generated by Algorithm 1 to accelerate the speed of our algorithm. We traverse B from front to end. In the kth iteration, we match each nonempty bin p
i to its corresponding kth qualified bin. If the kth qualified bin of p
i exists, we transfer the highest possible mass; otherwise, we do nothing. Since corresponding bins are already handled in the first step, we have no necessity to store the reference bin in the qualified bin set B.Algorithm3 gives the pseudocode of our algorithm. The proposed algorithm is named asLinearEMD, because its complexity is linear. The very important parameter B stores the qualified bins of original point (e.g., (
0,0
) for 2-dimensional histogram). Notice that B only needs to be computed once.Algorithm 3:
EM
D
+ between two histograms.
Input: P, Q: histogram; μ: threshold; B: qualified bins
of original point.
Output: Distance between P and Q
(1)cost
=
0;
(2)ω
=
0;
//  P
art 1, corresponding bin
(3)foreach  p
∈
P  do
(4)if  corresponding q exists  then
(5) mass
←
min
⁡
(
w
p
,
w
q
); // w
p: weight of bin p
(6) ω
←
ω
+
mass;
(7) w
p
←
w
p
-
mass;
(8) w
q
←
w
q
-
mass;
(9) if  w
p
=
0  then
(10)  remove p from P;
(11) if  w
q
=
0  then
(12)  remove q from Q;
//  Part 2, cross bin
(13)for  i
=
1 to B
.
size  do
(14)
foreach  p
∈
P  do
(15)if  ith qualified bin q exists  then
(16)mass
←
min
⁡
(
w
p
,
w
q
);
(17)ω
←
ω
+
mass;
(18)cost
←
cost
+
mass
×
distance
(
p
,
q
);
(19)w
p
←
w
p
-
mass, w
q
←
w
q
-
mass;
(20)if  w
p
=
0  then
(21)            removep from P;
(22)if  w
q
=
0  then
(23)            removeq from Q;
//  Part 3, compute the
EM
D
+
(24)EM
D
+
←
cost
+
(
max
⁡
(
w
(
P
)
,
w
(
Q
)
)
-
ω
)
×
μ;LinearEMD is composed of 3 parts. Part 1 focuses on the corresponding bins (lines 3~12). Since corresponding bin has no cost (d
i
j
=
0), we should move as much mass as possible from P to Q. Then we delete the empty bins (lines 9~12). Once part 1 is done, P and Q will shrink dramatically. The time complexity of part 1 is O
(
N
) where N
=
min
⁡
(
|
P
|
,
|
Q
|
). Part 2 adopts greedy idea to match cross bins (lines 13~23). The outer loop (line 13) guarantees that nearest bins get preferential treatment and each bin in P is shifted to the same position. Ideally, many bins in P or Q will be empty after several iterations. The speed will be faster and faster. Assume the average of bin numbers in each iteration is L (L
≪
N) and the size of B is K. K is a constant. The complexity of part 2 is O
(
K
L
)
=
O
(
L
). Part 3 which gets the result of EM
D
+ by a simple calculation has an O
(
1
) expensive (line 31). Therefore, the total cost is O
(
N
+
L
)
=
O
(
N
) (L
≪
N).


## 5. Experimental Evaluation

This section aims to experimentally test and verify the following things.(1)
LinearEMD is efficient enough for massive image retrieval.
(2)
Thresholding the ground distance can get a more effective result.
(3)
The result set obtained by LinearEMD is very close the exact value.To do that, we use image retrieval as an example to perform the experiments. Section5.1 gives the setup of our experiments and Section 5.2 describes some evaluation criteria. The experimental results are shown in Section 5.3.


### 5.1. Experimental Setup

Wang’s dataset [23] is a benchmark for image retrieval. This dataset contains 1,000 color images in 10 classes. Each class has 100 images. Since there are some ambiguous images in this dataset, we use a subset of Wang’s dataset filtered by [19] (we downloaded the filtered image list from http://www.seas.upenn.edu/ofirpele/FastEMD/). This subset is composed of 773 images. We exact the first 5 images from each classe as the queries. We perform our experiments on 5 different size RGB color histograms, that is, 8
*
8
*
8, 10
*
10
*
10, 12
*
12
*
12, 14
*
14
*
14, and 16
*
16
*
16, respectively. We use RGB color histogram because it is simple and has the ability to estimate the experimental goals. We do not claim that RGB color histogram is the best descriptor for image retrieval.We compareEM
D
+ with EMD
^ (in this section, when we say EMD
^, we mean EMD
^ with thresholded ground distance), EMD-L
1, and EMD. To compute EM
D
+, we use LinearEMD algorithm (we implement this algorithm in C++). For EMD-L
1, we use the authors implementation [10] (we downloaded the code from http://www.dabi.temple.edu/hbling/code_data.htm/). For EMD
^, we use FastEMD [19] (we downloaded the code from http://www.seas.upenn.edu/ofirpele/FastEMD/code/). For EMD, we also use FastEMD [19], since simplex is very expensive. The ground distance of EMD-L
1 is, of course, L
1-norm. For others, we choose both L
1- and L
2-norm. For EM
D
+ and EMD
^, the ground distance thresholds are 3.0, 3.5, 4.0, 4.5, and 5.0 corresponding to the size of RGB color histogram, respectively. We set the coefficient θ of EMD
^ to be equal to 1. Since EM
D
+ equals EMD
^, the value of EM
D
+ (computed by LinearEMD) is approximate value and the value of EMD
^ (computed by FastEMD) is the exact one.The experiments are run on a computer with one Intel Xeon E5645 2.40 GHz CPU, 8 GB memory. The OS is CentOS release 6.3 and gcc version is 4.4.6.


### 5.2. Evaluation Criteria

This section introduces some criteria to fairly evaluate the performance of these distance measures and algorithms.To evaluate the speed, we compare the average runtime of each query. The effectiveness is measured by the average number of correct retrieved images. We define the mean overlap rate (MOR) to calculate the common images retrieved by LinearEMD and FastEMD. MOR is the average Jaccard coefficient:(10)
MOR
=
100
%
n
∑
i
=
1
n
‍
|
R
A
i
∩
R
B
i
|
|
R
A
i
∪
R
B
i
|
,
where n is the number of query, R
A
i is the result set of the ith query retrieved by algorithm A, and R
B
i is the result set of the ith query retrieved by algorithm B.To measure how close the values computed by two methods are, we employ the symmetric mean absolute percentage error (SMAPE) [24]. SMAPE is defined as follows:
(11)
SMAPE
=
100
%
n
∑
i
=
1
n
‍
|
A
i
-
B
i
|
(
A
i
+
B
i
)
/
2
,
where n is the number of common retrieved images and A
i and B
i are the values of the ith image computed by algorithms A and B.


### 5.3.  Results

Firstly, we compare the speed of each method. Figure3 shows the average running time of each query. Obviously, our algorithm exceeds other methods in both ground distances. Our algorithm is faster by two to three orders of magnitude faster than EMD, two orders of magnitude faster than EMD
^ and an order of magnitude faster than EMD-L
1. The larger the size of RGB color histogram, the better our algorithm. When the size of RGB color histogram is 16
*
16
*
16, the LinearEMD is 300 to 500 times faster than the FastEMD and nearly 3000 to 4500 times faster than the EMD. Therefore, the LinearEMD is competent enough to handle massive image retrieval task.Figure 3
The average searching time of each query.The second thing of this section is to compare the retrieved images of each method. Figure4 shows the average number of correct retrieved images of each query. It is obvious that thresholding the ground distance of EMD (i.e., EM
D
+ and EMD
^) can obtain better retrieval result. That is, EM
D
+ and EMD
^ are more robust and have more discernibility. EM
D
+ and EMD
^ with L
1 norm as ground distance can achieve better results than L
2-norm as ground distance. However, for EMD, L
2 as ground distance retrieves more images. Generally, EM
D
+(L
1) (EM
D
+(L
p) denotes EM
D
+ with L
p-norm as the ground distance; EMD
^(L
p) and EMD(L
p) have similar meanings) gets the best results and EM
D
+(L
2) is better than EMD
^(L
2). In addition, Figure 4 illustrates the precision of EM
D
+ and EMD
^ will be improved while the size of histogram is increased. Although our method is the approximate algorithm, the precision of our method is no less, even higher, than the exact algorithm (FastEMD). Figure 5 shows an example of retrieved images. The images retrieved by EM
D
+(L
1), EM
D
+(L
2), and EMD
^(L
1) are all correct, yet the last image retrieved by EMD
^(L
2) is incorrect. However, there are only 2 correct images in the results retrieved by EMD-L
1, EMD(L
1), and EMD(L
2), respectively. This result shows the advantage of thresholding the ground distance once again.Figure 4
Results for image retrieval.An example of retrieved images. The image on the top is the query. (b)–(h) are the 5 nearest neighbors retrieved by different methods. The last two images in (c) look very similar, but they are not the same image. The first 4 images in (b)–(d) are the same.
(a)
Query
(b)
EMD
+
(
L
1
)
(c)
EMD
+
(
L
2
)
(d)
EMD
^
(
L
1
)
(e)
EMD
^
(
L
2
)
(f)
EMD-
L
1
(g)
EMD-
(
L
1
)
(h)
EMD-
(
L
2
)The third thing is to count the common images retrieved by our algorithm and the FastEMD. To this end, we employ the mean overlap rate (MOR). Figure6(a) shows the MOR of the two 50 nearest neighbors obtained by our algorithm and the FastEMD (we only count the correct retrieved images). The MOR(L
1)s (MOR(L
p) denotes the MOR value of LinearEMD and FastEMD with L
p-norm as the ground distance; SAMPE(L
p) has the same meaning.) are between 95% and 98% and the MOR(L
2)s are between 93% and 95%. Although the MOR(L
2)s are a little bit lower than the MOR(L
1)s, they are still very high. Figure 6(a) indicates that the results between our algorithm and the FastEMD share many common images.The MOR and SMAPE between our algorithm and FastEMD.
(a)
MOR of 50 nearest neighbors
(b)
SMAPE of 50 nearest neighborsFinally, we compute the SMAPE of our algorithm and the FastEMD. Like the MOR, we use the two 50 nearest neighbors and only count the correct retrieved images. The SMAPE results are shown in Figure6(b). Obviously, the SMAPE(L
1)s are lower than the SMAPE(L
2)s. In spite of this, the maximum value is lower than 4
%. Figure 6(b) illustrates that the values computed by our approximate algorithm highly approach the values computed by the FastEMD.In one word, our approximate algorithm is an optimal choice for comparing images.


## 5.1. Experimental Setup

Wang’s dataset [23] is a benchmark for image retrieval. This dataset contains 1,000 color images in 10 classes. Each class has 100 images. Since there are some ambiguous images in this dataset, we use a subset of Wang’s dataset filtered by [19] (we downloaded the filtered image list from http://www.seas.upenn.edu/ofirpele/FastEMD/). This subset is composed of 773 images. We exact the first 5 images from each classe as the queries. We perform our experiments on 5 different size RGB color histograms, that is, 8
*
8
*
8, 10
*
10
*
10, 12
*
12
*
12, 14
*
14
*
14, and 16
*
16
*
16, respectively. We use RGB color histogram because it is simple and has the ability to estimate the experimental goals. We do not claim that RGB color histogram is the best descriptor for image retrieval.We compareEM
D
+ with EMD
^ (in this section, when we say EMD
^, we mean EMD
^ with thresholded ground distance), EMD-L
1, and EMD. To compute EM
D
+, we use LinearEMD algorithm (we implement this algorithm in C++). For EMD-L
1, we use the authors implementation [10] (we downloaded the code from http://www.dabi.temple.edu/hbling/code_data.htm/). For EMD
^, we use FastEMD [19] (we downloaded the code from http://www.seas.upenn.edu/ofirpele/FastEMD/code/). For EMD, we also use FastEMD [19], since simplex is very expensive. The ground distance of EMD-L
1 is, of course, L
1-norm. For others, we choose both L
1- and L
2-norm. For EM
D
+ and EMD
^, the ground distance thresholds are 3.0, 3.5, 4.0, 4.5, and 5.0 corresponding to the size of RGB color histogram, respectively. We set the coefficient θ of EMD
^ to be equal to 1. Since EM
D
+ equals EMD
^, the value of EM
D
+ (computed by LinearEMD) is approximate value and the value of EMD
^ (computed by FastEMD) is the exact one.The experiments are run on a computer with one Intel Xeon E5645 2.40 GHz CPU, 8 GB memory. The OS is CentOS release 6.3 and gcc version is 4.4.6.


## 5.2. Evaluation Criteria

This section introduces some criteria to fairly evaluate the performance of these distance measures and algorithms.To evaluate the speed, we compare the average runtime of each query. The effectiveness is measured by the average number of correct retrieved images. We define the mean overlap rate (MOR) to calculate the common images retrieved by LinearEMD and FastEMD. MOR is the average Jaccard coefficient:(10)
MOR
=
100
%
n
∑
i
=
1
n
‍
|
R
A
i
∩
R
B
i
|
|
R
A
i
∪
R
B
i
|
,
where n is the number of query, R
A
i is the result set of the ith query retrieved by algorithm A, and R
B
i is the result set of the ith query retrieved by algorithm B.To measure how close the values computed by two methods are, we employ the symmetric mean absolute percentage error (SMAPE) [24]. SMAPE is defined as follows:
(11)
SMAPE
=
100
%
n
∑
i
=
1
n
‍
|
A
i
-
B
i
|
(
A
i
+
B
i
)
/
2
,
where n is the number of common retrieved images and A
i and B
i are the values of the ith image computed by algorithms A and B.


## 5.3.  Results

Firstly, we compare the speed of each method. Figure3 shows the average running time of each query. Obviously, our algorithm exceeds other methods in both ground distances. Our algorithm is faster by two to three orders of magnitude faster than EMD, two orders of magnitude faster than EMD
^ and an order of magnitude faster than EMD-L
1. The larger the size of RGB color histogram, the better our algorithm. When the size of RGB color histogram is 16
*
16
*
16, the LinearEMD is 300 to 500 times faster than the FastEMD and nearly 3000 to 4500 times faster than the EMD. Therefore, the LinearEMD is competent enough to handle massive image retrieval task.Figure 3
The average searching time of each query.The second thing of this section is to compare the retrieved images of each method. Figure4 shows the average number of correct retrieved images of each query. It is obvious that thresholding the ground distance of EMD (i.e., EM
D
+ and EMD
^) can obtain better retrieval result. That is, EM
D
+ and EMD
^ are more robust and have more discernibility. EM
D
+ and EMD
^ with L
1 norm as ground distance can achieve better results than L
2-norm as ground distance. However, for EMD, L
2 as ground distance retrieves more images. Generally, EM
D
+(L
1) (EM
D
+(L
p) denotes EM
D
+ with L
p-norm as the ground distance; EMD
^(L
p) and EMD(L
p) have similar meanings) gets the best results and EM
D
+(L
2) is better than EMD
^(L
2). In addition, Figure 4 illustrates the precision of EM
D
+ and EMD
^ will be improved while the size of histogram is increased. Although our method is the approximate algorithm, the precision of our method is no less, even higher, than the exact algorithm (FastEMD). Figure 5 shows an example of retrieved images. The images retrieved by EM
D
+(L
1), EM
D
+(L
2), and EMD
^(L
1) are all correct, yet the last image retrieved by EMD
^(L
2) is incorrect. However, there are only 2 correct images in the results retrieved by EMD-L
1, EMD(L
1), and EMD(L
2), respectively. This result shows the advantage of thresholding the ground distance once again.Figure 4
Results for image retrieval.An example of retrieved images. The image on the top is the query. (b)–(h) are the 5 nearest neighbors retrieved by different methods. The last two images in (c) look very similar, but they are not the same image. The first 4 images in (b)–(d) are the same.
(a)
Query
(b)
EMD
+
(
L
1
)
(c)
EMD
+
(
L
2
)
(d)
EMD
^
(
L
1
)
(e)
EMD
^
(
L
2
)
(f)
EMD-
L
1
(g)
EMD-
(
L
1
)
(h)
EMD-
(
L
2
)The third thing is to count the common images retrieved by our algorithm and the FastEMD. To this end, we employ the mean overlap rate (MOR). Figure6(a) shows the MOR of the two 50 nearest neighbors obtained by our algorithm and the FastEMD (we only count the correct retrieved images). The MOR(L
1)s (MOR(L
p) denotes the MOR value of LinearEMD and FastEMD with L
p-norm as the ground distance; SAMPE(L
p) has the same meaning.) are between 95% and 98% and the MOR(L
2)s are between 93% and 95%. Although the MOR(L
2)s are a little bit lower than the MOR(L
1)s, they are still very high. Figure 6(a) indicates that the results between our algorithm and the FastEMD share many common images.The MOR and SMAPE between our algorithm and FastEMD.
(a)
MOR of 50 nearest neighbors
(b)
SMAPE of 50 nearest neighborsFinally, we compute the SMAPE of our algorithm and the FastEMD. Like the MOR, we use the two 50 nearest neighbors and only count the correct retrieved images. The SMAPE results are shown in Figure6(b). Obviously, the SMAPE(L
1)s are lower than the SMAPE(L
2)s. In spite of this, the maximum value is lower than 4
%. Figure 6(b) illustrates that the values computed by our approximate algorithm highly approach the values computed by the FastEMD.In one word, our approximate algorithm is an optimal choice for comparing images.


## 6. Conclusion

In this paper, we have presented a new variant of EMD with thresholded ground distance, named asEM
D
+. We proved that EM
D
+ is a metric for any kind of histogram. To compute the value of EM
D
+, we have proposed a linear approximate algorithm. This algorithm is an effective and efficient method for image retrieval. Experimental results show that our method break through the comparing methods. And thus, our method is competent enough for massive image retrieval task. Our method may also be able to compare other histograms (e.g., HSV color histogram) and descriptors (e.g., SIFT).



---

*Source: 406358-2014-01-16.xml*