More than 1 year has passed since last update.

Richardson-Lucy deconvolutionをRGB画像で実装しよう

Last updated at 2023-05-13Posted at 2022-08-30

はじめに

Richardson-Lucy deconvolution（以下、RL法）は、画像等を鮮明にする手法である。一般的な撮像観測では、観測装置の回折等により光源が拡がって撮像される。この拡がり具合は、点拡がり関数(Point Spread Function, 以下、PSF)で表される。RL法は既知のPSFと撮像画像から、真の画像を推定する手法である。

本記事ではまず一般的なRL法のRGB画像での実装方法を紹介する。そして、RL法には画像の縁に歪みが出る性質があるが、その縁の歪みを軽減する方法についても紹介する。また、本記事での撮像画像は適当なPSFで畳み込むことで擬似的に作成する。

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

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

実装

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

各種インポート

import cv2
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import multivariate_normal
from skimage import restoration

実行時の各種バージョン

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

※scikit-imageのrestoration.richardson_lucyは、Version 0.19.1で一部変更されている。（変更箇所概要：Add small regularization to avoid zero-division in richardson_lucy (gh-5976))

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は何でも良いが、今回は多変量正規分布に従うPSFで実装する。多変量正規分布にはSciPyのmultivariate_normalを使用する。

一般的なRL法

PSF

ポイントとしては、PSFのサイズは必ず奇数にしてPSFの中心が画像の中央に来るようにしてください。RL法はPSFの畳み込み計算箇所と180°回転したPSFを使った畳み込み計算箇所の2つがあり、PSFのサイズが奇数でない場合は中心で一意に定義されないのでRL法が正しく行えない可能性がある。
(私自身この事は経験して気づいたのですが、プログラム上はPSFサイズが偶数でも画像サイズの半分+1で中心が定義されているようで畳み込み自体は行えてしまうが、その場合180°回転したPSFと中心が異なってしまうので注意してください。)

psf_px = 51
psf_sigma = 7
psf = get_psf(psf_px, psf_sigma)

fig, ax = plt.subplots()
ax.imshow(psf)
plt.show()

真の画像

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

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

# 真の画像の読み込み
img = cv2.imread('lena.png')

fig, ax = plt.subplots()
img_display = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
ax.imshow(img_display)
plt.show()

撮像画像

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

fig, ax = plt.subplots()
blur_img_display = cv2.cvtColor(blur_img, cv2.COLOR_BGR2RGB)
ax.imshow(blur_img_display )
plt.show()

RL法の結果

ポイントはrestoration.richardson_lucy()は、キーワード引数にclipがあり、この引数は-1より小さいあるいは1より大きい値が出た場合にそれぞれ-1, 1に置換が行われる。デフォルトではclip=Trueとなっていて、この機能を使うために0-1の範囲になるように255.で割っている。また、255.で割ることでnumpyのfloat64型に自動で変換されるので型変換も担っている。

# [0-1]float64画像に変換してから各チャンネルごとでRL法
iterations = 50
rl_img = np.zeros_like(blur_img, dtype=np.float64)
rl_img[..., 0] = restoration.richardson_lucy(blur_img[..., 0]/255., psf, iterations)
rl_img[..., 1] = restoration.richardson_lucy(blur_img[..., 1]/255., psf, iterations)
rl_img[..., 2] = restoration.richardson_lucy(blur_img[..., 2]/255., psf, iterations)

# [0-1]float64画像から[0-255]uint8画像に変換
rl_img *= 255.
rl_img = rl_img.astype(np.uint8)

fig, ax = plt.subplots()
rl_img_display = cv2.cvtColor(rl_img, cv2.COLOR_BGR2RGB)
ax.imshow(rl_img_display)
plt.show()

一般的なRL法では、エッジが目立ってしまう。そこで、エッジを目立たなくすために前処理として撮像画像の縁の画像値でpaddingをしてRL法を行う。

撮像画像の縁の画像値でpaddingの前処理を導入したRL法

撮像画像を縁の画像値でpaddingした撮像画像

