SourceURL:file:///home/pai4090/断层成像/2/医学断层成像实验crt.docx
- 平行投影解析算法获得Shepp-Logan模型的投影数据,并显示对应的sinogram图。
Shepp-Logan头模型参数
编号 | 椭圆对应密度 | 椭圆长轴a | 椭圆长轴b | 中心点x | 中心点y | 椭圆倾角 |
1 | 1 | 0.69 | 0.92 | 0 | 0 | 0 |
2 | -0.8 | 0.6624 | 0.874 | 0 | -0.0184 | 0 |
3 | -0.2 | 0.11 | 0.31 | 0.22 | 0 | -18 |
4 | -0.2 | 0.16 | 0.41 | -0.22 | 0 | 18 |
5 | 0.1 | 0.21 | 0.25 | 0 | 0.35 | 0 |
6 | 0.1 | 0.46 | 0.46 | 0 | -0.1 | 0 |
7 | 0.1 | 0.46 | 0.46 | 0 | -0.1 | 0 |
8 | 0.1 | 0.46 | 0.23 | -0.08 | -0.605 | 0 |
9 | 0.1 | 0.023 | 0.023 | 0 | -0.606 | 0 |
10 | 0.1 | 0.023 | 0.046 | 0.06 | -0.605 | 0 |
python代码:
import numpy as np
import matplotlib.pyplot as plt
import phantominator
from phantominator import shepp_logan
def forward_projection(theta_proj, N, N_d):
shep = np.array([[1, 0.69, 0.92, 0, 0, 0],
[-0.8, 0.6624, 0.8740, 0, -0.0184, 0],
[-0.2, 0.1100, 0.3100, 0.22, 0, -18],
[-0.2, 0.1600, 0.4100, -0.22, 0, 18],
[0.1, 0.2100, 0.2500, 0, 0.35, 0],
[0.1, 0.0460, 0.0460, 0, 0.1, 0],
[0.1, 0.0460, 0.0460, 0, 0.1, 0],
[0.1, 0.0460, 0.0230, -0.08, -0.605, 0],
[0.1, 0.0230, 0.0230, 0, -0.606, 0],
[0.1, 0.0230, 0.0460, 0.06, -0.605, 0]])
theta_num = len(theta_proj)
P = np.zeros((int(N_d), theta_num))
rho = shep[:, 0]
ae = 0.5 * N * shep[:, 1]
be = 0.5 * N * shep[:, 2]
xe = 0.5 * N * shep[:, 3]
ye = 0.5 * N * shep[:, 4]
alpha = shep[:, 5]
alpha = alpha * np.pi / 180
theta_proj = theta_proj * np.pi / 180
TT = np.arange(-(N_d - 1) / 2, (N_d - 1) / 2 + 1)
for k1 in range(theta_num):
P_theta = np.zeros(int(N_d))
for k2 in range(len(xe)):
a = (ae[k2] * np.cos(theta_proj[k1] - alpha[k2]))**2 + (be[k2] * np.sin(theta_proj[k1] - alpha[k2]))**2
temp = a - (TT - xe[k2] * np.cos(theta_proj[k1]) - ye[k2] * np.sin(theta_proj[k1]))**2
ind = temp > 0
P_theta[ind] += rho[k2] * (2 * ae[k2] * be[k2] * np.sqrt(temp[ind])) / a
P[:, k1] = P_theta
P_min = np.min(P)
P_max = np.max(P)
P = (P - P_min) / (P_max - P_min)
return P
N = 256
theta = np.arange(0, 180, 1)
I = shepp_logan(N) # Replace with your phantom generation code
N_d = 2 * np.ceil(np.linalg.norm(np.array(I.shape) - np.floor((np.array(I.shape) - 1) / 2) - 1)) + 3
P = forward_projection(theta, N, int(N_d))
plt.figure()
plt.imshow(I, cmap='gray')
plt.title('256x256 Head Phantom')
plt.figure()
plt.imshow(P, cmap='gray')
plt.colorbar()
plt.title('180° Parallel Beam Projection')
theta1 = np.arange(0, 45, 1)
N_d = 2 * np.ceil(np.linalg.norm(np.array(I.shape) - np.floor((np.array(I.shape) - 1) / 2) - 1)) + 3
P1 = forward_projection(theta1, N, int(N_d))
plt.figure()
plt.imshow(P1, cmap='viridis')
plt.colorbar()
plt.title('45° Parallel Beam Projection')
plt.show()
theta2 = np.arange(0, 300, 1)
N_d = 2 * np.ceil(np.linalg.norm(np.array(I.shape) - np.floor((np.array(I.shape) - 1) / 2) - 1)) + 3
P2 = forward_projection(theta2, N, int(N_d))
plt.figure()
plt.imshow(P2, cmap='viridis')
plt.colorbar()
plt.title('300° Parallel Beam Projection')
plt.show()
- shepp-logan头模型为点源。
点源参数
编号 | 椭圆对应密度 | 椭圆长轴a | 椭圆长轴b | 中心点x | 中心点y | 椭圆倾角 |
1 | 1 | 0.005 | 0.005 | 0.5 | 0 | 0 |
代码:
import numpy as np
import matplotlib.pyplot as plt
import phantominator
from phantominator import shepp_logan
def forward_projection(theta, N, N_d):
point = np.array([1, 0.005, 0.005, 0.5, 0, 0])
theta_num = len(theta)
P = np.zeros((N_d, theta_num))
rho = point[0]
ae = 0.5 * N * point[1]
be = 0.5 * N * point[2]
xe = 0.5 * N * point[3]
ye = 0.5 * N * point[4]
alpha = point[5]
alpha = alpha * np.pi / 180
theta = np.array(theta) * np.pi / 180
TT = np.arange(-(N_d - 1) / 2, (N_d - 1) / 2 + 1)
for k1 in range(theta_num):
P_theta = np.zeros(N_d)
for k2 in range(1): # 修改此处
a = (ae * np.cos(theta[k1] - alpha)) ** 2 + (be * np.sin(theta[k1] - alpha)) ** 2
temp = a - (TT - xe * np.cos(theta[k1]) - ye * np.sin(theta[k1])) ** 2
ind = temp > 0
P_theta[ind] += rho * (2 * ae * be * np.sqrt(temp[ind])) / a
P[:, k1] = P_theta
return P
N = 256
theta = np.arange(0, 360, 1)
I = shepp_logan(N) # 实现phantom函数生成头部图像的代码
N_d = int(2 * np.ceil(np.linalg.norm(np.array(I.shape) - np.floor((np.array(I.shape) - 1) / 2) - 1)) + 3)
P = forward_projection(theta, N, N_d)
plt.figure()
plt.imshow(P, cmap='gray')
plt.colorbar()
plt.title('point 360° Parallel Beam Projection')
plt.show()
theta1 = np.arange(0, 45, 1)
N_d = 2 * np.ceil(np.linalg.norm(np.array(I.shape) - np.floor((np.array(I.shape) - 1) / 2) - 1)) + 3
P1 = forward_projection(theta1, N, int(N_d))
plt.figure()
plt.imshow(P1, cmap='viridis')
plt.colorbar()
plt.title('point 45° Parallel Beam Projection')
plt.show()
theta2 = np.arange(0, 300, 1)
N_d = 2 * np.ceil(np.linalg.norm(np.array(I.shape) - np.floor((np.array(I.shape) - 1) / 2) - 1)) + 3
P2 = forward_projection(theta2, N, int(N_d))
plt.figure()
plt.imshow(P2, cmap='hot')
plt.colorbar()
plt.title('point 300° Parallel Beam Projection')
plt.show()