[PYTHON] 코드 효율성 향상 : 슬라이딩 윈도우의 표준 편차

코드 효율성 향상 : 슬라이딩 윈도우의 표준 편차

나는 픽셀의 이웃에 위치한 픽셀의 표준 편차를 이미지의 각 픽셀에 대해 계산하는 함수를 개선하려고 노력 중이다. 내 함수는 2 개의 임베디드 루프를 사용하여 행렬을 가로 지르며, 이는 내 프로그램의 병목입니다. 나는 numpy의 덕택으로 루프를 없애고 그것을 향상시킬 수있는 방법이있을 것이라고 생각하지만 진행 방법을 모른다. 어떤 조언을 환영합니다!

문안 인사

def sliding_std_dev(image_original,radius=5) :
    height, width = image_original.shape
    result = np.zeros_like(image_original) # initialize the output matrix
    hgt = range(radius,height-radius)
    wdt = range(radius,width-radius)
    for i in hgt:
        for j in wdt:
            result[i,j] = np.std(image_original[i-radius:i+radius,j-radius:j+radius])
    return result

해결법

1. ==============================

   1.시원한 속임수 : 윈도우의 값의 합과 제곱 값의 합계 만 주어진 표준 편차를 계산할 수 있습니다.

   시원한 속임수 : 윈도우의 값의 합과 제곱 값의 합계 만 주어진 표준 편차를 계산할 수 있습니다.

   따라서 데이터에 대한 일정한 필터를 사용하여 표준 편차를 매우 빠르게 계산할 수 있습니다.

   from scipy.ndimage.filters import uniform_filter
   
   def window_stdev(arr, radius):
       c1 = uniform_filter(arr, radius*2, mode='constant', origin=-radius)
       c2 = uniform_filter(arr*arr, radius*2, mode='constant', origin=-radius)
       return ((c2 - c1*c1)**.5)[:-radius*2+1,:-radius*2+1]
   

   이것은 원래 기능보다 훨씬 빠릅니다. 1024x1024 배열과 반경이 20 인 경우 이전 함수는 34.11 초 걸리고 새 함수는 0.11 초가 걸리고 속도는 300 배가됩니다.

   이것은 어떻게 수학적으로 작용합니까? 각 창에 대한 sqrt (평균 (x ^ 2) - 평균 (x) ^ 2)을 계산합니다. 표준 편차 sqrt (mean ((x - mean (x)) ^ 2))에서 다음과 같이이 양을 유도 할 수 있습니다.

   E를 기대 연산자 (기본적으로 mean ())라고하고, X는 데이터의 임의 변수입니다. 그때:

   E [(X-E [X]) ^ 2] = E [X ^ 2 - 2X * E [X] + E [X] ^ 2] E (X ^ 2) - E [2X * E [X]] + E [E [X] ^ 2] (기대 연산자의 선형성에 의함) = E [X ^ 2] - 2E [X] * E [X] + E [X] ^ 2 (다시 선형성으로 E [X]는 상수 임) = E [X ^ 2] - E [X] ^ 2

   이 기법을 사용하여 계산 된 양은 수학적으로 표준 편차와 동일하다는 것을 증명합니다.

