%---------------------------------------------------------------
% 3 dimensional optical flow based algorithm for image registration
% By: Catherine Wang for finished our publised codes as part of her
% post-doctoral work at Univ of Texas M.D. Anderson
% Cancer Center (supervisors by Dr. Rex Cheung and
% Dr. Lei Dong)
% By: Rex Cheung for correcting the first difference caluculations, and
% enhancement, and extensive documentation
% To remove the dependence on image processing toolbox by
% writing 1. zero-padding program and 2. smoothing
% program. This gives users better control of the
% parameters. These programs could be downloaded from my
% earlier postings.
% Please cite: Implementation and validation of a three-dimensional
% deformable registration algorithm for targeted prostate cancer
% radiotherapy. Wang H, Dong L, Lii MF, Lee AL, de Crevoisier R,
% Mohan R, Cox JD, Kuban DA, Cheung R. IJROBP. 2005.
% Call or email me if you have questions.
%--------------------------------------------------------------
%
%NOTE: THE PHYSICAL INTERPRETATION OF DEMONS ALGORITHM by Dr. Rex Cheung.
%NOTE: PLEASE READ THIS, IT IS IMPORTANT FOR UNDERSTANDING THIS PROGRAM
%The algorithm of Demons involves:
%1. find the difference at each voxel
%position using the fixed image as the reference, the delta=movingimage -
%fixedimage (3D, minus in matlab does it element by element.
%2. Once it is established the location where the local deformation occur (i.e. there is
%more intensity than expected by using the fixedimage reference), calculate
%the downhill grad in x-, y- and z-directions, get the resultant downhill
%gradient (this the the steepest slope to decent the hill).
%3. this direction is that of a 'stream' flowing prependicularly to the contours of the
%'encergy' (L2 normal in this case) hill, use a step factor to multiply
%this to get the actual actionable stepX (in x-, y- and z-directions) size.
%4. then simply compute the new voxel position by taking a step for each
%voxel.
%5. simply assign the intensity of the movingimage at (newX, newY,
%newZ) as that at this same coordinate on the fixedimage.
%6. repeat until reach the max number of allowed iterations or a stopping criterion.
%
%Intuitive interpretation of Demons: In the above algorithm, Demons could be seen as a
%force pushing the contour of iso-intensity out (like the riddles created
%by throwing a stone into the pond). Eventually, the energy (i.e. the square
%of the difference) is dissipated by repeatedly calculating the delta energy and smoothing
%the delta.
%
%In this interpretation, the symmetric force in symmetric force
%Demons implementation is more of a computational divice than having a
%physical interpretation.
%
%---------------------------------------------------------------
% -------Variables for distance traveled by the intensity shift ------
% NOTE: if du = 1, then set the intensity at (i+1, j, k) =
% intensity at (i, j, k), so the intensity flow du in this step
% u = total distance moved to, return by calling demons[u,v,w]=cal_demons(img11,img22,[],[],[],100);
% v = total distance moved to,return by calling demons[u,v,w]=cal_demons(img11,img22,[],[],[],100);
% w = total distance moved to,return by calling demons[u,v,w]=cal_demons(img11,img22,[],[],[],100);
% du = one step of u as in du(i,j,k)=di(i,j,k)*gx1(i,j,k)/tmp; where di is
% the difference img2-img1 at voxel (i,j,k), tmp is the squared
% quantity to scale each step lower, gx1 is the centered first
% difference at (i,j,k)
% dv
% dw
% ------------- coordinate variables --------------------------
% i = current coordinate in x dimension
% j
% k
% ni = new i = the calculated new corrdinate in x dimension see note on
% du, ni simply = 1+i or in general du+i, SO ni = new i in one
% step
% nj
% nk
% --------- variables used to calculate force --------------
% di = delta of images = moving - fixed, voxel by voxel
% di2 = square of di
% tmp=grad12(i,j,k)+stepfactor*di2(i,j,k); original Demons
% (grad of 1)^2 (delta of img)^2
% idea is that the highened intensity (delta of img) flow down
% the steepest slope (grad of 1) therefore is like pushing the
% iso-intensity contour outward until there is no more (delta of
% img) difference in intensities at the solution. NOTE: steepest
% decent is prependicular to contours.
% tmp1=(img(ny1,nx1,nz1)*dx2*dy2+img(ny1,nx2,nz1)*dx1*dy2+img(ny2,nx1,nz1)*dx2*dy1+img(ny2,nx2,nz1)*dx1*dy1);
% tmp2=(img(ny1,nx1,nz2)*dx2*dy2+img(ny1,nx2,nz2)*dx1*dy2+img(ny2,nx1,nz2)*dx2*dy1+img(ny2,nx2,nz2)*dx1*dy1);
% def(i,j,k)=tmp1*dz2+tmp2*dz1;
% img = img2, buffer of moving image, used for calculations without
% destroy the original copy of img2
% img1 = fixed image
% img2 = moving image
% img1_ori = a copy of img1
% img2_ori
% gx1 = gradient by first difference centered at (i,j,k) in x-dimension
% for image 1
% gy1
% gz1
% (grad1)2=gx1.*gx1+gy1.*gy1+gz1.*gz1; get the centered 2nd difference
%
% gx2 = gradient by first difference centered at (i,j,k) in x-dimension
% for image 2
% gy2
% gz2
%---------------------------------------------------------
%---------------------------------------------------------------
% demons3d - 3 dimensional optical flow based algorithm
% for image registration
% By: Catherine Wang,
% By: Rex Cheung, cellular modules, bug-fixes and enhancement.
% Please cite: Implementation and validation of a three-dimensional
% deformable registration algorithm for targeted prostate cancer
% radiotherapy.Wang H, Dong L, Lii MF, Lee AL, de Crevoisier R,
% Mohan R, Cox JD, Kuban DA, Cheung R. IJROBP. 2005.
% Call or email me if you have questions.
%--------------------------------------------------------------
%-----------------------------------------
% deform image using known field
% img: image to be deformed
% u,v,w: displacement in x,y,z direction.
% NOTE: u,v,w,x,y,z and img all have the
% same dimensions.
%----------------------------------------
function def=Deform3Dimage(img,u,v,w)
imgsize=size(img);
def=[];
for k=1:imgsize(3) %iterator through each voxel, and perform voxel-wise calculations
for i=1:imgsize(1) %note that all vectors have the same dimensions
for j=1:imgsize(2)
ni=i+v(i,j,k);
nj=j+u(i,j,k);
nk=k+w(i,j,k);
if ni<1 || ni>imgsize(1) || nj<1 || nj>imgsize(2) || nk<1 || nk>imgsize(3) %check out of bound, set value to 0
def(i,j,k)=0;
else
ny1=floor(ni); %new y1 = floor to be conservative
nx1=floor(nj);
nz1=floor(nk);
ny2=round(ny1+1); %new y2 = unit step beyond new y1
nx2=round(nx1+1);
nz2=round(nz1+1);
if isequal(ny1,ni)
ny2=ny1;
end
if isequal(nx1,nj)
nx2=nx1;
end
if isequal(nz1,nk)
nz2=nz1;
end
dy1=ni-ny1; %if current position is not the new position, i.e. there is a displcement occurred.
dx1=nj-nx1;
dz1=nk-nz1;
dy2=1-dy1; %delta y2 defined as 1- delta y1
dx2=1-dx1;
dz2=1-dz1;
tmp1=(img(ny1,nx1,nz1)*dx2*dy2+img(ny1,nx2,nz1)*dx1*dy2+img(ny2,nx1,nz1)*dx2*dy1+img(ny2,nx2,nz1)*dx1*dy1); %3D grad's
tmp2=(img(ny1,nx1,nz2)*dx2*dy2+img(ny1,nx2,nz2)*dx1*dy2+img(ny2,nx1,nz2)*dx2*dy1+img(ny2,nx2,nz2)*dx1*dy1);
def(i,j,k)=tmp1*dz2+tmp2*dz1; %from Thirion's demon pushing iso-intensity contour algorithm
end
end
end
end