はじめに

以前の記事では、PSF 自身に Richardson–Lucy（RL）法を適用し、復元後の PSF の広がり（enclosed 1σ）を見ることで、RL 法の復元精度を簡易的に評価する手法を試しました。

この考え方は Sakai et al. (2024) に基づいており、実際の天文データ解析でも PSF 劣化や復元精度の評価に用いられています。

本記事ではその延長として、PSF を真の分布とした点光源観測を模擬し、そこから得られる ポアソン統計に従う観測画像 に対して RL 法を適用した場合に、enclosed 1σ がどのように振る舞うかを確認したメモをまとめます。

厳密な理論検証を目的としたものではなく、挙動確認を目的とした簡単な検証という位置付けです。

手法の概要

行ったことは以下の通りです。

ガウシアン PSF を作成
総カウント数 N を指定してポアソンノイズを付加（観測を模擬）
ノイズなしの PSF をカーネルとして用い、ポアソンノイズを加えた観測画像に RL 法を適用
各反復ステップで enclosed 1σ（68% energy radius）を計算
この手順を複数回（30回）繰り返し、平均と不偏標準偏差を評価
総カウント数 N を変えて結果を比較

※本検証では分解能評価を簡単化するため、中心位置は固定しています。

ここで用いる enclosed 1σ は、ガウシアン分布の分散としての σ ではなく、PSF の総強度の 68% が内包される半径を指します。PSF 形状がガウシアンから外れたり、RL 法によって歪んだ場合でも一貫して評価できる量として、この定義を採用しています。

実装コード

以下が今回の検証に用いたコード一式です。ノイズなし（理想ケース）と、複数の総カウント数 N に対する結果を重ねて表示します。

実行コードは Google Colab 上
https://colab.research.google.com/drive/1y_mYCD5MQnbGCnc6gTGX5k4UFeg9JqDs?usp=sharing
からもお試しいただけます。

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

# -----------------------------
# PSF 作成
# -----------------------------
def create_gaussian_psf(size=201, sigma=0.5):
    x = np.linspace(-1, 1, size)
    y = np.linspace(-1, 1, size)
    xx, yy = np.meshgrid(x, y)
    r2 = xx**2 + yy**2
    psf = np.exp(-r2 / (2 * sigma**2))
    psf /= psf.sum()
    return psf

# -----------------------------
# Richardson–Lucy 法
# -----------------------------
def richardson_lucy_all_iter(image, psf, iterations=50):
    image = image.astype(float)
    psf = psf.astype(float)

    psf_mirror = np.flip(psf, (0, 1))
    im_deconv = np.full(image.shape, 1.0)
    im_deconv /= im_deconv.sum()

    ims = np.zeros((iterations + 1, *image.shape))
    ims[0] = im_deconv.copy()

    for i in range(1, iterations + 1):
        conv = convolve(im_deconv, psf, mode="same") + 1e-14
        im_deconv *= convolve(image / conv, psf_mirror, mode="same")
        ims[i] = im_deconv.copy()

    return ims

# -----------------------------
# enclosed energy 半径（68%）
# -----------------------------
def calculate_sigma(psf, threshold=0.68):
    size = psf.shape[0]
    center = size // 2
    y, x = np.indices(psf.shape)
    r = np.sqrt((x - center)**2 + (y - center)**2)

    order = np.argsort(r.flat)
    cumsum = np.cumsum(psf.flat[order])
    idx = np.searchsorted(cumsum, threshold)

    return r.flat[order][idx]

# -----------------------------
# ポアソンノイズ（観測模擬）
# -----------------------------
def add_poisson_noise(psf, total_counts):
    counts = np.random.poisson(psf * total_counts)
    return counts / counts.sum()

# -----------------------------
# パラメータ
# -----------------------------
psf_size   = 201
sigma_psf  = 0.5
iterations = 100
N_list     = [30, 50, 100]
n_real     = 30
np.random.seed(100)

# -----------------------------
# 真の PSF
# -----------------------------
psf = create_gaussian_psf(psf_size, sigma=sigma_psf)

# -----------------------------
# ノイズなし（理想ケース）
# -----------------------------
rl_ideal = richardson_lucy_all_iter(psf, psf, iterations)
sigmas_ideal = np.array([calculate_sigma(im) for im in rl_ideal])

# -----------------------------
# プロット準備
# -----------------------------
plt.figure(figsize=(8, 5))
x = np.arange(iterations + 1)

plt.plot(
    x,
    sigmas_ideal,
    color="black",
    lw=2.5,
    label="No noise (ideal)"
)

# -----------------------------
# Poisson ノイズあり（統計評価）
# -----------------------------
for i, N in enumerate(N_list):
    print(f"Poisson noise level: {i+1}/{len(N_list)} (N={N})")
    sigmas_all = []

    for _ in tqdm(range(n_real), desc=f"Realizations (N={N})", leave=False):
        psf_noisy = add_poisson_noise(psf, N)
        rl_poisson = richardson_lucy_all_iter(psf_noisy, psf, iterations)
        sigmas = np.array([calculate_sigma(im) for im in rl_poisson])
        sigmas_all.append(sigmas)

    sigmas_all = np.array(sigmas_all)

    mean = sigmas_all.mean(axis=0)
    std  = sigmas_all.std(axis=0, ddof=1)

    plt.plot(x, mean, lw=1.5, label=f"Poisson noise (N={N})")
    plt.fill_between(x, mean - std, mean + std, alpha=0.2)

# -----------------------------
# 表示
# -----------------------------
plt.xlabel("Iteration")
plt.ylabel("Enclosed 1σ (pixels)")
plt.title("Effective resolution limited by Poisson statistics")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

図では、黒線がノイズなしの理想ケース、色線が総カウント数 N=30, 50, 100 のポアソンノイズを含む場合の挙動を示しています。あみがけは ±不偏標準偏差で、統計が悪いほど分解能のばらつきが大きくなることが確認できます。

結果と簡単な考察

反復回数に対する enclosed 1σ の挙動を見ると、ノイズの有無によらず、反復とともに単調に鋭くなっていくことがわかります。一方で、統計が悪い場合ほど分解能のばらつきが大きくなることも確認できます。これは、ポアソン統計に由来する揺らぎによって、復元結果に試行ごとの差が生じるためだと考えられます。

興味深い点として、平均的にはノイズを入れたケースにおいて、ノイズなしで得られる分解能を下回る挙動は見られませんでした。このことから、ノイズなしで得られる enclosed 1σ は、RL 法における実効的な分解能の参考値として捉えられる可能性があります。実データ解析においても、反復回数や復元結果を考える際の一つの目安として使えるかもしれません。

まとめ

以前の記事で扱った PSF 自己評価の考え方をベースに、ポアソンノイズを含む観測を模擬した場合の挙動を簡単に確認しました。定量的な結論を与えるものではありませんが、RL 法を実データに適用する際に「統計がどこで効いてくるか」を把握するためのメモとしてまとめています。

反復回数の選定や、復元結果の解釈を考える際の参考になれば幸いです。

参考文献

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