图像重建技术与压缩感知

在医学成像领域,计算机断层扫描(CT)是一种常用的技术,它能够通过不同角度获取的一系列平行投影来重建图像。这种数据采集方法需要大量的投影数据,其数量通常与图像的线性尺寸(像素数)成正比。为了简化问题,假设图像是稀疏的,即只有物体边界上的像素具有非零值。这种数据可能对应于细胞材料等。然而,大多数图像在不同的基底上是稀疏的,例如Haar小波。因此,仅采集了图像线性尺寸的1/7的投影数据,这就需要利用样本的先验信息(其稀疏性),这是压缩感知的一个例子。

断层摄影投影操作是一种线性变换。除了对应于线性回归的数据保真项外,还对图像的L1范数进行惩罚,以考虑其稀疏性。由此产生的优化问题被称为Lasso。使用了Lasso类,它使用了坐标下降算法。重要的是,这种实现在稀疏矩阵上的计算效率比这里使用的投影操作更高。

L1惩罚的重建结果实现了零误差(所有像素都成功地被标记为0或1),即使在投影中添加了噪声。相比之下,L2惩罚(Ridge)为像素产生了大量的标记错误。在重建的图像上观察到重要的伪影,与L1惩罚相反。特别要注意的是,将角落的像素与中心盘分开的圆形伪影,这些像素对较少的投影做出了贡献。

代码实现

以下是一个使用Python语言和相关科学计算库实现的图像重建示例。该示例展示了如何生成合成的二进制数据,并使用Lasso优化算法进行图像重建。

import matplotlib.pyplot as plt import numpy as np from scipy import ndimage, sparse from sklearn.linear_model import Lasso, Ridge def _weights(x, dx=1, orig=0): x = np.ravel(x) floor_x = np.floor((x - orig) / dx).astype(np.int64) 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.0 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 = np.logical_and(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 def generate_synthetic_data(): ""“生成合成的二进制数据。”"" rs = np.random.RandomState(0) n_pts = 36 x, y = np.ogrid[0:l, 0:l] mask_outer = (x - l / 2.0)**2 + (y - l / 2.0)**2 < (l / 2.0)**2 mask = np.zeros((l, l)) points = l * rs.rand(2, n_pts) mask[(points[0]).astype(int), (points[1]).astype(int)] = 1 mask = ndimage.gaussian_filter(mask, sigma=l / n_pts) res = np.logical_and(mask > mask.mean(), mask_outer) return np.logical_xor(res, ndimage.binary_erosion(res)) # 生成合成图像和投影 l = 128 proj_operator = build_projection_operator(l, l // 7) data = generate_synthetic_data() proj = proj_operator @ data.ravel()[:, np.newaxis] proj += 0.15 * np.random.randn(*proj.shape) # 使用L2(Ridge)惩罚进行重建 rgr_ridge = Ridge(alpha=0.2) rgr_ridge.fit(proj_operator, proj.ravel()) rec_l2 = rgr_ridge.coef_.reshape(l, l) # 使用L1(Lasso)惩罚进行重建 rgr_lasso = Lasso(alpha=0.001) rgr_lasso.fit(proj_operator, proj.ravel()) rec_l1 = rgr_lasso.coef_.reshape(l, l) plt.figure(figsize=(8, 3.3)) plt.subplot(131) plt.imshow(data, cmap=plt.cm.gray, interpolation="nearest") plt.axis("off") plt.title("原始图像") plt.subplot(132) plt.imshow(rec_l2, cmap=plt.cm.gray, interpolation="nearest") plt.title("L2惩罚") plt.axis("off") plt.subplot(133) plt.imshow(rec_l1, cmap=plt.cm.gray, interpolation="nearest") plt.title("L1惩罚") plt.axis("off") plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1) plt.show()

上述代码首先定义了一组辅助函数,用于生成中心坐标、构建投影操作符以及生成合成的二进制数据。然后,代码展示了如何使用Lasso和Ridge算法对添加了噪声的投影数据进行重建,并比较了两种方法的效果。

沪ICP备2024098111号-1
上海秋旦网络科技中心:上海市奉贤区金大公路8218号1幢 联系电话:17898875485