PSFが自らの不定性を推定！Richardson-Lucy法の実践的な活用法

はじめに

デコンボリューション法は、ぼやけた画像を復元するための基本的な手法であり、天文学では観測装置による画像の劣化を補正するために広く利用されています。

本記事では、PSF（点拡がり関数）を用いたRL法の復元精度評価について解説します。この手法は、PSF自身を観測画像としてRL法を適用することで、不定性や復元精度を簡易的に見積もるものです。また、反復回数が復元結果に与える影響についても検討します。

なお、Sakai et al. (2024) では、この手法を天文学的観測データに適用し、光軸外でのPSF劣化や復元精度を評価しています。本記事では、その手法を簡潔に解説し、応用例も紹介します。

手法の概要

1. 基本的な評価手法

PSFは点光源の広がりを表す関数であり、これをRL法で復元した結果を比較することで、復元効果を評価します。手順は以下の通りです：

1. ガウシアンPSF画像を作成
   点光源を表現するPSFをシミュレーションで作成します。

2. PSF自身を観測画像としてRL法を適用
   作成したPSFをRL法に入力し、復元の変化を追跡します。

3. PSFの広がり（1σ）を比較
   RL法を適用してもPSFは完全な点には戻りません。特定の広がり（1σ）で収束する現象を評価することで、復元可能な範囲や不定性を定量化できます。これが本手法の鍵です！

この方法により、復元後の画像がどの程度元の状態に近づいたかを簡易的に評価できます（図1）。

注意：この手法はノイズや画像構造の影響を受けるため、実際の観測画像では結果が異なることがあります。評価結果を解釈する際には十分に注意してください。

2. 反復回数の重要性

RL法の反復回数は、復元精度に大きく影響を与えます。図2に示すように、反復回数に応じてPSFの広がり（σ）が変化し、最適な回数で収束します。このため、復元画像と同じ反復回数をPSF評価でも用いることが重要です。

実装と結果

以下に、PSFの作成、RL法の適用、および評価の結果を示します。
実装コードはGoogle Colabこちらと本記事の下部のコード全容に載せています。

実装コード（抜粋）

図1を作成するためのコードの抜粋

# ガウシアン PSF作成
psf = create_gaussian_psf(size=psf_size)  # 点拡がり関数の作成

# PSFにRL法を適用し、反復回数分の復元結果を取得
rl_all_with_psf = results = richardson_lucy_all_iter(psf, psf, iterations=iterations)

# 元のPSFと復元後のPSFの広がりを計算
sigma_psf = calculate_sigma(psf)  # 元のPSFの広がり（1σ）
sigma_rl_with_psf = calculate_sigma(rl_all_with_psf[iterations])  # RL法後の広がり

結果のグラフ

以下に、RL法適用前後のPSFの比較と反復回数に応じたσの変化を示します。

図1：PSFと復元後の比較

図1では、左側が元のPSF、右側がRL法後のPSFです。復元後のPSFが中心に集約され、劣化が補正されていることが確認できます。

図2：反復回数に応じたPSFのσの変化

図2は、横軸が反復回数、縦軸が復元されたPSFの広がり（σ）を示しています。反復回数が増えるにつれてσが収束することが確認できます。この結果は、観測画像に適用した反復回数と同じ段階で得られたPSFのσを評価に使用することが重要であることを示しています。

観測画像はノイズや構造の複雑性の影響を受けるため、シンプルなPSFと完全に一致するわけではありませんが、同じ反復段階の結果を比較することで、復元精度の評価において一貫性のある判断が可能になると考えられます。

応用例

本手法は、観測装置の性能評価や補正に応用できます。

応用例1: 適切な反復回数の選定

図2で示したように、RL法の反復回数が復元結果に大きく影響を与えます。こちらの記事では、観測画像のカイ二乗値を基に最適な反復回数を推定する方法が提案されています。この手法と本記事の評価手法を組み合わせることで、適切な反復回数に応じたRL後の不定性を評価することが可能になります。

応用例2: スムージング処理

復元画像にPSFの広がり（σ）を基準としたガウス関数でスムージングを適用することで、系統的な誤差を軽減できます。特に、光軸外のPSFが劣化している場合に有効です。

まとめ

本記事では、PSFを用いてRL法の復元画像の品質を評価する簡易的な手法を紹介しました。この方法は、特に光軸外でのPSF劣化におけるRL法後の復元精度の評価に有効であり、反復回数を考慮することで復元結果の一貫性と信頼性を向上させることが可能です。

また、Sakai et al. (2024) の研究では、この手法を天文学的観測データに適用し、実践的な利用方法を示しています。この研究は、RL法とPSF評価を組み合わせた解析が、観測装置の性能評価やデータ補正において有効であることを示唆しています。

今後、この手法はより複雑な観測データへの応用や、高度な統計手法との組み合わせによって、さらに幅広い分野での活用が期待されます。本記事で紹介した評価手法が、実践的な解析の第一歩として役立つことを願っています。

参考文献

• Sakai, Y. et al. 2024, ApJ 974 245, doi: 10.3847/1538-4357/ad739f

コード全容

コードはこちら

図1:PSFと復元後の比較

import numpy as np
from scipy.signal import convolve
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm


