Guest User

Optimal transport registration

a guest
May 7th, 2012
373
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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×