SHARE
TWEET

Optimal transport registration

a guest May 7th, 2012 271 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. function  [u dt cur energy] = match_pdfs_2d_haker( pdf1, pdf2)
  2.  
  3. % normalize pdfs
  4. cumsum1 = cumIntegrate(pdf1(:));
  5. cumsum2 = cumIntegrate(pdf2(:));
  6. pdf1 = pdf1/cumsum1(end);
  7. pdf2 = pdf2/cumsum2(end);
  8.  
  9. [M N] = size(pdf1);
  10.  
  11. % find an initial map, by doing 1D matchings in x and then in y.
  12. proj1 = zeros(1,N);
  13. proj2 = zeros(1,N);
  14. for i=1:N
  15.     cumsum1 = cumIntegrate(pdf1(:,i));
  16.     cumsum2 = cumIntegrate(pdf2(:,i));
  17.     proj1(i) = cumsum1(end);
  18.     proj2(i) = cumsum2(end);
  19. end;
  20.  
  21. a = find_warp(proj2, proj1);
  22. ga = gradientAccurate(a);
  23.  
  24. [X,Y] = meshgrid(1:N, 1:M);
  25. interpPDF2 = interp2(X,Y,pdf2, a, 1:M, 'spline',0.);
  26. b=zeros(M,N);
  27. for i=1:N
  28.     b(:,i) = find_warp( interpPDF2(:, i).*ga(i), pdf1(:,i));  
  29. end;
  30.  
  31. %u^0 = (a,b)
  32. u=reshape([repmat(a', [M 1]) b], M, N,2);
  33.  
  34. %let's do a constant number of iterations
  35. for t=1:3000
  36.      
  37.     %%%%%%%% compute grad^orthog Laplace^{-1} div(u^orthog) %%%%%%%%%
  38.     %%%%%%%% where orthog stands for a 90 degrees rotations,%%%%%%%%%
  39.     %%%%%%%% and Laplace^-1(g) solves Laplacian f = g       %%%%%%%%%
  40.  
  41.     [dudx dudy] = gradientAccurate(u);
  42.     divorthog = -dudy(:,:,1)+dudx(:,:,2); %div(u^orthog)
  43.     f=reshape(poicalc(divorthog(:),1,1,M,N), M, N); %Laplace^{-1} div(u^orthog)
  44.     [dfdx dfdy] = gradientAccurate(f); %grad Laplace^{-1} div(u^orthog)
  45.    
  46.     % 1./pdf1 * grad^orthog Laplace^{-1} div(u^orthog)
  47.     update1 = repmat(1./pdf1, [1 1 2]).*reshape([-dfdy dfdx], size(u));            
  48.  
  49.     %%%%%%%% upwind for Du %%%%%%%%%%%%%%%%%
  50.     [dudx dudy] = gradientAccurate(u);
  51.     dudxm = dudx;
  52.     dudym = dudy;
  53.     dudxp = dudx;
  54.     dudyp = dudy;
  55.     dudxm(:, 3:end-2, :) = (3*u(:, 3:end-2, :) - 4*u(:, 2:end-3, :) + u(:, 1:end-4, :))*0.5;
  56.     dudxp(:, 3:end-2, :) = (-3*u(:, 3:end-2, :) + 4*u(:, 4:end-1, :) - u(:, 5:end, :))*0.5;
  57.     dudym(3:end-2, :, :) = (3*u(3:end-2, :, :) - 4*u(2:end-3, :, :) + u(1:end-4, :, :))*0.5;
  58.     dudyp(3:end-2, :, :) = (-3*u(3:end-2, :, :) + 4*u(4:end-1, :, :) - u(5:end, :, :))*0.5;
  59.      
  60.     Dupp = abs(dudxp(:,:,1).*dudyp(:,:,2)-dudxp(:,:,2).*dudyp(:,:,1));
  61.     Dupm = abs(dudxp(:,:,1).*dudym(:,:,2)-dudxp(:,:,2).*dudym(:,:,1));
  62.     Dump = abs(dudxm(:,:,1).*dudyp(:,:,2)-dudxm(:,:,2).*dudyp(:,:,1));
  63.     Dumm = abs(dudxm(:,:,1).*dudym(:,:,2)-dudxm(:,:,2).*dudym(:,:,1));
  64.        
  65.     Du =  (update1(:, :,1)>0).*(update1(:, :,2)>0).*Dupp ...
  66.         + (update1(:, :,1)>0).*(update1(:, :,2)<0).*Dupm ...
  67.         + (update1(:, :,1)<0).*(update1(:, :,2)>0).*Dump ...
  68.         + (update1(:, :,1)<0).*(update1(:, :,2)<0).*Dumm;        
  69.        
  70.     dt(t) = 0.2*min(1./abs(update1(:))); %dt according to stability conditions
  71.    
  72.     % update : du/dt = 1./pdf1 * det(Jac(u)) * update1
  73.     u = u+dt(t).*repmat(Du, [1 1 2]).*update1;    
  74. end;
  75.  
  76.  
  77.  
  78. function y = find_warp(pdf1, pdf2) %1D matching
  79. cdf1 = cumIntegrate(pdf1);
  80. cdf2 = cumIntegrate(pdf2);
  81. cdf1 = cdf1/cdf1(end);
  82. cdf2 = cdf2/cdf2(end);
  83. N = length(cdf1);
  84. [ucdf1, uxx] = unique(cdf1, 'first');
  85. y = transpose(interp1(ucdf1,uxx,cdf2,'spline',N));
  86.    
  87.    
  88. function y = cumIntegrate(f)  % generic function to approximate the cumulative integral
  89.  y= cumtrapz(f);  % or y=cumsum, or y=intgrad1 with the external library
  90.  
  91. function varargout=gradientAccurate(f)
  92. if (length(f)==length(f(:)))
  93.     [varargout{1}] = gradient2(f);  %or my 5th order gradient2(f)  code
  94. else
  95.     [varargout{1} varargout{2}] = gradient2(f);  %or my 5th order gradient2(f) code
  96. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top