运用滤波反投影的方法对图像进行重建matlab仿真

1.算法描述

   直接由正弦图得到反投影图像，会存在严重的模糊，这是早期 CT 系统所存在的问题。傅立叶中心切片定理表明，投影的一维傅立叶变换是得到投影区域的二维傅立叶变换的切片。滤波反投影重建算法在反投影前将每一个采集投影角度下的投影进行卷积处理，从而改善点扩散函数引起的形状伪影，有效地改善了重建的图像质量。

   括号 [] 内部是一个一维傅里叶反变换，可以认为这是一个一维滤波器的传递函数。由于 ∣ ω ∣ |\omega|∣ω∣ 是一个不可积的斜坡函数（Slope function），可以通过对斜坡加窗进行限制，典型地如汉明窗（Hamming window）、韩窗（Hann window）。该式也可以使用空间卷积来实现：

   这表明，将对应的投影 g ( ρ , θ ) g(\rho, \theta)g(ρ,θ) 与斜坡滤波器传递函数 s ( ρ ) s(\rho)s(ρ) 的傅里叶反变换进行卷积，可以得到角度 θ \thetaθ 的各个反投影，整个反投影图像可以通过对所有反投影图像积分得到。

   滤波反投影 ，运用滤波反投影的方法对图像进行重建matlab仿真+程序操作视频 对比包括 Ram-Lak (Ramp),Shepp-Logan,Ram-Lak Cosine,Ram-Lak Hamming,Ram-Lak Hann。

Ram-Lak

   斜坡滤波器（甚至在被加窗时）在频率域的直流项为零，故每一幅反投影图像的均值将为零。这将意味着，每一幅反投影图像都将有正像素和负像素值，当所有的反投影图像相加形成最终的重建图像时，一些负像素值位置可能变成正像素，而平均值可能不为零，但是，典型地，最终的图像将还是有负像素值。Ram-Lak滤波函数如下所示： 

Shepp-Logan

   Shepp-Logan头模型是由L. A. Shepp和B. F. Logan于1974年首次提出的，它由10 个位置、大小、方向、密度各异的椭圆叠加而成，模拟一个脑部的断层。Shepp-Logan滤波函数如下所示： 

2.仿真效果预览
matlab2022a仿真结果如下：

3.MATLAB部分代码预览

    switch list1_val
      case 'None'
        filter_val = 'None';
      case 'Ram-Lak (Ramp)'
        filter_val = 'Ram-Lak';
      case 'Shepp-Logan'
        filter_val = 'Shepp-Logan';
      case 'Ram-Lak Cosine'
        filter_val = 'Cosine';
      case 'Ram-Lak Hamming'
        filter_val = 'Hamming';
      case 'Ram-Lak Hann'
        filter_val = 'Hann';
      case 'Special'
        filter_val = 'Special';
    end
    
    switch list2_val
      case 'Linear'
        interpol_val = 'linear';
      case 'Nearest'
        interpol_val = 'nearest';
      case 'Spline'
        interpol_val = 'spline';
      case 'Cubic'
        interpol_val = 'pchip';
    end
 
    theta = round(0 : 180/antal_projektioner : 179);
    del_sinogram = data.sinogram(:, theta+1);
 
    if strcmp(filter_val, 'Special') == 1
      if strcmp(file_name(end-4:end), '3.mat') == 1
        data.SB = data.SB + 0.3 * rand(size(data.SB));
      end
      
      if strcmp(file_name(end-4:end), 'f.mat') == 1
        filter_val = 'Hamming';
        rec_img = iradon(del_sinogram, theta, interpol_val, filter_val);
      else
        FFT_sinogram = fftshift(fft(del_sinogram),1);
        FFT_sinogram = FFT_sinogram .* data.SB(:, theta+1);
        del_sinogram = ifft(ifftshift(FFT_sinogram,1));
        rec_img = iradon(del_sinogram, theta, interpol_val);
      end
    else
      rec_img = iradon(del_sinogram, theta, interpol_val, filter_val);
    end
    
    figure, imagesc(rec_img), colormap(pink)
    title([data.txt ': ' num2str(antal_projektioner) ' projections, ' list1_val ' filter, ' list2_val ' interpolation'])
  end
A_023