B-spline Grid, Image and Point based Registration: [O_new,Spacing]=refine_grid(O_trans,Spacing,sizeI) - File Exchange

Code covered by the BSD License

Highlights from
B-spline Grid, Image and Point based Registration

manually_warp_images(varargin) MANUALLY_WARP_IMAGES, is a GUI which shows the controlpoint grid used
showcs3(varargin) Function showcs3(V) will display x,y,z cross sections of the volume V
A=imresize3d(V,scale,tsize,nt... This function resizes a 3D image volume to new dimensions
Dlocal=jacobiandet_transform_... Determinant of the Jacobian of the transformation grid
E=strain(Ux,Uy,Uz) Calculate the Eulerian strain from displacement images
F=backwards2forwards(B) This function will turn a backward transformation field in to a
I3=movepixels(I1,T,mode) This function movepixels, will (backwards) translate the pixels
I=affine_registration_image(p... This function affine_registration_image, uses affine transfomation of the
I=imgaussian(I,sigma,siz) IMGAUSSIAN filters an 1D, 2D color/greyscale or 3D image with an
Igrid=make_grid_image(Spacing... This function creates a uniform 2d or 3D image of grid
Iout=affine_transform(Iin,M,m... Function affine_transform, is a wrapper of the (mex) functions
Iout=affine_transform_2d_doub... Affine transformation function (Rotation, Translation, Resize)
Iout=movepixels_2d_double(Iin... This function movepixels, will translate the pixels of an image
M=make_transformation_matrix(... This function make_transformation_matrix.m creates an affine
O_error=bspline_registration_... Function Registration_image. This function will create
O_trans=bspline_grid_fitting(... calculate which is the closest point on the lattic to the top-left
O_trans=make_init_grid(Spacin... This function creates a uniform 2d or 3D b-spline control grid
W=bspline_coefficients(u,v,w) Calculate influences of all neighborh b-spline knots
[Fx,Fy]=backwards2forwards_2d... This function will turn a backwards transformation field into
[I,T]=bspline_transform(O,I,S... Function bspline_transform, is a wrapper of the mex
[Iout,Tx,Ty]=bspline_transfor... Bspline transformation grid function
[O_error, O_grad]=bspline_reg... Function registration_gradient. This function will calculate a registration
[O_error, O_grad]=jacobiandet... Convert Grid vector to grid matrix
[O_error, O_grad]=jacobiandet... Gradient and error of the Jacobian of the transformation grid
[O_new,Spacing]=refine_grid(O... Refine image transformation grid of 1D b-splines with use of spliting matrix
[O_penalty, O_grad]=penalties... This function calculates the 2D or 3D equivalent of the 2D
[O_trans,Spacing,Xreg]=point_... This function creates a 2D or 3D b-spline grid, which transform space to
[O_trans,Spacing]=MakeDiffeom... This function MakeDiffeomorphic will make the b-spline grid diffeomorphic
[O_trans2,Spacing2]=inversegr... Calculates the (numeric) inverse of the current b-spline grid
[Tlocal,Dlocal]=bspline_trans... If Spacing and points are not an integer, we can not use fast look up tables,
[V1,V2]=get_example_data get_example_data is used to make some rigid and nonrigid transformed
[e,egrad]=affine_registration... This function affine_registration_error, uses affine transfomation of the
[hist12, hist1, hist2]=mutual... This function makes a 2D joint histogram of 1D,2D...ND images
[s_Trans, s_Rotation, s_Scale... This function gives scaling values for different parameters, such as
[t,I]=image_difference(V,U,ty... This function gives a registration error and error image I between the two
[x,fval,exitflag,output,grad]... FMINLBFGS finds a local minimum of a function of several variables.
addpaths.m add all needed function paths
err=squared_difference_double...
image_interpolation(Iin,Tloca... This function is used to transform an 2D image, in a backwards way with an
image_registration(Imoving,Is... This function image_registration is the most easy way to register two
maxNumCompThreads(varargin) maxNumCompThreads returns the number of available CPU cores, works with
compile_c_files.m
example_2d_affine.m
example_2d_nonrigid_1.m
example_2d_nonrigid_2.m
example_2d_nonrigid_3.m
example_3d_affine.m
example_3d_nonrigid.m
View all files

from B-spline Grid, Image and Point based Registration by Dirk-Jan Kroon
B-spline registration of two 2D / 3D images or corrsp. points, affine and with smooth b-spline grid.

[O_new,Spacing]=refine_grid(O_trans,Spacing,sizeI)

function [O_new,Spacing]=refine_grid(O_trans,Spacing,sizeI)
% Refine image transformation grid of 1D b-splines with use of spliting matrix
%                 
% Msplit=(1/8)*[4 4 0 0;
%               1 6 1 0;
%               0 4 4 0;
%               0 1 6 1;
%               0 0 4 4];
%
%     [O_new,Spacing] = refine_grid(O_trans,Spacing,sizeI)
%
% Function is written by D.Kroon University of Twente (August 2010)

Spacing=Spacing/2;

if(ndims(O_trans)==3)
    % Refine B-spline grid in the x-direction
    O_newA=zeros(((size(O_trans,1)-2)*2-1)+2,size(O_trans,2),2);
    i=1:size(O_trans,1)-3; j=1:size(O_trans,2); h=1:2;
    [I,J,H]=ndgrid(i,j,h); 
    I=I(:); J=J(:); H=H(:);
    
    ind=sub2ind(size(O_trans),I,J,H);
    P0=O_trans(ind); 
    P1=O_trans(ind+1);
    P2=O_trans(ind+2); 
    P3=O_trans(ind+3);
    Pnew=split_knots(P0,P1,P2,P3);
     
    ind=sub2ind(size(O_newA),1+((I-1)*2),J,H);
    O_newA(ind)=Pnew(:,1);
    O_newA(ind+1)=Pnew(:,2);
    O_newA(ind+2)=Pnew(:,3);
    O_newA(ind+3)=Pnew(:,4);
    O_newA(ind+4)=Pnew(:,5);

    % Refine B-spline grid in the y-direction
    O_newB=zeros(size(O_newA,1),((size(O_newA,2)-2)*2-1)+2,2);
    i=1:size(O_newA,2)-3; j=1:size(O_newA,1); h=1:2;
    [J,I,H]=ndgrid(j,i,h); 
    I=I(:); J=J(:); H=H(:);
    
    ind=sub2ind(size(O_newA),J,I,H);
    P0=O_newA(ind);
    P1=O_newA(ind+size(O_newA,1)); 
    P2=O_newA(ind+size(O_newA,1)*2); 
    P3=O_newA(ind+size(O_newA,1)*3);
    Pnew=split_knots(P0,P1,P2,P3);
   
    ind=sub2ind(size(O_newB),J,1+((I-1)*2),H);
    O_newB(ind                 )=Pnew(:,1);
    O_newB(ind+  size(O_newB,1))=Pnew(:,2);
    O_newB(ind+2*size(O_newB,1))=Pnew(:,3);
    O_newB(ind+3*size(O_newB,1))=Pnew(:,4);
    O_newB(ind+4*size(O_newB,1))=Pnew(:,5);
    
    % Set the final refined matrix
    O_new=O_newB;
    
    % Make sure a new uniform grid will have the same dimensions (crop)
    dx=Spacing(1); dy=Spacing(2);
    X=ndgrid(-dx:dx:(sizeI(1)+(dx*2)),-dy:dy:(sizeI(2)+(dy*2)));
    O_new=O_new(1:size(X,1),1:size(X,2),1:2);
else
    % Refine B-spline grid in the x-direction
    O_newA=zeros(((size(O_trans,1)-2)*2-1)+2,size(O_trans,2),size(O_trans,3),3);
    k=1:size(O_trans,3); j=1:size(O_trans,2); i=1:size(O_trans,1)-3; h=1:3;
    [I,J,K,H]=ndgrid(i,j,k,h); 
    I=I(:); J=J(:); K=K(:); H=H(:);
    
    ind=sub2ind(size(O_trans),I,J,K,H);
    P0=O_trans(ind); 
    P1=O_trans(ind+1);
    P2=O_trans(ind+2); 
    P3=O_trans(ind+3);
    
    Pnew=split_knots(P0,P1,P2,P3);
    
    ind=sub2ind(size(O_newA),1+((I-1)*2),J,K,H);
    O_newA(ind)=Pnew(:,1);
    O_newA(ind+1)=Pnew(:,2);
    O_newA(ind+2)=Pnew(:,3);
    O_newA(ind+3)=Pnew(:,4);
    O_newA(ind+4)=Pnew(:,5);

    % Refine B-spline grid in the y-direction
    O_newB=zeros(size(O_newA,1),((size(O_newA,2)-2)*2-1)+2,size(O_newA,3),3);
    k=1:size(O_newA,3); j=1:size(O_newA,1); i=1:size(O_newA,2)-3; h=1:3;
    [J,I,K,H]=ndgrid(j,i,k,h); 
    I=I(:); J=J(:); K=K(:); H=H(:);
    
    ind=sub2ind(size(O_newA),J,I,K,H);
    P0=O_newA(ind);
    P1=O_newA(ind+size(O_newA,1)); 
    P2=O_newA(ind+size(O_newA,1)*2); 
    P3=O_newA(ind+size(O_newA,1)*3);
    Pnew=split_knots(P0,P1,P2,P3);
   
    ind=sub2ind(size(O_newB),J,1+((I-1)*2),K,H);
    O_newB(ind                 )=Pnew(:,1);
    O_newB(ind+  size(O_newB,1))=Pnew(:,2);
    O_newB(ind+2*size(O_newB,1))=Pnew(:,3);
    O_newB(ind+3*size(O_newB,1))=Pnew(:,4);
    O_newB(ind+4*size(O_newB,1))=Pnew(:,5);
    
    % Refine B-spline grid in the z-direction
    O_newC=zeros(size(O_newB,1),size(O_newB,2),((size(O_newB,3)-2)*2-1)+2,3);
    k=1:size(O_newB,2); j=1:size(O_newB,1); i=1:size(O_newB,3)-3; h=1:3;
    [J,K,I,H]=ndgrid(j,k,i,h); 
    I=I(:); J=J(:); K=K(:); H=H(:);
    
    ind=sub2ind(size(O_newB),J,K,I,H);
    a=size(O_newB,1)*size(O_newB,2);
    P0=O_newB(ind);
    P1=O_newB(ind+a); 
    P2=O_newB(ind+a*2); 
    P3=O_newB(ind+a*3);
    Pnew=split_knots(P0,P1,P2,P3);
    
    ind=sub2ind(size(O_newC),J,K,1+((I-1)*2),H);
    a=size(O_newC,1)*size(O_newC,2);
    O_newC(ind    )=Pnew(:,1);
    O_newC(ind+  a)=Pnew(:,2);
    O_newC(ind+2*a)=Pnew(:,3);
    O_newC(ind+3*a)=Pnew(:,4);
    O_newC(ind+4*a)=Pnew(:,5);
    
    % Set the final refined matrix
    O_new=O_newC;
    
    % Make sure a new uniform grid will have the same dimensions (crop)
    dx=Spacing(1); dy=Spacing(2); dz=Spacing(3);
    X=ndgrid(-dx:dx:(sizeI(1)+(dx*2)),-dy:dy:(sizeI(2)+(dy*2)),-dz:dz:(sizeI(3)+(dz*2)));
    O_new=O_new(1:size(X,1),1:size(X,2),1:size(X,3),1:3);
end

function Pnew=split_knots(P0,P1,P2,P3)
Pnew(:,1)=(4/8)*(P0+P1);
Pnew(:,2)=(1/8)*(P0+6*P1+P2);
Pnew(:,3)=(4/8)*(P1+P2);
Pnew(:,4)=(1/8)*(P1+6*P2+P3);
Pnew(:,5)=(4/8)*(P2+P3);

Highlights from B-spline Grid, Image and Point based Registration

Highlights from
B-spline Grid, Image and Point based Registration