【图像重建】基于ART算法和SIRT算法实现超声CT反演附MATLAB代码

✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

⛄ 内容介绍

计算机层析成像(即CT)不仅是放射诊断医学领域内的里程碑,而且在现代工业无损检测和勘探等领域中也是重要的研究手段. 滤波反投影法和代数重建法是CT中最为常见的算法,基于对这两种算法的深入分析。

⛄ 部分代码

function [A,b,x,theta,p,d] = paralleltomo(N,theta,p,d,isDisp,isMatrix)

%PARALLELTOMO  Creates a 2D parallel-beam tomography test problem

%

%   [A,b,x,theta,p,d] = paralleltomo(N)

%   [A,b,x,theta,p,d] = paralleltomo(N,theta)

%   [A,b,x,theta,p,d] = paralleltomo(N,theta,p)

%   [A,b,x,theta,p,d] = paralleltomo(N,theta,p,d)

%   [A,b,x,theta,p,d] = paralleltomo(N,theta,p,d,isDisp)

%   [A,b,x,theta,p,d] = paralleltomo(N,theta,p,d,isDisp,isMatrix)

%

% This function uses the "line model" to create a 2D X-ray tomography test

% problem with an N-times-N pixel domain, using p parallel rays for each

% angle in the vector theta.

%

% Input:

%   N         Scalar denoting the number of discretization intervals in

%             each dimesion, such that the domain consists of N^2 cells.

%   theta     Vector containing the projetion angles in degrees.

%             Default: theta = 0:1:179.

%   p         Number of rays for each angle. Default: p = round(sqrt(2)*N).

%   d         Scalar denoting the distance from the first ray to the last.

%             Default: d = p-1.

%   isDisp    If isDisp is non-zero it specifies the time in seconds

%             to pause in the display of the rays. If zero (the default),

%             no display is shown.

%   isMatrix  If non-zero, a sparse matrix is returned in A (default).

%             If zero, instead a function handle is returned.

%

% Output:

%   A         If input isMatrix is 1 (default): coefficient matrix with

%             N^2 columns and length(theta)*p rows.

%             If input isMatrix is 0: A function handle representing a

%             matrix-free version of A in which the forward and backward

%             operations are computed as A(x,'notransp') and A(y,'transp'),

%             respectively, for column vectors x and y of appropriate size.

%             The size of A can be retrieved using A([],'size'). The matrix

%             is never formed explicitly, thus saving memory.

%             elements in A.

%   b         Vector containing the rhs of the test problem.

%   x         Vector containing the exact solution, with elements

%             between 0 and 1.

%   theta     Vector containing the used angles in degrees.

%   p         The number of used rays for each angle.

%   d         The distance between the first and the last ray.

%

% See also: fancurvedtomo, fanlineartomo, seismictomo, seismicwavetomo,

%           sphericaltomo.

% Reference: A. C. Kak and M. Slaney, Principles of Computerized

% Tomographic Imaging, SIAM, Philadelphia, 2001.

% Code written by: Per Christian Hansen, Jakob Sauer Jorgensen, and 

% Maria Saxild-Hansen, 2010-2017.

% This file is part of the AIR Tools II package and is distributed under

% the 3-Clause BSD License. A separate license file should be provided as

% part of the package. 

% 

% Copyright 2017 Per Christian Hansen, Technical University of Denmark and

% Jakob Sauer Jorgensen, University of Manchester.

% Default value of the projection angles theta.

if nargin < 2 || isempty(theta)

    theta = 0:179;

end

% Make sure theta is double precison to prevent round-off issues caused by

% single input.

theta = double(theta);

% Default value of the number of rays.

if nargin < 3 || isempty(p)

    p = round(sqrt(2)*N);

end

% Default value of d.

if nargin < 4 || isempty(d)

    d = p-1;

end

% Default illustration: False

if nargin < 5 || isempty(isDisp)

    isDisp = 0;

end

% Default matrix or matrix-free function handle? Matrix

if nargin < 6 || isempty(isMatrix)

    isMatrix = 1;

end

% Construct either matrix or function handle

if isMatrix

    A = get_or_apply_system_matrix(N,theta,p,d,isDisp);

else

    A = @(U,TRANSP_FLAG) get_or_apply_system_matrix(N,theta,p,d,isDisp,...

        U,TRANSP_FLAG);

end

if nargout > 1

    % Create phantom head as a reshaped vector.

    x = phantomgallery('shepplogan',N);

    x = x(:);

    % Create rhs.

    if isMatrix

        b = A*x;

    else

        b = A(x,'notransp');

    end

end

function A = get_or_apply_system_matrix(N,theta,p,d,isDisp,u,transp_flag)

% Define the number of angles.

nA = length(theta);

% The starting values both the x and the y coordinates.

x0 = linspace(-d/2,d/2,p)';

y0 = zeros(p,1);

% The intersection lines.

x = (-N/2:N/2)';

y = x;

% Prepare for illustration

if isDisp

    AA = phantomgallery('smooth',N);

    figure

end

% Deduce whether to set up matrix or apply to input, depending on whether

% input u is given.

isMatrix = (nargin < 6);

if isMatrix

    % Initialize vectors that contains the row numbers, the column numbers

    % and the values for creating the matrix A effiecently.

    rows = zeros(2*N*nA*p,1);

    cols = rows;

    vals = rows;

    idxend = 0;

    II = 1:nA;

    JJ = 1:p;

