I have attached a 2D basis function warping routine that you may find useful.
Type "help basis2d" for the operating instructions. The code is fairly old,
but hopefully should work OK.
Regards,
-John
| I think you are going to have to use some specifically 2D method. John A
| might have some suggestions?
function [F1,DefX,DefY] = basis2d(G0,F0,params)
% Nonlinear warp of one image to another
% FORMAT [F1,DefX,DefY] = basis2d(G0,F0,params)
% G0 - target image
% F0 - object image
% params(1:2) - number of basis functions
% - suggest [8 8]
% params(3) - regularization
% - suggest 0.01
% - lower value - more regularization
% params(4) - # iterations
% - suggest 12
% F1 - image F0 warped to fit image G0
% DefX & DefY - the deformation field
tst = 1;
if nargin<3, params = [6 6 0.1 16]; end;
nx = params(1);
ny = params(2);
sd1 = params(3);
its = params(4);
% The (seperable) basis functions
BX = spm_dctmtx(size(G0,1), nx);
BY = spm_dctmtx(size(G0,2), ny);
% Inverse covariance describing the apriori distribution of the
% parameters
kx=(pi*((1:nx)'-1)/size(F0,1)).^2; ox=ones(nx,1);
ky=(pi*((1:ny)'-1)/size(F0,2)).^2; oy=ones(ny,1);
IC0 = kron(ky.*ky,ox) + kron(oy,kx.*kx) + 2*kron(ky,kx);
IC0 = IC0*(sd1.^(-2));
IC0 = diag([IC0 ; IC0 ; 0]);
% Smooth the images
G = spm_conv(G0,8);
F = spm_conv(F0,8);
% Starting estimates for parameters
x1 = zeros(nx*ny*2+1,1);
x1(nx*ny*2+1)=1;
for iter=1:its
% Generate deformation field from parameters
x11=reshape(x1(1:2*nx*ny),nx*ny,2);
DefX = ones(1,size(G,1))'*(1:size(G,2)) ...
+ BX*reshape(x11(:,1),nx,ny)*BY';
DefY = (1:size(G,1))'*ones(1,size(G,2)) ...
+ BX*reshape(x11(:,2),nx,ny)*BY';
% Sample the image according to the deformation field
F1 = slice2d(F,DefY,DefX);
if (tst)
F10 = slice2d(F0,DefY,DefX);
subplot(2,2,1);imagesc((G0)');title('G0');axis image xy off
subplot(2,2,2);imagesc((F10)');title('F1');axis image xy off;
subplot(2,2,3);imagesc((F0)');title('F0');axis image xy off
hold on;
plot(DefY(:,1:8:256),DefX(:,1:8:256),'r',DefY(1:8:256,:)',DefX(1:8:256,:)','r');
hold off;
subplot(2,2,4);imagesc((F10-G0)');title('F1-G0');axis image xy off
drawnow;
end
% Gradients of warped image
[dDX,dDY]=gradient(F);
DX = slice2d(dDX,DefY,DefX);
DY = slice2d(dDY,DefY,DefX);
sc = x1(nx*ny*2+1);
wt = F1~=0;
wt = conv2(wt,ones(3)/9,'same');
wt = wt>0.9;
G1=-G .*wt;
F1=F1.*wt + G1*sc;
DX=DX.*wt;
DY=DY.*wt;
% Genarate matrix Alpha and vector Beta, where:
% Alpha = A'*A;
% Beta = A'*b;
% A = [diag(DX(:))*kron(BY,BX) diag(DY(:))*kron(BY,BX) G(:)];
% b = F(:);
Alpha = zeros(2*nx*ny+1);
Beta = zeros(2*nx*ny+1,1);
for y=1:size(G,2)
BYBY = BY(y,:)'*BY(y,:);
BXDX = kron(DX(:,y),ones(1,nx)).*BX;
BXDY = kron(DY(:,y),ones(1,nx)).*BX;
A11= kron(BYBY, BXDX'*BXDX);
A12= kron(BYBY, BXDX'*BXDY);
A13= kron(BY(y,:)', BXDX'*G1(:,y));
A22= kron(BYBY, BXDY'*BXDY);
A23= kron(BY(y,:)', BXDY'*G1(:,y));
A33= G1(:,y)'*G1(:,y);
Alpha = Alpha + [A11 A12 A13
A12' A22 A23
A13' A23' A33];
B1 = kron(BY(y,:)', BXDX'*F1(:,y));
B2 = kron(BY(y,:)', BXDY'*F1(:,y));
B3 = G1(:,y)'*F1(:,y);
Beta = Beta + [B1 ; B2 ; B3];
end
% Note that certain elements of alpha and beta
% can be generated by:
% [A11,B1]=spm_kronutil(DX,F1,BX,BY);
% [A22,B2]=spm_kronutil(DY,F1,BX,BY);
d = F1;
ss = (d(:)'*d(:));
n = sum(wt(:));
tmp = diff(F1,1);
w1 = sqrt(ss/(2*tmp(:)'*tmp(:)));
tmp = diff(F1,2);
w2 = sqrt(ss/(2*tmp(:)'*tmp(:)));
mul = (2*pi*w1)^(-0.5) * (2*pi*w2)^(-0.5);
nu = (n-(2*nx*ny+1)) * mul;
Vr = ss/nu;
fprintf('iter # %.2d res = %g smo=(%g,%g: %g)\n', iter, Vr, w1,w2, mul);
% Solve for new value of x1, using Alpha and Beta
% using residual sum of squares to influence
% regularization
x1 = (Alpha + IC0*Vr)\(Alpha*x1 - Beta);
end
F1 = slice2d(F0,DefY,DefX);
|