壓縮感知
什麽是壓縮感知呢？
示例說明：使用L1先驗(Lasso)進行層析成像重建
此示例介紹了從一組沿不同角度獲取的平行投影來重建圖像的方法。這種數據集一般是通過CT檢查獲得的。CT即計算機斷層掃描或計算層析成像。
在樣本上沒有任何先驗信息的情況下，重建圖像所需的投影數約為圖像線性大小l(width x height)的量級(以像素為單位)。為簡單起見，我們在這裏考慮一個稀疏圖像，其中隻有圖像中對象邊界上的像素具有非零值。不過，大多數圖像在不同的基礎上都是稀疏的，例如Haar小波。因為數據的稀疏性，隻能得到l/7投影，因此像做圖像重建有必要使用樣本上可用的先驗信息(稀疏度)：這就是壓縮感知。
層析成像投影操作是線性變換。除了對應於線性回歸的數據保真項之外，我們還對圖像的L1範數進行了懲罰，以考慮其稀疏性。由此產生的優化問題稱為Lasso，它對應的scikit-learn類是sklearn.linear_model.Lasso，使用坐標下降算法求解實現。重要的是，與這裏使用的投影運算符相比，該實現在稀疏矩陣上的計算效率更高。
即使將噪聲添加到投影中，使用L1懲罰進行的重建也會產生零誤差的結果(所有像素均成功標記為0或1)。相比之下，L2懲罰(sklearn.linear_model.Ridge)會為像素產生大量錯誤標記。也就是說，L2與L1懲罰相比，重建的圖像上會觀察到明顯的偽影。
代碼實現[Python]
# -*- coding: utf-8 -*-
print(__doc__)
# Author: Emmanuelle Gouillart
# License: 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) / 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):
""" Compute the tomography design matrix.
Parameters
----------
l_x : int
linear size of image array
n_dir : int
number of angles at which projections are acquired.
Returns
-------
p : sparse matrix of shape (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 < 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():
""" Synthetic binary 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.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)
# Reconstruction with L2 (Ridge) penalization
rgr_ridge = Ridge(alpha=0.2)
rgr_ridge.fit(proj_operator, proj.ravel())
rec_l2 = rgr_ridge.coef_.reshape(l, l)
# 使用L1懲罰(Lasso)重構
# 最佳alpha值，采用交叉驗證(LassoCV)調優確定
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分9.761秒。
運行代碼輸出的圖片內容如下:
源碼下載
- Python版源碼文件: plot_tomography_l1_reconstruction.py
- Jupyter Notebook版源碼文件: plot_tomography_l1_reconstruction.ipynb