脳内デコンボリューションで解き明かす！PSFの影響か、真の構造の広がりか？

## はじめに

天文学のデータ解析では、観測される天体の構造は PSF（Point Spread Function）によってぼやけます。したがって、観測データの広がりが PSF の影響によるものなのか、それとも天体の本来の構造によるものなのかを見極めることは、解析において重要な課題となります。  

例えば、X線天文学では、PSF をシミュレーションや実測データをもとに見積もり、観測データの真の空間的スケールを推定します。その際、イメージデコンボリューション法を用いて直接真の天体像を推定することもありますが、ノイズの影響を受けやすいという課題があります。そのため、PSF の影響を直感的に理解し、**頭の中で逆畳み込みをイメージできる** ことが、より適切な解析手法の選択につながる場合があります。  

本記事では、PSFと真の構造をガウシアン分布とみなし、観測される広がりの変化を数学的に整理し、可視化を通じて直感的に理解することを目指します。

本記事の実装コードは Google Colab からも確認できます。

https://colab.research.google.com/drive/1aZlELJpp9IcXsQyMFIKpdRQGPHz7Mmil?usp=sharing

## 観測される広がりの数式

天文学や画像処理において、観測されるデータは真の空間構造とPSFの畳み込みとして得られます。まず、畳み込みの数式を導入し、その後、PSFと真の構造をガウシアン分布と仮定した場合の広がりを考えます。

### 畳み込みとは
観測される画像 $I(x)$ は、真の構造 $f(x)$ とPSF $p(x)$ の畳み込みで表されます。

```math
I(x) = (f * p)(x) = \int_{-\infty}^{\infty} f(t) p(x - t) dt
```

この畳み込みにより、PSFの影響で真の構造がぼやけます。

### PSFと真の構造がガウシアンの場合
PSFと真の構造をガウシアン分布と仮定すると、それぞれ次のように表されます。

- 真の構造
  ```math
  f(x) = \frac{1}{\sqrt{2\pi} \sigma_{\text{true}}} \exp \left(-\frac{x^2}{2\sigma_{\text{true}}^2} \right)
  ```

- PSF
  ```math
  p(x) = \frac{1}{\sqrt{2\pi} \sigma_{\text{psf}}} \exp \left(-\frac{x^2}{2\sigma_{\text{psf}}^2} \right)
  ```
  
ガウシアン同士の畳み込みもガウシアンになり、観測される広がりの標準偏差は
```math
\sigma_{\text{obs}} = \sqrt{\sigma_{\text{true}}^2 + \sigma_{\text{psf}}^2}
```
で表されます。

## 数式の可視化

$\sigma_{\text{obs}} = \sqrt{\sigma_{\text{true}}^2 + \sigma_{\text{psf}}^2}$の関係を可視化してみたいと思います。横軸に真の構造の標準偏差、縦軸に観測の標準偏差をとります。

<details><summary>数式可視化用のPythonコード</summary>


```python:PSFによる空間分解能の影響
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.size'] = 12

# パラメータ設定
sigma_true_array = np.linspace(0.01, 5, 100)  # x は真の構造のスケール
sigma_psf = 1  # PSFの空間分解能
sigma_obs_array = np.sqrt(sigma_true_array**2 + sigma_psf**2)

# プロット
plt.figure(figsize=(7, 7))
plt.plot(sigma_true_array, sigma_obs_array, label=fr"With PSF: $\sigma_{{\text{{obs}}}} = \sqrt{{\sigma_{{\text{{true}}}}^2 + \sigma_{{\text{{psf}}}}^2}}$, ($\sigma_{{\text{{psf}}}}={sigma_psf}$)", color='k')
plt.plot(sigma_true_array, sigma_true_array, label=fr"Without PSF: $\sigma_{{\text{{obs}}}}=\sigma_{{\text{{true}}}}$", color='red', linestyle="--")

plt.title("Observed Standard Deviation with PSF")
plt.xlabel(r"$\sigma_{\text{true}}$")
plt.ylabel(r"$\sigma_{\text{obs}}$")
plt.legend()
plt.grid(True)
plt.gca().set_aspect('equal')

plt.show()
```

