OpenCVの2Dガウシアンフィルタ：GaussianBlurの特徴、注意点と対策

cv2.GaussianBlurで2Dガウシアンフィルタが利用できる

長所は、シグマが0でも動くこと

import cv2 as cv

img=np.zeros((5,5))
img[2,2]=1
img
# array([[0., 0., 0., 0., 0.],
#        [0., 0., 0., 0., 0.],
#        [0., 0., 1., 0., 0.],
#        [0., 0., 0., 0., 0.],
#        [0., 0., 0., 0., 0.]])

cv.GaussianBlur(img, (3, 3), 0)
# array([[0.    , 0.    , 0.    , 0.    , 0.    ],
#        [0.    , 0.0625, 0.125 , 0.0625, 0.    ],
#        [0.    , 0.125 , 0.25  , 0.125 , 0.    ],
#        [0.    , 0.0625, 0.125 , 0.0625, 0.    ],
#        [0.    , 0.    , 0.    , 0.    , 0.    ]])

cv.GaussianBlur(img, (3, 3), 1)
# array([[0.        , 0.        , 0.        , 0.        , 0.        ],
#        [0.        , 0.07511361, 0.1238414 , 0.07511361, 0.        ],
#        [0.        , 0.1238414 , 0.20417996, 0.1238414 , 0.        ],
#        [0.        , 0.07511361, 0.1238414 , 0.07511361, 0.        ],
#        [0.        , 0.        , 0.        , 0.        , 0.        ]])

短所１．入力(img)が存在しないと、フィルタ自体を出力できない

※ここではど真ん中が１、それ以外が0の、1まわり大きい行列を先に用意することで解消していた

短所２．カーネルが偶数サイズだと動かない

※元ネタのgetGaussianKernelは偶数サイズでも動く

cv.GaussianBlur(img, (4, 4), 0)
# Traceback (most recent call last):
#   File "<stdin>", line 1, in <module>
# cv2.error: OpenCV(4.7.0) /home/conda/feedstock_root/build_artifacts/libopencv_1688234562235/work/modules/imgproc/src/smooth.dispatch.cpp:293: error: (-215:Assertion failed) ksize.width > 0 && ksize.width % 2 == 1 && ksize.height > 0 && ksize.height % 2 == 1 in function 'createGaussianKernels'

cv.getGaussianKernel(4, 0)
# array([[0.11165005],
#        [0.38834995],
#        [0.38834995],
#        [0.11165005]])

短所３．入力画像よりカーネルサイズが２小さくないと計が１になるよう正規化されない

※これは全体を自身のsumで割れば、正規化ぽく実装できるが、これは実際にはサイズを＋２大きくした時の正規化した行列と中身が異なるので、なんちゃって正規化である

import numpy as np
img.shape
# (5, 5)

np.sum(cv.GaussianBlur(img, (3, 3), 0))
# 1.0

cv.GaussianBlur(img, (5, 5), 0)
# array([[0.015625, 0.03125 , 0.046875, 0.03125 , 0.015625],
#        [0.03125 , 0.0625  , 0.09375 , 0.0625  , 0.03125 ],
#        [0.046875, 0.09375 , 0.140625, 0.09375 , 0.046875],
#        [0.03125 , 0.0625  , 0.09375 , 0.0625  , 0.03125 ],
#        [0.015625, 0.03125 , 0.046875, 0.03125 , 0.015625]])
np.sum(cv.GaussianBlur(img, (5, 5), 0))
# 1.265625

cv.GaussianBlur(img, (5, 5), 0) / np.sum(cv.GaussianBlur(img, (5, 5), 0))
# array([[0.01234568, 0.02469136, 0.03703704, 0.02469136, 0.01234568],
#        [0.02469136, 0.04938272, 0.07407407, 0.04938272, 0.02469136],
#        [0.03703704, 0.07407407, 0.11111111, 0.07407407, 0.03703704],
#        [0.02469136, 0.04938272, 0.07407407, 0.04938272, 0.02469136],
#        [0.01234568, 0.02469136, 0.03703704, 0.02469136, 0.01234568]])
np.sum(cv.GaussianBlur(img, (5, 5), 0) / np.sum(cv.GaussianBlur(img, (5, 5), 0)))
# 1.0

