TV正則化を用いたRichardson-Lucy deconvolutionを実装しよう

はじめに

Richardson-Lucy deconvolution（以下、RL法）[1, 2]は、画像等を鮮明にする手法である。一般的な撮像観測では、観測装置の回折等により光源が拡がって撮像される。この拡がり具合は、点拡がり関数（Point Spread Function, 以下、PSF）で表される。RL法は既知のPSFと撮像画像から、真の画像を推定する手法である。
RL法は、ノイズが多く入った画像ではノイズを増幅され、それを防ぐために正則化演算子を加えることがある。そのノイズの中でも、特にポアソンノイズの増幅抑制の効果のある正則化演算子の一つに、Total Variation (TV)正則化[3]がある。
そのTV正則化を用いたRL法だが、サンプルコードを見かけなかったので、本記事ではTV正則化を用いたRL法を実装する。そして、ポアソンノイズを加えたサンプル画像に対して、正則化項なしのRL法とTV正則化項を加えた場合の見た目の比較を行う。また、補足ではTV正則化の重みを適当に変化させた結果や定量評価を行なっているので興味のある方はそちらもご覧いただきたい。

RL法とは

RL法は、既知のPSFと撮像画像を用いて真の画像を推定する手法である。真の画像の分布、既知のPSF、撮像画像の関係をベイズの定理を用いて関係づける。そして、ベイズ推定を繰り返し計算することにより真の画像を推定する。理論式は

 W_{i,r+1} = W_{i,r}\sum_k \frac{P_{i,k}H_k}{\sum_j P_{j,k}W_{j,r}}\
\quad r = 0,\ 1,\  2,\  \dots

と表される。$j,k$はピクセル座標とする。$W$は真の画像、$H$は撮像画像、$r$は反復回数とする。$P_{j,k}$はPSFで、$W_j$が$H_k$で観測される確率を$P(H_k|W_j)=P_{j,k}$とする。式の導出は、私の記事で恐縮ですが、Richardson-Lucy deconvolutionの理論式を導出しようを参照されたい。

TV正則化を用いたRL法

Total Variation (TV) 正則化[3]を用いることで、画像中のエッジを保存し、均質な領域を平滑化することができる[4]。TV正則化を用いたRL法は、

 W_{i,r+1} = \frac{W_{i,r}}{1-\lambda_{TV}\textrm{div}\left( \frac{\nabla W_{i,r}}{|\nabla W_{i,r}|}\right)}\sum_k \frac{P_{i,k}H_k}{\sum_j P_{j,k}W_{j,r}}

と表される。$\lambda_{TV}$は、正則化の重みを表し、原論文[4]では、$\lambda_{TV}=0.002$を用いている。
ここで、$\lambda_{TV}$について補足すると、$\lambda_{TV}$が小さくなるほどTV正則化の影響が小さくなり、特に$\lambda_{TV}=0$のとき、RL法と同一の式となる。また、$\lambda_{TV}=1$のように大きい値を入れると、$1-\lambda_{TV}\textrm{div}\left( \frac{\nabla W_{i,r}}{|\nabla W_{i,r}|}\right)<0$となり、真の画像が負になってしまうので注意したい。

実装

Google Colabで作成した本記事のコードは、こちらにあります。

サンプル画像

本記事ので使用する真の画像は

よりダウンロードできる。

各種インポート

各種インポート

import cv2
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from scipy.signal import convolve

実行時の各種バージョン

Name	Version
opencv-python	4.6.0.66
numpy	1.21.6
matplotlib	3.2.2
scipy	1.7.3

PSF作成用の関数

PSF作成用の関数

def get_psf(px, sigma):
    half_px = px // 2
    x, y = np.mgrid[-half_px:half_px+1:, -half_px:half_px+1:]
    pos = np.dstack((x, y))
    mean = np.array([0, 0])  # 分布の平均
    cov = np.array([[sigma**2, 0], [0, sigma**2]])  # 分布の分散共分散行列
    rv = multivariate_normal(mean, cov)
    psf = rv.pdf(pos)
    return psf

