はじめに
一般的な撮像観測では、観測装置の回折等により光源が拡がって（畳み込まれて）撮像される。この拡がり具合は、点拡がり関数（Point Spread Function, 以下、PSF）で表される。Deconvolution手法の目標は、PSFで畳み込まれる前の真の画像を推定することである。
その手法には大きく二つある。まず、PSFが入手できるのであれば、Richardson-Lucy deconvolution（以下、RL法）[1, 2]やWiener deconvolutionなどの古典的に有名な手法がある。
一方、多くの場合PSFの入手は困難なため、PSFを用いないで撮像画像のみから真の画像を推定するBlind deconvolutionと呼ばれる手法もある。その一つにRL法に基づいたBlind deconvolution（以下、Blind-RL法）[3]がある。
そのBlind-RL法だが、サンプルコードを見かけなかったので、本記事で実装する。そして、Blind-RL法の注意点などを補足で紹介する。
RL法とは
RL法[1, 2]は、既知の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$は2次元のピクセル座標とする。$W$は真の画像、$H$は撮像画像、$r$は反復回数とする。$P_{j,k}$はPSFで、$W_j$が$H_k$で観測される確率を$P(H_k|W_j)=P_{j,k}$とする。式の導出は、私の記事で恐縮ですが、Richardson-Lucy deconvolutionの理論式を導出しようを参照されたい。
本記事でBlind-RL法[3]を説明するにあたり、RL法を畳み込みを使って、
\begin{aligned}
W_{r+1}(x)=W_r(x)\left\{\left[\frac{H(x)}{W_r(x) \otimes P(x)}\right] \otimes P(-x)\right\}\ \quad r = 0,\ 1,\ 2,\ \dots
\end{aligned}
と表せる。ここで、$x$は2次元のピクセル座標である。
Blind-RL法とは
Blind-RL法[3]は、基本的にRL法の概念が使われる。定性的に理解するために、まずは下図をご覧いただきたい。
数学的に撮像画像は、真の画像とPSFの畳み込みでかけ、任意の$f$と$g$は交換律より$f\otimes g = g\otimes f$が成り立つ。そのため、上図のように真の画像とPSFの対応が反対になっていても数学的に問題はない。このことを利用すると、真の画像が完全に分かっている場合は、RL法でPSFを推定することも可能である。
(また、物理的には、真の画像とPSFが完全に分からない場合は、上図右側のように真の画像が拡がっていてPSFが点が二つの可能性もありうることを意味している。)
Blind-RL法は、真の画像とPSFの初期値を与えて、反復回数毎にRL法により真の画像とPSFを交互に同時に推定していく手法である。数式で
\begin{aligned}
P_{r+1}(x)&=P_r(x)\left\{\left[\frac{H(x)}{P_r(x) \otimes W_r(x)}\right] \otimes W_r(-x)\right\} \\
W_{r+1}(x)&=W_r(x)\left\{\left[\frac{H(x)}{W_r(x) \otimes P_{r+1}(x)}\right] \otimes P_{r+1}(-x)\right\}\ \quad r = 0,\ 1,\ 2,\ \dots
\end{aligned}
と表される。$W$は真の画像、$H$は撮像画像、$r$は反復回数、$x$は2次元のピクセル座標とする。
実装
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作成用の関数
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は何でも良いが、本記事ではガウシアンフィルタで実装する。
Blind-RL法
def blind_richardson_lucy(image, psf_shape=(15, 15), iterations=50):
if image.shape < psf_shape:
print(f'Error: Specify the size of psf shape {psf_shape} to be less than or equal to that of image shape {image.shape}.')
if psf_shape[0] % 2 == 0 or psf_shape[1] % 2 == 0:
print('Error: The psf must be specified with odd size.')
# Set data type
float_type = np.promote_types(image.dtype, np.float32)
# Convert image size to odd size for convolution
_pad_h, _pad_w = (0, 0)
image = image.astype(float_type, copy=True)
if image.shape[0] % 2 == 0:
_pad_h = 1
if image.shape[1] % 2 == 0:
_pad_w = 1
image = np.pad(image, ((0, _pad_h), (0, _pad_w)), 'constant')
# Set initial estimated deconvolution image and PSF
im_deconv = np.full(image.shape, 0.5, dtype=float_type)
psf = np.full(psf_shape, 0.5, dtype=float_type)
# Regularization by psf region
_trim_h = range((image.shape[0]-psf_shape[0])//2, (image.shape[0]+psf_shape[0])//2)
_trim_w = range((image.shape[1]-psf_shape[1])//2, (image.shape[1]+psf_shape[1])//2)
_trim = np.ix_(_trim_h, _trim_w)
# Small regularization parameter used to avoid 0 divisions
eps = 1e-12
# Blind-RL deconvolution
for _ in range(iterations):
# Estimate psf
conv = convolve(im_deconv, psf, mode='same') + eps
relative_blur = image / conv
im_deconv_mirror = np.flip(im_deconv)
psf *= convolve(relative_blur, im_deconv_mirror, mode='same')[_trim]
# Estimate im_deconv
conv = convolve(im_deconv, psf, mode='same') + eps
relative_blur = image / conv
psf_mirror = np.flip(psf)
im_deconv *= convolve(relative_blur, psf_mirror, mode='same')
# Remove padded areas
im_deconv = im_deconv[:image.shape[0]-_pad_h, :image.shape[1]-_pad_w]
# Normalize psf
psf /= np.sum(psf)
return im_deconv, psf
入力
-
image: 撮像画像 -
psf_shape: PSFの縦横サイズ(サイズは中心が一意に決まるように必ず奇数とし、撮像画像よりも大きいサイズはPSFの推定が困難であるため指定しない。実際、本記事ではそのサイズを仕様上指定できない。) -
iteration: 反復回数
出力
-
im_deconv: 推定された真の画像 -
psf: 推定されたPSF
コードの説明
本記事のRL法のコードは、scikit-image version0.18.3, 0.19.2のrichardson_lucy()を参考にしている。一部変数名が重複してしまいバグを生みそうで申し訳ないが、同じような処理をする箇所のためあえて変数名を変更しなかった。
-
# Convert image size to odd size for convolutionからのブロックは、imageの画像サイズが奇数になるように調整している。psfのサイズに限らず、なぜimageも奇数になるようにしているかは、畳み込みの座標が一意に決まるようにしたいためである。コードの下にim_deconv_mirror = np.flip(im_deconv)やpsf_mirror = np.flip(psf)にあるように、180°回転した画像での畳み込みも計算される。そのため、回転前の畳み込みと回転後での畳み込みでピクセルの中心座標が異なってしまうのを回避するために、畳み込みで使う全ての画像は奇数にする必要がある。 -
# Set initial estimated deconvolution image and PSFのブロックは、im_deconv、psfの初期値を適当に全てフラットな値を入れている。「全部0.5で入れて画素値の合計がおかしくならないの？」と思われる方もいると思うが、RL法は分母分子で打ち消し合うので、適当な値を入れてもimageと合計が同じになるので問題はない。ただ、Blind-RL法では、psfもRL法の式で推定するので、imageの合計と同じ値になりpsf合計が1にならないので、コードの最後に規格化する必要がある。 -
# Regularization by psf regionのブロックは、Blind-RL法そのままでは、画像のサイズが撮像画像と同じ大きさになってしまい精度が落ちるため、psfを適当なサイズにトリミングするために使われる。なお、トリミングは画像の中心の座標からpsf_shapeで指定した範囲の領域を切り取っている。 -
im_deconv_mirrorは、im_deconv180度回転した画像で、畳み込み処理の$\otimes W_r(x)$用である。 - 1つ目の
convは、$P_r(x) \otimes W_r(x)$を計算している。ゼロ除算を避けるため、非常に小さい値(eps)を加えている。 - 1つ目の
relative_blurは、$\frac{H(x)}{P_r(x) \otimes W_r(x)}$を計算している。 -
psf *= convolve(relative_blur, im_deconv_mirror, mode='same')[_trim]は、$P_r(x){\left[\frac{H(x)}{P_r(x) \otimes W_r(x)}\right] \otimes W_r(-x)}$を計算している。_trimにより、psfが中心の座標からpsf_shapeのサイズになるようにトリミングしている。 -
# Estimate im_deconvのブロックも# Estimate psfのブロックとほとんど同様の処理である。 -
im_deconv = im_deconv[:image.shape[0]-_pad_h, :image.shape[1]-_pad_w]は、# Convert image size to odd size for convolutionのブロックで画像サイズが奇数になるように追加した結果、im_deconvも同じサイズの奇数になっているため、その領域を取り除いている。 -
psf /= np.sum(psf)は、psfを大きさが1になるように規格化している。なお、# Set initial estimated deconvolution image and PSFでも説明しているが、RL法の式から反復回数毎では、psfとim_deconvは、分母分子で打ち消し合うので規格化する必要がない。im_deconvを規格化せず、psfのみ規格化する理由は、RL法は撮像画像imageとim_deconvの合計が常に同じになる手法のため、im_deconvはすでに規格化されている。一方、psfもRL法を使って推定するため撮像画像(image)と合計が同じになってしまうため、確率の総和が1になるように規格化する必要がある。
結果
# 各パラメータ
true_img_fname = 'lena.png' # シミュレーションに使う真の画像
psf_px = 51 # PSF用のガウシアンフィルタのピクセルサイズ(中心座標が一意のピクセルとなる奇数推奨)
psf_sigma = 7 # PSF用のガウシアンフィルタのσ
iterations = 50 # RL法の反復回数
true_img = cv2.imread(true_img_fname)
true_img = cv2.cvtColor(true_img, cv2.COLOR_BGR2GRAY)
# PSF
true_psf = get_psf(psf_px, psf_sigma)
# 真の画像をPSFで畳み込む
blur_img = convolve(true_img, true_psf, 'full')
# Blind-RL法
psf_shape = (19, 19)
blind_rl_img, estimate_psf = blind_richardson_lucy(blur_img/255., psf_shape, iterations)
blind_rl_img[blind_rl_img > 1] = 1
blind_rl_img *= 255.
blind_rl_img = blind_rl_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_visible(False)
axs[1, 0].set_title('The true psf')
axs[1, 0].imshow(true_psf, cmap='gray')
axs[1, 1].set_title(f'The blur image by the true psf')
axs[1, 1].imshow(blur_img, cmap='gray')
axs[2, 0].set_title(f'Blind-RL deconvolution\niterations={iterations}, psf_shape={psf_shape}, estimate psf')
axs[2, 0].imshow(estimate_psf, cmap='gray')
axs[2, 1].set_title(f'Blind-RL deconvolution\niterations={iterations}, deconvolved image')
axs[2, 1].imshow(blind_rl_img, cmap='gray')
plt.show()
- 上段左が真の画像
- 中段左がガウシアンのPSF
- 中段右がガウシアンの中段左のPSFで畳み込んだ撮像画像
- 下段左がBlind-RL法の反復回数50回で推定されたPSF
- 下段右がBlind-RL法の反復回数50回で推定された真の画像
Blind-RL法により撮像画像のみから、下段右のような鮮明な画像を得ることができる。
補足
Blind-RL法は、既知のPSFを使用したRL法などと比べて非常に不安定である。そこで、補足では推定するPSFのサイズ、ブレ度合い、撮像画像が孤立しているかによる依存性について紹介する。
PSFの推定サイズによる依存性
真の画像とPSFで数学的に価値が同じなので、PSFを大きくすると真の画像の情報がPSFに含まれてしまう可能性がある。それを防ぐためには、推定するPSFをある程度程度小さくしておく必要があるが、そうすると本来の真のPSFよりも小さいサイズの場合は適当なPSFを表現できなくなり鮮明にならない。そこが、Blind-RL法の一つの課題である。
図で確認
以下の図の説明
1段目: 真の画像
2段目: 左は真のPSF、右は真のPSFよってたたみ込まれた撮像画像。
3段目: PSFのサイズを(9, 9)に指定して、Blind-RL法反復回数50回の結果。左が推定されたPSF、右が推定された真の画像。
4段目: PSFのサイズを(19, 19)に指定して、Blind-RL法反復回数50回の結果。左が推定されたPSF、右が推定された真の画像。
5段目: PSFのサイズを(29, 29)に指定して、Blind-RL法反復回数50回の結果。左が推定されたPSF、右が推定された真の画像。
PSFのサイズが(19, 19)の時、ある程度推定された真の画像が鮮明であることが確認できる。このように、PSFのサイズを大きくすると、PSFの表現力が増す分過剰な真の画像が得られてしまうので、PSFの事前知識により適切なPSFサイズにしておく必要がある。
PSFのブレの大きさによる影響
一応、ガウシアンフィルタのシグマサイズを変更してブレの大きさの違いによる影響をみてみる。PSFの推定サイズによる依存性で触れたように、適切なPSFサイズはブレ度合いによって異なる。
図で確認
この節の以下の図の説明
1段目: 真の画像
2段目: 左は真のPSF、右は真のPSFよってたたみ込まれた撮像画像。
3段目: PSFのサイズを(9, 9)に指定して、Blind-RL法反復回数50回の結果。左が推定されたPSF、右が推定された真の画像。
4段目: PSFのサイズを(19, 19)に指定して、Blind-RL法反復回数50回の結果。左が推定されたPSF、右が推定された真の画像。
5段目: PSFのサイズを(29, 29)に指定して、Blind-RL法反復回数50回の結果。左が推定されたPSF、右が推定された真の画像。
真のPSFの大きさ 小(psf_sigma=4)
真のPSFの大きさ 中(psf_sigma=7)
真のPSFの大きさ 大(psf_sigma=10)
撮像画像の畳み込みの種類による影響
推定されるPSFは非常に不安定なので、畳み込みの縁を真の画像と同じサイズにするconvolve(true_img, true_psf, 'same')の場合や適当にトリミングした画像では、画像内に外部の情報が含まれるのでPSFが歪んでしまうので注意したい。
図で確認
畳み込みの縁を真の画像と同じサイズになるconvolve(true_img, true_psf, 'same')の場合、画像の縁で特殊な畳み込み処理が行われる。そのため、推定されるPSFは歪でしまうことが確認できる。比較のために忠実な畳み込みconvolve(true_img, true_psf, 'full')も載せている。
この節の以下の図の説明
- 上段左が真の画像
- 中段左がガウシアンのPSF
- 中段右がガウシアンの中段左のPSFで畳み込んだ撮像画像
- 下段左がBlind-RL法の反復回数50回で推定されたPSF
- 下段右がBlind-RL法の反復回数50回で推定された真の画像
convolve(true_img, true_psf, 'full')の場合
convolve(true_img, true_psf, 'same')の場合
撮像画像の畳み込みの外からの漏れ込みによる影響
一般に孤立した画像のみが畳み込みまれて撮影されることは稀で、無限に広がる画像の一部がカメラなどに収まっているので撮像画像には外部からの漏れ込みも入り込む。そのような撮像画像に対してBlind-RL法をしたときの影響をみたいので、撮像画像をPSFで畳み込んだ後、適当に切り出した画像で実装する。比較のため、切り取った後の孤立した画像での忠実な畳み込みの結果も載せている。
下図より、畳み込み後に切り取った画像の場合、画像の外からの画像情報も入り込んでいるため、PSFの推定が上手くできなず歪んでしまうことが確認できる。このように、孤立した撮像画像に比べて外部からの漏れ込みがある撮像画像では、真の画像とPSFの推定は難易度が高い。
図で確認
まとめ
Blind-RL法[3]は、撮像画像のみから、PSFと真の画像の数学的な価値が同じであることを利用して、反復的にPSFと真の画像を交互にRL法で推定していく手法である。本記事では、公にその手法のコードを見かけなかったので実装しました。
補足では、推定するPSFのサイズ、ブレ度合い、撮像画像が孤立しているかによる依存性について紹介した。
また、天体画像では、構造がある程度まとまった孤立した領域を見つけることができるので、この手法は有用かもしれない。
参考文献
[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] DA Fish, AM Brinicombe, ER Pike, and JG Walker. Blind deconvolution by means of the richardson–lucy algorithm. JOSA A, Vol. 12, No. 1, pp. 58–65, 1995.