# 真の画像の読み込み
img = cv2.imread('lena.png')

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

# 撮像画像を縁の画素値でpadding
pad_px = 51
pad = (((pad_px, pad_px), (pad_px, pad_px), (0, 0)))
blur_pad_img = np.pad(blur_img, pad, 'edge')

fig, ax = plt.subplots()
blur_pad_img = cv2.cvtColor(blur_pad_img, cv2.COLOR_BGR2RGB)
ax.imshow(blur_pad_img)
plt.show()

撮像画像を縁の画像値でpaddingすることで、各縁のピクセル数を51ピクセル追加してもそこまで違和感を感じないと思う。よく見ると右肩が伸びているのが少し不自然だがそこまででもない。この画像で同様にRL法を行う。

paddingした撮像画像でRL法した結果

# [0-1]float64画像に変換してから各チャンネルごとでRL法
rl_pad_img = np.zeros_like(blur_pad_img, dtype=np.float64)
rl_pad_img[..., 0] = restoration.richardson_lucy(blur_pad_img[..., 0]/255., psf, iterations)
rl_pad_img[..., 1] = restoration.richardson_lucy(blur_pad_img[..., 1]/255., psf, iterations)
rl_pad_img[..., 2] = restoration.richardson_lucy(blur_pad_img[..., 2]/255., psf, iterations)

# [0-1]float64画像から[0-255]uint8画像に変換
rl_pad_img *= 255.
rl_pad_img = rl_pad_img.astype(np.uint8)

fig, ax = plt.subplots()
rl_pad_img_display = cv2.cvtColor(rl_pad_img, cv2.COLOR_BGR2RGB)
ax.imshow(rl_pad_img)
plt.show()

やはり画像の縁は歪んでしまうが、そもそも追加したものであるから追加した縁を取り除くとどうだろうか。

paddingした縁を取り除いた後の画像

rl_img = rl_pad_img[pad_px:-pad_px, pad_px:-pad_px]

fig, ax = plt.subplots()
rl_img_display = cv2.cvtColor(rl_img, cv2.COLOR_BGR2RGB)
ax.imshow(rl_img)
plt.show()

縁の歪みはほとんど目立たなくなるのがわかる。このように、適当に縁の画素値でpaddingすることでRL法の縁の歪みを目立たなくして見た目を良くすることができる。

コード全容

先ほどの撮像画像の縁の画素値をpaddingしたRL法のコードをまとめて関数化したので参考にしてください。RGB画像だけでなくグレースケール画像でのコードを用意している。

撮像画像の縁の画素値をpaddingしたRL法

RGB画像用

RGB画像用のRL法の関数

def rgb_img_RL(blur_img, psf, iterations, pad_px=51):
    """
    input
    blur_img: blurred 3-channel image (3D array)
    psf: Poinst spread function (2D array)
    iterations: Richardson-Lucy deconvolution iterations
    pad_px: number of pixels to pad `blur_img`
    output
    rl_img: image with discarded padded edge pixels after Richardson-Lucy deconvolution
    description
    This Richardson-Lucy deconvolutionis used for RGB 3-channel image 
    with edge distortion reduced by edge padding processing.
    """
    pad = (((pad_px, pad_px), (pad_px, pad_px), (0, 0)))
    blur_pad_img = np.pad(blur_img, pad, 'edge')

    rl_pad_img = np.zeros_like(blur_pad_img, dtype=np.float64)
    rl_pad_img[..., 0] = restoration.richardson_lucy(blur_pad_img[..., 0]/255., psf, iterations)
    rl_pad_img[..., 1] = restoration.richardson_lucy(blur_pad_img[..., 1]/255., psf, iterations)
    rl_pad_img[..., 2] = restoration.richardson_lucy(blur_pad_img[..., 2]/255., psf, iterations)

    rl_pad_img *= 255.
    rl_pad_img = rl_pad_img.astype(np.uint8)
    rl_img = rl_pad_img[pad_px:-pad_px, pad_px:-pad_px]

    return rl_img