def create_gaussian_psf(size=201, sigma=0.2):
    x = np.linspace(-1, 1, size)
    y = np.linspace(-1, 1, size)
    xx, yy = np.meshgrid(x, y)
    r_squared = xx**2 + yy**2
    psf = np.exp(-r_squared / (2 * sigma**2))
    psf /= psf.sum()  # PSFの規格化
    return psf

def richardson_lucy_all_iter(image, psf, iterations=50, normalize_initial=True):
    """
    Perform Richardson-Lucy deconvolution with an option to normalize the initial guess.

    Parameters:
        image (ndarray): Observed (blurred) image.
        psf (ndarray): Point Spread Function (PSF).
        iterations (int): Number of iterations.
        normalize_initial (bool): Whether to normalize the initial guess.

    Returns:
        im_deconvs (ndarray): Array of deconvolved images at each iteration.
    """
    # Promote data type to float32 or higher for computation
    float_type = np.promote_types(image.dtype, np.float32)
    image = image.astype(float_type, copy=False)
    psf = psf.astype(float_type, copy=False)

    # Mirror the PSF for convolution
    psf_mirror = np.flip(psf, (0, 1))

    # Initialize
    im_deconv = np.full(image.shape, 0.5, dtype=float_type)  # Initial estimate

    # Normalize the initial guess if requested
    if normalize_initial:
        im_deconv /= im_deconv.sum()

    # Store all iterations
    im_deconvs = np.empty((iterations + 1, *image.shape), dtype=float_type)
    im_deconvs[0] = im_deconv.copy()

    # Iterative Richardson-Lucy deconvolution
    for i in range(1, iterations + 1):
        conv = convolve(im_deconv, psf, mode="same") + 1e-14  # Avoid division by zero
        relative_blur = image / conv
        im_deconv *= convolve(relative_blur, psf_mirror, mode="same")

        # Save current iteration result
        im_deconvs[i] = im_deconv.copy()

    return im_deconvs

# シグマ計算関数
def calculate_sigma(psf, threshold=0.68):
    size = psf.shape[0]
    center = size // 2
    x, y = np.meshgrid(range(size), range(size))
    distances = np.sqrt((x - center)**2 + (y - center)**2)
    flat_psf = psf.flatten()
    flat_distances = distances.flatten()
    sorted_indices = np.argsort(flat_distances)
    sorted_psf = flat_psf[sorted_indices]
    sorted_distances = flat_distances[sorted_indices]

    cumsum = np.cumsum(sorted_psf)
    radius_index = np.searchsorted(cumsum, threshold)
    sigma = sorted_distances[radius_index]
    return sigma

# パラメータ
psf_size = 201
iterations = 100

# ガウシアン PSF作成
psf = create_gaussian_psf(size=psf_size)

# PSFにRL法を適用
rl_all_with_psf = results = richardson_lucy_all_iter(psf, psf, iterations=iterations)

# シグマ計算
sigma_psf = calculate_sigma(psf)
sigma_rl_with_psf = calculate_sigma(rl_all_with_psf[iterations])

# プロット
fig, axs = plt.subplots(1, 2, figsize=(13, 5))

# PSF
extent = (-psf_size / 2, psf_size / 2, -psf_size / 2, psf_size / 2)
im0 = axs[0].imshow(psf+1e-12, cmap="gray", norm=LogNorm(vmin=1e-5, vmax=1e-1), extent=extent)
circle_psf = plt.Circle((0, 0), sigma_psf, color="lime", fill=False, linestyle="--", linewidth=1.5, label="1σ")
axs[0].add_artist(circle_psf)
axs[0].set_title(f"PSF \n1σ: {sigma_psf:.2f} pixels")
axs[0].legend()
fig.colorbar(im0, ax=axs[0], orientation='vertical', label="Intensity")

# RL後のPSF
im1 = axs[1].imshow(rl_all_with_psf[iterations]+1e-12, cmap="gray", norm=LogNorm(vmin=1e-5, vmax=1e-1), extent=extent)
circle_rl = plt.Circle((0, 0), sigma_rl_with_psf, color="lime", fill=False, linestyle="--", linewidth=1.5, label="1σ")
axs[1].add_artist(circle_rl)
axs[1].set_title(f"RL with PSF (Iterations: {iterations})\n1σ: {sigma_rl_with_psf:.2f} pixels")
axs[1].legend()
fig.colorbar(im1, ax=axs[1], orientation='vertical', label="Intensity")

plt.tight_layout()
plt.show()

図2:反復回数に応じたPSFのσの変化

# シグマ値を保存する配列を初期化
sigmas = np.zeros(iterations+1)

# 各反復結果からシグマ値を計算して保存
for i in range(iterations+1):
    sigmas[i] = calculate_sigma(rl_all_with_psf[i])  # シグマ値を配列に格納

# プロット
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(np.arange(0, iterations + 1), sigmas, "o-", label="1σ")
ax.set_xlabel("Iteration")
ax.set_ylabel("Sigma (pixels)")
ax.set_title("Sigma Values Over Iterations")
ax.set_ylim(0, None)
ax.grid(True, linestyle="--", alpha=0.7)
ax.legend()

plt.show()