2. ==============================

   2.이미지 처리에서 이런 종류의 작업을 수행하는 데 가장 자주 사용되는 방법은 1984 년이 신문에서 소개 한 아이디어 인 합계 면적 테이블을 사용하는 것입니다. 아이디어는 창 위에 추가하여 수량을 계산하고 창을 이동하는 것입니다. 한 픽셀 씩 오른쪽으로 이동하면 새 창에 모든 항목을 추가 할 필요가 없으므로 합계에서 가장 왼쪽 열을 뺀 다음 가장 오른쪽에있는 새 열을 추가하면됩니다. 따라서 배열에서 두 차원 모두에 대해 누적 합계 배열을 만들면 두 개의 합계와 빼기가있는 창에 대해 합계를 얻을 수 있습니다. 배열과 그 사각형에 대한 합계 면적 테이블을 유지하는 경우, 그 둘로부터 분산을 얻는 것은 매우 쉽습니다. 구현 방법은 다음과 같습니다.

   이미지 처리에서 이런 종류의 작업을 수행하는 데 가장 자주 사용되는 방법은 1984 년이 신문에서 소개 한 아이디어 인 합계 면적 테이블을 사용하는 것입니다. 아이디어는 창 위에 추가하여 수량을 계산하고 창을 이동하는 것입니다. 한 픽셀 씩 오른쪽으로 이동하면 새 창에 모든 항목을 추가 할 필요가 없으므로 합계에서 가장 왼쪽 열을 뺀 다음 가장 오른쪽에있는 새 열을 추가하면됩니다. 따라서 배열에서 두 차원 모두에 대해 누적 합계 배열을 만들면 두 개의 합계와 빼기가있는 창에 대해 합계를 얻을 수 있습니다. 배열과 그 사각형에 대한 합계 면적 테이블을 유지하는 경우, 그 둘로부터 분산을 얻는 것은 매우 쉽습니다. 구현 방법은 다음과 같습니다.

   def windowed_sum(a, win):
       table = np.cumsum(np.cumsum(a, axis=0), axis=1)
       win_sum = np.empty(tuple(np.subtract(a.shape, win-1)))
       win_sum[0,0] = table[win-1, win-1]
       win_sum[0, 1:] = table[win-1, win:] - table[win-1, :-win]
       win_sum[1:, 0] = table[win:, win-1] - table[:-win, win-1]
       win_sum[1:, 1:] = (table[win:, win:] + table[:-win, :-win] -
                          table[win:, :-win] - table[:-win, win:])
       return win_sum
   
   def windowed_var(a, win):
       win_a = windowed_sum(a, win)
       win_a2 = windowed_sum(a*a, win)
       return (win_a2 - win_a * win_a / win/ win) / win / win
   

   이 작동 방식을 확인하려면 다음을 수행하십시오.

   >>> a = np.arange(25).reshape(5,5)
   >>> windowed_var(a, 3)
   array([[ 17.33333333,  17.33333333,  17.33333333],
          [ 17.33333333,  17.33333333,  17.33333333],
          [ 17.33333333,  17.33333333,  17.33333333]])
   >>> np.var(a[:3, :3])
   17.333333333333332
   >>> np.var(a[-3:, -3:])
   17.333333333333332
   

   이는 회선 기반 방법보다 몇 가지 노치를 빠르게 실행해야합니다.

