ASD-POCSによるスパースビュー再構成

導入
---

後輩氏「POCSって何ですか？」
後輩氏「圧縮センシングがやりたいです」
ワイ氏「ちょっと調べてみるわ」

ASD-POCSを書いてみた。

インポート
---
```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
```
シミュレーションデータ
---
```python

from skimage.io import imread
from skimage import data_dir
from skimage.transform import downscale_local_mean, radon, iradon, iradon_sart
from skimage.measure import compare_psnr

N_view = 15
N_sart = 20

theta = np.linspace(0, 180, N_view, endpoint=False)
im = imread(data_dir+"/phantom.png", as_grey=True)
im = downscale_local_mean(im, (2, 2))

sng = radon(im, theta=theta, circle=True)
fbp = iradon(sng, theta=theta, circle=True)

sart = np.zeros_like(im)
for i in range(N_sart): 
    sart= iradon_sart(sng, image=sart, theta=theta)    

fig, ax = plt.subplots(1, 3, figsize=(12, 4))
ax[0].imshow(im, cmap='gray', vmin=0, vmax=1)
ax[1].imshow(fbp, cmap='gray', vmin=0, vmax=1)
ax[2].imshow(sart, cmap='gray', vmin=0, vmax=1)

ax[0].set_title('真\nPSNR [dB]')
ax[1].set_title('FBP\n{:.2f}'.format(compare_psnr(im, fbp)))
ax[2].set_title('SART\n{:.2f}'.format(compare_psnr(im, sart)))

for a in ax:
    a.axis('off')

plt.savefig("result1.png", dpi=220)
```
![result1.png](https://qiita-image-store.s3.amazonaws.com/0/114365/7a61795b-4508-074e-b04d-c2cecff163b9.png)

TVの勾配
---
```python
def gradTV(im, eps=1e-8):
    """ grad of TV """
    x_minus_1 = np.roll(im, 1, axis=1)
    x_minus_1_y_plus_1 = np.roll(x_minus_1, -1, axis=0)
    y_minus_1 = np.roll(im, 1, axis=0)
    x_plus_1_y_minus_1 = np.roll(y_minus_1, -1, axis=1)
    x_plus_1 = np.roll(im, -1, axis=1)
    y_plus_1 = np.roll(im, -1, axis=0)
    
    v = (im - x_minus_1) / np.sqrt((im - x_minus_1) ** 2 + (x_minus_1_y_plus_1 - x_minus_1) ** 2 + eps)
    v += (im - y_minus_1) / np.sqrt((x_plus_1_y_minus_1 - y_minus_1) ** 2 + (im - y_minus_1) ** 2 + eps)
    v -= (x_plus_1 + y_plus_1 - im * 2) / np.sqrt((x_plus_1 - im) ** 2 + (y_plus_1 - im) ** 2 + eps)
    
    return v
```

ASD-POCS
---
```python
from skimage.measure import compare_psnr

N_iter = 1000
N_sart = 1

ng = 20
alpha = 0.2
r_max = 0.95
alpha_red = 0.95
eps = 0.

x = np.zeros_like(im)

for i in range(N_iter):
    x0 = x.copy()

    # SART
    for _ in range(N_sart): 
        x = iradon_sart(sng, image=x, theta=theta)
    x *= x > 0
    x_res = x.copy()
    
    # Adaptive steppest descent of TV minimization
    b = radon(x, theta=theta, circle=True)
    dd = np.sqrt(np.sum((b - sng) ** 2))
    dp = np.sqrt(np.sum((x - x0) ** 2)) 

    if i == 0:
        dtvg = alpha * dp

    x0 = x.copy()
    
    for _ in range(ng):
        dx = gradTV(x)
        dx = dx / (np.sqrt(np.sum(dx ** 2)))
#         dx = dx / (np.sqrt(np.sum(dx ** 2)) + 1e-8)
        x = x - dtvg * dx
        
    dg = np.sqrt(np.sum((x - x0) ** 2))
    
    if (dg > r_max * dp) and (dd > eps):
        dtvg = dtvg * alpha_red
    
    if (i + 1) % 10 == 0: 
        print(i + 1, compare_psnr(im, x_res), dtvg)
           
asd_pocs = x_res
```
結果
---
![result.png](https://qiita-image-store.s3.amazonaws.com/0/114365/1218d3c6-7379-da5d-38aa-6523e866e23b.png)

考察
---
* 元のASD-POCSはARTだが、書くのがめんどくさかったので、skimageのiradon_sartを使った
* ファントム下部の小さい3つのスポットは分解できていない

リンク
--
[ASD-POCSのJupyter notebook](https://github.com/kibo35/reconstruction/blob/master/ASD_POCS/ASD_POCS.ipynb)
[Radon transform -- skimage](http://scikit-image.org/docs/dev/auto_examples/transform/plot_radon_transform.html#sphx-glr-auto-examples-transform-plot-radon-transform-py)
[Sidky EY, Pan X, Phys Med Biol 53: 4777-4807, 2008](https://www.ncbi.nlm.nih.gov/pubmed/18701771)