使用するPSFは何でも良いが、本記事ではガウシアンフィルタで実装する。

TV正則化を用いたRL法

tv正則化を用いたRL法

def richardson_lucy_tv_reg(image, psf, iterations=50, clip=True, filter_epsilon=None, tv_lambda=0.002):
    float_type = np.promote_types(image.dtype, np.float32)
    image = image.astype(float_type, copy=False)
    psf = psf.astype(float_type, copy=False)
    im_deconv = np.full(image.shape, 0.5, dtype=float_type)
    psf_mirror = np.flip(psf)

    # Small regularization parameter used to avoid 0 divisions
    eps = 1e-12

    for _ in range(iterations):
        conv = convolve(im_deconv, psf, mode='same') + eps  # avoid 0 divisions
        if filter_epsilon:
            relative_blur = np.where(conv < filter_epsilon, 0, image / conv)
        else:
            relative_blur = image / conv
        
        # TV regularization
        grad_x = np.gradient(im_deconv, axis=1)
        grad_y = np.gradient(im_deconv, axis=0)
        grad_norm = np.sqrt(grad_x**2+grad_y**2) + eps  # avoid 0 divisions
        tv_reg = np.gradient(grad_x / grad_norm, axis=1) + np.gradient(grad_y / grad_norm, axis=0)

        im_deconv = im_deconv/(1-tv_lambda*tv_reg)*convolve(relative_blur, psf_mirror, mode='same')

    if clip:
        im_deconv[im_deconv > 1] = 1
        im_deconv[im_deconv < -1] = -1

    return im_deconv

コードの説明

RL法の実装部は、scikit-image version0.18.3, 0.19.2を参考にしているだけなので、ここではTV正則化の該当箇所のみの説明に留める。

TV正則化の該当箇所

# TV regularization
grad_x = np.gradient(im_deconv, axis=1)
grad_y = np.gradient(im_deconv, axis=0)
grad_norm = np.sqrt(grad_x**2+grad_y**2) + eps  # avoid 0 divisions
tv_reg = np.gradient(grad_x / grad_norm, axis=1) + np.gradient(grad_y / grad_norm, axis=0)

im_deconv = im_deconv/(1-tv_lambda*tv_reg)*convolve(relative_blur, psf_mirror, mode='same')

1. grad_x, grad_yで、真の画像の推定画像に対して、各軸の勾配を求めている。数式の$\nabla W_{i,r}$に対応している。
2. grad_normで、各軸の勾配を規格化する変数を用意している。ゼロ除算を避けるため、非常に小さい値(eps)を加えている。数式の$|\nabla W_{i,r}|$に対応している。
3. tv_regで、TV正則化演算子を求めている。数式の$\textrm{div}\left( \frac{\nabla W_{i,r}}{|\nabla W_{i,r}|}\right)$に対応している。
4. im_deconvで、TV正則化tv_reg演算子に重みをtv_lambdaでかけている。数式の$\frac{W_{i,r}}{1-\lambda_{TV}\textrm{div}\left( \frac{\nabla W_{i,r}}{|\nabla W_{i,r}|}\right)}\sum_k \frac{P_{i,k}H_k}{\sum_j P_{j,k}W_{j,r}}$に対応している。

結果

結果

# 各パラメータ
true_img_fname = 'lena.png'  # シミュレーションに使う真の画像
psf_px = 51  # PSF用のガウシアンフィルタのピクセルサイズ(中心座標が一意のピクセルとなる奇数推奨)
psf_sigma = 6  # PSF用のガウシアンフィルタのσ
iterations = 200  # RL法の反復回数
true_img = cv2.imread(true_img_fname)
true_img = cv2.cvtColor(true_img, cv2.COLOR_BGR2GRAY)

# PSF
psf = get_psf(psf_px, psf_sigma)

