Richardson-Lucy deconvolutionの特徴をシミュレーションで理解しよう（前編：シミュレーションの実装）

# はじめに
Richardson-Lucy deconvolution（以下、RL法）は、画像等を鮮明にする手法である。一般的な撮像観測では、観測装置の回折等により光源が拡がって撮像される。この拡がり具合は、点拡がり関数(Point Spread Function, 以下、PSF)で表される。RL法は既知のPSFと撮像画像から、真の画像を推定する手法である。
![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/1ce4c96a-b902-8d29-7d32-62842014d6bf.png)

# 本シリーズの目標
実際の撮像観測では真の画像は知り得ないので、RL法により推定された真の画像は、本当に真の画像に近づいているか分からない。そこで、本シリーズでは、**真の画像とPSFが完全にわかっていて、かつノイズが乗らない理想的な状況でのシミュレーションにより、RL法の特徴を理解する**。
本シリーズは3部構成（前編：シミュレーションの実装、中編：結果の解析、後編：補足）で、気長に記事にしていく予定である。（お楽しみに！）

# RL法とは
RL法は、既知のPSFと撮像画像を用いて真の画像を推定する手法である。真の画像の分布、既知のPSF、撮像画像の関係をベイズの定理を用いて関係づける。そして、ベイズ推定を繰り返し計算することにより真の画像を推定する。理論式は

