【图像重建】基于先验和运动的PRIMOR算法实现呼吸门控 CT图像重建附matlab代码

1 内容介绍

心胸小动物成像中呼吸门控的低剂量方案导致使用 Feldkamp-Davis-Kress (FDK) 方法重建的图像中的条纹伪影。我们提出了一种新颖的基于先验和基于运动的重建（PRIMOR）方法，该方法通过添加一个惩罚函数来改进基于先验的重建 (PBR)运动模型。先验图像是作为所有呼吸门的平均值生成的，用FDK重建。使用非刚性估计呼吸门之间的运动基于层次B样条的配准方法我们将 PRIMOR 与没有运动估计的等效 PBR 方法进行比较，使用高剂量数据。从这些用微型 CT 扫描仪获得的数据中，我们发现了不同的场景通过改变光子通量和投影数量来模拟。方法在对比度噪声比 (CNR)、均方误差 (MSE)、条纹伪影指标(SAI)、解误差范数 (SEN) 和呼吸运动校正。此外，要评估每种方法对肺研究量化的影响，我们计算了从分割每个图像获得的掩码的 Jaccard 相似性指数，与获得的相比较从高剂量重建。在所有情况下，这两种迭代方法都极大地改进了 FDK 重建。 PBR 容易出现条纹伪影并在骨骼中呈现模糊效果和肺组织同时使用少量投射和低剂量。采用 PBR作为参考，PRIMOR 将 CNR 提高了 33%，并降低了 MSE、SAI 和 SEN分别为 20%、4% 和 13%。 PRIMOR 还为呼吸运动提供了更好的补偿和更高的 Jaccard 相似指数。总之，提出的新方法小动物扫描仪中的低剂量呼吸门控显示图像改善质量，并允许减少剂量或减少之间的投影数量是以前 PBR 方法的两倍和三倍。

2 仿真代码

