導入
後輩氏「POCSって何ですか?」
後輩氏「圧縮センシングがやりたいです」
ワイ氏「ちょっと調べてみるわ」
ASD-POCSを書いてみた。
インポート
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
シミュレーションデータ
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)
TVの勾配
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
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
結果
考察
- 元のASD-POCSはARTだが、書くのがめんどくさかったので、skimageのiradon_sartを使った
- ファントム下部の小さい3つのスポットは分解できていない
リンク
ASD-POCSのJupyter notebook
Radon transform -- skimage
Sidky EY, Pan X, Phys Med Biol 53: 4777-4807, 2008