m基于小波变换的PET图像重建matlab仿真

1.算法仿真效果
matlab2022a仿真结果如下：

2.算法涉及理论知识概要

   PET测量数据（光子计数值）与生理模型紧密相关。然而，通常的重建方法多采用无生理意义的统计特性约束，从而忽略了反映生理特性的生理模型约束。本项目针对PET图像重建问题，提出利用H无穷粒子滤波混合算法融入数据统计特性模型、生理模型双重约束，并结合GPU完成重建过程的加速的方法。具体实现过程中以先验的共性生理特性作为约束，利用H无穷的鲁棒性处理不确定性，给出放射性药物浓度的粗估计，然后利用粒子滤波可以融合任何统计特性的能力，细调节估计值。以上迭代循环进行，最后给出放射性浓度分布与方差特性。由于采用了GPU加速策略，重建时间将显著缩短。该方法的实现将为PET图像重建提供新的技术手段。

   PET(Positron Emission Tomography)因其灵敏度高、特异性好、全身显像及安全性好等优势得到了越来越多研究者的关注,同时也是继x射线断层成像(x-ray computed tomography,X-CT)和核磁共振(magnetic resonance imaging,MRI)成像技术之后发展起来的一项需要高技术和高成本同时也是最先进和极其重要的医学领域成像技术。但是,在实际PET采样数据获取的过程中由于受到PET采样系统硬件和软件的限制以及伽马光子在体内传输过程中相互影响,会导致采集到的PET投影数据中往往含有大量噪声,从而会影响PET图像重建的质量。已知当重建出的图像中含有噪声或重建出的PET图像质量不高时会对后续医生的诊断有影响,则在PET成像技术的临床应用中如何通过设计重建算法重建出质量高的图像也是极其重要且非常有研究价值的一个步骤。

对应的代码简要说明：

这个部分对应的小波分解，通过一维快速变换得到最终的二维变换效果。

这个就是PET迭代重构。

最后的误差曲线。

3.MATLAB核心程序

I    = imread('shepp-logan.jpg');            
N    = 1;%分解级数
NUM  = 90;%迭代次数
if isrgb(I) == 1
   im = mat2gray(rgb2gray(I));
else
   im = mat2gray(I);
end
im2   = im;
[r,c] = size(im2);
 
 
figure;
subplot(131);
imshow(im,[]);
title('原始的图像');    
 
 
%得到行值与列值
[rows,cols]=size(im);                                                
 
 
for i=1:N      
    %先对行进行一维快速分解
    for j=1:cols/2^(i-1)             
        im(j,1:cols/2^(i-1)) = func_fast_wavelet_trans(im(j,1:cols/2^(i-1)),1);       
    end
    %对列进行一维快速分解
    for j=1:rows/2^(i-1)                                           
        im(1:rows/2^(i-1),j) = func_fast_wavelet_trans(im(1:rows/2^(i-1),j)',1);      
    end
end
%赋值给输出变量                                                          
im(abs(im)<0.1)=0;
subplot(132);
imshow(im)                                                      
title('变换之后的图像')
wl     = im;
im_tmp = im;
%计算概率矩阵A
%计算概率矩阵A
%计算概率矩阵A
for jj = 1:NUM
    im = im_tmp;
    jj
    % 进行重建
    % 进行重建
    for i=N:-1:1                                                         
       for j=rows/2^(i-1):-1:1                                   
            im(1:rows/2^(i-1),j) = func_image_recon(im(1:rows/2^(i-1),j)',1,jj,im2(1:rows/2^(i-1),j)',wl(1:rows/2^(i-1),j)');      
        end
       for j=cols/2^(i-1):-1:1                                        
            im(j,1:cols/2^(i-1)) = func_image_recon(im(j,1:cols/2^(i-1)),1,jj,im2(1:rows/2^(i-1),j)',wl(1:rows/2^(i-1),j)');       
        end
    end
 
    subplot(133);
    imshow(im)                                                      
    title('重建图像')
 
 
    f0 = round(255*im);
    f1 = round(255*im2);
 
    errors(jj) = abs(sum(sum(f1 - f0)))/(8*256*256);
    clear im        
end