はじめに
画像はしばしば「ぼやけ」ます。
たとえば天体画像では、本来は点のような天体であっても、望遠鏡や検出器の影響によって少し広がって見えます。こうしたぼやけを見ると、
元のシャープな画像に戻せないだろうか？
と思うのは自然なことです。
実際、そのような目的で発展してきたのがイメージデコンボリューション（Image Deconvolution）です。
ただし、ぼやけた画像だけを見ても、元の画像を完全に知ることはできません。そこで必要になるのが、「どのようにぼやけたのか」という情報です。
その手掛かりとしてよく使われるのが PSF（Point Spread Function）です。PSFとは、理想的には点として見える天体が、観測装置によって実際にはどのような形に広がって見えるかを表したものです。
たとえば、ある望遠鏡で点源を観測したときに、点ではなくガウシアンのような広がりを持って見えるのであれば、その広がり方がPSFです。
PSFが分かっていれば、
元画像 → PSFによってぼやける → 観測画像
という過程を数式で表すことができます。
PSFが画像全体で同じ形をしている場合、ノイズを無視すれば観測画像は
O(x,y)=I(x,y)*P(x,y)
と書けます。
ここで、$I(x,y)$ は真の画像、$P(x,y)$ は PSF、$O(x,y)$ は観測画像です。また、$*$ は畳み込みを表します。つまり、観測画像 $O(x,y)$ は、真の画像 $I(x,y)$ が PSF によってぼやけたものです。
デコンボリューションとは、このぼやけた画像から元の画像を推定する処理です。一見すると、PSFが分かっていれば元の画像を簡単に復元できそうに見えます。しかし実際には、観測にはノイズも含まれるため、単純にはうまくいきません。
今回は最も単純な方法として、フーリエ変換を使った逆畳み込みを実装しながら、なぜノイズが問題になるのかを見ていきます。
ただし、最初から2次元画像で考えると少し見通しが悪いので、まずは1次元の信号で考えてみます。
本記事の概要を簡単に示します。
実装コードは Google Colab からもお試しいただけます：
https://colab.research.google.com/drive/1NVHpDmfX03DjLB22B1QmRsKYDj58rv45?usp=sharing
1次元で考えてみる
真の信号を $I(x)$、PSF を $P(x)$、観測信号を $O(x)$ とします。ノイズを無視すれば、
O(x) = I(x) * P(x)
と書けます。
ここで大事なのは、畳み込みはフーリエ変換すると掛け算になるということです。
\tilde{O}(k) = \tilde{I}(k)\tilde{P}(k)
ここで、$\tilde{O}(k)$、$\tilde{I}(k)$、$\tilde{P}(k)$ は、それぞれ $O(x)$、$I(x)$、$P(x)$ をフーリエ変換したものです。
したがって、PSF が分かっていれば、
\tilde{I}(k) = \frac{\tilde{O}(k)}{\tilde{P}(k)}
として、元の信号を推定できます。
これが、フーリエ変換で見た一番単純なデコンボリューションです。
ノイズがない場合
まずはノイズがない場合で試してみます。ここでは、近い位置にある2つの鋭いピークを持つ信号を作り、それをガウシアン PSF でぼやけさせます。そのあと、フーリエ空間で割り算して、元の信号が戻るかを見てみます。
import numpy as np
import matplotlib.pyplot as plt
# -------------------------
# 1次元の真の信号を作ります
# -------------------------
N = 128
x = np.arange(N)
signal_true = np.zeros(N)
signal_true[60] = 0.5
signal_true[70] = 1.0
# -------------------------
# FFTで扱いやすい1次元ガウシアンPSFを作ります
# -------------------------
# FFTでは index=0 が原点として扱われます。
# そのため、PSFの中心を配列中央ではなく index=0 に置きます。
#
# np.minimum(x, N - x) にすると、
# 周期境界を考えたときに index=0 からの距離になります。
# これは FFT を使った畳み込みと相性がよい形です。
# -------------------------
sigma = 4.0
psf_x = np.minimum(x, N - x)
psf = np.exp(-psf_x**2 / (2 * sigma**2))
psf /= psf.sum()
# -------------------------
# FFTで畳み込みを計算し、信号をぼやけさせます
# -------------------------
F_true = np.fft.fft(signal_true)
F_psf = np.fft.fft(psf)
F_blur = F_true * F_psf
signal_blur = np.real(np.fft.ifft(F_blur))
# -------------------------
# フーリエ空間で割り算して戻します
# -------------------------
# ノイズがない理想的な場合は、基本的には F_psf で割れば戻ります。
# ただし、F_psf が数値的にかなり小さくなる周波数があるので、
# 非常に小さい tiny を足しています。
# -------------------------
tiny = 1e-20
F_deconv = F_blur / (F_psf + tiny)
signal_deconv = np.real(np.fft.ifft(F_deconv))
# -------------------------
# 表示します
# -------------------------
plt.figure(figsize=(7, 4))
plt.plot(signal_true, label="True signal")
plt.plot(signal_blur, label="Blurred signal")
plt.plot(signal_deconv, label="Deconvolved signal", color="gray", linewidth=2, linestyle="--")
plt.ylim(-0.1, 1.1)
plt.xlabel("x")
plt.ylabel("Value")
plt.legend(loc="upper left")
plt.tight_layout()
plt.show()
この例では、真の信号 signal_true を PSF でぼやけさせて signal_blur を作っています。そのあと、フーリエ空間で
\tilde{I} = \frac{\tilde{O}}{\tilde{P}}
を計算して、元の信号を戻しています。
ノイズがない理想的な場合には、元の鋭いピークがほぼそのまま戻ってきます。これは、ぼやけた信号を
\tilde{O} = \tilde{I}\tilde{P}
として作り、それをほぼ同じ $\tilde{P}$ で割っているからです。つまり、ノイズがなければ
\frac{\tilde{I}\tilde{P}}{\tilde{P}} = \tilde{I}
となります。
コードに tiny を加える理由
コード内の tiny = 1e-20 は、ノイズを抑えるためのものではありません。これは、$\tilde{P}$ が極端に小さい値になったときのゼロ割りや数値的不安定性を避けるための数値計算上の工夫です。また、この値は復元結果そのものを変えることを目的としているわけではなく、理論的な結果から大きくずれない範囲で十分小さな値を入れています。
ノイズがある場合
ここまで見ると、デコンボリューションは単に割り算すればよいように見えます。しかし、実際の観測画像にはノイズがあります。
ノイズを $N$ とすると、観測信号は
O = I * P + N
と書けます。
これをフーリエ変換すると、
\tilde{O}
=
\tilde{I}\tilde{P}
+
\tilde{N}
です。
ここで、ノイズがない場合と同じように $\tilde{P}$ で割ると、
\frac{\tilde{O}}{\tilde{P}}
=
\tilde{I}
+
\frac{\tilde{N}}{\tilde{P}}
となります。
つまり、欲しい成分 $\tilde{I}$ だけでなく、
\frac{\tilde{N}}{\tilde{P}}
という項も出てきます。ここが、単純な逆畳み込みの難しいところです。
PSF は信号をなめらかにする働きを持つため、多くの場合、高周波成分では $\tilde{P}$ が小さくなります。そのため、$\tilde{P}$ が小さいところでは、ノイズ成分 $\tilde{N}$ が大きく増幅されてしまいます。
実際に、先ほどの1次元信号にノイズを足してみます。
# -------------------------
# ぼやけた信号にノイズを足します
# -------------------------
rng = np.random.default_rng(0)
noise_level = 0.001
signal_blur_noisy = signal_blur + noise_level * rng.normal(size=signal_blur.shape)
# -------------------------
# ノイズあり信号をフーリエ空間で割り算して戻します
# -------------------------
# ここでも tiny は非常に小さいままにします。
# そのため、ほぼ単純な割り算になり、
# ノイズが増幅される様子を見ることができます。
# -------------------------
F_blur_noisy = np.fft.fft(signal_blur_noisy)
F_deconv_noisy = F_blur_noisy / (F_psf + tiny)
signal_deconv_noisy = np.real(np.fft.ifft(F_deconv_noisy))
# -------------------------
# 表示します
# -------------------------
plt.figure(figsize=(7, 4))
plt.plot(signal_true, label="True signal")
plt.plot(signal_blur_noisy, label="Blurred + noise")
plt.plot(signal_deconv_noisy, label="Simple deconvolution", color="gray", alpha=0.3, linestyle="--")
plt.ylim(-0.1, 1.1)
plt.xlabel("x")
plt.ylabel("Value")
plt.legend(loc="upper left")
plt.tight_layout()
plt.show()
ノイズを少し足しただけでも、単純な逆畳み込みでは復元信号が大きく乱れることがあります。
式で見ると自然です。
\frac{\tilde{O}}{\tilde{P}}
=
\tilde{I}
+
\frac{\tilde{N}}{\tilde{P}}
となるため、信号だけでなくノイズも $\tilde{P}$ で割られてしまうからです。特に、$\tilde{P}$ が小さい周波数では、ノイズが大きく増幅されます。
つまり、単純なフーリエ逆畳み込みでは、ぼやけを戻しているつもりでも、同時にノイズも強くしてしまいます。
少しだけ安定化してみる
単純な逆畳み込みは
\tilde{I}
=
\frac{\tilde{O}}{\tilde{P}}
でした。
しかし、$\tilde{P}$ が小さい周波数ではノイズが大きく増幅されます。そこで、完全に割り算するのではなく、$\tilde{P}$ が小さい周波数を少し抑えてみます。
たとえば、単純な逆畳み込みに
\frac{|\tilde{P}|^2}
{|\tilde{P}|^2 + \epsilon}
という係数を掛けます。つまり、
\tilde{I}
=
\frac{\tilde{O}}{\tilde{P}}
\frac{|\tilde{P}|^2}
{|\tilde{P}|^2 + \epsilon}
とします。
ここで、$\epsilon$ は小さな正の値です。この係数は、$|\tilde{P}|^2$ が十分大きいところでは
\frac{|\tilde{P}|^2}
{|\tilde{P}|^2 + \epsilon}
\simeq 1
となります。つまり、その周波数では普通の逆畳み込みに近くなります。
一方で、$|\tilde{P}|^2$ が非常に小さいところでは
\frac{|\tilde{P}|^2}
{|\tilde{P}|^2 + \epsilon}
\simeq 0
となります。つまり、ノイズが大きく増幅されそうな周波数成分を抑えることができます。
見方を変えると、PSFによってほとんど失われてしまった周波数成分については、無理に復元しようとしないという考え方です。
ガウシアンPSFの場合、一般に高周波成分ほど強く抑えられるため、$|\tilde{P}|$ が小さくなるのは主に高周波側です。そのため、この安定化は「高周波の復元を少し諦める代わりに、ノイズの増幅を防ぐ」と考えることができます。
ここまでの議論では、
\tilde{I}
=
\frac{\tilde{O}}{\tilde{P}}
\frac{|\tilde{P}|^2}
{|\tilde{P}|^2 + \epsilon}
という形で考えてきました。
一方で、実際にコードへ実装するときには、少し変形した
\tilde{I}
=
\frac{\tilde{O}\tilde{P}^{*}}
{|\tilde{P}|^2 + \epsilon}
という形で書くことがよくあります。ここで、$\tilde{P}^{*}$ は $\tilde{P}$ の複素共役です。こちらの形は、NumPy などで実装するときにそのままコードへ書きやすいため、実際にはこちらを使うことが多いです。
念のため、両者が同じであることを確認します。ここで、
|\tilde{P}|^2
=
\tilde{P}\tilde{P}^{*}
なので、
\frac{\tilde{O}}{\tilde{P}}
\frac{|\tilde{P}|^2}
{|\tilde{P}|^2 + \epsilon}
=
\frac{\tilde{O}}{\tilde{P}}
\frac{\tilde{P}\tilde{P}^{*}}
{|\tilde{P}|^2 + \epsilon}
=
\frac{\tilde{O}\tilde{P}^{*}}
{|\tilde{P}|^2 + \epsilon}
となり、両者が同じであることが分かります。つまり、新しい考え方を導入しているわけではなく、先ほどの式を書き換えただけです。
次に、$\epsilon$ の値をいくつか変えて、安定化の強さを比較してみます。
# -------------------------
# eps_stable を変えて、安定化の強さを比較します
# -------------------------
eps_list = [1e-4, 1e-3, 1e-2]
plt.figure(figsize=(7, 4))
plt.plot(signal_true, label="True signal")
plt.plot(signal_blur_noisy, label="Blurred + noise")
plt.plot(signal_deconv_noisy, label="Simple deconvolution", color="gray", alpha=0.3, linestyle="--")
for eps_stable in eps_list:
F_deconv_stable = (
F_blur_noisy * np.conj(F_psf)
/ (np.abs(F_psf)**2 + eps_stable)
)
signal_deconv_stable = np.real(np.fft.ifft(F_deconv_stable))
plt.plot(
signal_deconv_stable,
label=f"Stabilized, eps={eps_stable:g}",
alpha=0.7,
)
plt.ylim(-0.1, 1.1)
plt.xlabel("x")
plt.ylabel("Value")
plt.legend(loc="upper left")
plt.tight_layout()
plt.show()
この方法では、単純な割り算に比べてノイズの増幅は抑えられます。ただし、これは完全な復元ではありません。
$\epsilon$ が小さいと、単純な逆畳み込みに近くなり、細かい構造は戻りやすくなります。その一方で、ノイズも残りやすくなります。逆に $\epsilon$ を大きくすると、ノイズは抑えられますが、ピークもなまりやすくなります。
つまり、デコンボリューションでは「細かい構造を戻したい」という気持ちと、「ノイズを増やしたくない」という気持ちの間でバランスを取る必要があります。
なお、ここで使った式は、単純な逆畳み込みを安定化したものです。Wiener filter もこの考え方に近いですが、Wiener filter では $\epsilon$ に相当する部分を、信号とノイズのパワースペクトルに基づいて決めます。
その意味で、ここでの方法は Wiener filter そのものというより、Wiener filter につながる簡単な安定化だと見ると分かりやすいです。
2次元画像で考えてみる
ここまで1次元信号で見てきましたが、2次元画像でも考え方は同じです。
真の画像を $I(x,y)$、PSF を $P(x,y)$、観測画像を $O(x,y)$ とします。ノイズを無視すれば、
O(x,y) = I(x,y) * P(x,y)
です。
これを2次元フーリエ変換すると、
\tilde{O}(k_x,k_y)
=
\tilde{I}(k_x,k_y)\tilde{P}(k_x,k_y)
になります。
したがって、PSF が分かっていれば、
\tilde{I}(k_x,k_y)
=
\frac{\tilde{O}(k_x,k_y)}
{\tilde{P}(k_x,k_y)}
として、元の画像を推定できます。
実際に2次元画像でも試してみます。
# -------------------------
# 2次元の真の画像を作ります
# -------------------------
N = 128
image_true = np.zeros((N, N))
image_true[64, 60] = 0.5
image_true[64, 70] = 1.0
# -------------------------
# FFTで扱いやすい2次元ガウシアンPSFを作ります
# -------------------------
# FFTでは (0, 0) が原点として扱われます。
# そのため、PSFの中心を画像中央ではなく (0, 0) に置きます。
# -------------------------
y, x = np.indices((N, N))
dx = np.minimum(x, N - x)
dy = np.minimum(y, N - y)
sigma = 4.0
psf = np.exp(-(dx**2 + dy**2) / (2 * sigma**2))
psf /= psf.sum()
# -------------------------
# FFTで畳み込みを計算し、画像をぼやけさせます
# -------------------------
F_true = np.fft.fft2(image_true)
F_psf = np.fft.fft2(psf)
F_blur = F_true * F_psf
image_blur = np.real(np.fft.ifft2(F_blur))
# -------------------------
# フーリエ空間で割り算して戻します
# -------------------------
tiny = 1e-20
F_deconv = F_blur / (F_psf + tiny)
image_deconv = np.real(np.fft.ifft2(F_deconv))
# -------------------------
# 表示します
# -------------------------
plt.figure(figsize=(10, 3))
plt.subplot(1, 3, 1)
plt.imshow(image_true, origin="lower", cmap="gray")
plt.title("True image")
plt.colorbar(fraction=0.046, pad=0.04)
plt.subplot(1, 3, 2)
plt.imshow(image_blur, origin="lower", cmap="gray")
plt.title("Blurred image")
plt.colorbar(fraction=0.046, pad=0.04)
plt.subplot(1, 3, 3)
plt.imshow(image_deconv, origin="lower", cmap="gray")
plt.title("Deconvolved image")
plt.colorbar(fraction=0.046, pad=0.04)
plt.tight_layout()
plt.show()
2次元でも、やっていることは1次元の場合と同じです。ノイズがない場合には、
O = I * P
をフーリエ変換して、
\tilde{O} = \tilde{I}\tilde{P}
と見ます。そして、$\tilde{P}$ で割ることで、
\tilde{I} = \frac{\tilde{O}}{\tilde{P}}
を計算しています。
このように、ノイズがない理想的な場合には、フーリエ空間で割り算するだけで元の画像に近いものが戻ってきます。
2次元画像でもノイズを入れてみる
次に、2次元画像にノイズを足してみます。
# -------------------------
# ぼやけた画像にノイズを足します
# -------------------------
rng = np.random.default_rng(0)
noise_level = 0.001
image_blur_noisy = image_blur + noise_level * rng.normal(size=image_blur.shape)
# -------------------------
# ノイズあり画像を単純に逆畳み込みします
# -------------------------
F_blur_noisy = np.fft.fft2(image_blur_noisy)
F_deconv_noisy = F_blur_noisy / (F_psf + tiny)
image_deconv_noisy = np.real(np.fft.ifft2(F_deconv_noisy))
# -------------------------
# eps_stable を変えて、安定化の強さを比較します
# -------------------------
eps_list = [1e-5, 1e-4, 1e-3, 1e-2]
image_deconv_stable_list = []
for eps_stable in eps_list:
F_deconv_stable = (
F_blur_noisy * np.conj(F_psf)
/ (np.abs(F_psf)**2 + eps_stable)
)
image_deconv_stable = np.real(np.fft.ifft2(F_deconv_stable))
image_deconv_stable_list.append(image_deconv_stable)
# -------------------------
# 表示します
# -------------------------
plt.figure(figsize=(12, 6))
# -------------------------
# 上段：真の画像から単純な逆畳み込みまで
# -------------------------
plt.subplot(2, 4, 1)
plt.imshow(image_true, origin="lower", cmap="gray")
plt.title("True image")
plt.colorbar(fraction=0.046, pad=0.04)
plt.subplot(2, 4, 2)
plt.imshow(image_blur, origin="lower", cmap="gray")
plt.title("Blurred")
plt.colorbar(fraction=0.046, pad=0.04)
plt.subplot(2, 4, 3)
plt.imshow(image_blur_noisy, origin="lower", cmap="gray")
plt.title("Blurred + noise")
plt.colorbar(fraction=0.046, pad=0.04)
plt.subplot(2, 4, 4)
plt.imshow(image_deconv_noisy, origin="lower", cmap="gray")
plt.title("Simple deconv.")
plt.colorbar(fraction=0.046, pad=0.04)
# -------------------------
# 下段：安定化した逆畳み込みの比較
# -------------------------
for i, (eps_stable, image_deconv_stable) in enumerate(
zip(eps_list, image_deconv_stable_list)
):
plt.subplot(2, 4, 5 + i)
plt.imshow(image_deconv_stable, origin="lower", cmap="gray")
plt.title(f"eps={eps_stable:g}")
plt.colorbar(fraction=0.046, pad=0.04)
plt.tight_layout()
plt.show()
上段では、真の画像、ぼやけた画像、ノイズを加えた画像、そして単純な逆畳み込みの結果を示しています。ノイズを少し足しただけでも、単純な逆畳み込みでは復元画像が大きく乱れることが分かります。
下段では、安定化した式
\tilde{I}
=
\frac{\tilde{O}\tilde{P}^{*}}
{|\tilde{P}|^2 + \epsilon}
を使い、$\epsilon$ の値を変えた結果を並べています。
$\epsilon$ が小さいと単純な逆畳み込みに近くなりますが、ノイズも残りやすくなります。一方で、$\epsilon$ を大きくするとノイズは抑えられますが、画像もなまりやすくなります。
このように、2次元画像でも、1次元信号の場合と同じく、単純な逆畳み込みはノイズに弱く、安定化の強さによって復元結果が変わります。
補足：PSF の中心位置について
この記事のコードでは、PSF の中心を配列の中央ではなく、1次元では index=0、2次元では (0, 0) に置いています。
これは、FFT を使った畳み込みでは、配列の index=0 が原点として扱われるためです。
そのため、1次元では
psf_x = np.minimum(x, N - x)
2次元では
dx = np.minimum(x, N - x)
dy = np.minimum(y, N - y)
として、周期境界を考えたときの原点からの距離を使っています。
このあたりは少し実装寄りの話なので、最初は「FFTで畳み込みをするときは、PSFの中心位置に注意が必要」くらいに思っておけば十分です。
まとめ
この記事では、フーリエ変換を使って画像のぼやけを戻してみました。
PSF が画像全体で同じ形をしていて、ノイズを無視できる場合、観測画像は
O = I * P
と書けます。
畳み込みはフーリエ空間では掛け算になるので、
\tilde{O} = \tilde{I}\tilde{P}
です。
そのため、ノイズがない理想的な場合には、
\tilde{I}
=
\frac{\tilde{O}}{\tilde{P}}
として元の画像を推定できます。
しかし、ノイズがある場合には、
O = I * P + N
なので、
\frac{\tilde{O}}{\tilde{P}}
=
\tilde{I}
+
\frac{\tilde{N}}{\tilde{P}}
となります。つまり、ノイズ成分も $\tilde{P}$ で割られてしまいます。
PSF は画像をなめらかにするため、高周波側では $\tilde{P}$ が小さくなります。そのため、単純なフーリエ逆畳み込みでは、ノイズが大きく増幅されることがあります。
これを少し抑えるために、単純な逆畳み込みに
\frac{|\tilde{P}|^2}
{|\tilde{P}|^2 + \epsilon}
という係数を掛ける方法を試しました。これは
\tilde{I}
=
\frac{\tilde{O}\tilde{P}^{*}}
{|\tilde{P}|^2 + \epsilon}
と書くこともできます。
この方法では、$\tilde{P}$ が小さい周波数成分を抑えることで、ノイズの増幅をある程度防ぐことができます。ただし、これは完全な復元ではありません。
細かい構造を戻そうとするとノイズも増えやすくなり、ノイズを抑えようとすると細かい構造も戻りにくくなります。
単純なフーリエ逆畳み込みは、考え方としてはとても分かりやすいです。しかし、実際の画像解析ではノイズや PSF の扱いが重要になります。
そのため、実データでは Wiener filter や Richardson-Lucy 法のような、より安定したデコンボリューション手法が使われます。
関連記事