```python:PSFによる中心部の強度の影響
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.size'] = 12

# パラメータ設定
sigma_true_array = np.linspace(0.01, 5, 100)  # x は真の構造のスケール
sigma_psf = 1  # PSFの空間分解能
sigma_obs_array = np.sqrt(sigma_true_array**2 + sigma_psf**2)

# 中心強度の計算
I_true_values = 1 / (np.sqrt(2 * np.pi) * sigma_true_array)
I_obs_values = 1 / (np.sqrt(2 * np.pi) * sigma_obs_array)

# プロット
plt.figure(figsize=(7, 7))
plt.plot(sigma_true_array, I_obs_values, label=fr"With PSF: $\sigma_{{\text{{obs}}}} = \sqrt{{\sigma_{{\text{{true}}}}^2 + \sigma_{{\text{{psf}}}}^2}}$, ($\sigma_{{\text{{psf}}}}={sigma_psf}$)", color='k')
plt.plot(sigma_true_array, I_true_values, label=fr"Without PSF: $\sigma_{{\text{{obs}}}}=\sigma_{{\text{{true}}}}$", color='red', linestyle="--")

plt.title("Observed Peak Intensity with PSF")
plt.xlabel(r"$\sigma_{\text{true}}$")
plt.ylabel(r"$I_{\text{obs}}(0)$")
plt.legend()
plt.grid(True)
plt.yscale('log')

plt.show()
```

</details>


**出力結果**
- <font color="red">**PSFがある場合**</font>（下図は、$\sigma_{\text{psf}}=1$に固定して表示したもの）
- **PSFがない場合**

| 観測される広がりの変化 | 中心部の強度 |
|:-:|:-:|
| ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/bce4baf8-fd58-44a2-9efd-d70d88631b8d.png) | ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/0ce4a3e1-ac7f-4cae-aa6c-5badaff2af7f.png)|


この可視化から、以下のことが分かります。
- 真の構造がPSFより大きいと、PSFの影響は小さい
- 真の構造がPSFより小さいと、観測される広がりが大きくなる

## 実際のガウス分布で確認
ここでは、実際のガウス分布を用いて、PSFの影響を直感的に捉える練習をします。

<details><summary>PSFと真の空間構造のガウス分布の確認コード</summary>

```python:実際のガウス分布で確認
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.size'] = 12

# 1次元ガウシアン関数
def gaussian_1d(x, sigma):
    return np.exp(-x**2 / (2 * sigma**2)) / (np.sqrt(2 * np.pi) * sigma)

# パラメータ設定
sigma_true = 3  # 真の構造の広がり
sigma_psf = 1  # PSFの広がり
sigma_obs = np.sqrt(sigma_true**2 + sigma_psf**2)  # 観測される広がり

# x軸の範囲
x_values = np.linspace(-20, 20, 500)

# 各ガウシアンを計算
y_true = gaussian_1d(x_values, sigma_true)
y_true /= np.sum(y_true)
y_psf = gaussian_1d(x_values, sigma_psf)
y_psf /= np.sum(y_psf)

# 畳み込み演算
y_obs = np.convolve(y_true, y_psf, mode='same')

# プロット
plt.figure(figsize=(1.618*5, 5))
plt.plot(x_values, y_true, label=fr"$\sigma_{{\text{{true}}}} = {sigma_true}$", color='k')
plt.plot(x_values, y_psf, label=fr"$\sigma_{{\text{{psf}}}} = {sigma_psf}$", color='gray', linestyle=":")
plt.plot(x_values, y_obs, label=fr"$\sigma_{{\text{{obs}}}} = \sqrt{{\sigma_{{\text{{true}}}}^2 + \sigma_{{\text{{psf}}}}^2}}={sigma_obs:.2f}$", color='red', linestyle="--")

plt.xlabel("x")
plt.ylabel("Intensity")
plt.title("1D Gaussian Convolution: Effect of PSF")
plt.legend()
plt.grid(True)

plt.show()
```

