• API
• FAQ
• Tools
• Archive
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);
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.
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)
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 %%%%%%%%%%%%%%%%%
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.