非刚性图像配准 matlab简单示例 demons算法

2011-05-25 17:21

非刚性图像配准 matlab简单示例 demons算法，

% Clean clc; clear all; close all;

% Compile the mex files %compile_c_files

% Read two images I1=im2double(imread('ssftrinew1.png'));
I2=im2double(imread('ssftri.png'));

% Set static and moving image S=I2; M=I1;

% Alpha (noise) constant alpha=2.5;

% Velocity field smoothing kernel Hsmooth=fspecial('gaussian',[60 60],10);

% The transformation fields Tx=zeros(size(M)); Ty=zeros(size(M)); Tz=zeros(size(M));

[Sy,Sx] = gradient(S); for itt=1:200 % Difference image between moving and static image Idiff=M-S;

% Default demon force, (Thirion 1998) %Ux = -(Idiff.*Sx)./((Sx.^2+Sy.^2)+Idiff.^2); %Uy = -(Idiff.*Sy)./((Sx.^2+Sy.^2)+Idiff.^2);

% Extended demon force. With forces from the gradients from both % moving as static image. (Cachier 1999, He Wang 2005) [My,Mx] = gradient(M); Ux = -Idiff.* ((Sx./((Sx.^2+Sy.^2)+alpha^2*Idiff.^2))+(Mx./((Mx.^2+My.^2)+alpha^2*Idiff.^2)));
Uy = -Idiff.*
((Sy./((Sx.^2+Sy.^2)+alpha^2*Idiff.^2))+(My./((Mx.^2+My.^2)+alpha^2*Idiff.^2))); % When divided by zero Ux(isnan(Ux))=0; Uy(isnan(Uy))=0;

% Smooth the transformation field Uxs=3*imfilter(Ux,Hsmooth); Uys=3*imfilter(Uy,Hsmooth);

% Add the new transformation field to the total transformation field. Tx=Tx+Uxs; Ty=Ty+Uys; %M=movepixels(I1,Tx,Ty,Tz,0); M=movepixels_2d_double(I1,Tx,Ty,0); end gridelment=gridshow(); gridelment=movepixels_2d_double(im2double(gridelment),Tx,Ty,0); subplot(1,3,1), imshow(I1,[]); title('image 1'); subplot(1,3,2), imshow(I2,[]); title('image 2'); subplot(1,3,3), imshow(M,[]); title('Registered image 1'); figure,subplot(131),imshow(I1),subplot(132),imshow(abs(I2-M)),subplot(133),imshow(abs(I2-I1)) figure,imshow(gridelment)

function gridelment=gridshow() gridelment=ones(256,256)*255; for i=1:5:256 gridelment(i,:)=0; end for j=1:5:256 gridelment(:,j)=0; end gridelment=uint8(gridelment); imshow(gridelment);

function Iout=movepixels_2d_double(Iin,Tx,Ty,mode) % This function movepixels, will translate the pixels of an image % according to x and y translation images (bilinear interpolated). % % Iout = movepixels_2d_double(I,Tx,Ty,mode); % % Inputs; % Tx, Ty: The transformation images, describing the % (backwards) translation of every pixel in x and y direction. % mode: If 0: linear interpolation and outside pixels set to nearest pixel % 1: linear interpolation and outside pixels set to zero % (cubic interpolation only supported by compiled mex file) % 2: cubic interpolation and outsite pixels set to nearest pixel % 3: cubic interpolation and outside pixels set to zero % % Outputs, % Iout : The transformed image % % Function is written by D.Kroon University of Twente (February 2009) % Make all x,y indices [x,y]=ndgrid(0:size(Iin,1)-1,0:size(Iin,2)-1);

% Calculate the Transformed coordinates Tlocalx = x+Tx; Tlocaly = y+Ty;

% 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;

% limit indexes to boundaries check_xBas0=(xBas0<0)|(xBas0>(size(Iin,1)-1)); check_yBas0=(yBas0<0)|(yBas0>(size(Iin,2)-1)); xBas0(check_xBas0)=0;
yBas0(check_yBas0)=0;
check_xBas1=(xBas1<0)|(xBas1>(size(Iin,1)-1)); check_yBas1=(yBas1<0)|(yBas1>(size(Iin,2)-1)); xBas1(check_xBas1)=0;
yBas1(check_yBas1)=0;

Iout=zeros(size(Iin)); for i=1:size(Iin,3); Iin_one=Iin(:,:,i); % Get the intensities intensity_xyz0=Iin_one(1+xBas0+yBas0*size(Iin,1)); intensity_xyz1=Iin_one(1+xBas0+yBas1*size(Iin,1)); intensity_xyz2=Iin_one(1+xBas1+yBas0*size(Iin,1)); intensity_xyz3=Iin_one(1+xBas1+yBas1*size(Iin,1)); % Make pixels before outside Ibuffer mode if(mode==1||mode==3) intensity_xyz0(check_xBas0|check_yBas0)=0; intensity_xyz1(check_xBas0|check_yBas1)=0; intensity_xyz2(check_xBas1|check_yBas0)=0; intensity_xyz3(check_xBas1|check_yBas1)=0; end Iout_one=intensity_xyz0.*perc0+intensity_xyz1.*perc1+intensity_xyz2.*perc2+intensity_xyz3.*perc3; Iout(:,:,i)=reshape(Iout_one, [size(Iin,1) size(Iin,2)]); end