3. ==============================

   3.첫째,이 작업을 수행하는 방법은 여러 가지가 있습니다.

   첫째,이 작업을 수행하는 방법은 여러 가지가 있습니다.

   속도면에서 가장 효율적인 것은 아니지만 scipy.ndimage.generic_filter를 사용하면 움직이는 창 위에 임의의 파이썬 함수를 쉽게 적용 할 수 있습니다.

   빠른 예를 들면 다음과 같습니다.

   result = scipy.ndimage.generic_filter(data, np.std, size=2*radius)
   

   경계 조건은 모드 kwarg로 제어 할 수 있습니다.

   이것을 수행하는 또 다른 방법은 효과적으로 다양한 움직이는 창을보기 위해 다양한 스트라이핑 트릭을 사용하고 마지막 축을 따라 np.std를 적용하는 것입니다. (참고 : 이것은 이전 답변 중 하나에서 가져온 것입니다 : https://stackoverflow.com/a/4947453/325565)

   def strided_sliding_std_dev(data, radius=5):
       windowed = rolling_window(data, (2*radius, 2*radius))
       shape = windowed.shape
       windowed = windowed.reshape(shape[0], shape[1], -1)
       return windowed.std(axis=-1)
   
   def rolling_window(a, window):
       """Takes a numpy array *a* and a sequence of (or single) *window* lengths
       and returns a view of *a* that represents a moving window."""
       if not hasattr(window, '__iter__'):
           return rolling_window_lastaxis(a, window)
       for i, win in enumerate(window):
           if win > 1:
               a = a.swapaxes(i, -1)
               a = rolling_window_lastaxis(a, win)
               a = a.swapaxes(-2, i)
       return a
   
   def rolling_window_lastaxis(a, window):
       """Directly taken from Erik Rigtorp's post to numpy-discussion.
       <http://www.mail-archive.com/numpy-discussion@scipy.org/msg29450.html>"""
       if window < 1:
          raise ValueError, "`window` must be at least 1."
       if window > a.shape[-1]:
          raise ValueError, "`window` is too long."
       shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
       strides = a.strides + (a.strides[-1],)
       return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
   

   언뜻보기에 여기서 무슨 일이 벌어지고 있는지 이해하는 것은 약간 어렵습니다. 내 대답 중 하나를 연결하지 말고 설명을 다시 입력하고 싶지 않으므로 여기를 살펴보십시오. https://stackoverflow.com/a/4924433/325565 전에 "striding"트릭.

   우리가 타이밍을 반지름이 5 인 100x100 무작위 수레 배열과 비교하면 원래 또는 generic_filter 버전보다 ~ 10 배 빠릅니다. 그러나이 버전에서는 경계 조건에 유연성이 없습니다. (이것은 현재하고있는 것과 동일 합니다만, generic_filter 버전은 속도면에서 많은 유연성을 제공합니다.)

   # Your original function with nested loops
   In [21]: %timeit sliding_std_dev(data)
   1 loops, best of 3: 237 ms per loop
   
   # Using scipy.ndimage.generic_filter
   In [22]: %timeit ndimage_std_dev(data)
   1 loops, best of 3: 244 ms per loop
   
   # The "stride-tricks" version above
   In [23]: %timeit strided_sliding_std_dev(data)
   100 loops, best of 3: 15.4 ms per loop
   
   # Ophion's version that uses `np.take`
   In [24]: %timeit new_std_dev(data)
   100 loops, best of 3: 19.3 ms per loop
   

   "stride-tricks"버전의 단점은 "보통"번쩍이는 롤링 윈도우 트릭과 달리이 버전은 복사본을 만들고 원래 배열보다 훨씬 크기 때문입니다. 큰 배열에서 이것을 사용하면 메모리 문제가 발생합니다! (사이드 노트에서, 그것은 메모리 사용과 속도면에서 @ Ophion의 대답과 기본적으로 동일합니다. 같은 일을하는 다른 방법 일뿐입니다.)

4. ==============================

   4.먼저 인덱스를 얻은 다음 np.take를 사용하여 새 배열을 만들 수 있습니다.

   먼저 인덱스를 얻은 다음 np.take를 사용하여 새 배열을 만들 수 있습니다.

   def new_std_dev(image_original,radius=5):
       cols,rows=image_original.shape
   
       #First obtain the indices for the top left position
       diameter=np.arange(radius*2)
       x,y=np.meshgrid(diameter,diameter)
       index=np.ravel_multi_index((y,x),(cols,rows)).ravel()
   
       #Cast this in two dimesions and take the stdev
       index=index+np.arange(rows-radius*2)[:,None]+np.arange(cols-radius*2)[:,None,None]*(rows)
       data=np.std(np.take(image_original,index),-1)
   
       #Add the zeros back to the output array
       top=np.zeros((radius,rows-radius*2))
       sides=np.zeros((cols,radius))
   
       data=np.vstack((top,data,top))
       data=np.hstack((sides,data,sides))
       return data
   

   먼저 임의의 데이터를 생성하고 타이밍을 확인하십시오.

   a=np.random.rand(50,20)
   
   print np.allclose(new_std_dev(a),sliding_std_dev(a))
   True
   
   %timeit sliding_std_dev(a)
   100 loops, best of 3: 18 ms per loop
   
   %timeit new_std_dev(a)
   1000 loops, best of 3: 472 us per loop
   

   큰 배열의 경우 메모리가 충분하면 항상 더 빠릅니다.

   a=np.random.rand(200,200)
   
   print np.allclose(new_std_dev(a),sliding_std_dev(a))
   True
   
   %timeit sliding_std_dev(a)
   1 loops, best of 3: 1.58 s per loop
   
   %timeit new_std_dev(a)
   10 loops, best of 3: 52.3 ms per loop
   

   원래 함수는 아주 작은 배열 일수록 빠르며, hgt * wdt> 50 일 때 break even point처럼 보입니다. 뭔가 함수를 메모하는 사각형 프레임을 가져 와서 인덱스를 샘플링하지, 오른쪽 하단의 색인에있는 표준 dev에 배치됩니다. 이것은 의도적입니까?

from https://stackoverflow.com/questions/18419871/improving-code-efficiency-standard-deviation-on-sliding-windows by cc-by-sa and MIT license