本示例描述了从一组沿不同角度获取的平行投影来重建图像的方法。本数据集是通过**“计算机断层扫描”**(computed tomography,即CT)获取的。 在样本上没有任何先验信息的情况下,重建图像所需的投影数量约为图像的线性大小l
(以像素为单位)。为简单起见,我们在这里考虑一个稀疏图像,其中只有对象边界上的像素具有非零值,例如,这样的数据可以对应于蜂窝材料,但是大多数图像在不同的基础上都是稀疏的,例如Haar小波, 因为仅获取
l/7
投影,因此必须使用样本上可用的先验信息(其稀疏性):这是压缩感知的一个示例。
断层投影算子(tomography projection operation)是一种线性变换。除了对应于线性回归的数据保真度项外,我们还惩罚图像的L1范数以解释其稀疏性,由此产生的优化问题称为
Lasso
。我们使用
sklearn.linear_model.Lasso
类,该类使用坐标下降算法,与这里使用的投影算子(projection operator)相比,这种实现在稀疏矩阵上的计算效率更高。
即使在投影中加入了噪声,用L1惩罚法重建得到的结果误差为零(所有像素都成功地标记为0或1)。相比之下,L2惩罚(
sklearn.linear_model.Ridge
)会产生大量像素标记错误,在重建图像上观察到重要的伪影(artifacts)。特别要注意的是,圆角伪影将角落中的像素分开,与中央磁盘(central disk)相比,圆角伪影贡献的投影更少。
print(__doc__)# 作者: Emmanuelle Gouillart # 许可证: BSD 3 clause
import numpy as np
from scipy import sparse
from scipy import ndimage
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge
import matplotlib.pyplot as plt
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) / dxreturn 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 - centerreturn X, Y
def build_projection_operator(l_x, n_dir):""" 计算断层扫描设计矩阵。
参数
----------
l_x : int
图像阵列的线性大小
n_dir : int
获取投影的角度数。
返回
-------
p : 形状为(n_dir l_x, l_x**2)的稀疏矩阵
"""
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 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.) ** 2 + (y - l / 2.) ** 2 mask = np.zeros((l, l))
points = l * rs.rand(2, n_pts)
mask[(points[0]).astype(np.int), (points[1]).astype(np.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(岭)惩罚重建
rgr_ridge = Ridge(alpha=0.2)
rgr_ridge.fit(proj_operator, proj.ravel())
rec_l2 = rgr_ridge.coef_.reshape(l, l)# L1(Lasso)惩罚重建# 通过LassoCV交叉验证确定alpha的最佳值
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('original image')
plt.subplot(132)
plt.imshow(rec_l2, cmap=plt.cm.gray, interpolation='nearest')
plt.title('L2 penalization')
plt.axis('off')
plt.subplot(133)
plt.imshow(rec_l1, cmap=plt.cm.gray, interpolation='nearest')
plt.title('L1 penalization')
plt.axis('off')
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
right=1)
plt.show()
脚本的总运行时间: ( 0 分 8.929 秒)
估计的内存使用量: 40 MB
下载python源代码:plot_tomography_l1_reconstruction.py
下载Jupyter notebook源代码:plot_tomography_l1_reconstruction.ipynb
由Sphinx-Gallery生成的画廊
☆☆☆为方便大家查阅,小编已将scikit-learn学习路线专栏文章统一整理到公众号底部菜单栏,同步更新中,关注公众号,点击左下方“系列文章”,如图:
欢迎大家和我一起沿着scikit-learn文档这条路线,一起巩固机器学习算法基础。(添加微信:mthler,备注:sklearn学习,一起进【sklearn机器学习进步群】开启打怪升级的学习之旅。)