Code covered by the BSD License  

Highlights from
3D Cone beam CT (CBCT) projection backprojection FDK MLEM reconstruction MATLAB codes for students

image thumbnail

3D Cone beam CT (CBCT) projection backprojection FDK MLEM reconstruction MATLAB codes for students

by

 

10 Mar 2012 (Updated )

3D Cone beam CT (CBCT) projection backprojection MLEM FDK reconstruction MATLAB source codes

Demo3_MLEM.m
addpath('bin');

%% Parameter setting %%

param.nx = 64;
param.ny = 64;
param.nz = 64;

param.sx = 64; % mm
param.sy = 64; % mm
param.sz = 64; % mm

%The real detector panel pixel density (number of pixels)
param.nu = 128;		
param.nv = 100;

% Detector setting, according to Varian Trilogy OBI (real size)
param.su = 128;	% mm
param.sv = 100;	% mm

% X-ray source and detector setting
param.DSD = 1000;    %  Distance source to detector 
param.DSO = 600;	%  X-ray source to object axis distance

% angle setting
param.dir = -1;   % gantry rotating direction
param.dang = 1;
param.deg = 0:param.dang:360;
param.deg = param.deg*param.dir;
param.nProj = length(param.deg);


% filter='ram-lak','cosine', 'hamming', 'hann' 
param.filter='ram-lak'; % high pass filter

param.dx = param.sx/param.nx;
param.dy = param.sy/param.ny;
param.dz = param.sz/param.nz;
param.du = param.su/param.nu;
param.dv = param.sv/param.nv;

param.off_u = 0; param.off_v = 0; % detector rotation shift (real size)

% % % Geometry calculation % % %
param.xs = [-(param.nx-1)/2:1:(param.nx-1)/2]*param.dx;
param.ys = [-(param.ny-1)/2:1:(param.ny-1)/2]*param.dy;
param.zs = [-(param.nz-1)/2:1:(param.nz-1)/2]*param.dz;

param.us = (-(param.nu-1)/2:1:(param.nu-1)/2)*param.du + param.off_u;
param.vs = (-(param.nv-1)/2:1:(param.nv-1)/2)*param.dv + param.off_v;

param.interptype = 'linear'; % 'linear', 'nearest', 'cubic'

%% Measurement read

load Measurement.mat % proj read
    

%% Recon case 2 -  Iterative reconstruction: MLEM

img = ones(param.nx, param.ny, param.nz,'single');
Norimg = CTbackprojection(ones(param.nu, param.nv, param.nProj, 'single'), param);

for iter = 1:50
    tic;
    proj_ratio = proj./CTprojection(img,param);
    proj_ratio(isnan(proj_ratio)) = 0;
    proj_ratio(isinf(proj_ratio)) = 0;
    
    img_ratio = CTbackprojection(proj_ratio, param)./Norimg;
    img_ratio(isnan(img_ratio)) = 0;
    img_ratio(isinf(img_ratio)) = 0;
    
    img = img.*img_ratio;
    toc;
    figure(3); imagesc(max(img(:,:,round(end/2)),0)); axis off; axis equal; colormap gray; colorbar;
    title(['Iteration - ',num2str(iter)]);
    pause(0.1);
end









Contact us