This may be of interest to some of the methodology geeks among you...
> Recently I read your paper 'A fast diffeomorphic image
> registration algorithm'. I am interested in the inversion and composition
> in both small deformation setting and diffeomorphic setting.
>
> Now I am trying to reproduce the result shown in Figure 1 in your paper.
> But I don't quite understand how to achieve the composition of x+u(x) and
> x-u(x). Would you please give me some help?
If you have MATLAB, then the attached source code may be helpful. The
functions can be compiled (in MATLAB) by:
mex dartel2.c optimizer2d.c diffeo2d.c -O
The composition function in diffeo2d.c actually does the composition. You
can also compose with Jacobian tensor fields (composition_jacobian) and
Jacobian determinant fields (composition_detjac).
Here is some example MATLAB code to illustrate the functions using randomly
generated deformations....
% Generate Some random small displacement
n = 64;
lam = 10;
L = spdiags(repmat([-1 2.01 -1],[n,1]),[-1 0 1],n,n);
L(1,n) = -1;
L(n,1) = -1;
O = speye(n);
L = kron(O,L)+kron(L,O);
L = L'*L*lam;
Ux = reshape(L\randn(n^2,1),n,n);Ux=Ux-mean(Ux(:));
Uy = reshape(L\randn(n^2,1),n,n);Uy=Uy-mean(Uy(:));
% Identity transform
[X,Y]=ndgrid(1:n,1:n);
% Generate small deformation forward and inverse
Def1=cat(3,X+Ux,Y+Uy);
Def2=cat(3,X-Ux,Y-Uy);
% Compose them
Comp1 = dartel2('comp',Def1,Def2);
Comp2 = dartel2('comp',Def2,Def1);
% Function for plotting deformations
showdef = inline(['plot(Y(:,:,1),Y(:,:,2),''k'','...
'Y(:,:,1)'',Y(:,:,2)'',''b'');'], 'Y');
% Plot the deformations
figure
subplot(2,2,1); showdef(Def1); axis image; axis([1 n 1 n]);
subplot(2,2,2); showdef(Def2); axis image; axis([1 n 1 n]);
subplot(2,2,3); showdef(Comp1); axis image; axis([1 n 1 n]);
subplot(2,2,4); showdef(Comp2); axis image; axis([1 n 1 n]);
% DARTEL-like deformations and compositions
V =cat(3, Ux, Uy);
Y1=dartel2('Exp', V,8);
Y2=dartel2('Exp',-V,8);
C1=dartel2('comp',Y1,Y2);
C2=dartel2('comp',Y2,Y1);
% Plot them
figure
subplot(2,2,1); showdef(Y1); axis image; axis([1 n 1 n]);
subplot(2,2,2); showdef(Y2); axis image; axis([1 n 1 n]);
subplot(2,2,3); showdef(C1); axis image; axis([1 n 1 n]);
subplot(2,2,4); showdef(C2); axis image; axis([1 n 1 n]);
Best regards,
-John
|