Cách nhanh hơn để tính tổng khác biệt bình phương giữa một hình ảnh (M, N) và một mẫu (3, 3) để khớp mẫu?

Question

Tôi đang triển khai thuật toán cho Tổng hợp kết cấu như được nêu here. Đối với điều này, tôi cần tính tổng của các khác biệt bình phương, một số liệu để ước tính lỗi giữa template và các vị trí khác nhau trên image. Tôi có triển khai làm việc chậm tại chỗ như sau:

total_weight = valid_mask.sum() 
for i in xrange(input_image.shape[0]): 
    for j in xrange(input_image.shape[1]): 
     sample = image[i:i + window, j:j + window] 
     dist = (template - sample) ** 2 
     ssd[i, j] = (dist * valid_mask).sum()/total_weight 

Ở đây, total_weight chỉ để chuẩn hóa. Một số pixel có cường độ không xác định, vì vậy tôi sử dụng valid_mask để che khuất chúng. Vòng lặp lồng nhau này nằm bên trong của 2 vòng, vì vậy đó là 4 vòng lồng nhau rõ ràng là một kẻ giết người hiệu suất!

Có cách nào tôi có thể làm cho nó nhanh hơn trong NumPy hoặc Python, thay thế cho vòng lặp lồng nhau này không? Vectorization là có thể? Tôi cần phải làm việc trên (3, 3) một phần của số image bằng (3, 3) của số template.

Tôi sau đó sẽ thực hiện điều này trong Cython, vì vậy nhanh hơn tôi có thể làm cho nó hoạt động chỉ bằng cách sử dụng NumPy, tốt hơn.

Bạn có thể tìm mã hoàn chỉnh here. Dòng 62 - 67 được trích dẫn ở đây.

Cảm ơn,
Chintak

Answer 1

Điều này về cơ bản là cải tiến so với câu trả lời của Warren Weckesser. Con đường để đi là rõ ràng với một cửa sổ đa chiều cửa sổ của mảng ban đầu, nhưng bạn muốn giữ cho rằng xem từ kích hoạt một bản sao. Nếu bạn mở rộng sum((a-b)**2), bạn có thể biến nó thành sum(a**2) + sum(b**2) - 2*sum(a*b) và thao tác nhân này sau đó bạn có thể thực hiện với toán tử đại số tuyến tính, với cải thiện đáng kể về hiệu suất và bộ nhớ sử dụng:

def sumsqdiff3(input_image, template): 
    window_size = template.shape 
    y = as_strided(input_image, 
        shape=(input_image.shape[0] - window_size[0] + 1, 
          input_image.shape[1] - window_size[1] + 1,) + 
          window_size, 
        strides=input_image.strides * 2) 
    ssd = np.einsum('ijkl,kl->ij', y, template) 
    ssd *= - 2 
    ssd += np.einsum('ijkl, ijkl->ij', y, y) 
    ssd += np.einsum('ij, ij', template, template) 

    return ssd 

In [288]: img = np.random.rand(500, 500) 

In [289]: template = np.random.rand(3, 3) 

In [290]: %timeit a = sumsqdiff2(img, template) # Warren's function 
10 loops, best of 3: 59.4 ms per loop 

In [291]: %timeit b = sumsqdiff3(img, template) 
100 loops, best of 3: 18.2 ms per loop 

In [292]: np.allclose(a, b) 
Out[292]: True 

Tôi đã để nguyên thông số valid_mask vì tôi không hiểu hoàn toàn cách bạn sử dụng nó. Về nguyên tắc, chỉ cần zeroing các giá trị tương ứng trong template và/hoặc input_image cũng nên thực hiện thủ thuật tương tự.

Answer 2

Tôi nghĩ bạn đã làm một công việc tốt thực hiện thuật toán của bạn. Vectorization là một tùy chọn nhưng tôi khuyến khích bạn sử dụng trình biên dịch tối ưu hóa Numba để biên dịch cú pháp Python thành Mã máy. Hiệu ứng Numba khá ấn tượng. Numba vs. Cython: Take 2 là phần giới thiệu ngắn gọn về Numba và so sánh hiệu suất.

Answer 3

Có thể đáng giá khi bạn kiểm tra cách nó hoạt động nếu bạn sắp xếp lại thuật toán để thực hiện tính toán hàng khôn ngoan. Ý tưởng là bạn có thể đang sử dụng bộ nhớ cache CPU tốt hơn nếu đọc bộ nhớ liên tục.

Pseudo code:

for template_row in template: 
    for row in image: 
    for col in image: 
     # find distance template_row to sample_row 
     # add sum to ssd[row - template_row, col] 

mã thực tế (sau khi Warren):