```math
 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
```
と表される。$j,k$はピクセル座標とする。$W$は真の画像、$H$は撮像画像、$r$は反復回数とする。$P_{j,k}$はPSFで、$W_j$が$H_k$で観測される確率を$P(H_k|W_j)=P_{j,k}$とする。式の導出は、私の記事で恐縮ですが、[Richardson-Lucy deconvolutionの理論式を導出しよう](https://qiita.com/yusuke_s_yusuke/items/1f442023d9d1326da2a1)を参照されたい。


# シミュレーションの実装
Google Colabで作成した本記事のコードは、[こちら](https://colab.research.google.com/drive/19MpY6Ar3lOIP8CpQnVR_n8CckpW4ydLd?usp=sharing)にあります。


## 各種インポート
```Python:各種インポート
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.cm as cm
from scipy.stats import norm  # おまけの項で使います
from scipy.stats import multivariate_normal
from scipy.signal import convolve
from skimage import restoration
```
実行時の各種バージョン

| Name | Version |
|:-:|:-:|
| numpy | 1.21.5 |
| matplotlib | 3.2.2 |
| scipy | 1.4.1 |
| scikit-image | 0.18.3 |

※[scikit-image](https://scikit-image.org/)の[``restoration.richardson_lucy``](https://scikit-image.org/docs/stable/api/skimage.restoration.html#richardson-lucy)は、Version 0.19.1で一部変更されている。（変更箇所概要：[Add small regularization to avoid zero-division in richardson_lucy (gh-5976)](https://scikit-image.org/docs/dev/release_notes.html#bugfixes))

## PSF作成用の関数
```Python: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](https://scipy.org/)の[``multivariate_normal``](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.multivariate_normal.html#scipy.stats.multivariate_normal)を使用する。
多変量正規分布の性質については、[多変量正規分布の性質を図で解説する](https://qiita.com/yutera12/items/383847ca9f474cb05ed8)や[多変量正規分布をPythonでplotして理解する](https://qiita.com/g-k/items/698c7f9e4a213d73197b)を参照されたい。

## RL法のシミュレーション
```Python:RL法のシミュレーション
"""シミュレーション"""
# 真の画像
true_px = 128
true_img = np.zeros([true_px, true_px])
true_img[64, 52] = 0.5
true_img[64, 76] = 0.5

# PSF
psf_px = 127
psf_sigma = 17
psf = get_psf(psf_px, psf_sigma)

# 真の画像をPSFで畳み込む（撮像画像）
obtained_img = convolve(true_img, psf, "same")

# RL法
iterations = 50
rl_img = restoration.richardson_lucy(obtained_img, psf, iterations)

"""シミュレーションの結果を図に表示"""
fig, (ax1, ax2, ax3, ax4) = plt.subplots(figsize=(18, 7), ncols=4)
cbar_min = 1e-5
cbar_max = 1e-2

# 真の画像
ax1.set_title("The true image ({0}×{0} pixel)\n".format(true_px))
cbar1 = ax1.imshow(true_img, cmap=cm.CMRmap, vmin=0, vmax=1)
ax1.axis("off")
fig.colorbar(cbar1, ax=ax1, orientation="horizontal")

# PSF
ax2.set_title("PSF σ={0} ({1}×{1} pixel)\n".format(psf_sigma, psf_px))
psf_display = np.where(psf < cbar_min, cbar_min, psf)
ax2.imshow(psf_display, cmap=cm.CMRmap,
           norm=colors.LogNorm(vmin=cbar_min, vmax=cbar_max))
ax2.axis("off")

# 撮像画像
ax3.set_title("The obtained image ({0}×{0} pixel)\n".format(true_px))
obtained_img_display = np.where(obtained_img < cbar_min, cbar_min, obtained_img)
ax3.imshow(obtained_img_display, cmap=cm.CMRmap,
           norm=colors.LogNorm(vmin=cbar_min, vmax=cbar_max))
ax3.axis("off")

# RL法による生成画像
ax4.set_title("Richardson-Lucy deconvolution\n image iterations={0}".format(iterations))
rl_img_display = np.where(rl_img < cbar_min, cbar_min, rl_img)
cbar2 = ax4.imshow(rl_img_display, cmap=cm.CMRmap,
                   norm=colors.LogNorm(vmin=cbar_min, vmax=cbar_max))
ax4.axis("off")
fig.colorbar(cbar2, ax=[ax2, ax3, ax4], orientation="horizontal", aspect=70)

plt.show()
```
* 出力結果
![Unknown-118.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/44e5dd72-ca69-34b9-3706-c164b6dff9eb.png)

* コードの概念図
![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/568d500a-dfb0-8896-eeaf-fdced2cd740a.png)

* コードの説明
    1. 真の画像とPSFを用意する。上の例では、真の画像はピクセル座標``[64, 52], [64, 76]``に``0.5``の値を入れたものを用いる。PSFは多変量正規分布の分散共分散行列の標準偏差σを、``psf_sigma=17``として用いる。<br>※PSFのピクセルサイズ``psf_px``は、奇数が望ましい。奇数だとPSFの中心のピクセル座標が一意に決まるので、真の画像をPSFで畳み込むときやRL法が上手くいく。（偶数のピクセルサイズでも畳み込み自体は実行できるが、意図せぬ（ある法則に従った）ピクセル座標を中心とした畳み込みが行われるので、注意したい。）
    1. 真の画像をPSFで畳み込み、撮像画像を作成する。
    1. 撮像画像と2と同じPSFで、RL法を行う。<br>※``***_display = np.where(*** < cbar_min, cbar_min, ***)``は、対数のカラーバーでは``0``や非常に小さい値は上手く表示できないため、カラーバーの最小値より小さい値をカラーバーの最小値に置換している。

## 反復回数を変化させる
* 反復回数5, 10, 15, 20回の結果
![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/3c364e4a-a80b-f3b1-c89b-020709e9c452.png)
* 反復回数25, 50, 75, 100回の結果
![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/3d494bdd-8eec-3eca-a6af-cc6e8375bca7.png)


このように、反復回数が大きくなるにつれて、真の画像に近づくことが分かる。また、反復回数が5~20回では、25~100回に比べて急激に画像が変化していることが分かる。


## PSFを変化させる
* ``psf_sigma=21, 17, 13, 9``で、反復回数5, 10, 15, 20回の結果
![図2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/32965845-beca-c39f-f80e-f78241348c1d.png)


* ``psf_sigma=21, 17, 13, 9``で、反復回数25, 50, 75, 100回の結果
![図3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/627d1b0a-b73a-58e6-5d9f-3d5682ab3b44.png)

このように、PSFが小さいほど、同じ反復回数で比べたときに真の画像に近いことが分かる。

## 色々なサンプルでシミュレーション
RL法の特徴を理解するために、偏りのある点源の場合、点源が適当に散らばった場合、拡がった構造の場合、偏りのある拡がった構造の場合でシミュレーションする。

### 偏りのある点源の場合
```Python:偏りのある点源の場合（真の画像のコード部分）
"""シミュレーション"""
# 真の画像
true_px = 128
true_img = np.zeros([true_px, true_px])
true_img[64, 52] = 0.2
true_img[64, 76] = 0.8
"""以下省略"""
```
[RL法のシミュレーション](#rl法のシミュレーション)のコードの、真の画像を変更している。

* 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/aaad1214-2148-bf5a-e1c6-36d3be4d7b55.png)
* ``psf_sigma=21, 17, 13, 9``で、反復回数5, 10, 15, 20回の結果
![図4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/0db3d07b-3cd1-c801-f71f-0a6e7b133be5.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数25, 50, 75, 100回の結果
![図5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/2224a839-37ef-26ea-0f6e-255b90654b1a.png)

このように、大きい値（明るい）の点源によって小さい値（暗い）の点源が見えにくくなることが分かる。ただ、PSFが小さい場合は反復回数を増やすことで、それぞれの点源に見分けられている。

### 適当に散らばった点源の場合
```Python:適当に散らばった点源の場合（真の画像のコード部分）
"""シミュレーション"""
# 真の画像
true_px = 128
true_img = np.zeros([true_px, true_px])
true_img[77, 44] = 0.1
true_img[34, 60] = 0.1
true_img[77, 49] = 0.1
true_img[61, 93] = 0.1
true_img[57, 40] = 0.1
true_img[38, 95] = 0.1
true_img[89, 48] = 0.1
true_img[36, 70] = 0.1
true_img[30, 45] = 0.1
true_img[42, 77] = 0.1
"""以下省略"""
```
[RL法のシミュレーション](#rl法のシミュレーション)のコードの、真の画像と真の画像のカラーバーの最大値を``0.1``に変更している。
※適当に散らばった点源のピクセル座標はPythonの疑似乱数（[random](https://docs.python.org/ja/3/library/random.html)）により生成しているが、本題とは関係が薄いためコードを割愛する。

* 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/615804b6-91d5-2f4a-9c9d-0c2d3cfde218.png)
* ``psf_sigma=21, 17, 13, 9``で、反復回数5, 10, 15, 20回の結果
![図6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/f3a37875-705e-9e82-af1d-56e8935f871e.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数25, 50, 75, 100回の結果
![図7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/55c56621-b31d-67f3-4e2c-73a8aa1857ce.png)

このように、PSFが小さいほど見分けられる構造が増えることが分かる。また、RL法を行った方が見分けられる構造が増えるとことも分かる。

### 拡がった構造の場合
```Python:拡がった構造の場合（真の画像のコード部分）
"""シミュレーション"""
# 真の画像
true_px = 128
true_half_px = true_px // 2
x, y = np.mgrid[-true_half_px:true_half_px:, -true_half_px:true_half_px:]
pos = np.dstack((x, y))
rv1 = multivariate_normal([0, -8], [[4, 0], [0, 18]])
rv2 = multivariate_normal([0, 8], [[4, 0], [0, 18]])
true_img = 0.5*rv1.pdf(pos) + 0.5*rv2.pdf(pos)
"""以下省略"""
```
[RL法のシミュレーション](#rl法のシミュレーション)のコードの、真の画像と真の画像のカラーバーの表示形式を対数に変えている。

* 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/5a7ac37e-827e-9e50-6b6b-a736a24e8120.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数5, 10, 15, 20回の結果
![図8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/0dda801f-4a16-5d56-b6c9-b398fcd11f84.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数25, 50, 75, 100回の結果
![図9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/b8dbea0d-b8eb-d95d-f016-7d1a3f445188.png)

このように、点源のときと同様PSFが小さいほど、同じ反復回数で比べたときに真の画像に近いことが分かる。
また、PSFよりも小さい領域では[PSFを変化させる](#psfを変化させる)と似た形がRL法で得られることが分かる。このことから、PSFの拡がりよりも小さい領域において、**RL法では点源なのか拡がった構造なのか見分けることが困難である**と考えられる。

### 偏りのある拡がった構造の場合
```Python:偏りのある拡がった構造の場合（真の画像のコード部分）
"""シミュレーション"""
# 真の画像
true_px = 128
true_half_px = true_px // 2
x, y = np.mgrid[-true_half_px:true_half_px:, -true_half_px:true_half_px:]
pos = np.dstack((x, y))
rv1 = multivariate_normal([0, -8], [[4, 0], [0, 18]])
rv2 = multivariate_normal([0, 8], [[4, 0], [0, 18]])
true_img = 0.2*rv1.pdf(pos) + 0.8*rv2.pdf(pos)
"""以下省略"""
```
[RL法のシミュレーション](#rl法のシミュレーション)のコードの、真の画像と真の画像のカラーバーの表示形式を対数に変えている。

* 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/7dfbc90e-2901-18bd-3ea5-236af97acf6c.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数5, 10, 15, 20回の結果
![図10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/0064e0c9-20b3-da28-b462-e02c169bf3fd.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数25, 50, 75, 100回の結果
![図11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/9a9ec47f-0361-4300-cf58-5ead16b59989.png)

このように、点源のときと同様大きい値（明るい）の構造によって小さい値（暗い）の構造が見えにくくなることが分かる。ただ、PSFが小さい場合は反復回数を増やすことで、それぞれの点源に見分けられている。
また、PSFよりも小さい領域では[偏りのある点源の場合](#偏りのある点源の場合)と似た形がRL法で得られることが分かる。このことから、PSFの拡がりよりも小さい領域において、**RL法では点源なのか拡がった構造なのか見分けることが困難である**と考えられる。

# 1次元でのRL法のシミュレーションの実装（おまけ）
1次元でのRL法のシミュレーションの実装についても簡単に紹介する。
基本的なコードは2次元の場合と同じである。ただ、1次元のPSFは正規分布（SciPyの[``norm``](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy-stats-norm)）を使い、2次元のときのような高級な多変量正規分布は使わずに実装する。

## 点源の場合
```Python:RL法のシミュレーション（1次元バージョン）
"""シミュレーション（1次元バージョン）"""
# 真の画像
true_px = 128
true_img = np.zeros(true_px)
true_img[52] = 0.5
true_img[76] = 0.5

# PSF
psf_px = 127
mean = psf_px // 2
psf_sigma = 13
psf = norm.pdf(np.arange(psf_px), loc=mean, scale=psf_sigma)

# 真の画像をPSFで畳み込む（撮像画像）
obtained_img = convolve(true_img, psf, "same")

# RL法
iterations = 50
rl_img = restoration.richardson_lucy(obtained_img, psf, iterations)

"""シミュレーションの結果を図に表示"""
fig, (ax1, ax2, ax3, ax4) = plt.subplots(figsize=(18, 4), ncols=4)

# 真の画像
ax1.set_title("The true image ({0}×{0} pixel)\n".format(true_px))
ax1.plot(true_img)
ax1.set_xlabel("pixel")
ax1.set_ylabel("value")
ax1.set_ylim(0, 1)

# PSF
ax2.set_title("PSF σ={0} ({1}×{1} pixel)\n".format(psf_sigma, psf_px))
ax2.plot(np.arange(-63, 64, 1), psf)
ax2.set_xlabel("pixel")
ax2.set_ylabel("probability")
ax2.set_ylim(bottom=0)

# 撮像画像
ax3.set_title("The obtained image ({0}×{0} pixel)\n".format(true_px))
ax3.plot(obtained_img)
ax3.set_xlabel("pixel")
ax3.set_ylabel("value")
ax3.set_ylim(bottom=0)

# RL法による生成画像
ax4.set_title("Richardson-Lucy deconvolution\n image iterations={0}".format(iterations))
ax4.plot(rl_img)
ax4.set_xlabel("pixel")
ax4.set_ylabel("value")
ax4.set_ylim(bottom=0)

fig.tight_layout()
plt.show()
```

* 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/e6191f19-0124-9d33-d45e-ff7a36389c1d.png)

* コードの概念図
![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/5c574c1b-2c65-7b91-0a48-d66b2e7f228c.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数10~100回の結果<br>次図を作成したコードは、本題とは関係が薄いためここでは割愛してGoogle Colab内[こちら](https://colab.research.google.com/drive/19MpY6Ar3lOIP8CpQnVR_n8CckpW4ydLd?usp=sharing)で紹介する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/5c1f65ca-006e-72ba-324b-0464555f5db8.png)

## 偏りのある点源の場合
```Python:偏りのある点源の場合（真の画像のコード部分）
"""シミュレーション（1次元バージョン）"""
# 真の画像
true_px = 128
true_img = np.zeros(true_px)
true_img[52] = 0.2
true_img[76] = 0.8
"""以下省略"""
```

* 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/d56fdc21-78aa-59b7-ea1e-2f9a57850eb3.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数10~100回の結果<br>次図を作成したコードは、本題とは関係が薄いためここでは割愛してGoogle Colab内[こちら](https://colab.research.google.com/drive/19MpY6Ar3lOIP8CpQnVR_n8CckpW4ydLd?usp=sharing)で紹介する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/86838107-22b5-5331-65ed-4770deab5b8a.png)

## 適当に散らばった点源の場合

```Python:適当に散らばった点源の場合（真の画像のコード部分）
"""シミュレーション（1次元バージョン）"""
# 真の画像
true_px = 128
true_img = np.zeros(true_px)
true_img[98] = 0.2
true_img[107] = 0.2
true_img[10] = 0.2
true_img[66] = 0.2
true_img[124] = 0.2
"""以下省略"""
```

※適当に散らばった点源のピクセル座標はPythonの疑似乱数（[random](https://docs.python.org/ja/3/library/random.html)）により生成しているが、本題とは関係が薄いためコードを割愛する。

* 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/76661150-d782-8c2f-1fae-520a15f1a515.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数10~100回の結果<br>次図を作成したコードは、本題とは関係が薄いためここでは割愛してGoogle Colab内[こちら](https://colab.research.google.com/drive/19MpY6Ar3lOIP8CpQnVR_n8CckpW4ydLd?usp=sharing)で紹介する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/1b9181c8-0b7e-6179-d845-18b3385bc16b.png)

## 拡がった構造の場合

```Python:拡がった構造の場合（真の画像のコード部分）
"""シミュレーション（1次元バージョン）"""
# 真の画像
true_px = 128
signal_1 = norm.pdf(np.arange(true_px), loc=52, scale=6)
signal_2 = norm.pdf(np.arange(true_px), loc=76, scale=6)
true_img = 0.5*signal_1 + 0.5*signal_2
"""以下省略"""
```

* 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/69749017-7323-7945-68ec-0b804b4c1362.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数10~100回の結果<br>次図を作成したコードは、本題とは関係が薄いためここでは割愛してGoogle Colab内[こちら](https://colab.research.google.com/drive/19MpY6Ar3lOIP8CpQnVR_n8CckpW4ydLd?usp=sharing)で紹介する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/4d2e13fc-4e7f-706b-26a5-d5bd04ef89cb.png)

## 偏りのある拡がった構造の場合

```Python:偏りのある拡がった構造の場合（真の画像のコード部分）
"""シミュレーション（1次元バージョン）"""
# 真の画像
true_px = 128
signal_1 = norm.pdf(np.arange(true_px), loc=52, scale=6)
signal_2 = norm.pdf(np.arange(true_px), loc=76, scale=6)
true_img = 0.2*signal_1 + 0.8*signal_2
"""以下省略"""
```

* 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/4c4e732b-0342-9d80-92c2-d5f8b6d0edda.png)

* ``psf_sigma=21, 17, 13, 9``で、反復回数10~100回の結果<br>次図を作成したコードは、本題とは関係が薄いためここでは割愛してGoogle Colab内[こちら](https://colab.research.google.com/drive/19MpY6Ar3lOIP8CpQnVR_n8CckpW4ydLd?usp=sharing)で紹介する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1599258/f76002fb-95ad-0969-e359-68cba67ad271.png)


# まとめ
**真の画像とPSFが完全にわかっていて、かつノイズが乗らない理想的な状況でのシミュレーション**により、わかったことを以下にまとめる。

* RL法により得られた画像は、撮像画像と比べて真の画像に近いことがわかった。
* RL法の反復回数が5~20回で得られる画像は、25回以上で得られる画像よりも急激に画像が変化していた。
* RL法の反復回数が多くなるにつれて、真の画像により近い画像が得られた。
* PSFが小さいほど、RL法を行ったときに見分けられる部分が増えた（真の画像に近い画像が得られた）。
* 点源と拡がった構造のシミュレーションの結果、PSFの拡がりよりも小さい領域では似たようなイメージが得られた。<br>このことから、**RL法はPSFよりもある程度小さい領域では、そこに点源なのか拡がった構造が存在するのか見極めることが困難である**と考えられる。

中編ではRL法の結果をもう少し深掘りする予定である。

# 参考文献
* Richardson, William Hadley. 1972, JOSA, 62, 55–59
* Lucy, L. B. 1974, Astronomical Journal, 79, 745–754

# 関連資料
* [Richardson-Lucy deconvolutionの理論式を導出しよう](https://qiita.com/yusuke_s_yusuke/items/1f442023d9d1326da2a1)
* [Richardson-Lucy deconvolutionの特徴をシミュレーションで理解しよう（中編：結果の解析）](https://qiita.com/yusuke_s_yusuke/items/186d2b2af04391bcd7e0)