function dt = symearthmover2( v1,w1,v2,w2,thres  )
% dt = symearthmover2( v1,w1,v2,w2,thres  )
%  Calculates a distance between true weighted fiber directions and estimated weighted fiber directions
%   using an antipodally symmetric Earth Mover's Distance
% 
% Inputs:
%  - v1, 3xK0 matrix: true diffusion-weighted fiber directions
%  - w1, K0 x 1 vector: true fiber weights
%  - v2, 3xK matrix: estimated diffusion-weighted fiber directions
%  - w2, K x1 vector: estimated fiber weights
%  - thres, small scalar: all weights below threshold and corresponding fibers
%     will be dropped before computing the distance
if (sum(w1)==0);
  dt=2;
  return;
end
if (sum(w2)==0);
  dt=2;
  return;
end
w1=reshape(w1,[],1);
w2=reshape(w2,[],1);
w1=w1/sum(w1);
w2=w2/sum(w2);
v1n = v1(:,w1 > thres);
v2n = v2(:,w2 > thres);
w1n = w1(w1 > thres);
w2n = w2(w2 > thres);
w1n=w1n/sum(w1n);
w2n=w2n/sum(w2n);
d1 = reshape(pdist2(v1n',v2n'),[],1);
d2 = reshape(pdist2(-v1n',v2n'),[],1);
d12 = [d1,d2];
dd = min(d12')';
k1 = size(v1n,2);
k2 = size(v2n,2);
a1 = repmat(eye(k1),1,k2);
a2 = kron(eye(k2),ones(1,k1));
aeq= [a1; a2];
beq = [w1n;w2n];
options=optimset('MaxIter',2000,'Display','off');
[~,dt]=linprog(dd,[],[],aeq,beq,dd*0,dd*0+1,dd*0,options);
end