def sumsqdiffr(input_image, template, valid_mask=None): 
    if valid_mask is None: 
     valid_mask = np.ones_like(template) 
    total_weight = valid_mask.sum() 
    window_size = template.shape 
    ssd = np.zeros((input_image.shape[0] - window_size[0] + 1, 
        input_image.shape[1] - window_size[1] + 1)) 

    for tr in xrange(template.shape[0]): 
     for i in xrange(tr, ssd.shape[0] + tr): 
      for j in xrange(ssd.shape[1]): 
       sample = input_image[i, j:j + window_size[1]] 
       dist = (template[tr] - sample) ** 2 
       ssd[i - tr, j] += (dist * valid_mask[tr]).sum() 
    return ssd 

Đó là chậm hơn hai lần hơn việc thực hiện ban đầu.

(Nếu ai đó muốn soi sáng cho tôi dù toàn bộ ý tưởng là sai hoặc những gì đang gây ra điều này, tôi muốn được vui mừng khi đạt được một số hiểu biết trong nó)

Answer 4

Bạn có thể làm một số điều tuyệt vời với các as_strided chức năng kết hợp với phát sóng của numpy. Dưới đây là hai phiên bản chức năng của bạn:

import numpy as np 
from numpy.lib.stride_tricks import as_strided 


def sumsqdiff(input_image, template, valid_mask=None): 
    if valid_mask is None: 
     valid_mask = np.ones_like(template) 
    total_weight = valid_mask.sum() 
    window_size = template.shape 
    ssd = np.empty((input_image.shape[0] - window_size[0] + 1, 
        input_image.shape[1] - window_size[1] + 1)) 
    for i in xrange(ssd.shape[0]): 
     for j in xrange(ssd.shape[1]): 
      sample = input_image[i:i + window_size[0], j:j + window_size[1]] 
      dist = (template - sample) ** 2 
      ssd[i, j] = (dist * valid_mask).sum() 
    return ssd 


def sumsqdiff2(input_image, template, valid_mask=None): 
    if valid_mask is None: 
     valid_mask = np.ones_like(template) 
    total_weight = valid_mask.sum() 
    window_size = template.shape 

    # Create a 4-D array y, such that y[i,j,:,:] is the 2-D window 
    #  input_image[i:i+window_size[0], j:j+window_size[1]] 
    y = as_strided(input_image, 
        shape=(input_image.shape[0] - window_size[0] + 1, 
          input_image.shape[1] - window_size[1] + 1,) + 
          window_size, 
        strides=input_image.strides * 2) 

    # Compute the sum of squared differences using broadcasting. 
    ssd = ((y - template) ** 2 * valid_mask).sum(axis=-1).sum(axis=-1) 
    return ssd 

Đây là phiên ipython để so sánh chúng.

Mẫu mà tôi sẽ sử dụng cho bản demo:

In [72]: template 
Out[72]: 
array([[-1, 1, -1], 
     [ 1, 2, 1], 
     [-1, 1, -1]]) 

Một đầu vào nhỏ vì vậy chúng tôi có thể kiểm tra kết quả:

In [73]: x 
Out[73]: 
array([[ 0., 1., 2., 3., 4., 5., 6.], 
     [ 7., 8., 9., 10., 11., 12., 13.], 
     [ 14., 15., 16., 17., 18., 19., 20.], 
     [ 21., 22., 23., 24., 25., 26., 27.], 
     [ 28., 29., 30., 31., 32., 33., 34.]]) 

Áp dụng hai chức năng để x và kiểm tra xem chúng tôi nhận cùng một kết quả:

In [74]: sumsqdiff(x, template) 
Out[74]: 
array([[ 856., 1005., 1172., 1357., 1560.], 
     [ 2277., 2552., 2845., 3156., 3485.], 
     [ 4580., 4981., 5400., 5837., 6292.]]) 

In [75]: sumsqdiff2(x, template) 
Out[75]: 
array([[ 856., 1005., 1172., 1357., 1560.], 
     [ 2277., 2552., 2845., 3156., 3485.], 
     [ 4580., 4981., 5400., 5837., 6292.]]) 

Bây giờ tạo một đầu vào lớn hơn "hình ảnh":

In [76]: z = np.random.randn(500, 500) 

và kiểm tra việc thực hiện:

In [77]: %timeit sumsqdiff(z, template) 
1 loops, best of 3: 3.55 s per loop 

In [78]: %timeit sumsqdiff2(z, template) 
10 loops, best of 3: 33 ms per loop 

Không quá tồi tàn. :)

Hai nhược điểm:

• Việc tính toán trong sumsqdiff2 sẽ tạo ra một mảng tạm rằng, đối với một mẫu 3x3, sẽ là 9 lần so với kích thước của input_image. (Nói chung nó sẽ là template.size lần kích thước của input_image.)
• Những "thủ đoạn sải chân" này sẽ không giúp bạn khi bạn Cythonize mã. Khi chuyển đổi sang Cython, bạn thường kết thúc việc đưa trở lại các vòng lặp mà bạn đã loại bỏ khi vector hóa với vón cục.