【图像重建】基于布雷格曼迭代（bregman alteration）算法集合ART算法实现CT图像重建附matlab代码

1 简介

Fluorescence diffuse optical tomography (fDOT) is a noninvasive imaging technique that makes it possible to quantify the spatial distribution of fluorescent tracers in small animals. fDOT image reconstruction is commonly performed by means of iterative methods such as the algebraic reconstruction technique (ART). The useful results yielded by more advanced l1-regularized techniques for signal recovery and image reconstruction, together with the recent publication of Split Bregman (SB) procedure, led us to propose a new approach to the fDOT inverse problem, namely, ART-SB. This method alternates a cost-efficient reconstruction step (ART iteration) with a denoising filtering step based on minimization of total variation of the image using the SB method, which can be solved efficiently and quickly. We applied this method to simulated and experimental fDOT data and found that ART-SB provides substantial benefits over conventional ART.

2 部分代码

% Demo_ART_SB_Reconstruction.m% % Demo for a novel iterative reconstruction method that alternates a% computationally efficient linear solver (ART) with a fast denoising step% based on the Split Bregman formulation and that is used to reconstruct% fluorescence diffuse optical tomography data, as in the paper J% Chamorro-Servent, J F P J Abascal, J Aguirre, S Arridge, T Correia, J% Ripoll, M Desco, J J Vaquero. Use of Split Bregman denoising for% iterative reconstruction in fluorescence diffuse optical tomography. J% Biomed Opt, 18(7):076016, 2013. http://dx.doi.org/10.1117/1.JBO.18.7.076016       % % This work proposes to combine computationally efficient solver and% denoising steps with the potential to handle large scale problems in a% fast an efficient manner. Here we combine Algebraic Reconstruction% Technique (ART) and Split Bregman denoising but these methods can be% substituted by the method of choice. %% Our particular choices are explained as follows: %    ART has been widely used in tomography to handle large data sets, as% it can work with one data point (projection) at a time. % %    The Split Bregman formulation allows to solve the total variation% denoising problem in a fast and computationally efficient way. At each% iteration the solution is given by analytical formulas. A linear system  % can be solved in the Fourier domain, keeping the size of the problem (the% Hessian) equal to the number of voxels in the image, or by using% Gauss-Seidel method, which exploits the block diagonal structure of the% Hessian. %% In this demo we solve the image reconstruction problem of fuorescence% diffuse optical tomography but it can be applied to other imaging% modalities. %% % If you use this code, please cite Chamorro-Servent et al. Use of Split% Bregman denoising for iterative reconstruction in fluorescence diffuse% optical tomography. J Biomed Opt, 18(7):076016, 2013.% http://dx.doi.org/10.1117/1.JBO.18.7.076016        %% Judit Chamorro-Servent, Juan FPJ Abascal, Juan Aguirre% Departamento de Bioingenieria e Ingenieria Aeroespacial% Universidad Carlos III de Madrid, Madrid, Spain% juchamser@gmail.com, juanabascal78@gmail.com, desco@hggm.es % READ SIMULATED DATA % Load data, Jacobian matrix and target imageload('DataRed','data','JacMatrix','uTarget');%% Replace the previous line by the following for a larger matrix size% load('Data','data','JacMatrix','uTarget');% but it requires to download data set from the following link: % https://www.dropbox.com/s/461ub2ps0ur8uvd/Data.mat?dl=0% Actually the reduced Jacobian matrix (from DataRed), 700x4000, has been% obtained from a larger matrix (from Data), 6561x4000, by random selection% of rows, providing similar results!  [nr nc]     = size(JacMatrix);% Display target image in the discretized domainN           = [20 20 10];   x           = linspace(-10,10,N(1));y           = linspace(-10,10,N(2));z           = linspace(0,10,N(3));[X,Y,Z]     = meshgrid(x,y,z);% Uncomment below to get a plot of all slices% h           = Plot2DMapsGridSolution(reshape(uTarget,N),X,Y,Z,3); colorbar;rand('seed',0);% Add Gaussian noisedata       = data + 0.01*max(data)*randn(length(data),1);% Randomized index of Jacobian rows for ART reconstruction% indRand     = randperm(length(data));% indRand     = 1:nr; % Uncomment to remove randomized index% -------------------------------------------------------------------------% ART%% It requires selecting few parameters: %    numIterART = number of iterations. The more iterations the better fit% of the data (in this case it converges in 30 iterations) %    relaxParam = relaxation parameter (between 0 and 1). Choose 0.1 to% remove noise and 0.9 to get a closer fit of the datanumIterART  = 30;   % relaxParam  = 0.1;u0          = zeros(nc,1);fprintf('Solving ART reconstruction ... (it takes around 1 s)\n');tic; % uART = ARTReconstruction(JacMatrix,data,relaxParam,numIterART,u0); % The following version runs 3-4 times faster; must try with larger% matrices to assess the differenceuART = ARTReconstruction_Fast(JacMatrix,data,relaxParam,numIterART,u0); tocuART     = reshape(uART,N);%  h = Plot2DMapsGridSolution(uART,X,Y,Z,3); colorbar;       % -------------------------------------------------------------------------% ART-SB  %% It requres selecting the following paramters for ART and Split Bregman % denoising:   %    ART (see previous section). Here we provide high relaxation for a% better fit of the data, as noise will be removed in the denoising step) %%    SPLIT BREGMAN DENOISING:% Split Bregman parameters (mu, lambda, alpha) allow to tune the algorithm% if needed but it is generally quite robust and selecting all parameters% to 1 (or as suggested in the paper) works fine. Then, the only parameters% that need to be chosen are the number of iterations %%    mu     = weight of the fidelity term. (Generally 0.1 to 1 works fine;%    see the paper for more details). The larger the more weight to the%    noisy image%    lamdba = weight of the functionals that impose TV constraints. The% larger the higher the regularization. (Usually chosen larger than mu and% any value around 1 works fine) %    alpha  = the weight of the functional that imposes TV sparsity% (L1-norm). No need to tune this parameter, value 1 should work%    nInner = Number of inner iterations that impose TV constraints. The% higher the number of iterations the more is imposed the TV constraint%    nOuter = Number of outer iterations that impose data fidelity% constraint. Choose this larger than 1 (2 works fine), as TV may provide% an output image with lower constrast and the Bregman iteration can% correct for that. numIter     = 30;numIterART  = 10;      relaxParam  = 0.9;uARTSB      = zeros(N);mu          = 0.3;lambda      = 2*mu;alpha       = 1;nInner      = 5;nOuter      = 2; % indRand = randperm(length(data));fprintf('Solving ART-SB reconstruction ... (it takes around 12 s)\n');h = waitbar(0,'Solving ART-SB reconstruction') ;ticfor it = 1:numIter    % ART reconstruction step: Iterative linear solver    sol = ARTReconstruction_Fast(JacMatrix,data,relaxParam,numIterART,uARTSB(:));     %sol = ARTReconstruction(JacMatrix,data,relaxParam,numIterART,uARTSB(:));     solGrid     = reshape(sol,N);       % SB denoising step    uARTSB      = TV_SB_denoising_3D(solGrid,mu,lambda,alpha,nInner,nOuter);        % Uncomment below to do 2D slice-by-slice smoothing instead. It takes    % similar time but it could be parallelized, which can be faster in some    % applications for large scale problems (now SB denoising takes less    % than a second), by changing the for to a parfor loop  %     for iz = 1:N(3)%         uARTSB(:,:,iz) = TV_SB_denoising_2D(solGrid(:,:,iz),mu,lambda,nInner,nOuter);%     end    % Compute solution error norm    err(it) = norm(uARTSB(:)-uTarget(:))/norm(uTarget(:));        waitbar(it/numIter);end % ittocclose(h);% Display resultsfigure; plot(err); ylabel('Solution error'); xlabel('Number of iterations');title('Convergence of ART-SB');% Target imagePlot2DMapsGridSolution(reshape(uTarget,N),X,Y,Z,3); set(gcf,'name','TARGET','numbertitle','off'); colormap gray;% Reconstructed images0Plot2DMapsGridSolution(uART,X,Y,Z,3); set(gcf,'name','ART reconstruction','numbertitle','off') colormap gray;Plot2DMapsGridSolution(uARTSB,X,Y,Z,3); set(gcf,'name','ART-SB reconstruction','numbertitle','off') colormap gray;  % -------------------------------------------------------------------------    %

3 仿真结果

4 参考文献

[1]马敏, 王涛, 范广永. 基于分割Bregman迭代的内-外置电极电容层析成像[J]. 科学技术与工程, 2018, 018(032):212-216.

[2] Chamorro-Servent J , Abascal J F , Aguirre J , et al. Use of Split Bregman denoising for iterative reconstruction in fluorescence diffuse optical tomography[J]. Journal of Biomedical Optics, 2013, 18.

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

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