文章来源于网络收集而来,版权归原创者所有,如有侵权请及时联系!
生成投影
代码
def _weights(x, dx=1, orig=0):
x = np.ravel(x)
floor_x = np.floor((x - orig) / dx)
alpha = (x - orig - floor_x * dx) / dx
return np.hstack((floor_x, floor_x + 1)), np.hstack((1 - alpha, alpha))
def _generate_center_coordinates(l_x):
X, Y = np.mgrid[:l_x, :l_x].astype(np.float64)
center = l_x / 2.
X += 0.5 - center
Y += 0.5 - center
return X, Y
def build_projection_operator(l_x, n_dir):
X, Y = _generate_center_coordinates(l_x)
angles = np.linspace(0, np.pi, n_dir, endpoint=False)
data_inds, weights, camera_inds = [], [], []
data_unravel_indices = np.arange(l_x ** 2)
data_unravel_indices = np.hstack((data_unravel_indices,
data_unravel_indices))
for i, angle in enumerate(angles):
Xrot = np.cos(angle) * X - np.sin(angle) * Y
inds, w = _weights(Xrot, dx=1, orig=X.min())
mask = (inds >= 0) & (inds < l_x)
weights += list(w[mask])
camera_inds += list(inds[mask] + i * l_x)
data_inds += list(data_unravel_indices[mask])
proj_operator = sparse.coo_matrix((weights, (camera_inds, data_inds)))
return proj_operator
投影运算符
l = 128
proj_operator = build_projection_operator(l, l // 7)
proj_operator
'''
<2304x16384 sparse matrix of type '<class 'numpy.float64'>'
with 555378 stored elements in COOrdinate format>
'''
维度:角度( l // 7
),位置( l
),每个图像( l x l
)
proj_t = np.reshape(proj_operator.todense().A, (l//7,l,l,l))
第一个坐标指的是线的角度,第二个坐标指代线的位置。
索引为 3 的角度的直线:
plt.imshow(proj_t[3,0], cmap='gray');
plt.imshow(proj_t[3,1], cmap='gray');
plt.imshow(proj_t[3,2], cmap='gray');
plt.imshow(proj_t[3,40], cmap='gray');
垂直位置 40 处的其他直线:
plt.imshow(proj_t[4,40], cmap='gray');
plt.imshow(proj_t[15,40], cmap='gray');
plt.imshow(proj_t[17,40], cmap='gray');
X 射线和数据之间的交点
接下来,我们想看看直线如何与我们的数据相交。 请记住,这就是数据的样子:
plt.figure(figsize=(5,5))
plt.imshow(data, cmap=plt.cm.gray)
plt.axis('off')
plt.savefig("images/data.png")
proj = proj_operator @ data.ravel()[:, np.newaxis]
角度为 17,位置为 40 的穿过数据的 X 射线:
plt.figure(figsize=(5,5))
plt.imshow(data + proj_t[17,40], cmap=plt.cm.gray)
plt.axis('off')
plt.savefig("images/data_xray.png")
它们相交的地方。
both = data + proj_t[17,40]
plt.imshow((both > 1.1).astype(int), cmap=plt.cm.gray);
那条 X 射线的强度:
np.resize(proj, (l//7,l))[17,40]
# 6.4384498372605989
角度为 3,位置为 14 的穿过数据的 X 射线:
plt.imshow(data + proj_t[3,14], cmap=plt.cm.gray);
它们相交的地方。
both = data + proj_t[3,14]
plt.imshow((both > 1.1).astype(int), cmap=plt.cm.gray);
CT 扫描的测量结果在这里是一个小数字:
np.resize(proj, (l//7,l))[3,14]
# 2.1374953737965541
proj += 0.15 * np.random.randn(*proj.shape)
关于 *args
a = [1,2,3]
b = [4,5,6]
c = list(zip(a, b))
c
# [(1, 4), (2, 5), (3, 6)]
list(zip(*c))
# [(1, 2, 3), (4, 5, 6)]
投影(CT 读取)
plt.figure(figsize=(7,7))
plt.imshow(np.resize(proj, (l//7,l)), cmap='gray')
plt.axis('off')
plt.savefig("images/proj.png")
如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
发布评论