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

filtering_fft(proj,param )
function [ proj ] = filtering_fft(proj,param )
%FILTERING Summary of this function goes here
%   Detailed explanation goes here

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

[uu,vv] = meshgrid(-us,vs);

w = (param.DSD)./sqrt((param.DSD)^2+uu.^2 + vv.^2);

for i=1:param.nProj
    proj(:,:,i) = proj(:,:,i).*w';
end

d = 1;

[filt, filt_len] = Filter(param.filter, param.nu, d);

fproj = (zeros(filt_len,param.nv*param.nProj,'single'));
fproj(filt_len/2-param.nu/2+1:filt_len/2+param.nu/2,:) = (reshape(proj,param.nu,param.nv*param.nProj));


fproj = fft(fproj);

filt = repmat(filt',[1 param.nv*param.nProj]);

fproj = fproj.*filt;

fproj = (real(ifft(fproj)));

proj = reshape(fproj(end/2-param.nu/2+1:end/2+param.nu/2,:)/2/param.du*(2*pi/param.nProj)/2*(param.DSD/param.DSO),param.nu,param.nv,param.nProj);
   


end

function [filt, order] = Filter(filter, len, d)

order = max(64,2^nextpow2(2*len));

filt = 2*( 0:(order/2) )./order;
w = 2*pi*(0:size(filt,2)-1)/order;   % frequency axis up to Nyquist 

switch lower(filter)
    case 'ram-lak'
        % Do nothing
    case 'shepp-logan'
        % be careful not to divide by 0:
        filt(2:end) = filt(2:end) .* (sin(w(2:end)/(2*d))./(w(2:end)/(2*d)));
    case 'cosine'
        filt(2:end) = filt(2:end) .* cos(w(2:end)/(2*d));
    case 'hamming'  
        filt(2:end) = filt(2:end) .* (.54 + .46 * cos(w(2:end)/d));
    case 'hann'
        filt(2:end) = filt(2:end) .*(1+cos(w(2:end)./d)) / 2;
    otherwise
        filter
        error('Invalid filter selected.');
end

filt(w>pi*d) = 0;                      % Crop the frequency response
filt = [filt , filt(end-1:-1:2)];    % Symmetry of the filter
return

end

Contact us