【Python】OpenCVのGaussianBlurについて

## はじめに

OpenCVの`GaussianBlur`を使用していたのですが、自分が想定していたのとは異なる数値結果が出力されたので、どのような仕様になっているのか少し調べてみました。


## Gaussian Blurとは

`cv.GaussianBlur` は、OpenCVライブラリで提供されている関数で、画像にガウシアン（ガウス分布に基づく）ブラー（ぼかし）効果を適用するために使います。詳しくは[公式ドキュメント](https://docs.opencv.org/4.x/d4/d13/tutorial_py_filtering.html)等を参照して下さい。

下図左が原画像で、右がブラー後の画像です。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/d041f3c1-11ce-bba1-11e9-95cde70ea1e6.png" width=200px>

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/7f3382eb-4016-4277-b09f-2336d4d9201d.png" width=200px>

画像引用元：https://github.com/opencv/opencv/blob/master/samples/data/baboon.jpg

ブラー後の輝度は、ガウス分布に基づくカーネルと元画像を畳み込むことで計算されます。カーネルを画像にすると次のようになります。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/bbf93d59-755b-25ff-1c33-952a21b685df.png" width=200px>


ガウシアンブラーは、画像処理の基本的な処理の一つであり、多くの画像処理タスクで広く使われています。例えば以下のような場面で使われます。

- ノイズ除去
    ガウシアンブラーは、画像からノイズを除去するのに効果的です。

- 画像の前処理：
    コンピュータビジョンや画像処理のアプリケーションにおいて、画像のディテールを減らすためにガウシアンブラーを適用することがあります。

## 環境

- opencv-python==4.7.0.72

## 使い方

基本的な使い方を説明します。

```python
import cv2 as cv

# 画像を読み込む
img = cv.imread('path/to/image.jpg')

# ガウシアンブラーを適用
blurred_img = cv.GaussianBlur(img, ksize=(31, 31), sigmaX=0)
```

- カーネルサイズ (ksize) は、カーネルの大きさです。カーネルが大きいほど、畳み込む範囲が大きくなり、ブラー効果は強くなります。ksize は奇数である必要があります。 ksize を0に設定すると、標準偏差から自動的にカーネルサイズを計算します。
- 標準偏差 (sigmaX, sigmaY) を変更することで、ブラーの「滑らかさ」を調整できます。標準偏差が大きいほど、ブラーは滑らかになりますが、細部が失われやすくなります。sigmaX と sigmaY の両方を0に設定すると、カーネルサイズから標準偏差を自動的に計算します。

このサンプルコードでは、31x31のカーネルサイズと、X方向の標準偏差を0（OpenCVが自動的に計算）でガウシアンブラーを適用しています。

## 自動計算について

この記事の本題。自動計算方法が思っていた処理と微妙に異なるところがあったので、それについて説明します。

### カーネルサイズから標準偏差を自動計算する方法

[ドキュメント](https://docs.opencv.org/4.8.0/d4/d86/group__imgproc__filter.html#gac05a120c1ae92a6060dd0db190a61afa)に計算方法が書いてあります。

```math
\sigma = 0.3 \times ((\text{ksize} - 1) \times 0.5 - 1) + 0.8
```

ただし、この計算方法はksizeが13以上のときのみ当てはまり、ksizeが3, 5, 7, 9, 11のときは予め用意された数値で計算されます。

カーネルは`cv.getGaussianKernel`関数で求めることができます。この関数では1次元のカーネルが計算されます。2次元にしたい場合は、外積を計算すれば2次元にすることができます。

```python
cv.getGaussianKernel(ksize=3, sigma=0)
# array([[0.25],
#        [0.5 ],
#        [0.25]])

cv.getGaussianKernel(ksize=5, sigma=0)
# array([[0.0625],
#        [0.25  ],
#        [0.375 ],
#        [0.25  ],
#        [0.0625]])

cv.getGaussianKernel(ksize=7, sigma=0)
# array([[0.03125 ],
#        [0.109375],
#        [0.21875 ],
#        [0.28125 ],
#        [0.21875 ],
#        [0.109375],
#        [0.03125 ]])

cv.getGaussianKernel(ksize=9, sigma=0)
# array([[0.015625  ],
#        [0.05078125],
#        [0.1171875 ],
#        [0.19921875],
#        [0.234375  ],
#        [0.19921875],
#        [0.1171875 ],
#        [0.05078125],
#        [0.015625  ]])

cv.getGaussianKernel(ksize=11, sigma=0)
# array([[0.00881223],
#        [0.02714358],
#        [0.06511406],
#        [0.12164907],
#        [0.17699836],
#        [0.20056541],
#        [0.17699836],
#        [0.12164907],
#        [0.06511406],
#        [0.02714358],
#        [0.00881223]])
```

`getGaussianKernel`関数はksizeが正の値でないとエラーになるので注意が必要です。

ksize が13以上の場合は次の関数で同じカーネルを得られます（精度の問題か、全桁は一致しません）。

```python
def get_gaussian_kernel(ksize, sigma):
    if ksize <= 0:
        raise ValueError("ksize must be > 0.")

    if not isinstance(ksize, int):
        raise ValueError("ksize must be int.")

    if ksize % 2 == 0:
        raise ValueError("ksize must be odd.")

    if sigma < 0:
        raise ValueError("sigma must be >= 0.")

    if sigma == 0:
        sigma = 0.3 * ((ksize - 1) * 0.5 - 1) + 0.8

    x = np.arange(ksize)
    g = np.exp(-(x - (ksize - 1) / 2)**2 / (2 * sigma**2))
    g = g / g.sum()
    return g.reshape(-1, 1)
```


### 標準偏差からカーネルサイズを自動計算する方法

[ソースコード](https://github.com/opencv/opencv/blob/4.x/modules/imgproc/src/smooth.dispatch.cpp)の`createGaussianKernels`関数を見ると、カーネルサイズの計算式は、次のようになっています。

```math
\begin{array}{ll}
(\sigma \times D \times 2 + 1) | 1
\end{array}
```

$\sigma$は標準偏差、$D$は画像が8ビットであれば3、それ以外であれば4です。$|$はビットごとのOR演算子で$|1$は偶数であれば1を足す、奇数であれば何もしません。

これを使用して再現すると次のようになると思います。

```python
def gaussian_blur(img, ksize, sigma):
    if ksize <= 0:
        if img.dtype == 'uint8':
            D = 3
        else:
            D = 4
        new_ksize = round(sigma * D * 2 + 1)
        if new_ksize % 2 == 0:
            new_ksize = new_ksize + 1
    else:
        new_ksize = ksize
    print("ksize", new_ksize)

    kernel = cv.getGaussianKernel(
        ksize=new_ksize, sigma=sigma, ktype=cv.CV_32F)
    dst = cv.sepFilter2D(img, ddepth=-1, kernelX=kernel, kernelY=kernel)
    
    return dst
```

あとはこれを`cv.GaussianBlur`と比較して一致すれば万々歳...だったのですが、なぜか一致しません。

```python
path = "baboon.jpg"
im = cv.imread(path)

k = 0
sigma = 5

dst1 = cv.GaussianBlur(im, (k, k), sigma, cv.BORDER_REPLICATE)

dst2 = gaussian_blur(im, ksize=k, sigma=sigma)

(dst1 - dst2).mean()  # 9.298937479654947
```

左が`dst1`で右が`dst2`です。見た目はほぼ同じです。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/01a3ce88-e317-4ee5-4d7e-571e1af11b18.png" width=200px>

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/c6ad6527-8ae3-fe2e-6be9-24f060662ea1.png" width=200px>

次のようにするとかなり違うことがわかります。

```python
img = np.zeros((256, 256), dtype=np.uint8)
img[:, 0] = 255

k = 0
sigma = 40

dst1 = cv.GaussianBlur(im, (k, k), sigma, cv.BORDER_REPLICATE)

dst2 = gaussian_blur(im, ksize=k, sigma=sigma)
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/a97854c6-e950-c880-8ffe-676fe5a65f63.png" width=200px>

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/204b5dd2-f588-f502-eaa4-de58ac492805.png" width=200px>

型が異なる場合でもまた違います。

```python
img = np.zeros((256, 256), dtype=np.float64)
img[:, 0] = 255

k = 0
sigma = 40

dst1 = cv.GaussianBlur(im, (k, k), sigma, cv.BORDER_REPLICATE)

dst2 = gaussian_blur(im, ksize=k, sigma=sigma)
```

3枚目は1枚目と2枚目の差です。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/ffb4d3df-f064-6ae8-636f-97834d85a715.png" width=200px>

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/69deefdd-618b-acc5-e83a-7877441ea5ff.png" width=200px>

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/ae52c669-e5ec-63bc-0268-a4bd13109bb2.png" width=200px>

`GaussianBlur`は内部で`sepFilter2D`を使用していると思うのですが、処理が違うようです。これ以上はよく分かりませんでした。

## おわりに

当初書く予定の記事とは違ったものになってしまいました。

OpenCVの処理を再現する場合は、使用する関数はもちろん引数も完全に一致させておかないと、想定と異なる処理になっているかもしれないので注意が必要です。