返回介绍

生成投影

发布于 2025-01-01 12:38:40 字数 4844 浏览 0 评论 0 收藏 0

代码

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 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。
列表为空,暂无数据
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
    原文