Code for Earth Mover Distance computation is available at
http://vision.stanford.edu/~rubner/emd/default.htm

We used the code for emd.c on that page with no changes.
We modified emd.h as follows:
/* DEFINITIONS */
//#define MAX_SIG_SIZE   100
#define MAX_SIG_SIZE   600
//#define MAX_ITERATIONS 500
#define MAX_ITERATIONS 5000
typedef struct {
    float X,Y;
} feature_t;

Our main program for calling emd.c was as follows, 
modified from the example provided on the web page above.

Input data for the main program consisted of 3 text files, 
lt5.txt, no5.txt, and gr5.txt.  Each was a text file of 500 rows containing 3 columns with
the responses to the optimal grating stimulus, the F1 response to the contrast modulated 
noise stimulus, and root mean sqaure response (all in spikes/sec)
of the 500 most responsive neurons in each of the 3 conditions 
ltmd (lt) =before recovery
nois (no) =after recovery viewing the contrast modulated noise stimulus
grat (gr) =after recovery viewing the grating stimuli
The F1 noise responses were all multiplied by 3.3416 to equalize the mean to that
of the optimal grating responses so that the two conditions would be weighted
equally for determination of the 500 mmost responsive neurons and the Flow vectors.

The identical code was run twice, once for ltmd->nois (as shown below) and a
second time for ltmd->grat data.

The output of interest was the set of Flow vectors identifying the corresponding points
in the ltmd->grat and ltmd->nois distributions (in files names "lt2no_points" and "lt2gr_points").  
These Flow vectors were displayed as heat maps using the Matlab script below.

//
//  main.c
//  ex1
//
//  Created by Michael Stryker on 5/6/14.
//  Copyright (c) 2014 Michael Stryker. All rights reserved.
//

/* example1.c */

#include <stdio.h>
#include <math.h>

#include "emd.h"

float dist(feature_t *F1, feature_t *F2)
{
    float dX = F1->X - F2->X, dY = F1->Y - F2->Y;
    return sqrt(dX*dX + dY*dY);
}

#define MPSSIZE 500

int main(int argc, const char * argv[])
{
    int i,k,n1,n2;
    FILE * file1, * file2, * fpout;
    
    feature_t f1[MPSSIZE] = {0};
    feature_t f2[MPSSIZE] = {0};
    int flosiz = (2*MPSSIZE)-1;
    flow_t flo[flosiz];
    float w1[MPSSIZE] , w2[MPSSIZE];
    float junk;
    
    for (i=0;i<MPSSIZE;i++) {
        w1[i]=1.;
        w2[i]=1.;
    }
    fpout = fopen("/Users/stryker/Documents/workspace/lt2no.csv",  "w");
   
    k=0;
    file1 = fopen("/Users/stryker/Documents/workspace/lt5.txt","r");
    while (fscanf(file1,"%f,%f,%f", &f1[k].X, &f1[k].Y,&junk) == 3) {
        k++;
        if (k >= 500) break;
    }
    n1 = k;
    printf("n1=%i\n",n1);
    fclose(file1);

    fprintf (fpout, "\n%s %i\n","/Users/stryker/Documents/workspace/lt5.csv",n1);
    for(i=0; i<k; i++) {
     if ((i % 1) == 0) {
        fprintf(fpout, "%i %f %f\n", i, f1[i].X, f1[i].Y);
     }
    }

    k=0;
    file2 = fopen("/Users/stryker/Documents/workspace/no5.txt","r");
    while (fscanf(file2,"%f,%f,%f", &f2[k].X, &f2[k].Y, &junk) == 3) {
        k++;
        if (k >= 500) break;
    }
    n2 = k;
    printf("n2=%i\n",n2);
    fclose(file2);
    
    fprintf (fpout, "\n%s %i\n","/Users/stryker/Documents/workspace/no5.csv",n2);
    for(i=0; i<k; i++) {
        if ((i % 1) == 0) {
            fprintf(fpout,"%i %f %f\n", i, f2[i].X, f2[i].Y);
        }
    }
    signature_t s1 = { n1, f1, w1};
    signature_t s2 = { n2, f2, w2};
    float       e;
    flosiz = n1+n2-1;
    
    e = emd(&s1, &s2, dist, &flo[0], &flosiz);
  
    for (i=0;i<flosiz;i++) {
        fprintf(fpout, "%i -> %i flo[%i]=%f\n",flo[i].from,flo[i].to,i,flo[i].amount);
    }
    fclose(fpout);
    
    fpout = fopen("/Users/stryker/Documents/workspace/lt2no_points",  "w");
    for (i=0;i<flosiz;i++) {
        fprintf(fpout, "%i %i %f %f %f %f\n",flo[i].from, flo[i].to,
                s1.Features[flo[i].from].X,s1.Features[flo[i].from].Y,
                s2.Features[flo[i].to].X,s2.Features[flo[i].to].Y);
    }
    fclose(fpout);
    printf("emd=%f\n", e);
   return 0;
}





%% Matlab script for producing heat maps to show changes from ltmd to noise or gratings 
%% using corresponding points files from earth mover distance calculation

pathname = 'C:\Users\stryker\Dropbox\workspace\';
filename = 'lt2gr_points';
pointsfile = fullfile(pathname,filename)
fp=fopen(pointsfile,'r');
b=fscanf(fp,'%i %i %f %f %f %f',[6,inf]);
fclose(fp);
hx=zeros(20);
hy=zeros(20);

for i=1:500
    % plot at the midpoint of the Flow vector and truncate rare data that would be outside the plotting area
    x = min(20, 1 + fix((b(5,i) + b(3,i))/2));  
    y = min(20, 1 + fix((b(6,i) + b(4,i))/2));  
    deltax = b(5,i) - b(3,i);
    deltay = b(6,i) - b(4,i);
    hx(x,y) = hx(x,y) + deltax; %accumulate x- and y- components of Flow vectors 
    hy(x,y) = hy(x,y) + deltay;
end
hxgr=transpose(hx(1:19,1:19));  %remove truncated data; transpose to match abscissa and ordinate of parts ABC of figure
figure('Name',['hxx' filename]);
image(hxgr);
set(gca,'YDir','normal'); 
hygr=transpose(hy(1:19,1:19));
figure('Name',['hyy' filename]);
image(hygr);
set(gca,'YDir','normal'); 

filename = 'lt2no_points';
pointsfile = fullfile(pathname,filename)
fp=fopen(pointsfile,'r');
b=fscanf(fp,'%i %i %f %f %f %f',[6,inf]);
fclose(fp);
hx=zeros(20);
hy=zeros(20);

for i=1:500
    x = min(20, 1 + fix((b(5,i) + b(3,i))/2));
    y = min(20, 1 + fix((b(6,i) + b(4,i))/2));
    deltax = b(5,i) - b(3,i);
    deltay = b(6,i) - b(4,i);
    hx(x,y) = hx(x,y) + deltax;
    hy(x,y) = hy(x,y) + deltay;
end
hxno=transpose(hx(1:19,1:19));
figure('Name',['hxx' filename]);
image(hxno);
set(gca,'YDir','normal'); 
hyno=transpose(hy(1:19,1:19));
figure('Name',['hyy' filename]);
image(hyno);
set(gca,'YDir','normal'); 


