Code covered by the BSD License

B-spline Grid, Image and Point based Registration

Dirk-Jan Kroon (view profile)

26 May 2008 (Updated )

B-spline registration of two 2D / 3D images or corrsp. points, affine and with smooth b-spline grid.

[Fx,Fy]=backwards2forwards_2d_double(Bx,By,H)
```function [Fx,Fy]=backwards2forwards_2d_double(Bx,By,H)
% This function will turn a backwards transformation field into
% a forwards transformation field.
%
% [Fx,Fy]=backwards2forwards_2d_double(Bx,By,H);
%
% inputs,
%   Bx,By : The backward transformation fields
%   H : The Splatting kernel
% outputs,
%   Fx,Fy : The forward transformation fields
%
% Function is written by D.Kroon University of Twente (Februari 2009)

% Gaussian kernel
hx_center=-floor(size(H,1)/2)-1;
hy_center=-floor(size(H,2)/2)-1;

Fx=zeros(size(Bx)); Fy=zeros(size(Bx));
Num=zeros(size(Bx));
[x,y]=ndgrid(1:size(Bx,1),1:size(Bx,2));
x=x(:); y=y(:);
nx=x+Bx(:); ny=y+By(:);
valx=-Bx(:); valy=-By(:);
for i=1:length(nx);
Tlocalx=nx(i); Tlocaly=ny(i);
% All the neighborh pixels involved in linear interpolation.
xBas0=floor(Tlocalx);  yBas0=floor(Tlocaly);
xBas1=xBas0+1; yBas1=yBas0+1;
% Linear interpolation constants (percentages)
xCom=Tlocalx-xBas0; yCom=Tlocaly-yBas0;
perc0=(1-xCom).*(1-yCom); perc1=(1-xCom).*yCom; perc2=xCom.*(1-yCom); perc3=xCom.*yCom;
for iHx=1:size(H,1)
for iHy=1:size(H,2)
for t=1:4
switch(t),
case(1), tx=xBas0+iHx+hx_center; ty=yBas0+iHy+hy_center; perc=perc0*H(iHx,iHy);
case(2), tx=xBas0+iHx+hx_center; ty=yBas1+iHy+hy_center; perc=perc1*H(iHx,iHy);
case(3), tx=xBas1+iHx+hx_center; ty=yBas0+iHy+hy_center; perc=perc2*H(iHx,iHy);
case(4), tx=xBas1+iHx+hx_center; ty=yBas1+iHy+hy_center; perc=perc3*H(iHx,iHy);
end
if(~((tx<1)||(tx>size(Bx,1))||(ty<1)||(ty>size(Bx,2)))),
Fx(tx,ty)=Fx(tx,ty)+valx(i)*perc; Fy(tx,ty)=Fy(tx,ty)+valy(i)*perc;
Num(tx,ty)=Num(tx,ty)+perc;
end
end
end
end
end
Fx=Fx./(Num+eps);
Fy=Fy./(Num+eps);
```