</details>


- 真の空間スケール **<<** PSFの空間分解能: $\sigma_{\text{psf}}=1, \sigma_{\text{true}}=0.5$
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/1f256884-9ab3-406d-9bd7-81f04f56c5db.png)
- 真の空間スケール **=** PSFの空間分解能: $\sigma_{\text{psf}}=1, \sigma_{\text{true}}=1$
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/8adec0c0-ee9e-480f-855f-3161b42c7ff1.png)
- 真の空間スケール **>>** PSFの空間分解能: $\sigma_{\text{psf}}=1, \sigma_{\text{true}}=2$
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/f850be08-e2e9-4bf0-a95f-302aaaad6a34.png)

パラメータを変えることで、PSFの影響を直感的に理解できます。実際に試してみると、ルートの項が効いてくるため、PSFと真の構造の広がりの関係が、直感的な予想と少し異なると感じる読者もいるかもしれません。

## GIFアニメーションで確認
PSFの影響がどのように変化するかをアニメーションでも確認します。

<details><summary>GIFアニメーション作成コード</summary>

```python:GIFアニメーションで確認
import imageio
import os

sigma_psf = 1

# sigma_true の範囲
sigma_true_values = np.arange(0.2, 5.2, 0.2)  # 5 から 20 まで変化

# x軸の範囲
x_values = np.linspace(-20, 20, 500)

# 画像ファイルリスト
image_files = []

for sigma_true in sigma_true_values:
    # 観測される広がりを計算
    sigma_obs = np.sqrt(sigma_true**2 + sigma_psf**2)

    # 各ガウシアンを計算
    y_true = gaussian_1d(x_values, sigma_true)
    y_true /= np.sum(y_true)
    y_psf = gaussian_1d(x_values, sigma_psf)
    y_psf /= np.sum(y_psf)

    # 畳み込み演算
    y_obs = np.convolve(y_true, y_psf, mode='same')

    # プロット
    plt.figure(figsize=(1.618 * 5, 5))
    plt.plot(x_values, y_true, label=fr"$\sigma_{{\mathrm{{true}}}} = {sigma_true:.1f}$", color='k')
    plt.plot(x_values, y_psf, label=fr"$\sigma_{{\mathrm{{psf}}}} = {sigma_psf}$", color='gray', linestyle=":")
    plt.plot(x_values, y_obs, label=fr"$\sigma_{{\mathrm{{obs}}}} = \sqrt{{\sigma_{{\mathrm{{true}}}}^2 + \sigma_{{\mathrm{{psf}}}}^2}}={sigma_obs:.2f}$", color='red', linestyle="--")

    plt.xlabel("x")
    plt.ylabel("Intensity")
    plt.title("1D Gaussian Convolution: Effect of PSF")
    plt.legend()
    plt.grid(True)

    # 画像保存
    image_path = os.path.join("./", f"frame_{sigma_true:.2f}.png")
    plt.savefig(image_path)
    image_files.append(image_path)
    plt.close()

# GIF アニメーションを作成
imageio.mimsave("./gaussian_convolution.gif", [imageio.imread(img) for img in image_files], fps=3)
```

</details>

![gaussian_convolution.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/81245e37-9d94-4503-afec-3b8027bdf543.gif)

このGIFから、PSFによる広がりの変化を視覚的に確認できます。小さい構造ほどぼやけやすく、大きな構造では影響が小さいことが分かります。

## まとめ
本記事では、PSF（Point Spread Function）が観測データに与える影響について、ガウシアン分布を用いて説明し、可視化を通じて直感的に理解できるようにしました。  

結果として、**小さい構造ほどPSFの影響を強く受け、十分に大きな構造では影響が小さくなる** ことが確認できました。PSFの影響を正しく理解することで、観測データの解釈や分析の精度を高めることができます。  

PSFの影響を意識することで、より正確なデータ解析が可能になります。本記事が、その理解の一助となれば幸いです。