mri.m

采用网格化GRIDDING算法对螺旋轨迹采集的MRI原始数据进行图像重建。

global dcf_jackson dcf_meyer dcf_pipe dcf_qian dcf method;
load DCFs.mat;
dcf=dcf_jackson;
method=0;
msgbox('DCF Data Import Finished');

% --------------------------------------------------------------------
function help_menu_Callback(hObject, eventdata, handles)
% hObject    handle to help_menu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function version_Callback(hObject, eventdata, handles)
% hObject    handle to version (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
helpdlg('Designed by Gary Wong, 2004 June, HUST','Copy Right');
 

% --------------------------------------------------------------------
function import_results_Callback(hObject, eventdata, handles)
% hObject    handle to import_results (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------

function large_matrix_fast() 
% ------------------------ larger matrix gridding method (fast) -------------------------%
global fid dcf;
process=waitbar(0,'Please wait......');
%-----find the minimun upper limit of the submatrix----%
 M=max(fid.kx)*2;
 H=512;
 while H<M
     H=H*2;
 end
%------initialize the parameters------%
submatrix=zeros(H);
ws=dcf.*fid.c;
% --------subsample to a larger matrix----------%
for n=1:fid.len
    x=round(fid.kx(n)+H/2);
    y=round(fid.ky(n)+H/2);
    submatrix(x,y)=submatrix(x,y)+ws(n);
    waitbar(n/fid.len*0.8,process,'mapping onto a larger matrix......');
end
% -------------FFT the spiral signal---------------%
waitbar(0.8,process,'doing fast fourier transform......')
I_fft=fft2(submatrix);
waitbar(0.87,process,'exchanging the pixel position......');
I_fft=fftshift(I_fft);
waitbar(0.95,process,'chop the appropriate pixels to display......');
% --------chop the large matrix to an appropriate one and display--------%
for  idx_x=1:512
     for idx_y=1:512
         I(idx_x,idx_y)=I_fft(idx_x+H/2-256,idx_y+H/2-256);
     end
end
 waitbar(1,process,'ready to display the image......');
 MRImage=abs(I);
 close(process);
 imshow(MRImage,[]);
 title('MRI reconstruction using larger matrix'); 

 
 
function large_matrix() 
% ------------------------ larger matrix gridding method (precise) -------------------------%
global fid dcf;
process=waitbar(0,'Please wait......');
%-----find the minimun upper limit of the submatrix----%
 M=max(fid.kx)*2;
 H=512;
 while H<M
     H=H*2;
 end
%------initialize the parameters------%
submatrix=zeros(H);
weight=zeros(H);
ws=dcf.*fid.c;
% --------subsample to a larger matrix----------%
for n=1:fid.len
    x=round(fid.kx(n)+H/2);
    y=round(fid.ky(n)+H/2);
    submatrix(x,y)=submatrix(x,y)+ws(n);
    weight(x,y)=weight(x,y)+1;
    waitbar(n/fid.len*0.2,process);
end
% -----------nomalize the energy----------------%
 [w.x,w.y,w.v]=find(weight);
 for n=1:length(w.v)
     submatrix(w.x(n),w.y(n))=submatrix(w.x(n),w.y(n))/w.v(n); 
     waitbar(n/length(w.v)*0.7+0.2,process);
 end
% -------------FFT the spiral signal---------------%
I_fft=fft2(submatrix);
I_fft=fftshift(I_fft);
waitbar(0.95,process);
% --------chop the large matrix to an appropriate one and display--------%
for  idx_x=1:512
     for idx_y=1:512
         I(idx_x,idx_y)=I_fft(idx_x+H/2-256,idx_y+H/2-256);
     end
end
 waitbar(1,process);
 MRImage=abs(I);
 close(process);
 imshow(MRImage,[]);
 title('MRI reconstruction using larger weighted matrix'); 
 
%-------------------------- Jackson Double Size (data driven) -------------------------% 
function jackson_data()
global fid dcf fov;
process=waitbar(0,'Please wait......');
%----initialize the parameters-----%
% fov=25;
count=200;
dx=0.01*fov/512;
dk=(1.0/dx)/1024;
k_mri=zeros(1024);
k_weight=zeros(1024);
%----generate the convolution kernel------%
kernel=KaiserBessel(2,count,18.5547); 
kernel=kernel';
kernel=cat(1,1,kernel);
%----normalize the k-space coordination----%
norm.kx=fid.kx/dk;
norm.ky=fid.ky/dk;
kx=floor(norm.kx);
ky=floor(norm.ky);
index.xmin=floor(norm.kx-0.0001-2)-kx+1;
index.ymin=floor(norm.ky-0.0001-2)-ky+1;
%-------data driven interpolation-------%
for num=1:fid.len
        for idx_x=index.xmin(num):2
            icx=round(abs(kx(num)+idx_x-norm.kx(num))*count/2);
        for idx_y=index.ymin(num):2
            icy=round(abs(ky(num)+idx_y-norm.ky(num))*count/2);
			wfn=kernel(icx+1)*kernel(icy+1)*dcf(num);
			k_mri(kx(num)+idx_x+512,512-ky(num)-idx_y)=k_mri(kx(num)+idx_x+512,512-ky(num)-idx_y)+fid.c(num)*wfn;
			k_weight(kx(num)+idx_x+512,512-ky(num)-idx_y)=k_weight(kx(num)+idx_x+512,512-ky(num)-idx_y)+wfn;  
        end
        end
     waitbar(num/fid.len*0.8,process);   
end

%----normalize grid data to 1 scale----%
[w.x,w.y,w.c]=find(k_weight);
for num=1:length(w.c)
    k_mri(w.x(num),w.y(num))=k_mri(w.x(num),w.y(num))/w.c(num);
end
waitbar(0.9,process);
%-----FFT the signal------%
I=fft2(k_mri);
I=fftshift(I);
waitbar(0.95,process);
%-----chop the central 512x512 point----%
for  idx_x=256:767
    for idx_y=256:767
        head.complex(idx_x-255,768-idx_y)=I(idx_x,idx_y);
    end
end
waitbar(1,process);
%-----display the image-----%
head.modulus=abs(head.complex);
imshow(head.modulus,[]);
close(process);
title('MRI reconstruction using gridding method of data-driven interpolation');

%-------------------------- Jackson Double Size (grid driven) -------------------------%
function jackson_grid()
global fid dcf fov;
process=waitbar(0,'Please wait......');
%----initialize the parameters-----%
% fov=25;
count=200;
dx=0.01*fov/512;
dk=(1/dx)/1024;
k_mri=zeros(1024);
k_weight=zeros(1024);
%----generate the convolution kernel------%
c=KaiserBessel(2,count,18.5547); 
c=c';
c=cat(1,1,c);
%----normalize the k-space coordination----%
norm.kx=fid.kx/dk;
norm.ky=fid.ky/dk;
%-------grid driven interpolation-------%
for num=1:fid.len
    for i=1:1024 		
      kxi=i-512;
       dux=abs(norm.kx(num)-kxi);
  			if dux<=2
			    for j=1:1024
				    kyj=j-512;
					duy=abs(norm.ky(num)-kyj);
                 	if duy<=2
					    icx=floor(dux/0.01+0.5);
						icy=floor(duy/0.01+0.5);
						wccf=c(icx+1)*c(icy+1)*dcf(num);
						k_mri(i,j)=k_mri(i,j)+wccf*fid.c(num);
                        k_weight(i,j)=k_weight(i,j)+wccf;
                    end
                end
            end
     end
     waitbar(num/fid.len*0.8,process)
end
%----normalize grid data to 1 scale----%
[w.x,w.y,w.c]=find(k_weight);
for num=1:length(w.c)
    k_mri(w.x(num),w.y(num))=k_mri(w.x(num),w.y(num))/w.c(num);
end
waitbar(0.9,process);
%-----FFT the signal------%
I=fft2(k_mri);
%-----exchange position-----%
I=fftshift(I);
waitbar(0.95,process);
%-----chop the central 512x512 point----%
for  idx_x=256:767
    for idx_y=256:767
        head.complex(idx_x-255,768-idx_y)=I(idx_x,idx_y);
    end
end
waitbar(1,process);
%-----display the image-----%
head.modulus=abs(head.complex);
close(process);
imshow(head.modulus,[]);
title('MRI reconstruction using Jackson double-sized gridding method');


%----------------------------------- Direct Fourier Transform-------------------------------%
function direct()
global fid dcf fov;
process=waitbar(0,'Please wait......');
pix=512;
 mr=zeros(pix);
 ws=dcf.*fid.c;
%  fov=25;
 dx=0.01*fov/pix;
 dy=0.01*fov/pix;
   for idx_x=1:pix
     xi=(idx_x-pix/2)*dx;
     for idx_y=1:pix
          yj=(idx_y-pix/2)*dy;
           for num=1:fid.len
                temp=i*2*pi*(fid.kx(num)*xi+fid.ky(num)*yj);
           mr(idx_x,idx_y)=mr(idx_x,idx_y)+exp(temp)*ws(num);
           end
           waitbar(idx_x*idx_y/(512*512),process);
     end
  end
 I=abs(mr);
 close(process);
 imshow(I,[]);
 title('MRI reconstruction using direct fourier transform');
 
function k=KaiserBessel(s,n,beta)
% Kaiser-Bessel function established for MRI reconstruction
% W = KaiserBessel(S,N,BETA) returns the BETA-valued N-point Kaiser function for a S-multiple N*N cartesian grid 
kdu=s/(n-1);
bes = abs(besseli(0,beta));
	for ii=1:n
		xx=beta*sqrt(1-(ii*kdu/s)*(ii*kdu/s));
		k(ii)=besseli(0,xx)/bes;
    end