소개

입력 신호에 대하여, 어느 필터 사이즈로 max/min을 취해 가는 처리가 필요한 경우가 있다. 제일의 유스 케이스는, 화상 처리에 있어서의 dilation(팽창)이나 erosion(수축) 처리이다. 개인적으로는, 옛날, 어느 폭을 가지는 non-maximum suppression 처리를 실시하기 위해서 구현한 적이 있다.
이 처리는, 나이브에 실장하면 신호 길이×필터 사이즈의 오더의 계산량이 필요하지만, 신호 길이의 오더의 계산량으로 억제할 수 있는 「van Herk/Gil-Werman 알고리즘」이라고 하는 것이 존재해 , 과연-이라는 느낌이므로 소개한다.

또한, 긴 신호에 대해 어느 필터 사이즈의 min/max를 계속 걸어 갔을 때에, 계산량이 필터 사이즈에 의존하지 않게 할 수 있다고 하는 것이므로, 1회의 min/max가 정수 오더라고 하는 것은 아니.

그건 그렇고,이 알고리즘은 문헌 1

van Herk/Gil-Werman 알고리즘

2
상기 사이트를 참고로 van Herk/Gil-Werman 알고리즘을 설명한다. 도면도 상기보다 인용하고 있지만, 첨자가 잘못되어 있는 부분을 수정하였다.

이 알고리즘에 대한 입력은
$I$: 입력 신호
$L$: 필터 크기(홀수)
된다. 출력은 입력 신호와 동일한 길이의 min/max 필터링 된 신호 $ I '$입니다.
이하에서는 max의 경우에서의 설명을 실시한다.

h tp // w w.ぇp 어쨌든. 이 m/g등 ys인가ぇー도 rp호우gy. HTML

이 알고리즘에서 입력 신호는 필터 크기 $ L $마다 처리됩니다. 여기서, 상술 한 바와 같이 제 1 처리 단위로 설명한다.
이 $L$ 마다의 처리에서는, 우선 윈도우의 중심으로부터 전후$L-1$까지 max를 취해 간 길이 $2L-1$의 배열 $A$를 준비한다. 이 때문에, 전후에 $(L-1)/2$의 마진이 필요하게 된다.
즉, 윈도우의 중심은 $A[L-1] = I[L-1]$이며,
$A[k] = max(I[k], A[k+1])$($k$ = $L-2, L-3, ..., 0$)
$A[k] = max(A[k-1], I[k])$ ($k$ = $L, L+1, ..., 2L-2$)
그리고 순차적으로 계산할 수 있다.



이 배열 $A$가 가능하면 출력 $I'$를 구하기 쉽고, $I'[J] = max(A[J - (L-1)/2], A[J + (L-1)/2])$가 된다. 이것은, 윈도우폭$L$내의 max를,$I[L-1]$의 전후로 분할해 행한 후, 통합하고 있을 것이다.

구현

파이썬에서 구현할 의미는 없지만 알고리즘 확인까지.

import numpy as np
import time


def main():
    L = 7
    sample_cnt = L * 50000 - 1
    I = np.random.random(sample_cnt)
    A = np.zeros(L * 2 - 1, dtype=np.float)
    N = (sample_cnt - L + 1) // L
    result1 = np.zeros(sample_cnt, dtype=np.float)
    result2 = np.zeros(sample_cnt, dtype=np.float)

    # van Herk/Gil-Wermanアルゴリズム
    start_time = time.time()

    for i in range(N):
        center = (i + 1) * L - 1
        start = center - L // 2
        A[L - 1] = I[center]

        for j in range(1, L):
            A[L - 1 - j] = max(A[L - j], I[center - j])
            A[L - 1 + j] = max(A[L - 2 + j], I[center + j])

        for j in range(L):
            result1[start + j] = max(A[j], A[j + L - 1])

    print("{} [s] elapsed".format(time.time() - start_time))

    # naive実装
    start_time = time.time()

    for i in range(L // 2, sample_cnt - L // 2):
        result2[i] = I[i - L // 2: i + L // 2 + 1].max()

    print("{} [s] elapsed".format(time.time() - start_time))


if __name__ == '__main__':
    main()

0.6796669960021973 [s] elapsed
0.9065239429473877 [s] elapsed

나이브 구현보다 빠릅니다.

오치↓

from scipy.ndimage import grey_dilation

# 略
start_time = time.time()
result3 = grey_dilation(I, size=(L,))
print("{} [s] elapsed".format(time.time() - start_time))

그렇다면,

0.009512186050415039 [s] elapsed

된다.



M. Herk, "A Fast Algorithm for Local Minimum and Maximum Filters on Rectangular and Octagonal Kernels,"in PRL, Vol. 13, No. 7, pp. 517-521, 1992.

J. Gil and M. Werman, "Computing 2-D Min, Median and Max Filters,"in TPAMI, Vol. 15, No. 5, pp. 504-507, 1993.  ↩