编程记录——研究一下python对shepp

参考博客：
CT典型数据——shepp_logan体模数据的生成 python版本
Python实现离散Radon变换
声明：除了对上述第二个链接的代码有少量修改、以及展示的运行结果外，其余工作都是对上述两个链接内容的整理总结。

生成一个shepp_logan体模数据

import phantominator
from phantominator import shepp_logan
import matplotlib.pyplot as plt

ph = shepp_logan(512)
plt.imshow(ph, cmap='gray')
plt.show()
print(ph.shape)

运行结果：

尝试对该体模数据进行radon变换：

from scipy import ndimage
import numpy as np
import matplotlib.pyplot as plt
import imageio
from cv2 import cv2

def DiscreteRadonTransform(image, views, detectors):
    res = np.zeros((detectors, views), dtype='float64')
    for s in range(views):
        # 将image逆时针旋转-s*180/views度的结果
        rotation = ndimage.rotate(image, -s*180/views, reshape=False)
        # sum()得到一个vector，其中元素是rotation中每列元素之和
        # 这个vector代表一个view的projection，然后存入res的列
        res[:,s] = sum(rotation)
    return res

# 256个view，512个detectors
# 编程时发现，好像detector的数量就必须和image的size保持一致
# 另一种写法可能更合理：
# radon = DiscreteRadonTransform(ph, 256, len(ph[0]))
# 这样让detectors数量能和图像的size保持一致，避免犯错
radon = DiscreteRadonTransform(ph, 256, 512)
radon1 = DiscreteRadonTransform(ph, 512, 512)
print(radon.shape)
print(radon1.shape)

#绘制原始图像和对应的sinogram图
plt.subplot(2, 2, 1)
plt.imshow(ph, cmap='gray')
plt.subplot(2, 2, 2)
plt.imshow(radon, cmap='gray')
plt.subplot(2, 2, 3)
plt.imshow(radon1, cmap='gray')
plt.show()

运行结果：

获得的sinogram，每一列是一个view的projection。