# Chapter 6. 連続画像

## 線形でShift-invariantな変換(transformation)
**Def. 変換**: 画像$f(\vec r)$→画像$g(\vec r)$の写像$\phi$

この章で主に扱う対象: 変換の中でも特に
* 線形性($\phi[a_1f_1 + a_2f_2] = a_1\phi[f_1] + a_2\phi[f_2]$)
* Shift-invariance($\phi[f(\vec r+\vec r_0)](\vec r) = \phi[f(\vec r)](\vec r+\vec r_0)$

を満たすもの

例:
* ぼかし
* x方向の微分

現実ではこれは近似的にしか成り立たないことが多い。例えば、
* 出力(明るさ)に上限がある
* 明るさは負にならない
* 周辺視野は光量が下がる

## 畳み込み
:::warning
**Def. 畳み込み**
$$(f\otimes h)(\vec r) := \iint_{-\infty}^{\infty}f(\vec r')h(\vec r-\vec r')d\vec r'$$

(注: この資料では、$d\vec r$とは$dxdy$の略記であるとする。また、今後無限積分を省略することがある。)
:::
**Prop**. (ある固定された関数$h(\vec r)$との)畳み込みは**線形**で**Shift-invariant**である。

証明: 線形性は積分の線形性、Shift-invarianceは変数変換からわかる

## デルタ超関数

**Def. (ディラックの)デルタ超関数**

任意の画像$f$について、
$$f=f\otimes \delta$$

が成り立つような超関数$\delta$を **(ディラックの)デルタ超関数** と呼ぶ。

注: $\delta$は普通の意味での関数ではなく、超関数である。

## 超関数とは

(この本よりも)猿渡先生のスライドに詳しい。

超関数$h$は、急減少関数(任意の指数関数よりも早く減衰する関数)である任意のテスト関数$t(x)$について、
$$\int_{-\infty}^{\infty}h(x)t(x)dx:=h[t]$$

の値が定まっているようなものである。(この式の左辺は形式的に積分のように書いているだけであり、実際にルベーグ積分を表すわけではない。)

注: 上では一次元で定義したが、二次元になっても同じで、$x$が$\vec r$になるだけである。

## Riemann-Lebesgueの補助定理


## 極限としてのデルタ超関数

* $\displaystyle \delta(x) = \lim_{\epsilon\to0}\begin{cases}0&(|x|>\epsilon)\\1/2\epsilon&(|x|\le \epsilon)\end{cases}$

* $\displaystyle \delta(x) = \lim_{\epsilon\to0}\frac{1}{\sqrt{2\pi}\epsilon}\exp\left(-\frac{1}{2}\left(\frac{x}{\epsilon}\right)^2\right)$
* $\displaystyle \delta(x) = \lim_{k\to\infty}\frac{\sin kx}{\pi x}$


## 畳み込みの普遍性

**Thm.** 全ての線形でShift-invariantな変換$\phi$はある固定された超関数$h$との畳み込みで表せる。

証明: $\phi[\delta] =: h$とすれば良い。
$$\begin{align}
\phi[f]&=\phi[f\otimes \delta]\\
&= \phi\left[\iint f(\vec r')\delta(\vec r - \vec r')d\vec r'\right]&&\text{(デルタ超関数の定義)}\\
&= \iint f(\vec r')\phi[\delta(\vec r - \vec r')]d\vec r' &&\text{(線形性)}\\
&= \iint f(\vec r')h(\vec r - \vec r')d\vec r' &&\text{(Shift-invariance)}\\
&= f\otimes h\\
\end{align}$$

注: $\phi[\delta]$は、テスト関数$t(\vec r)$に対して$\phi[\delta][t]=\delta[\phi[t]]$と定義される超関数である。
注: $h$は信号処理でインパルス応答と呼ばれているものである。ここでは*拡散関数*(Point-spread function)と呼ぶことにする。

## 畳み込みの性質

:::warning
**Prop**. 畳み込みについて次の性質が成り立つ。
* **可換性**: $f\otimes g = g\otimes f$
* **結合性**: $(a\otimes b) \otimes c = a\otimes (b\otimes c)$
:::

証明: 可換性は変数変換、結合性は積分順序の交換からわかる。

注: これは、複数の線形でShift-invariantな変換は一つにまとめることができることを意味している。

## 畳み込みとフーリエ変換

**Thm**. 全ての線形でShift-invariantな変換$\phi$は、フーリエ基底$e^{i\boldsymbol k\vec r}$を固有関数として持つ。すなわち、ある超関数$\lambda(\boldsymbol k)$が存在して、$$\phi[e^{i\boldsymbol k\vec r}]=\lambda(\boldsymbol k) e^{i\boldsymbol k\vec r}$$を満たす。(ただし、$\boldsymbol k$は任意の2次元横ベクトルである。)

証明: (本では省略されているが、)
$\phi$の拡散関数を$h$とすると、
$$\begin{align}
\phi[e^{i\boldsymbol k\vec r}]&=h\otimes(e^{i\boldsymbol k\vec r})\\
&=\iint h(\vec r')e^{i\boldsymbol k(\vec r-\vec r')}d\vec r'\\
&=\left(\iint h(\vec r')e^{-i\boldsymbol k\vec r'}d\vec r'\right)e^{i\boldsymbol k\vec r}\\
\end{align}$$

あとは$\lambda(\boldsymbol k)=\iint h(\vec r')e^{-i\boldsymbol k\vec r'}d\vec r'$とすれば良い。

**Thm**. フーリエ基底は完全系である。すなわち、任意の画像$f(\vec r)$は、
$$f(\vec r)=\frac{1}{4\pi^2}\iint A(\boldsymbol k)e^{i\boldsymbol k\vec r}d\boldsymbol k$$

という形で表せる。ただし、$A(\boldsymbol k)$は超関数である。

証明: $A(\boldsymbol k)=\iint f(\vec r)e^{-i\boldsymbol k\vec r}d\vec r$とし、代入して計算する。
$$\begin{align}\frac{1}{4\pi^2}\iint \iint f(\vec r')e^{-i\boldsymbol k\vec r'}d\vec r'e^{i\boldsymbol k\vec r}d\boldsymbol k &= \frac{1}{4\pi^2}\iint \iint f(\vec r')e^{i\boldsymbol k(\vec r - \vec r')}d\vec r'd\boldsymbol k
\\
&=\frac{1}{4\pi^2}\iint f(\vec r')\left[\frac{\sin k_x(x-x')\sin k_y(y-y')}{(x-x')(y-y')}\right]_{k_x, k_y=-\infty}^{\infty}d\vec r'\\
&=\frac{1}{4\pi^2}\iint f(\vec r')(2\pi\delta(x'))(2\pi\delta(y'))d\vec r'\\
&=\iint f(\vec r')\delta(\vec r')d\vec r'\\
&=f(\vec r)d\vec r\\
\end{align}$$

:::warning
$$f(\vec r)\mapsto A(\boldsymbol k)=\iint f(\vec r)e^{-i\boldsymbol k\vec r}d\vec r\tag{フーリエ変換}$$

$$\displaystyle A(\boldsymbol k)\mapsto f(\vec r)=\frac{1}{4\pi^2}\iint A(\boldsymbol k)e^{i\boldsymbol k\vec r}d\boldsymbol k\tag{逆フーリエ変換}$$
:::
注: 上の証明がフーリエ変換して逆変換すると元に戻ることの証明となっている。また、フーリエ変換と逆フーリエ変換がほぼ対称な関係になっていることに注意すれば、逆フーリエ変換のフーリエ変換が元に戻ることもほぼ自明である。

**Cor. 畳み込みのフーリエ変換はフーリエ変換の積**
$f(\vec r)$のフーリエ変換を$A(\boldsymbol k)$, $h(\vec r)$のフーリエ変換を$\lambda(\boldsymbol k)$とすると、$(h\otimes f)$のフーリエ変換は$\lambda(\boldsymbol k)A(\boldsymbol k)$である。

証明: 上の2つの定理から明らか。

## 微分とフーリエ変換
:::warning

**Prop**. $f(\vec r)$のフーリエ変換が$F(\boldsymbol k)$とすると、$\displaystyle \frac{\partial f}{\partial \vec v}(\vec r)$のフーリエ変換は$i\boldsymbol k\vec vF(\boldsymbol k)$である。

:::
証明:
$$ f(\vec r) = \frac{1}{4\pi^2}\iint F(\boldsymbol k)e^{i\boldsymbol k\vec r}d\boldsymbol k$$より、
$$\begin{align}\frac{\partial f}{\partial \vec v}(\vec r) &= \frac{1}{4\pi^2}\iint F(\boldsymbol k)\frac{\partial}{\partial \vec v}e^{i\boldsymbol k\vec r}d\boldsymbol k\\
&= \frac{1}{4\pi^2}\iint F(\boldsymbol k)i\boldsymbol k\vec ve^{i\boldsymbol k\vec r}d\boldsymbol k\\
\end{align}$$
(ここで微分と積分を交換していることに注意。変数が違うので問題ない。)

両辺をフーリエ変換すれば、(右辺は逆フーリエ変換のフーリエ変換が元に戻ることを用いて)
$$\iint \frac{\partial f}{\partial \vec v}(\vec r)e^{-i\boldsymbol k\vec r}d\vec r = F(\boldsymbol k)i\boldsymbol k\vec v$$

## ラプラシアンとフーリエ変換
:::warning
**Prop**. $f(\vec r)$のフーリエ変換が$F(\boldsymbol k)$とすると、$\displaystyle \Delta f := \frac{\partial^2 f}{\partial x^2 }+\frac{\partial^2 f}{\partial y^2 }$のフーリエ変換は$-|\boldsymbol k|^2F(\boldsymbol k)$である。
:::
証明: 上で証明した微分のフーリエ変換を合計4回使えば出る。

## 回転対称な変換
例:
* ラプラシアン
* 勾配の大きさの二乗
* 二階変分の二乗$\displaystyle \left(\frac{\partial^2 f}{\partial x^2}\right)^2+2\left(\frac{\partial^2 f}{\partial x\partial y}\right)\left(\frac{\partial^2 f}{\partial y\partial x}\right)+\left(\frac{\partial^2 f}{\partial y^2}\right)^2$

## ハンケル変換
今、回転対称な変換についての二次元フーリエ変換をよりきれいにできないか考える。

フーリエ変換
$$F(\boldsymbol k)=\iint f(\vec r)e^{-i\boldsymbol k\vec r}d\vec r$$

を、回転座標系に変数変換する。すなわち、
$$\vec r = r\left[\begin{matrix}\cos \theta\\\sin\theta\end{matrix}\right], \boldsymbol k = k\left[\begin{matrix}\cos\alpha & \sin\alpha\end{matrix}\right]$$
とする。$f(\vec r)$は回転対称であるので、$\bar f(r)$と書け、それをフーリエ変換した$F(\boldsymbol k)$も回転対称なので$\bar F(k)$と書ける。ここで、
$$\boldsymbol k\vec r = kr\cos(\theta-\alpha)$$

$$d\vec r = rdrd\theta$$

に注意して、上の積分を書き直すと、
$$\bar F(k) = \int_{r=0}^{\infty}\int_{\theta=0}^{2\pi} r\bar f(r)e^{-irk\cos(\theta-\alpha)}drd\theta$$

$\theta$に関する積分だけを取り出して、$\theta-\alpha =: \theta'$と置き換えると、
$$\begin{align}\int_0^{2\pi}e^{-irk\cos(\theta-\alpha)}d\theta
&= \int_0^{2\pi}e^{-irk\cos\theta'}d\theta'&&\text{(変数変換)}\\
&= 2\int_0^{\pi}\cos(rk\cos\theta')d\theta'&&\text{(周期性と対称性を利用)}\\
&= 2\pi J_0(rk)\\
\end{align}$$

ただし、$J_n(z)$は第一種ベッセル関数を表す。

よって、ハンケル変換
:::warning
$$\bar F(k) = 2\pi\int_{0}^{\infty}rJ_0(rk)\bar f(r)dr$$
:::
が導かれる。同様に、逆ハンケル変換
:::warning
$$\bar f(r) = \frac{1}{2\pi}\int_{0}^{\infty}kJ_0(rk)\bar F(k)dk$$
:::

も導ける。

注: $F(\boldsymbol k)$は一般に複素数だが、$\bar F(k)$は実数である。これは、回転対称性が線対称性も含むからだ。

▲$r J_0(r)$のグラフ

### [Bessel関数に関する補足](http://math.univ-lyon1.fr/~iohara/Doc/Jap/Bessel.pdf)
Bessel関数$J_n(z)$は次の母関数のローラン展開の係数として定義される:
$$e^{-\frac{1}{2}z\left(t - \frac{1}{t}\right)} = \sum_{n\in \mathbb Z}J_n(z)t^n$$

留数定理より、
$$J_0(z)=\frac{1}{2\pi i}\oint_C (e^{-\frac{1}{2}z\left(t - \frac{1}{t}\right)}/t) dt$$

ただし、$C$は原点を周る閉路。これに$t=e^{i\beta}$で置換すると、

$$\begin{align}J_0(z)&=\frac{1}{2\pi i}\int_0^{2\pi} (e^{-\frac{1}{2}z\left(e^{i\beta} - e^{-i\beta}\right)}/e^{i\beta}) ie^{i\beta}d\beta\\
&=\frac{1}{2\pi}\int_0^{2\pi} e^{-\frac{1}{2}z\left(e^{i\beta} - e^{-i\beta}\right)}d\beta\\
&=\frac{1}{2\pi}\int_0^{2\pi} e^{iz\sin \beta}d\beta\\
&=\frac{1}{2\pi}\int_0^{2\pi} \cos(z\sin \beta)d\beta\\
&=\frac{1}{2\pi}\int_0^{2\pi} \cos(z\cos \beta')d\beta'\\
\end{align}$$



▲[太鼓の振動(波動方程式の円境界条件を満たす解)](http://www.natural-science.or.jp/article/20101108184753.php)

## Hankel変換の応用

画像に対する「理想的な」ローパスフィルターを考えたい

信号処理で言えば周波数特性が矩形になる理想的なフィルタに相当するもの

**ただし、回転対称性を課す**

すなわち、Hankel変換が
$$\bar H(k) = \begin{cases}1&(k < k_c)\\0&(k\ge k_c)\end{cases}$$

となるような回転対称フィルタ$h(r)$を求めたい

逆Hankel変換を用いて、

$$\begin{align}\bar h(r) &= \frac{1}{2\pi}\int_{0}^{\infty}kJ_0(rk)\bar H(k)dk\\
&=\frac{1}{2\pi}\int_{0}^{k_c}kJ_0(rk)dk\\
\end{align}$$

ここで、$rk=:z$と置き換える。

$$\begin{align}\bar h(r) &=\frac{1}{2\pi r^2}\int_{0}^{rk_c}zJ_0(z)dz\\
\end{align}$$

ここで、Bessel関数の性質
$$\frac{d}{dz}(zJ_1(z))=zJ_0(z)$$

を用いると、

$$\begin{align}\bar h(r) &=\frac{k_c}{2\pi r}J_1(rk_c) = \frac{k_c^2}{2\pi }\frac{J_1(rk_c)}{rk_c}\\
\end{align}$$

ここで、$\displaystyle \lim_{z\to 0}\frac{J_1(z)}{z} = \frac{1}{2}$である。


注: 理想フィルタは上のように波打つので、現実的にはGaussianのようなものを使うことが多い。周波数特性がそんなに急峻に落ちてくれなくても実用的には問題ない。

## ベッセル関数についての補足2

$\displaystyle\frac{d}{dz}(zJ_1(z))=zJ_0(z)$を証明する。

ベッセル関数の定義
$$e^{-\frac{1}{2}z\left(t - \frac{1}{t}\right)} = \sum_{n\in \mathbb Z}J_n(z)t^n$$

の両辺を$z$で微分して、

$$J_{n-1}(z)-J_{n+1}(z)=2\frac{d}{dz}J_n(z)$$

また、定義式の両辺を$t$で微分して、
$$J_{n-1}(z)+J_{n+1}(z)=\frac{2n}{z}J_n(z)$$

これらの両辺を足して、
$$J_{n-1}(z)=\frac{d}{dz}J_n(z) + \frac{n}{z}J_n(z)$$

両辺に$z^n$を掛けると、
$$\begin{align}z^nJ_{n-1}(z) &= z^n\frac{d}{dz}J_n(z) + nz^{n-1}J_n(z) \\&= \frac{d}{dz}(z^nJ_n(z))\end{align}$$

これに$n=1$を代入したものが最初の式である。
なお、途中で出てきた2つの漸化式を引くことで、

$$\frac{d}{dz}(z^{-n}J_{n+1}(z))=-z^{-n}J_{n}(z)$$

も導ける。

## ぼかし
ぼかしは様々な種類があるが、代表的なものには
* ガウシアンぼかし
* $\displaystyle h(\vec r) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}{2}\frac{|\vec r|^2}{\sigma^2}\right)$
* これのフーリエ変換は(ハンケル変換を用いることもできるが、xとyに変数分離可能なのでそうやって求めたほうが簡単で、)
$\displaystyle H(\boldsymbol k) = \exp\left(-\frac{1}{2}\sigma^2|\boldsymbol k|^2\right)$
* $\sigma$が大きいほど、低い周波数しか通さない
$\sigma$と「通す周波数の高さに比例する量」の積は一定(cf. 不確定性原理)
* ピンぼけ
* ピンぼけは、一点に集まる光がある範囲(レンズの大きさによる)に一様に広まると考える
$\displaystyle h(\vec r) = \begin{cases}1/\pi R^2&(|\vec r| < R)\\0&(|\vec r|\ge R)\end{cases}$
* そのハンケル変換は、(理想フィルタと同じように求まり、)
$\displaystyle \bar H(k) = 2\frac{J_1(Rk)}{Rk}$
これも$R$が大きいほど低い周波数しか通さなくなる
全く通さない周波数もある
* モーションブラー(被写体ブレ)
* これは回転対称じゃない
簡単なため、$x$軸に沿って$2l$だけ動いたとすると、拡散関数が線分
$\displaystyle h(\vec r) = \begin{cases}\delta(y)/2l&(|x| < l)\\0&(|x|\ge l)\end{cases}$
となる。
* これのフーリエ変換は(信号処理でもおなじみの)
$\displaystyle H(\boldsymbol k) = \frac{\sin(lk_x)}{lk_x}$

## 復元

得られた画像は、ほしい画像にある既知のフィルタ$H(\boldsymbol k)$が掛けられたものだった(e.g. ピンぼけしていた)として、ほしい画像を復元したい。

理想的には、フーリエ変換が$H'(\boldsymbol k) = 1/H(\boldsymbol k)$なフィルタを作って得られた画像にかければよい。(言い換えれば、$h(\vec r) \otimes h'(\vec r) = \delta(\vec r)$となるフィルタ$h'(\vec r)$を求めて画像にかければ良い)

しかし、これを愚直に行うと次のような問題に当たる。
* $H(\boldsymbol k)$が0(全く通さない波数)のとき、その逆数が発散する
* $H(\boldsymbol k)$がほとんど0のとき、その逆数が大きくなり、ノイズに弱くなる

そのため、低すぎる$H(\boldsymbol k)$の値に対して発散しないように対策を行う必要がある。例えば、
* $|H'(\boldsymbol k)| = 1/\max(|H(\boldsymbol k)|, \epsilon)$
* $H'(\boldsymbol k) = H(\boldsymbol k)/(H(\boldsymbol k)^2+ \epsilon^2)$
など

## 相関とPower spectrum

**Def. 相関**

画像$a(\vec r)$と$b(\vec r)$の相関$a*b$を、
$$(a*b)(\vec r) := \iint a(\vec r'-\vec r)b(\vec r') d\vec r$$

と定義する。

注: 相関は、$a$を反転させて$b$と畳み込んだものと等しい。
注: $a*a$を$a$の**自己相関**という。自己相関は点対称である。
注: コーシー・シュワルツの定理から、自己相関は原点$(a*a)(\vec 0)$で最大値をとる。同様に、$b$が$a$をずらしたものであれば、$a*b$はそのずれ量のところで最大値を取る。ただし、画像が周期的な場合などは、他のところでも同じ最大値を取る可能性がある。
<!--
**Thm. Parsevalの定理**

画像$a(\vec r)$と画像$b(\vec r)$の内積を
$$\left<a, b\right> =\iint a(\vec r)b(\vec r)d\vec r$$

と定義する。$A(\boldsymbol k), B(\boldsymbol k)$をそれぞれのフーリエ変換とすると、それらの内積
$$\left<A, B\right> =\iint A^*(\boldsymbol k)B(\boldsymbol k)d\vec k$$

は元の内積に等しい。すなわち、
:::warning
$$\left<a, b\right> = \left<A, B\right>$$
:::
が成り立つ。

証明: 略(信号処理の授業でやった)
-->
**Thm. 自己相関とPower spectrum**

自己相関$(a*a)(\vec r)$のフーリエ変換とPower spectrum$|A(\boldsymbol k)|^2=A^*(\boldsymbol k)A(\boldsymbol k)$は等しい。

証明: 略(信号処理の授業でやった)

注: パワースペクトルは画像が平行移動したりしても変わらない

**Prop. White noiseのフーリエ変換は、位相が一様ランダムで大きさが一定となる。**

証明: 略(信号処理の授業でやった)

## Wiener Filter

**Thm. Wiener Filter**

Power spectrumが$\Phi_{bb}(\boldsymbol k)$の(未知)画像$b(\vec r)$に、Power spectrumが$\Phi_{nn}(\boldsymbol k)$のノイズ$n(\vec r)$が乗った画像$i(\vec r)=b(\vec r)+n(\vec r)$が得られているとする。このとき、期待二乗誤差を最小にするようなフィルタ$h(\vec r)$
$$\arg\min_{h(\vec r)}\mathbb E\left[\iint |h(\vec r)\otimes (b(\vec r) + n(\vec r)) - b(\vec r)|^2d\vec r\right]$$

のフーリエ変換$H(\boldsymbol k)$は、次の式で書ける:
:::warning
$$H(\boldsymbol k) = \frac{\Phi_{bb}(\boldsymbol k)}{\Phi_{bb}(\boldsymbol k)+\Phi_{nn}(\boldsymbol k)}$$
:::

証明: これも信号処理の授業でやった。変分法を使うと求まる。(「非因果性Wienerフィルタ」と対応する。)

**Cor.**
既知のフィルタ$p(\vec r)$がかかった信号$p(\vec r)\otimes b(\vec r)$との期待二乗誤差を最小にするフィルタ
$$\arg\min_{h(\vec r)}\mathbb E\left[\iint |h(\vec r)\otimes (b(\vec r) + n(\vec r)) - \vec p(\vec r)\otimes b(\vec r)|^2d\vec r\right]$$

のフーリエ変換$H(\boldsymbol k)$は、
$$H(\boldsymbol k) = \frac{\Phi_{bb}(\boldsymbol k)}{\Phi_{bb}(\boldsymbol k)+\Phi_{nn}(\boldsymbol k)}P(\boldsymbol k)$$

で与えられる。

## Image Model
Wiener filterを適用するには、現実の「正しい」画像のPower spectrumを知る必要がある。

* 高さ$H$, 幅$W$の長方形のフーリエ変換は、(おなじみのように)
$$F(\boldsymbol k) = \frac{\sin(k_x W)}{k_x}\frac{\sin(k_y H)}{k_y}$$

に比例する。振動成分を無視して、振幅が$1/k_xk_y$に比例するので、大雑把に$|k|^{-2}$に比例すると言える。
* 半径$R$の円のフーリエ変換は、上で見たように
$$\bar F(k) = 2R^2\frac{J_1(kR)}{kR}$$

に比例する。$z$が大きいときには$\displaystyle J_1(z)\approx\sqrt{\frac{2}{\pi z}}\sin(z-\pi/4)$という性質を使うと、振幅が$k^{-3/2}$に比例している。

このように、現実の画像は低い周波数にパワーが集中する傾向にある。

さらに、光学系の性能によって高い周波数ではこれ以上に早く減衰することもある。(例:望遠鏡では波長とレンズ口径によって決まる絶対的なカットオフ周波数があり、それ以上の周波数は全く入ってこれない。顕微鏡でも似たような限界が存在する。)

## Image reproduction

自然の画像はダイナミックレンジ(暗いところと明るいところの明るさの比)が広いが、それを表示する機材はそこまでダイナミックレンジが広くないことが多い。例えば、
* 写真フィルム: 1:100強
* 新聞紙: 1:10程度

そこに広いダイナミックレンジを無理やり「押し込める」ための手法の例:
* 自然の画像$I(\vec r)$を$I(\vec r)^\gamma$のように指数$0<\gamma<1$を用いて調整する
(人間の感覚が暗い領域で明るさの変化に敏感なため)
* 弱いハイパスフィルターを通して低い周波数に集中しているパワーを分散させる

# Chapter 7. 離散画像
## 画像サイズが有限な場合

**Prop. δ超関数列のフーリエ変換はδ超関数列になる**
$$f(\vec r) = \sum_{n, m\in \mathbb Z}\delta\left(\vec r-\left[\begin{matrix}nW\\mH\end{matrix}\right]\right)$$

のフーリエ変換は、
$$F(\boldsymbol k) = \sum_{n, m\in \mathbb Z}\delta\left(\boldsymbol k-\left[\begin{matrix}\frac{2\pi}{W}k_x&\frac{2\pi}{H}k_y\end{matrix}\right]\right)$$

証明: 計測通論Cでやった(ものの2次元版)ので省略

注: 有限な画像サイズの場合、フーリエ変換は全ての$\boldsymbol k$について求める必要がない。
∵画像を周期拡張して考えると、周期的な信号のフーリエ変換は離散的になるから
(周期的な画像は、ある元画像とδ超関数列の畳み込みで表せるので、フーリエ変換するとδ超関数列と元画像フーリエ変換の積で表せる)

## 画像が離散な場合

上の話の逆フーリエ変換版を考えると、離散画像のフーリエ変換は周波数空間で周期的になることがわかる。そのため、離散画像の場合はフーリエ変換を一周期分だけ覚えておけばいい。

## サンプリング定理

**Thm. Sampling Theorem**
もし元画像$I(\vec r)$がナイキスト周波数$\frac{\pi}{d}$以上の周波数成分を含まないなら、間隔$d$でサンプリングされた画像から元画像が完全に復元できる。

証明: 計測通論Cでやった。
間隔$d$でサンプリングされた画像をフーリエ変換すると、周期$\frac{2\pi}{d}$の周波数特性が得られる。元画像がナイキスト周波数以上の周波数成分を含まないなら、「安全に」$\pm\frac{\pi}{d}$以下の周波数成分のみを残して他をカットすることができ、それを逆フーリエ変換すると元画像が再現できる。

注: ナイキスト周波数以上の周波数成分を含む場合、ナイキスト周波数で折り返したような画像ができてしまう(エイリアシング)。そのため、サンプリングをする時は予めナイキスト周波数以上の周波数をカットしておかなければいけない。

注: この復元過程は、sincフィルタを畳み込むことに相当する。

注: ナイキスト周波数$K$以下ということは、
$$\iint I''(\vec r)^2d\vec r \le K^4 \iint I(\vec r)^2d\vec r$$

を意味する。すなわち、画像があまり激しく変化できない。

注: 二次元の場合、画像をサンプリングする時の素子の並べ方は普通は四角いが、六角形に並べると実はより効率が良くなる。四角く並べるとナイキスト周波数の限界の線が周波数空間で四角く書けるわけだが、六角形に並べるとそれが六角形で書け、結果として最も低い方向のナイキスト周波数を同じに保った時の素子の数を減らせる。

## 離散フーリエ変換

以上のことを使って、離散フーリエ変換を定義できる。

**Def. 離散フーリエ変換**
サイズが$W\times H$の離散画像$f[\vec r]$の離散フーリエ変換$F[\boldsymbol k]$は、

$$F[\boldsymbol k] = \sum_{\vec r} f[\vec r] e^{-2\pi i \left(\frac{xk_x}{W} + \frac{yk_y}{H}\right)}$$

離散逆フーリエ変換$f[\vec r]$は、

$$f[\vec r] = \frac{1}{WH}\sum_{\boldsymbol k} F[\boldsymbol k] e^{2\pi i \left(\frac{xk_x}{W} + \frac{yk_y}{H}\right)}$$

ただし、和は全ての範囲についての和を意味する。

注: 画像を周期拡張することを考えているが、普通は左端と右端は同じ色ではないので、不連続的になって、大きなノイズが乗る。これを回避するための方法は次のようなものがある。
* 周期拡張する代わりに、鏡像を周りに拡張する。これによって周期が2倍になるが、偶関数となるためフーリエ変換の虚数部分がなくなる(離散コサイン変換)。
* 周辺で値が小さくなるように窓関数を掛けてからフーリエ変換する。
例えば、Hanning窓
$$\cos^2(\pi \frac{x}{W})\cos^2(\pi \frac{y}{H})$$

などが考えられる。
Hanning窓を掛けた場合、周波数領域では
が畳み込まれる

注: 離散フーリエ変換は、高速フーリエ変換(FFT)というアルゴリズムを用いて$O(WH(\log W +\log H))$の時間で求められる。

**Prop. 正規分布に従うランダムノイズの離散フーリエ変換**
$f[\vec r]$が平均0、標準偏差$\sigma$の正規分布に従う独立な乱数の場合、$F[\boldsymbol k]$は平均0、標準偏差$\sqrt{WH}\sigma$の正規分布に独立に従う乱数となる。

## 循環畳み込み
離散フーリエ変換して積を求めて離散逆フーリエ変換することで、畳込みを求めることができる。

ただし、周期拡張して考えているので、反対側の端から色が「漏れ出す」ことに注意(循環畳み込みになってしまう)。これを防ぐには、
* 周囲から漏れ出した部分を使わないように、畳み込んだ後真ん中だけ切り出して使う
* 画像の周りに畳み込むFilterのSupportの分だけ空白を追加してから畳み込む(何の色を「空白」として追加すればいいかという問題はある)

という方法が考えられる。