2d optical flow demon for mono-modal image registration

by

 

InSPIRE: Matlab wrap component

opticalflow2d()
%Int J Radiat Oncol Biol Phys. 2005 Mar 1;61(3):725-35.
%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.
%Please do email or call me to discuss your comments.
%This was written by my post-doc Dr. H. Wang in 2004.
%For Matlab Central posting, I made it usable for png files. 
%Mainline, PA, Februry 29, 2012.


function opticalflow2d()
clear all;

figure
colormap(gray);
%fid1=fopen('img7_256.raw','r');
%fid2=fopen('img4_256.raw','r');
%ssize=256;
%ct1=[];
%ct2=[];
%clims=[700,1200];
%for i=1:12
%    refimg=fread(fid1,[ssize ssize],'uint16');
%    ct1=refimg';    %refimg(1:2:ssize,1:2:ssize)';
%    defimg=fread(fid2,[ssize ssize],'uint16');
%    ct2=defimg';    %defimg(1:2:ssize,1:2:ssize)';
%end

%imwrite(ct1, 'ctone.png', 'png');
%imwrite(ct2, 'cttwo.png', 'png');

ct1 = imread('prostatemask1.20.png');
ct2 = imread('prostatemask2.20.png');
ssize = 256;
ct1 = double(ct1);
ct2 = double(ct2);

ct2_ori=ct2;
subplot(221)
imagesc(ct1);
subplot(222);
imagesc(ct2);
alpha=20;    
% havg=fspecial('gaussian',7,2);
 havg=[1 0; 0 1];

% 
%     h = fspecial('gaussian',3,1);

    ex1=(circshift(ct1,[0,-1])-circshift(ct1,[0,1]))/2;
    ey1=(circshift(ct1,[-1,0])-circshift(ct1,[1,0]))/2;
au=zeros(ssize,ssize);
av=zeros(ssize,ssize);
for n=1:10
alpha=alpha+20;
    
    ex2=(circshift(ct2,[0,-1])-circshift(ct2,[0,1]))/2;
    ey2=(circshift(ct2,[-1,0])-circshift(ct2,[1,0]))/2;
    ex=(ex1+ex2)/2;
    ey=(ey1+ey2)/2;
    et=ct2-ct1;
u=zeros(ssize,ssize);
v=zeros(ssize,ssize);
uu=u;
vv=v;
    u=uu-ex.*(ex.*uu+ey.*vv+et)./(alpha*alpha+ex.*ex+ey.*ey);
    v=vv-ey.*(ex.*uu+ey.*vv+et)./(alpha*alpha+ex.*ex+ey.*ey);    
        %uu=smooth(u,0.1);
        %vv=smooth(v,0.1);
        windowSize = 5;
        uu=filter(ones(1,windowSize)/windowSize,1,u); %running average filter
        vv=filter(ones(1,windowSize)/windowSize,1,v);
        %uu=u;
        %vv=v;
%au=au+uu;
%av=av+vv;

def=zeros(ssize,ssize);
def1=zeros(ssize,ssize);
for i=80:180 %i=1:ssize
    for j=80:180%1:ssize
        ii=round(i+v(i,j));
        jj=round(j+u(i,j));
        if ii>0 & ii<ssize+1 & jj>0 & jj<ssize+1
            def1(i,j)=ct2(ii,jj);
        end
    end
end
ct2=def1;
subplot(223);
%imshow(def1-ct1,[]);

imagesc(et);
subplot(224)
imagesc(def1);
    pause(1);
    n
imwrite(ct2,'output.png','png');
end


Contact us