Question

引言

我们将使用生成的数据（不是真正的 CT 扫描）。 生成数据涉及一些有趣的 numpy 和线性代数，我们稍后会再回过头来看。

代码来自 Scikit-Learn 示例 压缩感知：使用 L1 先验的层析成像重建（Lasso） 。

生成图像

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 < (l / 2) ** 2
    mx,my = rs.randint(0, l, (2,n_pts))
    mask = np.zeros((l, l))
    mask[mx,my] = 1
    mask = ndimage.gaussian_filter(mask, sigma=l / n_pts)
    res = (mask > mask.mean()) & mask_outer
    return res ^ ndimage.binary_erosion(res)

l = 128
data = generate_synthetic_data()

plt.figure(figsize=(5,5))
plt.imshow(data, cmap=plt.cm.gray);

generate_synthetic_data 在做什么

l=8; n_pts=5
rs = np.random.RandomState(0)

x, y = np.ogrid[0:l, 0:l]; x,y

'''
(array([[0],
        [1],
        [2],
        [3],
        [4],
        [5],
        [6],
        [7]]), array([[0, 1, 2, 3, 4, 5, 6, 7]]))
'''

x + y

'''
array([[ 0,  1,  2,  3,  4,  5,  6,  7],
       [ 1,  2,  3,  4,  5,  6,  7,  8],
       [ 2,  3,  4,  5,  6,  7,  8,  9],
       [ 3,  4,  5,  6,  7,  8,  9, 10],
       [ 4,  5,  6,  7,  8,  9, 10, 11],
       [ 5,  6,  7,  8,  9, 10, 11, 12],
       [ 6,  7,  8,  9, 10, 11, 12, 13],
       [ 7,  8,  9, 10, 11, 12, 13, 14]])
'''

(x - l/2) ** 2 

'''
array([[ 16.],
       [  9.],
       [  4.],
       [  1.],
       [  0.],
       [  1.],
       [  4.],
       [  9.]])
'''

(x - l/2) ** 2 + (y - l/2) ** 2

'''
array([[ 32.,  25.,  20.,  17.,  16.,  17.,  20.,  25.],
       [ 25.,  18.,  13.,  10.,   9.,  10.,  13.,  18.],
       [ 20.,  13.,   8.,   5.,   4.,   5.,   8.,  13.],
       [ 17.,  10.,   5.,   2.,   1.,   2.,   5.,  10.],
       [ 16.,   9.,   4.,   1.,   0.,   1.,   4.,   9.],
       [ 17.,  10.,   5.,   2.,   1.,   2.,   5.,  10.],
       [ 20.,  13.,   8.,   5.,   4.,   5.,   8.,  13.],
       [ 25.,  18.,  13.,  10.,   9.,  10.,  13.,  18.]])
'''

mask_outer = (x - l/2) ** 2 + (y - l/2) ** 2 < (l/2) ** 2; mask_outer

'''
array([[False, False, False, False, False, False, False, False],
       [False, False,  True,  True,  True,  True,  True, False],
       [False,  True,  True,  True,  True,  True,  True,  True],
       [False,  True,  True,  True,  True,  True,  True,  True],
       [False,  True,  True,  True,  True,  True,  True,  True],
       [False,  True,  True,  True,  True,  True,  True,  True],
       [False,  True,  True,  True,  True,  True,  True,  True],
       [False, False,  True,  True,  True,  True,  True, False]], dtype=bool)
'''

plt.imshow(mask_outer, cmap='gray')

# <matplotlib.image.AxesImage at 0x7efcd9303278>

mask = np.zeros((l, l))
mx,my = rs.randint(0, l, (2,n_pts))
mask[mx,my] = 1; mask

'''
array([[ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.]])
'''

plt.imshow(mask, cmap='gray')

# <matplotlib.image.AxesImage at 0x7efcd9293940>

mask = ndimage.gaussian_filter(mask, sigma=l / n_pts)

plt.imshow(mask, cmap='gray')

# <matplotlib.image.AxesImage at 0x7efcd922c0b8>

res = np.logical_and(mask > mask.mean(), mask_outer)
plt.imshow(res, cmap='gray');

plt.imshow(ndimage.binary_erosion(res), cmap='gray');

plt.imshow(res ^ ndimage.binary_erosion(res), cmap='gray');