<span style="color:#333333"><span style="background-color:rgba(0, 0, 0, 0.03)"><code><span style="color:#afafaf">%</span> Demo_PBR_PRIMOR_CT.m</code><code><span style="color:#afafaf">%</span> </code><code><span style="color:#afafaf">%</span> Demo <span style="color:#ca7d37">for</span> reconstructing respiratory gated-CT data with a novel prior- and</code><code><span style="color:#afafaf">%</span> motion-based reconstruction (PRIMOR) method proposed <span style="color:#ca7d37">in</span> the publication </code><code><span style="color:#afafaf">%</span> JFPJ Abascal, M Abella, E Marinetto, J Pascau, M Desco. A novel prior-</code><code><span style="color:#afafaf">%</span> and motion-based compressed sensing method <span style="color:#ca7d37">for</span> small-animal respiratory</code><code><span style="color:#afafaf">%</span> gated CT. Plos One, 2016 (<span style="color:#ca7d37">in</span> press). </code><code><span style="color:#afafaf">%</span></code><code>​</code><code>% Faster implementations: </code><code><span style="color:#afafaf">%</span> </code><code><span style="color:#afafaf">%</span> Change forward operator, G, and backprojection operator, G<span style="color:#dd1144">', with your</span></code><code><span style="color:#afafaf">%</span><span style="color:#dd1144"> own GPU implementations to make the algorithm up to orders of magnitude</span></code><code><span style="color:#afafaf">%</span><span style="color:#dd1144"> faster.  </span></code><code><span style="color:#afafaf">%</span></code><code><span style="color:#afafaf">%</span></code><code><span style="color:#afafaf">%</span><span style="color:#dd1144"> If you use this code, please reference the paper JFPJ Abascal et al. A</span></code><code><span style="color:#afafaf">%</span><span style="color:#dd1144"> novel prior- and motion-based compressed sensing method for small-animal</span></code><code><span style="color:#afafaf">%</span><span style="color:#dd1144"> respiratory gated CT. Plos One, 2016 (in press). </span></code><code><span style="color:#afafaf">%</span> </code><code><span style="color:#afafaf">%</span><span style="color:#dd1144"> Juan FPJ Abascal, Monica Abella</span></code><code><span style="color:#afafaf">%</span><span style="color:#dd1144"> Departamento de Bioingenieria e Ingenieria Aeroespacial</span></code><code><span style="color:#afafaf">%</span><span style="color:#dd1144"> Universidad Carlos III de Madrid, Madrid, Spain</span></code><code><span style="color:#afafaf">%</span><span style="color:#dd1144"> mabella@hggm.es, juanabascal78@gmail.com, desco@hggm.es</span></code><code>​</code><code>%<span style="color:#dd1144"> -------------------------------------------------------------------------</span></code><code><span style="color:#afafaf">%</span><span style="color:#dd1144"> READ Data and forward operator G (mapping 2D image to the data)</span></code><code>​</code><code>%<span style="color:#dd1144"> These have been computed using Fessler'</span>s IRT software (J.A. Fessler,</code><code><span style="color:#afafaf">%</span> Image reconstruction toolbox (IRT), 2011, retrieved from</code><code><span style="color:#afafaf">%</span> http://www.eecs.umich.edu/~fessler/code/index.html. </code><code><span style="color:#afafaf">%</span></code><code><span style="color:#afafaf">%</span> The following lines provides the parameters used to generate the matrix</code><code><span style="color:#afafaf">%</span> operator, G, using IRT:</code><code><span style="color:#afafaf">%</span> n_m         = 1;</code><code><span style="color:#afafaf">%</span> numPerProj  = 350;</code><code><span style="color:#afafaf">%</span> N           = [350 350];</code><code><span style="color:#afafaf">%</span> n_x         = N(1);</code><code><span style="color:#afafaf">%</span> numProj     = [360]</code><code><span style="color:#afafaf">%</span> totalAngle  = [360];</code><code><span style="color:#afafaf">%</span> binning     = 4;</code><code><span style="color:#afafaf">%</span> n_ang = numProj;</code><code><span style="color:#afafaf">%</span> det_z_count = 1; detector_count= n_x; pixel_count=n_x*n_x;</code><code><span style="color:#afafaf">%</span> ds          = 50e-4*binning; % pixel sixe= min pixel*binning</code><code><span style="color:#afafaf">%</span> dx          = ds/1.6;</code><code><span style="color:#afafaf">%</span> cg          = sino_geom(<span style="color:#dd1144">'fan'</span>,<span style="color:#dd1144">'nb'</span>,n_x, <span style="color:#dd1144">'na'</span>, n_ang, ...</code><code><span style="color:#afafaf">%</span>     <span style="color:#dd1144">'ds'</span>,ds, <span style="color:#dd1144">'orbit'</span>,-totalAngle, ...</code><code><span style="color:#afafaf">%</span>     <span style="color:#dd1144">'offset_s'</span>, 0.0, ...</code><code><span style="color:#afafaf">%</span>     <span style="color:#dd1144">'dsd'</span>, 35.2, <span style="color:#dd1144">'dod'</span>, 13.2, <span style="color:#dd1144">'dfs'</span>, inf);</code><code><span style="color:#afafaf">%</span> ig          = image_geom(<span style="color:#dd1144">'nx'</span>, n_x, <span style="color:#dd1144">'ny'</span>, n_x,<span style="color:#dd1144">'dx'</span>, dx);</code><code><span style="color:#afafaf">%</span> G           = Gtomo2_dscmex(cg, ig);</code><code><span style="color:#afafaf">%</span> geom        = fbp2(cg, ig);</code><code>​</code><code>% FORWARD OPERATOR</code><code>if 0</code><code>    % It requires IRT software</code><code>    % Load Forward Operator precomputed in Linux system. If it does not</code><code>    % work, run the code above </code><code>    load('G_linux','geom','G');                                         </code><code>end</code><code>​</code><code>% READ SIMULATED DATA </code><code><span style="color:#afafaf">%</span></code><code><span style="color:#afafaf">%</span> HIGH DOSE DATA (ideal data). </code><code><span style="color:#afafaf">%</span> Target image acquired from a high dose protocol (dosis four <span style="color:#ca7d37">times</span></code><code><span style="color:#afafaf">%</span> the dosis used <span style="color:#ca7d37">for</span> static imaging)</code><code>load('image_HighDose','uTarget');</code><code>N           = size(uTarget);</code><code>​</code><code>% Display four respiratory gates</code><code>figure;</code><code>for it = 1:4</code><code>    imagesc(uTarget(:,:,it)); colormap gray; axis image; colorbar;</code><code>    title(['High dose, gate ' num2str(it)]);</code><code>    pause(0.3);</code><code>end</code><code>​</code><code>% STANDARD DOSE DATA: Standard dose <span style="color:#ca7d37">for</span> static imaging leads to irregularly</code><code><span style="color:#afafaf">%</span> distributed </code><code><span style="color:#afafaf">%</span> Data <span style="color:#ca7d37">for</span> different scenarios can be found at http://dx.doi.org/10.5281/zenodo.15685</code><code>nameData    = 'data_I0By6_120p';         % Dose reduction by 6</code><code>​</code><code>load(nameData,'dataAll');</code><code>Nd          = size(dataAll);</code><code>​</code><code>% Undersampling pattern</code><code>RAll        = dataAll>0;</code><code><span style="color:#afafaf">%</span> -------------------------------------------------------------------------</code><code><span style="color:#afafaf">%</span> PRIOR image</code><code><span style="color:#afafaf">%</span></code><code><span style="color:#afafaf">%</span> Prior image as the average image</code><code>Rsum        = sum(RAll,3);</code><code>dataAllAv   = sum(dataAll,3);</code><code><span style="color:#afafaf">dataAllAv(Rsum></span>0) = dataAllAv(Rsum>0)./Rsum(Rsum>0);</code><code>if 0 </code><code>    % To compute the reference image download the IRT code and run this</code><code>    % part</code><code>    uref     = fbp2(dataAllAv, geom); uref(uref<0) = 0;</code><code>    gaussianFilter = fspecial('gaussian', [5, 5], 3); % [7 7], 5 stronger filter if the prior edges is a problem</code><code>    uref_aux = imfilter(uref, gaussianFilter, 'symmetric', 'conv');</code><code>    uref        = uref_aux;    </code><code>else</code><code>    % Load precomputed prior image</code><code>    load('RefImage_I0By6_120p','uref');</code><code>end</code><code><span style="color:#afafaf">%</span> -------------------------------------------------------------------------</code><code><span style="color:#afafaf">%</span> FBP reconstruction using IRT software</code><code>frameThis   = 1;</code><code>if 0 </code><code>    % FDK reconstruction, it requires IRT code </code><code>    im          = fbp2(dataAll(:,:,frameThis), geom); im(im<0) = 0;</code><code>else</code><code>    % Load precomputed FDK reconstruction    </code><code>    load('FDK_I0By6_120p','im');</code><code>end</code><code>figure;</code><code>subplot(2,2,1);</code><code>imagesc(uTarget(:,:,frameThis)); title('Target gate 1 (high dose image)'); colorbar;</code><code>subplot(2,2,2);</code><code>imagesc(dataAll(:,:,frameThis)); title('Low dose data (120 proj., dose by 6)'); colorbar;</code><code>subplot(2,2,3);</code><code>imagesc(im*prod(Nd(1:2))/nnz(RAll(:,:,frameThis))); title('FDK reconstruction'); colorbar; axis image;</code><code>subplot(2,2,4);</code><code>imagesc(uref(:,:,frameThis)); title('Prior image: sum of four gates'); colorbar; colormap gray;</code><code>​</code><code>% -------------------------------------------------------------------------</code><code><span style="color:#afafaf">%</span> Create the support (circle)</code><code>[n_x_u,n_y_u,n_z_u] = size(uTarget);</code><code>[X,Y]       = meshgrid(linspace(1,n_x_u,n_x_u),linspace(1,n_x_u,n_x_u));</code><code>X           = X-n_x_u/2;</code><code>Y           = Y-n_x_u/2;</code><code>indCir      = find(sqrt(X.^2+Y.^2)<=n_x_u/2);</code><code>mask        = zeros(N(1:2));</code><code>mask(indCir)= 1;</code><code>​</code><code>% Apply mask to the prior image</code><code>uref        = uref.*mask;</code><code><span style="color:#afafaf">%</span> -------------------------------------------------------------------------</code><code><span style="color:#afafaf">%</span> PBR method</code><code>if 0</code><code>    % To run PBR method, run this part</code><code>    matlabpool(4); % comment if matlabpool not available</code><code>    mu          = 2;</code><code>    lambda      = 1;</code><code>    nBreg       = 25;</code><code>    alpha       = 0.4;</code><code>    beta        = 0.2;</code><code>    [uPbr,errPbr] = PBR_CT(G,dataAll,RAll,N,uref,mu,lambda,alpha,beta,nBreg,uTarget);    </code><code>    matlabpool close; </code><code>else</code><code>    % Load PBR result already computed</code><code>    load('PBR_Rec_I0By6_120p','uPbr','errPbr');</code><code>end</code><code><span style="color:#afafaf">%</span> -------------------------------------------------------------------------</code><code><span style="color:#afafaf">%</span> PRIMOR method</code><code>​</code><code>% Registration step</code><code>if 0</code><code>    % To compute the registration, download the FFD registration software</code><code>    % and run this part (it takes ~2min)</code><code>    </code><code>    % Compute the registration step</code><code>    matlabpool(8); % comment if matlabpool not available</code><code>    </code><code>    % Parameters for the registration step</code><code>    Options.Penalty     = 1e-4;</code><code>    Options.Registration= 'NonRigid';</code><code>    Options.MaxRef      = 3;</code><code>    Options.Similarity  = 'sd';</code><code>    </code><code>    % Estimate the temporal operator, computed by the registration of</code><code>    % consecutive gates from previous reconstruction (we used PBR)</code><code>    u0          = uPbr;</code><code>    [GridAll,SpacingAll,GridDAll,SpacingDAll] = ComputeSplineRegTwoDirectionsPoolDouble(u0,Options);</code><code>    TParameters.GridAll         = GridAll;</code><code>    TParameters.SpacingAll      = SpacingAll;</code><code>    TParameters.GridDAll        = GridDAll;</code><code>    TParameters.SpacingDAll     = SpacingDAll;</code><code>    matlabpool close; </code><code>else</code><code>    % Load registration result already computed</code><code>    load('TParameters_I0By6_120p','TParameters');</code><code>end</code><code>​</code><code>% Reconstruction step</code><code>mu          = 2;</code><code>lambda      = 1;</code><code>alpha       = 0.4;</code><code>beta        = 0.2;</code><code>gamma       = 0.5;</code><code>nBreg       = 50;</code><code>if 0</code><code>    % To run PRIMOR method download FFD registration software and IRT</code><code>    % software and execute this part  (it takes ~ 1h 20min)</code><code>    </code><code>    matlabpool(4); % comment if matlabpool not available               </code><code>​</code><code>    [uPrimor,errPrimor] = PRIMOR_CT(TParameters,G,dataAll,RAll,N,uref,mu,lambda,gamma,alpha,beta,nBreg,uTarget);</code><code>    % Reconstructed image and auxiliary variables are displayed for TV and</code><code>    % prior terms, for some iteration numbers. The number of nonzero</code><code>    % coefficients on the respective transformed domains are given as a</code><code>    % precentage</code><code>    matlabpool close;</code><code>else</code><code>    % Load PRIMOR result already computed</code><code>   load('PRIMOR_Rec_I0By6_120p','uPrimor','errPrimor');</code><code>end</code><code><span style="color:#afafaf">%</span> -------------------------------------------------------------------------</code><code><span style="color:#afafaf">%</span> Images <span style="color:#ca7d37">for</span> ideal image (high dose) and respiratory gated data with</code><code><span style="color:#afafaf">%</span> six-fold dose reduction reconstructed with FDK, PBR and PRIMOR methods</code><code>figure;</code><code>subplot(2,2,1);</code><code>imagesc(uTarget(:,:,frameThis)); axis image; axis off; colormap gray;</code><code>title('FDK, high dose');</code><code>ca = caxis;</code><code>subplot(2,2,2);</code><code>imagesc(im*prod(Nd(1:2))/nnz(RAll(:,:,frameThis))); axis image; </code><code>axis off; colormap gray;</code><code>caxis(ca); title('FDK, low dose');</code><code>subplot(2,2,3);</code><code>imagesc(uPbr(:,:,frameThis)); axis image; axis off; colormap gray;</code><code>caxis(ca); title('PBR, low dose');</code><code>subplot(2,2,4);</code><code>imagesc(uPrimor(:,:,frameThis)); axis image; axis off; colormap gray;</code><code>caxis(ca); title('PRIMOR, low dose');</code><code>​</code><code>% Zoom image</code><code>xZoom        = 120:280;</code><code>yZoom        = 70:240;</code><code>​</code><code>figure;</code><code>subplot(2,2,1);</code><code>imagesc(uTarget(xZoom,yZoom,frameThis)); axis image; axis off; colormap gray;</code><code>title('FDK, high dose');</code><code>ca = caxis;</code><code>subplot(2,2,2);</code><code>imagesc(im(xZoom,yZoom,1)*prod(Nd(1:2))/nnz(RAll(:,:,frameThis))); axis image; axis off; colormap gray;</code><code>caxis(ca); title('FDK, low dose');</code><code>subplot(2,2,3);</code><code>imagesc(uPbr(xZoom,yZoom,frameThis)); axis image; axis off; colormap gray;</code><code>caxis(ca); title('PBR, low dose');</code><code>subplot(2,2,4);</code><code>imagesc(uPrimor(xZoom,yZoom,frameThis)); axis image; axis off; colormap gray;</code><code>caxis(ca); title('PRIMOR, low dose');</code><code>​</code><code>% Convergence: solution error vs iteration number</code><code>figure; plot(mean(errPbr,2));</code><code>hold on; </code><code>plot(mean(errPrimor,2),'r');</code><code>legend('PBR','PRIMOR'); </code><code>xlabel('Iteration number'); ylabel('Solution error');</code><code>​</code><code>% -------------------------------------------------------------------------</code><code><span style="color:#afafaf">%</span></code></span></span>

3 运行结果

4 参考文献

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

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