Code for the Earth Movers Distance (EMD)

Introduction
This is an implementation of the Earth Movers Distance, as described in . The EMD computes the distance between two distributions, which are represented by signatures. The signatures are sets of weighted features that capture the distributions. The features can be of any type and in any number of dimensions, and are defined by the user.

The EMD is defined as the minimum amount of work needed to change one signature into the other. The notion of "work" is based on the user-defined ground distance which is the distance between two features. The size of the two signatures can be different. Also, the sum of weights of one signature can be different than the sum of weights of the other (partial match). Because of this, the EMD is normalized by the smaller sum.

The code is implemented in C, and is based on the solution for the Transportation problem as described in 

Please let me know of any bugs you find, or any questions, comments, suggestions, and criticisms you have. If you find this code useful for your work, I would like very much to hear from you. Once you do, I'll inform you of any improvements, etc. Also, an acknowledgment in any publication describing work that uses this code would be greatly appreciated.

Usage
To use the code, perform the following steps:
1. download the files `emd.h` and `emd.c` (check the log of changes for latest changes).
2. In `emd.h`, modify the line
```typedef int feature_t;
```

to reflect your feature data type. Structures may be used too, for example

```typedef struct {
int X,Y,Z;
} feature_t;
```
3. To compute an EMD, call:
```float emd(signature_t *Signature1, signature_t *Signature2,
float (*func)(feature_t *, feature_t *),
flow_t *Flow, int *FlowSize);
```

where

`Signature1, Signature2`:
Pointers to the two signatures that their distance we want to compute.
`Dist`:
Pointer to the ground distance function. i.e. the function that computes the distance between two features.
`Flow`:
(Optional) Pointer to a vector of flow_t (defined in emd.h) where the resulting flow will be stored. Flow must have n1+n2-1 elements, where n1 and n2 are the sizes of the two signatures respectively. If `NULL`, the flow is not returned.
`FlowSize`:
(Optional) In case `Flow` is not `NULL`, `FlowSize` should point to a integer where the number of flow elements (always less or equal to n1+n2-1) will be written.
4. Compile `emd.c` and link it to your code.

The signature data type `signature_t` is defined in `emd.h` as:

```typedef struct
{
int n;                /* Number of features in the distribution */
feature_t *Features;  /* Pointer to the features vector */
float *Weights;       /* Pointer to the weights of the features */
} signature_t;
```

The feature data type `feature_t` is defined in `emd.h` and should be modified by the user to reflect his feature type.

In special cases, the user might want to change some of the values in the definitions in `emd.h`:

`#define MAX_SIG_SIZE 100`
The maximum allowed number of features in a signature
`#define MAX_ITERATIONS 100`
Maximum number of iterations. For ordinary problems, 100 should be more than enough.
`#define INFINITY 1e20`
`INFINITY` should be much larger than any other value in the problem
`#define EPSILON 1e-6`
`EPSILON` determines the accuracy of the solution.

Examples
1. `example1.c`
In this example the features are in three-dimensional space, and the ground distance is the Euclidean distance.
To try this example, you need to modify `feature_t` in `emd.h` to be

typedef struct { int X,Y,Z; } feature_t;
2. `example2.c`
Here, instead of providing a function to compute the ground distances, they are given in a predefined matrix. This is done by defining `feature_t` as `int` (the default), and setting the features in each signature to have consecutive numbers. The ground distance function uses these numbers as indexes to the predefined cost matrix.
The resulting flow is returned from the `emd` function by passing as the last parameters a vector of `flow_t`, and a pointer to `int` where the actual number of flow elements will be written.
Also, in this example the total weights of the two signatures are different. The first signature has a total weight of 1 while the second signature has a total weight of 0.9.

Log of changes
Date Decription
05/10/98 Changed from absolute error check to relative error check + Fixed a bug that caused the "Unexpected error in findBasicVariables" error message (thanks Mark Ruzon)
03/04/98 Fixed a bug that sometime caused a crash when the second signature had only one feature (thanks Mark Ruzon)
03/31/98 Fixed a bug in the case where *both* signatures have only one feature (thanks Rajesh Batra)

References
 Y. Rubner, C. Tomasi, and L. J. Guibas. A Metric for Distributions with Applications to Image Databases. Proceedings of the 1998 IEEE International Conference on Computer Vision, Bombay, India, January 1998, pp. 59-66.

 F. S. Hillier and G. J. Lieberman. Introduction to Mathematical Programming McGraw-Hill, 1990.