【图像重建】基于先验和运动的重建 (PRIMOR)实现口腔CT图像重建附matlab代码和论文

1 简介

Low-dose protocols for respiratory gating in cardiothoracic small-animal imaging lead to

streak artifacts in the images reconstructed with a Feldkamp-Davis-Kress (FDK) method.

We propose a novel prior- and motion-based reconstruction (PRIMOR) method, which

improves prior-based reconstruction (PBR) by adding a penalty function that includes a

model of motion. The prior image is generated as the average of all the respiratory gates,

reconstructed with FDK. Motion between respiratory gates is estimated using a nonrigid

registration method based on hierarchical B-splines. We compare PRIMOR with an equiva

lent PBR method without motion estimation using as reference the reconstruction of high

dose data. From these data acquired with a micro-CT scanner, different scenarios were

simulated by changing photon flux and number of projections. Methods were evaluated in

terms of contrast-to-noise-ratio (CNR), mean square error (MSE), streak artefact indicator

(SAI), solution error norm (SEN), and correction of respiratory motion. Also, to evaluate the

effect of each method on lung studies quantification, we have computed the Jaccard similar

ity index of the mask obtained from segmenting each image as compared to those obtained

from the high dose reconstruction. Both iterative methods greatly improved FDK reconstruc

tion in all cases. PBR was prone to streak artifacts and presented blurring effects in bone

and lung tissues when using both a low number of projections and low dose. Adopting PBR

as a reference, PRIMOR increased CNR up to 33% and decreased MSE, SAI and SEN up

to 20%, 4% and 13%, respectively. PRIMOR also presented better compensation for respi

ratory motion and higher Jaccard similarity index. In conclusion, the new method proposed

for low-dose respiratory gating in small-animal scanners shows an improvement in image

quality and allows a reduction of dose or a reduction of the number of projections between

two and three times with respect to previous PBR approaches.

2 部分代码

