参考博客:
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变换
尝试对该体模数据进行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。