PET图像重建OSEM算法

 一、代码

代码是一个用于图像重建的函数，特别是用于正电子发射断层扫描（PET）或X射线计算机断层扫描（CT）的有序子集期望最大化（OSEM）算法。

from skimage.io import imread  # 导入图像读取函数
from skimage import data_dir  # 导入数据目录，用于获取示例数据集
from skimage.transform import radon, rescale, iradon  # 导入Radon变换、重缩放和反Radon变换函数
from skimage.draw import disk  # 导入绘制圆盘的函数
# from skimage.draw import  circle  # 导入绘制圆的函数（当前未使用）
from scipy.ndimage import gaussian_filter  # 导入高斯滤波函数
import numpy as np  # 导入NumPy库，用于数学运算
import matplotlib.pyplot as plt  # 导入Matplotlib的pyplot模块，用于绘图

# 定义OSEM算法函数
def osem(sinogram, nsub=1, niter=1, sfwhm=1.333, obj=True):
    # 获取正弦图的形状
    shape = sinogram.shape
    # 断言正弦图必须是正方形
    assert shape[0] == shape[1], "osem sino must have same size"
    # 创建角度数组，用于Radon变换,shape[0]等于400，即图片长度180度分为400份
    theta = np.linspace(0., 180., shape[0], endpoint=False)


    # 初始化重建图像为零矩阵
    recon = np.zeros(shape)
    # 创建一个圆盘，用于初始化重建图像的中心
    rr, cc = disk((shape[0] / 2, shape[1] / 2), shape[0] / 2 - 1)
    recon[rr, cc] = 1


    # 创建归一化矩阵，此例中nview等于400，norm为全一矩阵
    nview = len(theta)
    norm = np.ones(shape)
    wgts = []
    # 对每个子集进行迭代
    for sub in range(nsub):
        # 计算子集视图的索引
        views = range(sub, nview, nsub)
        # 计算权重
        wgt = iradon(norm[:, views], theta=theta[views], filter_name=None, circle=True)
        wgts.append(wgt)

    # 迭代过程 对应ELEM算法的数学公式进行迭代
    objs = []
    for iter in range(niter):
        # 随机排列子集的顺序
        order = np.random.permutation(range(nsub))
        for sub in order:
            # 计算子集视图的索引
            views = range(sub, nview, nsub)
            # 计算当前重建图像的Radon变换
            fp = radon(recon, theta=theta[views], circle=True)
            # 计算正弦图与Radon变换的比率，1e-6:输出科学计数法表示的小数1乘以10的负6次方。
            sinogram_test = sinogram[:, views]
            ratio = sinogram[:, views] / (fp + 1e-6)
            # 计算反投影
            bp = iradon(ratio, theta=theta[views], filter_name=None, circle=True)
            # 更新重建图像
            recon *= bp / (wgts[sub] + 1e-6)

            # 如果obj为True，计算目标函数，评估重建图像的质量
            if obj:
                fp = radon(recon, theta=theta, circle=True)
                ndx = np.where(fp > 0)
                val = -(sinogram[ndx] * np.log(fp[ndx]) - fp[ndx]).sum() / 1e6
                objs.append(val)

    # 使用反Radon变换进行滤波反投影重建
    fbp = iradon(sinogram, theta=theta, circle=True)
    # 如果sfwhm大于0，应用高斯滤波
    if sfwhm > 0:
        fbp = gaussian_filter(fbp, sfwhm / 2.355)
        recon = gaussian_filter(recon, sfwhm / 2.355)
    # 以下代码块被注释掉了，用于绘制收敛图和重建图像
    # plt.figure()
    # plt.suptitle('Convergence')
    # plt.plot(objs, '.-')
    # plt.xlabel('sub iteration')
    # plt.ylabel('negative log likelihood [1e+06]')
    # plt.savefig('objective.png')
    #
    # plt.figure()
    # plt.imshow(recon)
    # plt.axis('off')

    # 返回重建图像、滤波反投影重建图像和Radon变换图像
    return recon, fbp, fp
请注意，这段代码中有一些部分被注释掉了，这些部分通常用于调试和可视化结果。在实际使用时，可以根据需要取消注释。

二、代码与osem数学表达式的对应关系 

三、生成的结果

recon[rr, cc] = 1rr是对应圆盘的横坐标，cc是对应圆盘的纵坐标，recon[rr, cc] = 1是把圆盘全部元素置一

创建的圆盘