Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

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

導入

後輩氏「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)

result1.png

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

結果

result.png

考察

  • 元の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

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away