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 6load(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.