else

    % If transp_flag is 'size', only return size of operator, otherwise set

    % proper size of output.

    switch transp_flag

        case 'size'

            A = [p*nA,N^2];

            return

        case 'notransp' % Forward project.

            if length(u) ~= N^2

                error('Incorrect length of input vector')

            end

            A = zeros(p*nA,1);

        case 'transp' % Backproject

            if length(u) ~= p*nA

                error('Incorrect length of input vector')

            end

            A = zeros(N^2,1);

    end

    % If u is a Cartesian unit vector then we only want to compute a single

    % row of A; otherwise we want to multiply with A or A'.

    if strcmpi(transp_flag,'transp') && nnz(u) == 1 && sum(u) == 1

        % Want to compute a single row of A, stored as a column vector.

        ell = find(u==1);

        JJ = mod(ell,p);  if JJ==0, JJ = p; end

        II = (ell-JJ)/p + 1;

    else

        % Want to multiply with A or A'.

        II = 1:nA;

        JJ = 1:p;

    end

end

% Loop over the chosen angles.

for i = II

    % Illustration of the domain.

    if isDisp

        clf

        pause(isDisp)

        imagesc((-N/2+.5):(N/2-0.5),(-N/2+.5):(N/2-0.5),flipud(AA))

        colormap gray

        axis xy

        hold on

        axis equal

        axis(0.7*[-N N -N N])

    end

    % All the starting points for the current angle.

    x0theta = cosd(theta(i))*x0-sind(theta(i))*y0;

    y0theta = sind(theta(i))*x0+cosd(theta(i))*y0;

    % The direction vector for all rays corresponding to the current angle.

    a = -sind(theta(i));

    b = cosd(theta(i));

    % Loop over the rays.

    for j = JJ

        % Use the parametrisation of line to get the y-coordinates of

        % intersections with x = constant.

        tx = (x - x0theta(j))/a;

        yx = b*tx + y0theta(j);

        % Use the parametrisation of line to get the x-coordinates of

        % intersections with y = constant.

        ty = (y - y0theta(j))/b;

        xy = a*ty + x0theta(j);

        % Illustration of the rays.

        if isDisp

            plot(x,yx,'-','color',[220 0 0]/255,'linewidth',1.5)

            plot(xy,y,'-','color',[220 0 0]/255,'linewidth',1.5)

            set(gca,'Xtick',[],'Ytick',[])

            pause(isDisp)

        end

        % Collect the intersection times and coordinates.

        t = [tx; ty];

        xxy = [x; xy];

        yxy = [yx; y];

        % Sort the coordinates according to intersection time.

        [~,I] = sort(t);

        xxy = xxy(I);

        yxy = yxy(I);

        % Skip the points outside the box.

        I = (xxy >= -N/2 & xxy <= N/2 & yxy >= -N/2 & yxy <= N/2);

        xxy = xxy(I);

        yxy = yxy(I);

        % Skip double points.

        I = (abs(diff(xxy)) <= 1e-10 & abs(diff(yxy)) <= 1e-10);

        xxy(I) = [];

        yxy(I) = [];

        % Calculate the length within cell and determines the number of

        % cells which is hit.

        aval = sqrt(diff(xxy).^2 + diff(yxy).^2);

        col = [];

        % Store the values inside the box.

        if numel(aval) > 0

            % If the ray is on the boundary of the box in the top or to the

            % right the ray does not by definition lie with in a valid cell.

            if ~((b == 0 && abs(y0theta(j) - N/2) < 1e-15) || ...

                    (a == 0 && abs(x0theta(j) - N/2) < 1e-15)       )

                % Calculates the midpoints of the line within the cells.

                xm = 0.5*(xxy(1:end-1)+xxy(2:end)) + N/2;

                ym = 0.5*(yxy(1:end-1)+yxy(2:end)) + N/2;

                % Translate the midpoint coordinates to index.

                col = floor(xm)*N + (N - floor(ym));

            end

        end

        if ~isempty(col)

            if isMatrix

                % Create the indices to store the values to vector for

                % later creation of A matrix.

                idxstart = idxend + 1;

                idxend = idxstart + numel(col) - 1;

                idx = idxstart:idxend;

                % Store row numbers, column numbers and values.

                rows(idx) = (i-1)*p + j;

                cols(idx) = col;

                vals(idx) = aval;

            else

                % If any nonzero elements, apply forward or back operator

                if strcmp(transp_flag,'notransp')

                    % Insert inner product with u into w.

                    A( (i-1)*p+j ) = aval'*u(col);

                else  % Adjoint operator.

                    A(col) = A(col) + u( (i-1)*p+j )*aval;

                end

            end

        end

    end % end j

end % end i

if isMatrix

    % Truncate excess zeros.

    rows = rows(1:idxend);

    cols = cols(1:idxend);

    vals = vals(1:idxend);

    % Create sparse matrix A from the stored values.

    A = sparse(rows,cols,vals,p*nA,N^2);

end

⛄ 运行结果

⛄ 参考文献

[1] 徐茂林. 透射及散射断层成像中若干反演算法的研究[D]. 浙江大学, 2004.

[2] 龙芸, 程久龙. 基于SIRT算法的声波CT反演及工程应用[C]// 2016中国地球科学联合学术年会论文集(三十)——专题54:煤炭资源与矿山地球物理. 中国地球物理学会; 中国地质学会; 中国地震学会; 中国力学学会; 中国岩石力学与工程学会, 2016.

[3] 赵祥模, 李娜, 关可,等. 基于LTI射线追踪的改进的ART算法[J]. 中国图象图形学报A, 2009.

[4] 颜华, 孙灿烽, 王伊凡. 基于改进SIRT算法的声学CT温度场重建[J]. 沈阳工业大学学报, 2021.

[5] 张顺利, 张定华, 王凯,等. 一种基于ART算法的快速图像重建技术[J]. 核电子学与探测技术, 2007(3):5.

⛳️ 完整代码

❤️部分理论引用网络文献，若有侵权联系博主删除

❤️ 关注我领取海量matlab电子书和数学建模资料