# 真の画像をPSFで畳み込む
blur_img = cv2.filter2D(true_img, -1, psf)

# 真の画像をPSFで畳み込んだ画像に基づいて、ポアソンノイズを入れた画像の作成
np.random.seed(2023)
poisson_func = lambda x: np.random.poisson(x)
blurNoisy_img = poisson_func(blur_img)
blurNoisy_img[blurNoisy_img > 255] = 255
blurNoisy_img = blurNoisy_img.astype(np.uint8)

# 標準的なRL法(tv_lambda=0)
rl_img = richardson_lucy_tv_reg(blurNoisy_img/255., psf, iterations, tv_lambda=0)
rl_img *= 255.
rl_img = rl_img.astype(np.uint8)

# TV正則化を用いたRL法(tv_lambda=0.002)
rl_tv_img = richardson_lucy_tv_reg(blurNoisy_img/255., psf, iterations, tv_lambda=0.002)
rl_tv_img *= 255.
rl_tv_img = rl_tv_img.astype(np.uint8)

fig, axs = plt.subplots(2, 2, figsize=(10, 11))
axs[0, 0].set_title('The true image')
axs[0, 0].imshow(true_img, cmap='gray')
axs[0, 1].set_title(f'The blur and poisson noise image')
axs[0, 1].imshow(blurNoisy_img, cmap='gray')
axs[1, 0].set_title(f'Richardson-Lucy deconvolution\niterations={iterations}, TV_reg=0')
axs[1, 0].imshow(rl_img, cmap='gray')
axs[1, 1].set_title(f'Richardson-Lucy deconvolution\niterations={iterations}, TV_reg=0.002')
axs[1, 1].imshow(rl_tv_img, cmap='gray')

plt.show()

• 上段左が真の画像
• 上段右が真の画像をガウシアンフィルタで畳み込んだ後、その画素値に基づいてポアソンノイズを入れた画像(撮像画像)
• 下段左が撮像画像をTV正則化0のRL法(標準的なRL法)の反復回数200回の結果
• 下段右が撮像画像をTV正則化0.002のRL法の反復回数200回の結果

このようにTV正則化を用いた下段右のRL法の方が、ポアソンノイズによる影響が軽減されているのが確認できる。

補足

TV正則化の重み

TV正則化の重み

原論文[4]では、TV正則化の重みが0.002としていたが、別の重みでの結果も紹介する。

TV正則化の重み表示コード

TV正則化の重み表示コード

# 各パラメータ
true_img_fname = 'lena.png'  # シミュレーションに使う真の画像
psf_px = 51  # PSF用のガウシアンフィルタのピクセルサイズ(中心座標が一意のピクセルとなる奇数推奨)
psf_sigma = 6  # PSF用のガウシアンフィルタのσ
iterations = 200  # RL法の反復回数
true_img = cv2.imread(true_img_fname)
true_img = cv2.cvtColor(true_img, cv2.COLOR_BGR2GRAY)

# PSF
psf = get_psf(psf_px, psf_sigma)

# 真の画像をPSFで畳み込む
blur_img = cv2.filter2D(true_img, -1, psf)

# 真の画像をPSFで畳み込んだ画像に基づいて、ポアソンノイズを入れた画像(撮像画像)
np.random.seed(2023)
poisson_func = lambda x: np.random.poisson(x)
blurNoisy_img = poisson_func(blur_img)
blurNoisy_img[blurNoisy_img > 255] = 255
blurNoisy_img = blurNoisy_img.astype(np.uint8)

# 標準的なRL法(tv_lambda=0)
rl_img = richardson_lucy_tv_reg(blurNoisy_img/255., psf, iterations, tv_lambda=0)
rl_img *= 255.
rl_img = rl_img.astype(np.uint8)

# TV正則化を用いたRL法(tv_lambda=0.0002)
rl_tv00002_img = richardson_lucy_tv_reg(blurNoisy_img/255., psf, iterations, tv_lambda=0.0002)
rl_tv00002_img *= 255.
rl_tv00002_img = rl_tv00002_img.astype(np.uint8)

