毕业设计&课设-与CT重建和模拟相关的Matlab代码.zip

function bp = ctBackProject( sinogram, thetas, dSize, cx, cy, Nx, Ny, ... dx, dy, varargin ) % Written by Nicholas Dwork % sinogram: 2D array, each row corresponds to a theta, each column % to a detector % thetas: a 1D array of thetas corresponding to each sinogram row % dSize: the size of each detector in meters % cx: center x position of reconstructed region % cy: center y position of reconstructed region % xFOV: horizontal field of view in meters % yFOV: vertical field of view in meters % Nx: number of pixels horizontally in reconstruction % Ny: number of pixels vertically in reconstruction % dx: horizontal size of pixel in meters % dy: vertical size of pixel in meters % type (optional): 'iso' or 'fast' % 'iso' uses a rotation that's an isometry % 'fast' faster implementation defaultType = 'fast'; expectedTypes = { 'iso', 'fast' }; p = inputParser; p.addRequired('sinogram', @(x) ndims(x)==2); p.addRequired('thetas', @(x) ndims(x)==2); p.addRequired('dSize',@isnumeric); p.addRequired('cx',@isnumeric); p.addRequired('cy',@isnumeric); p.addRequired('Nx',@isnumeric); p.addRequired('Ny',@isnumeric); p.addRequired('dx',@isnumeric); p.addRequired('dy',@isnumeric); p.addOptional('type',defaultType, ... @(x) any(validatestring(x,expectedTypes)) ); p.parse( sinogram, thetas, dSize, cx, cy, Nx, Ny, ... dx, dy, varargin{:} ); inputs = p.Results; type = inputs.type; if strcmp(type, 'fast') bp = ctBackProjectFast( sinogram, thetas, dSize, cx, cy, Nx, Ny, ... dx, dy ); elseif strcmp(type, 'iso') bp = ctBackProjectIso( sinogram, thetas, dSize, cx, cy, Nx, Ny, ... dx, dy ); end end function bp = ctBackProjectIso( sinogram, thetas, dSize, cx, cy, Nx, Ny, ... dx, dy ) if dx ~= dy err = MException('ctCodes:ctBackProjectIso', ... 'dx and dy must be equal for this version of back projection'); throw(err); end [nThetas, nDetectors] = size(sinogram); dOffset = 0; % detector center offset dLocs = ( [0:nDetectors-1] - floor(0.5*nDetectors) ) * dSize - dOffset; % create the adjoint of the interpolation matrix if mod( Nx, 2 )==0 locs = ( ([0:Nx-1]) - 0.5*Nx + 0.5 ) * dx; else locs = ( ([0:Nx-1]) - floor(0.5*Nx) ) * dx; end B = zeros(nDetectors, Nx); for i=1:Nx tmp = zeros(1,Nx); tmp(i) = 1; interped = interp1( locs, tmp, dLocs,'linear',0 ); B(:,i) = interped; end Bt = B'; bp = zeros(Ny,Nx); parfor thIndx = 1:nThetas theta = thetas( thIndx ); line = sinogram(thIndx,:); tmp = Bt * line'; smeared = ones(Ny,1) * tmp'; bp = bp + isoRot( smeared, -theta ); if mod(thIndx,5)==0 disp([ 'ctBackProject Theta Indx / Theta: ', ... num2str(thIndx), '/', num2str(theta) ]); end end end function bp = ctBackProjectFast( sinogram, thetas, dSize, cx, cy, Nx, Ny, ... dx, dy ) [nThetas, nDetectors] = size(sinogram); dOffset = 0; % detector center offset dLocs = ( [0:nDetectors-1] - floor(0.5*nDetectors) ) * dSize - dOffset; % Make arrays of x and y positions of each pixel if mod( Nx, 2 )==0 lineXs = ( ([0:Nx-1]) - 0.5*Nx + 0.5 ) * dx + cx; else lineXs = ( ([0:Nx-1]) - floor(0.5*Nx) ) * dx + cx; end if mod( Ny, 2 )==0 lineYs = ( ([0:Ny-1]) - 0.5*Ny + 0.5 ) * dy + cy; else lineYs = ( ([0:Ny-1]) - floor(0.5*Ny) ) * dy + cy; end xs = ones(Ny,1) * lineXs; ys = lineYs' * ones(1,Nx); xs=xs(:) + cx; ys=ys(:) + cy; angles = atan2(ys,xs); pixDs = sqrt( xs.*xs + ys.*ys ); bp = zeros(Ny,Nx); parfor thIndx = 1:nThetas theta = thetas( thIndx ); projections = pixDs .* cos( angles - theta ); interped = interp1( dLocs, sinogram(thIndx,:), projections, ... 'linear', 0); bp = bp + reshape(interped,Ny,Nx); if mod(thIndx,10)==0 disp([ 'ctBackProject Theta Indx / Theta: ', ... num2str(thIndx), '/', num2str(theta) ]); end end end