pad_pxは撮像画像やPSFのサイズに合わせて適宜調整してください。このpad_pxは大きければ大きいほど画像の歪みが軽減されると予想されるが、計算時間がかかるため筆者は目安としてPSFのサイズ程度を使うようにしている。理由としては、RL法の計算上ではPSFのサイズ以上に拡がるが、現実的にはPSFのサイズ以上に拡がる量はほとんどないためである。

RGB画像用の関数を使った結果

# PSF
psf_px = 51
psf_sigma = 7
psf = get_psf(psf_px, psf_sigma)

# 真の画像の読み込み
img = cv2.imread('lena.png')

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

# 画像の縁を歪みを目立たなくするRL法
iterations = 50
rl_img = rgb_img_RL(blur_img, psf, iterations, pad_px=51)

fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(18, 6))
ax1.set_title('The true image')
img_display = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
ax1.imshow(img_display)
ax1.axis('off')
ax2.set_title('The obtained image')
blur_img_display = cv2.cvtColor(blur_img, cv2.COLOR_BGR2RGB)
ax2.imshow(blur_img_display)
ax2.axis('off')
rl_img_display = cv2.cvtColor(rl_img, cv2.COLOR_BGR2RGB)
ax3.set_title(f'Richardson-Lucy deconvolution iterations={iterations}')
ax3.imshow(rl_img_display)
ax3.axis('off')
plt.show()

グレースケール画像用

グレースケール画像用のRL法の関数

def gray_img_RL(blur_img, psf, iterations, pad_px=51):
    """
    input
    blur_img: blurred 1-channel image (2D array)
    psf: Poinst spread function (2D array)
    iterations: Richardson-Lucy deconvolution iterations
    pad_px: number of pixels to pad `blur_img`
    output
    rl_img: image with discarded padded edge pixels after Richardson-Lucy deconvolution
    description
    This Richardson-Lucy deconvolutionis used for gray 1-channel image 
    with edge distortion reduced by edge padding processing.
    """
    pad = (((pad_px, pad_px), (pad_px, pad_px)))
    blur_pad_img = np.pad(blur_img, pad, 'edge')

    rl_pad_img = np.zeros_like(blur_pad_img, dtype=np.float64)
    rl_pad_img = restoration.richardson_lucy(blur_pad_img/255., psf, iterations)

    rl_pad_img *= 255.
    rl_pad_img = rl_pad_img.astype(np.uint8)
    rl_img = rl_pad_img[pad_px:-pad_px, pad_px:-pad_px]

    return rl_img

グレースケール画像用の関数を使った結果

# PSF
psf_px = 51
psf_sigma = 7
psf = get_psf(psf_px, psf_sigma)

# 真の画像の読み込み
img = cv2.imread('lena.png', cv2.IMREAD_GRAYSCALE)

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

# 画像の縁を歪みを目立たなくするRL法
iterations = 50
rl_img = gray_img_RL(blur_img, psf, iterations, pad_px=51)

fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(18, 6))
ax1.set_title('The true image')
ax1.imshow(img, vmin=0, vmax=255, cmap='gray')
ax1.axis('off')
ax2.set_title('The obtained image')
ax2.imshow(blur_img, vmin=0, vmax=255, cmap='gray')
ax2.axis('off')
ax3.set_title(f'Richardson-Lucy deconvolution iterations={iterations}')
ax3.imshow(rl_img, vmin=0, vmax=255, cmap='gray')
ax3.axis('off')
plt.show()

まとめ

RL法をRGB画像で実装する方法と、RL法の問題点の1つである画像の縁の歪みを軽減させる方法の1例を紹介しました。ただ、本手法は撮像画像の縁の画素値をpaddingで補完しているだけで、本当は存在しない画像情報を使ってのRL法なので真に正しい手法ではないことは抑えておきましょう。

参考文献

Richardson, William Hadley. 1972, JOSA, 62, 55–59
Lucy, L. B. 1974, Astronomical Journal, 79, 745–754