# TV正則化を用いたRL法(tv_lambda=0.002)
rl_tv0002_img = richardson_lucy_tv_reg(blurNoisy_img/255., psf, iterations, tv_lambda=0.002)
rl_tv0002_img *= 255.
rl_tv0002_img = rl_tv0002_img.astype(np.uint8)

# TV正則化を用いたRL法(tv_lambda=0.02)
rl_tv002_img = richardson_lucy_tv_reg(blurNoisy_img/255., psf, iterations, tv_lambda=0.02)
rl_tv002_img *= 255.
rl_tv002_img = rl_tv002_img.astype(np.uint8)


fig, axs = plt.subplots(3, 2, figsize=(10, 16))
axs[0, 0].set_title('The true image')
axs[0, 0].imshow(true_img, cmap='gray')
axs[0, 1].set_title(f'The blur and poisson noise image')
axs[0, 1].imshow(blurNoisy_img, cmap='gray')
axs[1, 0].set_title(f'Richardson-Lucy deconvolution\niterations={iterations}, TV_reg=0')
axs[1, 0].imshow(rl_img, cmap='gray')
axs[1, 1].set_title(f'Richardson-Lucy deconvolution\niterations={iterations}, TV_reg=0.0002')
axs[1, 1].imshow(rl_tv00002_img, cmap='gray')
axs[2, 0].set_title(f'Richardson-Lucy deconvolution\niterations={iterations}, TV_reg=0.002')
axs[2, 0].imshow(rl_tv0002_img, cmap='gray')
axs[2, 1].set_title(f'Richardson-Lucy deconvolution\niterations={iterations}, TV_reg=0.02')
axs[2, 1].imshow(rl_tv002_img, cmap='gray')

plt.show()

TV正則化が大きくなるほどポアソンノイズが軽減されるが、TV正則化の重みが0.2のように大きい場合、微細な構造が平滑化されてしまうので注意したい。
また、この画像のTV正則化のパラメータ毎の結果を見る限りは、TV正則化の重みが0.002がノイズも軽減されバランスの取れた鮮明な画像になっている。

PSFの広がり具合を変えた結果

PSFの広がり具合を変えた結果

PSFの広がり具合によるTV正則化による効果を見たいので、psf_sigmaを5, 6, 7, 8に変更した結果を紹介する。

psf_sigma = 5の結果

psf_sigma = 6の結果

psf_sigma = 7の結果

psf_sigma = 8の結果

どの広がりのPSFにおいても、TV正則化を用いた方がノイズを軽減されていることがわかる。

定量評価

定量評価

標準的なRL法とTV正則化を加えたRL法の反復回数200回の結果について、どの結果が鮮明かを定量的に評価したい。そこで、真の画像との類似度比較をPeak signal-to-noise ratio (PSNR), Structural similarity index measure (SSIM)を用いて行う。PSNRとSSIMの違いは、画質の客観評価（SNR, PSNR, SSIM）についてによくまとまっているので興味のある方はご覧いただきたい。
2つの指標の違いを簡単に説明すると、以下の通りである。

• PSNRは画像全体のみならず局所領域などの劣化によって数値が大きく影響される場合があり、必ずしも人間の知覚と一致しない。
• SSIMは、輝度、コントラスト、構造の3種類の項を組み合わせるため、より人間の知覚に近いような指標である。

2つの指標はいずれも高いほど2つの画像間の類似度が高い傾向にあり、すなわち鮮明になっていることを評価できる。

PSNR, SSIMの計算は、OpenCVのcv2.PSNR(), cv2.quality.QualitySSIM_compute()を用いる。また、RL法の結果の画像の縁に歪みがあるため、画像の縁それぞれ50ピクセルを除いた領域で計算している。

PSNRによる評価(数字が大きい方が真の画像に近い)