img=np.zeros((7,7))
img[3,3]=1
img.shape
# (7, 7)
cv.GaussianBlur(img, (5, 5), 0)
# array([[0.        , 0.        , 0.        , 0.        , 0.        , 0.        , 0.        ],
#        [0.        , 0.00390625, 0.015625  , 0.0234375 , 0.015625  , 0.00390625, 0.        ],
#        [0.        , 0.015625  , 0.0625    , 0.09375   , 0.0625    , 0.015625  , 0.        ],
#        [0.        , 0.0234375 , 0.09375   , 0.140625  , 0.09375   , 0.0234375 , 0.        ],
#        [0.        , 0.015625  , 0.0625    , 0.09375   , 0.0625    , 0.015625  , 0.        ],
#        [0.        , 0.00390625, 0.015625  , 0.0234375 , 0.015625  , 0.00390625, 0.        ],
#        [0.        , 0.        , 0.        , 0.        , 0.        , 0.        , 0.        ]])
np.sum(cv.GaussianBlur(img, (5, 5), 0))
# 1.0

https://docs.opencv.org/4.x/d4/d13/tutorial_py_filtering.html
https://stackoverflow.com/questions/61394826/how-do-i-get-to-show-gaussian-kernel-for-2d-opencv

この短所を3つとも同時解消するには、numpyで自作するしかない

import numpy as np
import math

def gaussian(px, sigma):
    if sigma == 0:
        n = px - 1
        coeffs = np.array([math.comb(n, k) for k in range(n + 1)], dtype=float)
        kernel = np.outer(coeffs, coeffs)
        kernel /= kernel.sum()
        return kernel
    else:
        x, y = np.mgrid[-half_px:half_px+1, -half_px:half_px+1]
        squared_distance = x**2 + y**2
        psf = np.exp(-squared_distance / (2 * sigma**2))
        psf /= psf.sum()
        return psf

gaussian(3, 0)
# array([[0.0625, 0.125 , 0.0625],
#        [0.125 , 0.25  , 0.125 ],
#        [0.0625, 0.125 , 0.0625]])

gaussian(3, 1)
# array([[0.07511361, 0.1238414 , 0.07511361],
#        [0.1238414 , 0.20417996, 0.1238414 ],
#        [0.07511361, 0.1238414 , 0.07511361]])

gaussian(4, 0)
# array([[0.015625, 0.046875, 0.046875, 0.015625],
#        [0.046875, 0.140625, 0.140625, 0.046875],
#        [0.046875, 0.140625, 0.140625, 0.046875],
#        [0.015625, 0.046875, 0.046875, 0.015625]])

gaussian(5, 0)
# array([[0.00390625, 0.015625  , 0.0234375 , 0.015625  , 0.00390625],
#        [0.015625  , 0.0625    , 0.09375   , 0.0625    , 0.015625  ],
#        [0.0234375 , 0.09375   , 0.140625  , 0.09375   , 0.0234375 ],
#        [0.015625  , 0.0625    , 0.09375   , 0.0625    , 0.015625  ],
#        [0.00390625, 0.015625  , 0.0234375 , 0.015625  , 0.00390625]])

np.sum(gaussian(3, 0))
# 1.0
np.sum(gaussian(4, 0))
# 1.0
np.sum(gaussian(5, 0))
# 1.0

ちなみにmath.comb(import math)も使いたくないならこれもnumpyで書けばよい

import numpy as np
def binom_coeff_row(n):
    k = np.arange(n+1, dtype=float)
    coeffs = np.empty(n+1, dtype=float)
    coeffs[0] = 1.0
    if n > 0:
        coeffs[1:] = np.cumprod((n - np.arange(n)) / (np.arange(1, n+1)))
    return coeffs

def gaussian(px, sigma):
    if sigma == 0:
        n = px - 1
        coeffs = binom_coeff_row(n)
        kernel = np.outer(coeffs, coeffs)
        kernel /= kernel.sum()
        return kernel
    else:
        x, y = np.mgrid[-half_px:half_px+1, -half_px:half_px+1]
        squared_distance = x**2 + y**2
        psf = np.exp(-squared_distance / (2 * sigma**2))
        psf /= psf.sum()
        return psf