​%% The following lines provides the parameters used to generate the matrix% operator, G, using IRT:% n_m         = 1;% numPerProj  = 350;% N           = [350 350];% n_x         = N(1);% numProj     = [360]% totalAngle  = [360];% binning     = 4;% n_ang = numProj;% det_z_count = 1; detector_count= n_x; pixel_count=n_x*n_x;% ds          = 50e-4*binning; % pixel sixe= min pixel*binning% dx          = ds/1.6;% cg          = sino_geom('fan','nb',n_x, 'na', n_ang, ...%     'ds',ds, 'orbit',-totalAngle, ...%     'offset_s', 0.0, ...%     'dsd', 35.2, 'dod', 13.2, 'dfs', inf);% ig          = image_geom('nx', n_x, 'ny', n_x,'dx', dx);% G           = Gtomo2_dscmex(cg, ig);% geom        = fbp2(cg, ig);​% FORWARD OPERATORif 0    % It requires IRT software    % Load Forward Operator precomputed in Linux system. If it does not    % work, run the code above     load('G_linux','geom','G');                                         end​% READ SIMULATED DATA %% HIGH DOSE DATA (ideal data). % Target image acquired from a high dose protocol (dosis four times% the dosis used for static imaging)load('image_HighDose','uTarget');N           = size(uTarget);​% Display four respiratory gatesfigure;for it = 1:4    imagesc(uTarget(:,:,it)); colormap gray; axis image; colorbar;    title(['High dose, gate ' num2str(it)]);    pause(0.3);end​% STANDARD DOSE DATA: Standard dose for static imaging leads to irregularly% distributed % Data for different scenarios can be found at http://dx.doi.org/10.5281/zenodo.15685nameData    = 'data_I0By6_120p';         % Dose reduction by 6​load(nameData,'dataAll');Nd          = size(dataAll);​% Undersampling patternRAll        = dataAll>0;% -------------------------------------------------------------------------% PRIOR image%% Prior image as the average imageRsum        = sum(RAll,3);dataAllAv   = sum(dataAll,3);dataAllAv(Rsum>0) = dataAllAv(Rsum>0)./Rsum(Rsum>0);if 0     % To compute the reference image download the IRT code and run this    % part    uref     = fbp2(dataAllAv, geom); uref(uref<0) = 0;    gaussianFilter = fspecial('gaussian', [5, 5], 3); % [7 7], 5 stronger filter if the prior edges is a problem    uref_aux = imfilter(uref, gaussianFilter, 'symmetric', 'conv');    uref        = uref_aux;    else    % Load precomputed prior image    load('RefImage_I0By6_120p','uref');end% -------------------------------------------------------------------------% FBP reconstruction using IRT softwareframeThis   = 1;if 0     % FDK reconstruction, it requires IRT code     im          = fbp2(dataAll(:,:,frameThis), geom); im(im<0) = 0;else    % Load precomputed FDK reconstruction        load('FDK_I0By6_120p','im');endfigure;subplot(2,2,1);imagesc(uTarget(:,:,frameThis)); title('Target gate 1 (high dose image)'); colorbar;subplot(2,2,2);imagesc(dataAll(:,:,frameThis)); title('Low dose data (120 proj., dose by 6)'); colorbar;subplot(2,2,3);imagesc(im*prod(Nd(1:2))/nnz(RAll(:,:,frameThis))); title('FDK reconstruction'); colorbar; axis image;subplot(2,2,4);imagesc(uref(:,:,frameThis)); title('Prior image: sum of four gates'); colorbar; colormap gray;​% -------------------------------------------------------------------------% Create the support (circle)[n_x_u,n_y_u,n_z_u] = size(uTarget);[X,Y]       = meshgrid(linspace(1,n_x_u,n_x_u),linspace(1,n_x_u,n_x_u));X           = X-n_x_u/2;Y           = Y-n_x_u/2;indCir      = find(sqrt(X.^2+Y.^2)<=n_x_u/2);mask        = zeros(N(1:2));mask(indCir)= 1;​% Apply mask to the prior imageuref        = uref.*mask;% -------------------------------------------------------------------------% PBR methodif 0    % To run PBR method, run this part    matlabpool(4); % comment if matlabpool not available    mu          = 2;    lambda      = 1;    nBreg       = 25;    alpha       = 0.4;    beta        = 0.2;    [uPbr,errPbr] = PBR_CT(G,dataAll,RAll,N,uref,mu,lambda,alpha,beta,nBreg,uTarget);        matlabpool close; else    % Load PBR result already computed    load('PBR_Rec_I0By6_120p','uPbr','errPbr');end% -------------------------------------------------------------------------% PRIMOR method​% Registration stepif 0    % To compute the registration, download the FFD registration software    % and run this part (it takes ~2min)        % Compute the registration step    matlabpool(8); % comment if matlabpool not available        % Parameters for the registration step    Options.Penalty     = 1e-4;    Options.Registration= 'NonRigid';    Options.MaxRef      = 3;    Options.Similarity  = 'sd';        % Estimate the temporal operator, computed by the registration of    % consecutive gates from previous reconstruction (we used PBR)    u0          = uPbr;    [GridAll,SpacingAll,GridDAll,SpacingDAll] = ComputeSplineRegTwoDirectionsPoolDouble(u0,Options);    TParameters.GridAll         = GridAll;    TParameters.SpacingAll      = SpacingAll;    TParameters.GridDAll        = GridDAll;    TParameters.SpacingDAll     = SpacingDAll;    matlabpool close; else    % Load registration result already computed    load('TParameters_I0By6_120p','TParameters');end​% Reconstruction stepmu          = 2;lambda      = 1;alpha       = 0.4;beta        = 0.2;gamma       = 0.5;nBreg       = 50;if 0    % To run PRIMOR method download FFD registration software and IRT    % software and execute this part  (it takes ~ 1h 20min)        matlabpool(4); % comment if matlabpool not available               ​    [uPrimor,errPrimor] = PRIMOR_CT(TParameters,G,dataAll,RAll,N,uref,mu,lambda,gamma,alpha,beta,nBreg,uTarget);    % Reconstructed image and auxiliary variables are displayed for TV and    % prior terms, for some iteration numbers. The number of nonzero    % coefficients on the respective transformed domains are given as a    % precentage    matlabpool close;else    % Load PRIMOR result already computed   load('PRIMOR_Rec_I0By6_120p','uPrimor','errPrimor');end% -------------------------------------------------------------------------% Images for ideal image (high dose) and respiratory gated data with% six-fold dose reduction reconstructed with FDK, PBR and PRIMOR methodsfigure;subplot(2,2,1);imagesc(uTarget(:,:,frameThis)); axis image; axis off; colormap gray;title('FDK, high dose');ca = caxis;subplot(2,2,2);imagesc(im*prod(Nd(1:2))/nnz(RAll(:,:,frameThis))); axis image; axis off; colormap gray;caxis(ca); title('FDK, low dose');subplot(2,2,3);imagesc(uPbr(:,:,frameThis)); axis image; axis off; colormap gray;caxis(ca); title('PBR, low dose');subplot(2,2,4);imagesc(uPrimor(:,:,frameThis)); axis image; axis off; colormap gray;caxis(ca); title('PRIMOR, low dose');​% Zoom imagexZoom        = 120:280;yZoom        = 70:240;​figure;subplot(2,2,1);imagesc(uTarget(xZoom,yZoom,frameThis)); axis image; axis off; colormap gray;title('FDK, high dose');ca = caxis;subplot(2,2,2);imagesc(im(xZoom,yZoom,1)*prod(Nd(1:2))/nnz(RAll(:,:,frameThis))); axis image; axis off; colormap gray;caxis(ca); title('FDK, low dose');subplot(2,2,3);imagesc(uPbr(xZoom,yZoom,frameThis)); axis image; axis off; colormap gray;caxis(ca); title('PBR, low dose');subplot(2,2,4);imagesc(uPrimor(xZoom,yZoom,frameThis)); axis image; axis off; colormap gray;caxis(ca); title('PRIMOR, low dose');​% Convergence: solution error vs iteration numberfigure; plot(mean(errPbr,2));hold on; plot(mean(errPrimor,2),'r');legend('PBR','PRIMOR'); xlabel('Iteration number'); ylabel('Solution error');​% -------------------------------------------------------------------------%

3 仿真结果

4 参考文献

[1] Abascal J ,  Monica A ,  Eugenio M , et al. A Novel Prior- and Motion-Based Compressed Sensing Method for Small-Animal Respiratory Gated CT[J]. Plos One, 2016, 11(3):e0149841.

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

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