psf_sigma	tv_lambda=0(標準的なRL法)	tv_lambda=0.0002	tv_lambda=0.002	tv_lambda=0.02
5	23.26	23.49	24.46	24.47
6	23.20	23.35	23.94	23.77
7	22.90	23.00	23.37	23.19
8	22.41	22.49	22.85	22.73

psf_sigma=5以外では、tv_lambda=0.002の画像が一番PSNRが大きくなり、鮮明になっていることがわかる。

SSIMによる評価(数字が大きい方が真の画像に近い)

psf_sigma	tv_lambda=0(標準的なRL法)	tv_lambda=0.0002	tv_lambda=0.002	tv_lambda=0.02
5	0.52	0.54	0.63	0.62
6	0.55	0.56	0.62	0.60
7	0.56	0.57	0.60	0.58
8	0.56	0.56	0.59	0.57

tv_lambda=0.002の画像が一番PSNRが大きくなり、鮮明になっていることがわかる。

PSNR, SSIMから、今回のサンプル画像のケースでは、tv_lambda=0.002が最も真の画像に近い傾向があることがわかる。

TV正則化の説明

TV正則化の説明

TV正則化

\textrm{div}\left( \frac{\nabla W_{i,r}}{|\nabla W_{i,r}|}\right)

とは何かについて、可視化して簡単に説明する。

TV正則化の説明画像のコード

TV正則化の説明画像

def func(x, y):
    z = 2*abs((x*np.sin(0.1*x) + y*np.cos(0.1*y)))
    return z

# 適当な画像
y_shape, x_shape = (64, 64)
true_img = np.empty((y_shape, x_shape))
for y in range(y_shape):
    for x in range(x_shape):
        true_img[y, x] = func(x, y)

# ポアソンノイズを追加
poisson_func = lambda x: np.random.poisson(x)
noisy_img = poisson_func(true_img)

# ポアソンノイズのない画像のTV regularization
grad_x = np.gradient(true_img, axis=1)
grad_y = np.gradient(true_img, axis=0)
grad_norm = np.sqrt(grad_x**2+grad_y**2) + 1e-12  # avoid 0 divisions
tv_reg_true_img = np.gradient(grad_x / grad_norm, axis=1) + np.gradient(grad_y / grad_norm, axis=0)

# ポアソンノイズを加えた画像のTV regularization
grad_x = np.gradient(noisy_img, axis=1)
grad_y = np.gradient(noisy_img, axis=0)
grad_norm = np.sqrt(grad_x**2+grad_y**2) + 1e-12  # avoid 0 divisions
tv_reg_noisy_img = np.gradient(grad_x / grad_norm, axis=1) + np.gradient(grad_y / grad_norm, axis=0)

fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].set_title(f'True image')
axs[0, 0].imshow(true_img, cmap='CMRmap')
axs[0, 1].set_title('TV regularization of True image')
axs[0, 1].imshow(tv_reg_true_img, vmin=-np.max(np.abs(tv_reg_true_img)), vmax=np.max(np.abs(tv_reg_true_img)),cmap='bwr')
axs[1, 0].set_title(f'Noisy image')
axs[1, 0].imshow(noisy_img, cmap='CMRmap')
axs[1, 1].set_title('TV regularization of Noisy image')
axs[1, 1].imshow(tv_reg_noisy_img, vmin=-np.max(np.abs(tv_reg_noisy_img)), vmax=np.max(np.abs(tv_reg_noisy_img)),cmap='bwr')

plt.show()

• 上段左が真の画像(青から赤にかけて値が大きい)
• 上段右が真の画像のTV正則化の結果(青=負, 赤=正を表す)
• 下段左が真の画像の画素値に基づいたポアソンノイズを追加した画像
• 下段右が真の画像のTV正則化の結果

