编程记录——研究一下python对shepp_logan体模数据实现radon变换

参考博客:
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。

  • 1
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
好的,以下是基于MATLAB的Shepp-Logan实验代码。首先,我们需要生成Shepp-Logan: ```matlab % Shepp-Logan参数设置 N = 256; % 图像大小 a = 1; % 长轴大小 b = 0.6; % 短轴大小 c = 0.2; % 密度 phi = [0, 0, 0, 0, 0, 0, 0]; % 旋转角度 theta = [0, pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi]; % 旋转方向 % 生成Shepp-Logan [X,Y] = meshgrid(-1:2/(N-1):1); m = zeros(N,N); for i = 1:length(theta) t = ((X*cos(theta(i))+Y*sin(theta(i))).^2)/a^2 + ((Y*cos(theta(i))-X*sin(theta(i))).^2)/b^2; idx = t <= 1; m(idx) = m(idx) + c*((-1)^i)*((phi(i)*pi/180)*ones(size(t(idx)))); end ``` 接下来,我们将生成的Shepp-Logan进行投影,拟CT扫描: ```matlab % 投影角度设置 theta = 0:179; % 投影得到sinogram sinogram = radon(m, theta); % 显示sinogram figure; imagesc(theta, linspace(-1, 1, size(sinogram, 1)), sinogram); colormap(gray); xlabel('Projection angle (degrees)'); ylabel('Projection position (normalized)'); title('Sinogram of Shepp-Logan phantom'); ``` 最后,我们可以将投影结果进行反投影,重建原始图像: ```matlab % 反投影重建图像 reconstructed = iradon(sinogram, theta, 'linear', 'Ram-Lak', 1, N); % 显示重建图像 figure; imagesc(reconstructed); colormap(gray); axis square; title('Reconstructed image of Shepp-Logan phantom'); ``` 完整的代码如下所示: ```matlab % Shepp-Logan参数设置 N = 256; % 图像大小 a = 1; % 长轴大小 b = 0.6; % 短轴大小 c = 0.2; % 密度 phi = [0, 0, 0, 0, 0, 0, 0]; % 旋转角度 theta = [0, pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi]; % 旋转方向 % 生成Shepp-Logan [X,Y] = meshgrid(-1:2/(N-1):1); m = zeros(N,N); for i = 1:length(theta) t = ((X*cos(theta(i))+Y*sin(theta(i))).^2)/a^2 + ((Y*cos(theta(i))-X*sin(theta(i))).^2)/b^2; idx = t <= 1; m(idx) = m(idx) + c*((-1)^i)*((phi(i)*pi/180)*ones(size(t(idx)))); end % 投影角度设置 theta = 0:179; % 投影得到sinogram sinogram = radon(m, theta); % 显示sinogram figure; imagesc(theta, linspace(-1, 1, size(sinogram, 1)), sinogram); colormap(gray); xlabel('Projection angle (degrees)'); ylabel('Projection position (normalized)'); title('Sinogram of Shepp-Logan phantom'); % 反投影重建图像 reconstructed = iradon(sinogram, theta, 'linear', 'Ram-Lak', 1, N); % 显示重建图像 figure; imagesc(reconstructed); colormap(gray); axis square; title('Reconstructed image of Shepp-Logan phantom'); ``` 希望对您有所帮助!

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值