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_conv(proj,param )
function [ proj ] = filtering_conv(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

filt_len = param.nu*2;
h = ramp_flat(filt_len);

for iview = 1:param.nProj
    for iv = 1:param.nv
        proj(:,iv,iview) = conv(proj(:,iv,iview),h,'same');    
    end
end

proj = proj/param.du*(2*pi/param.nProj)/2*(param.DSD/param.DSO);

end

function [h, nn] = ramp_flat(n)
nn = [-(n/2):(n/2-1)]';
h = zeros(size(nn),'single');
h(n/2+1) = 1 / 4; 
odd = mod(nn,2) == 1;
h(odd) = -1 ./ (pi * nn(odd)).^2;
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