右図の赤は周辺の傾きが増加した領域、青は周辺の領域の傾きが減少した領域を表している。
TV正則化を用いたRL法では、赤い領域のRL法の反復回数後の値の更新を大きくし、反対に青い領域の反復回数後の値の更新を小さくしている。
このTV正則化によってRL法の更新時に領域毎にペナルティを付けることができ、良い感じのTV正則化の重みのときに、ポアソンノイズによる影響を抑えた鮮明な画像を生成することを可能にする。

一応、本記事のサンプル画像でも同様にTV正則化の結果を見てみましょう。

サンプル画像のTV正則化のコード

サンプル画像のTV正則化

true_img_fname = 'lena.png'  # シミュレーションに使う真の画像
psf_px = 51  # PSF用のガウシアンフィルタのピクセルサイズ(中心座標が一意のピクセルとなる奇数推奨)
psf_sigma = 6  # PSF用のガウシアンフィルタのσ
iterations = 20  # RL法の反復回数
true_img = cv2.imread(true_img_fname)
true_img = cv2.cvtColor(true_img, cv2.COLOR_BGR2GRAY)

# PSF
psf = get_psf(psf_px, psf_sigma)

# 真の画像をPSFで畳み込む
blur_img = cv2.filter2D(true_img, -1, psf)

# 真の画像をPSFで畳み込んだ画像に基づいて、ポアソンノイズを入れた画像(撮像画像)
np.random.seed(2023)
poisson_func = lambda x: np.random.poisson(x)
blurNoisy_img = poisson_func(blur_img)
blurNoisy_img[blurNoisy_img > 255] = 255
blurNoisy_img = blurNoisy_img.astype(np.uint8)

# TV正則化を用いたRL法(tv_lambda=0.002)
rl_tv_img = richardson_lucy_tv_reg(blurNoisy_img/255., psf, iterations, tv_lambda=0.002)
rl_tv_img *= 255.
rl_tv_img = rl_tv_img.astype(np.uint8)

# TV regularization
grad_x = np.gradient(rl_tv_img, axis=1)
grad_y = np.gradient(rl_tv_img, axis=0)
grad_norm = np.sqrt(grad_x**2+grad_y**2) + 1e-12  # avoid 0 divisions
tv_reg = np.gradient(grad_x / grad_norm, axis=1) + np.gradient(grad_y / grad_norm, axis=0)

fig, axs = plt.subplots(1, 2, figsize=(10, 11))
axs[0].set_title(f'Richardson-Lucy deconvolution\niterations={iterations}, TV_reg=0.002')
axs[0].imshow(rl_tv_img, cmap='gray')
axs[1].set_title('TV regularization of left image')
axs[1].imshow(tv_reg, vmin=-np.max(np.abs(tv_reg)), vmax=np.max(np.abs(tv_reg)),cmap='bwr')

plt.show()

良い感じに輪郭が赤くなっていていそうですね。

まとめ

RL法はノイズを増幅してしまうという課題があり、特にポアソンノイズの場合、TV正則化を用いたRL法が有効であることが知られている。公にTV正則化を用いたRL法のコードを見かけなかったので、本記事で紹介しました。
また、補足ではTV正則化のパラメータを適当に変化させた結果や定量評価を行い、TV正則化の重みが0.002で真の画像に最も近い画像が得られることを確認しました。

参考資料

参考文献

[1] William Hadley Richardson. Bayesian-based iterative method of image restoration. JoSA, Vol. 62, No. 1, pp. 55–59, 1972.
[2] Leon B Lucy. An iterative technique for the rectification of observed distributions. The astronomical journal, Vol. 79, p. 745, 1974.
[3] Leonid I Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: nonlinear phenomena, Vol. 60, No. 1-4, pp.259–268, 1992.
[4] Nicolas Dey, Laure Blanc-Feraud, Christophe Zimmer, Pascal Roux, Zvi Kam, Jean-Christophe Olivo-Marin, and Josiane Zerubia. Richardson–lucy algorithm with total variation regularization for 3d confocal microscope deconvolution. Microscopy research and technique, Vol. 69, No. 4, pp. 260–266, 2006