Как я могу эффективно обработать пустой массив в блоках, похожих на функцию Matlab blkproc (blockproc) - PullRequest
25 голосов
/ 22 февраля 2011

Я ищу хороший подход для эффективного разделения изображения на маленькие области, обработки каждой области отдельно, а затем повторной сборки результатов каждого процесса в одно обработанное изображение.Для этого у Matlab был инструмент под названием blkproc (в новых версиях Matlab вместо blockproc).

В идеальном мире функция или класс также поддерживают перекрытие между делениями в матрице ввода.В справке Matlab blkproc определяется как:

B = blkproc (A, [mn], [mborder nborder], fun, ...)

  • A - это ваша матрица ввода,
  • [mn] - размер блока.
  • [mborder, nborder] - размер области границы (необязательно).
  • .функция, применяемая к каждому блоку

Я объединил подход, но он кажется мне неуклюжим, и я уверен, что есть гораздо лучший способ.На риск моего собственного смущения, вот мой код:


import numpy as np

def segmented_process(M, blk_size=(16,16), overlap=(0,0), fun=None):
    rows = []
    for i in range(0, M.shape[0], blk_size[0]):
        cols = []
        for j in range(0, M.shape[1], blk_size[1]):
            cols.append(fun(M[i:i+blk_size[0], j:j+blk_size[1]]))
        rows.append(np.concatenate(cols, axis=1))
    return np.concatenate(rows, axis=0)

R = np.random.rand(128,128)
passthrough = lambda(x):x
Rprime = segmented_process(R, blk_size=(16,16), 
                           overlap=(0,0), 
                           fun=passthrough)

np.all(R==Rprime)

Ответы [ 6 ]

21 голосов
/ 22 февраля 2011

Вот несколько примеров другого (безциклового) способа работы с блоками:

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

A= np.arange(36).reshape(6, 6)
print A
#[[ 0  1  2  3  4  5]
# [ 6  7  8  9 10 11]
# ...
# [30 31 32 33 34 35]]

# 2x2 block view
B= ast(A, shape= (3, 3, 2, 2), strides= (48, 8, 24, 4))
print B[1, 1]
#[[14 15]
# [20 21]]

# for preserving original shape
B[:, :]= np.dot(B[:, :], np.array([[0, 1], [1, 0]]))
print A
#[[ 1  0  3  2  5  4]
# [ 7  6  9  8 11 10]
# ...
# [31 30 33 32 35 34]]
print B[1, 1]
#[[15 14]
# [21 20]]

# for reducing shape, processing in 3D is enough
C= B.reshape(3, 3, -1)
print C.sum(-1)
#[[ 14  22  30]
# [ 62  70  78]
# [110 118 126]]

Так что попытка просто скопировать функциональность matlab в numpy не всегда является лучшим способомспособ продолжить.Иногда требуется нестандартное мышление.

Предостережение :
В общем случае реализации, основанные на уловках шага , могут (но в этом нет необходимости) понести некоторые штрафы за производительность.Так что будьте готовы всеми способами измерить вашу производительность.В любом случае целесообразно сначала проверить, все ли готовые функции (или достаточно похожие, чтобы их можно было легко адаптировать) были реализованы в numpy или scipy.

Обновление :
Обратите внимание, что с strides здесь не связано действительное magic, поэтому я предоставлю простую функцию для получения block_view любого подходящего 2D numpy массива.Итак, поехали:

from numpy.lib.stride_tricks import as_strided as ast

def block_view(A, block= (3, 3)):
    """Provide a 2D block view to 2D array. No error checking made.
    Therefore meaningful (as implemented) only for blocks strictly
    compatible with the shape of A."""
    # simple shape and strides computations may seem at first strange
    # unless one is able to recognize the 'tuple additions' involved ;-)
    shape= (A.shape[0]/ block[0], A.shape[1]/ block[1])+ block
    strides= (block[0]* A.strides[0], block[1]* A.strides[1])+ A.strides
    return ast(A, shape= shape, strides= strides)

if __name__ == '__main__':
    from numpy import arange
    A= arange(144).reshape(12, 12)
    print block_view(A)[0, 0]
    #[[ 0  1  2]
    # [12 13 14]
    # [24 25 26]]
    print block_view(A, (2, 6))[0, 0]
    #[[ 0  1  2  3  4  5]
    # [12 13 14 15 16 17]]
    print block_view(A, (3, 12))[0, 0]
    #[[ 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 35]]
11 голосов
/ 22 февраля 2011

Процесс по частям / просмотрам.Конкатенация очень дорогая.

for x in xrange(0, 160, 16):
    for y in xrange(0, 160, 16):
        view = A[x:x+16, y:y+16]
        view[:,:] = fun(view)
7 голосов
/ 02 марта 2011

Я взял оба входа, а также мой оригинальный подход и сравнил результаты. Как правильно указывает @eat, результаты зависят от характера ваших входных данных. Удивительно, но сцепление превосходит обработку представления в нескольких случаях. У каждого метода есть приятное место. Вот мой контрольный код:

import numpy as np
from itertools import product

def segment_and_concatenate(M, fun=None, blk_size=(16,16), overlap=(0,0)):
    # truncate M to a multiple of blk_size
    M = M[:M.shape[0]-M.shape[0]%blk_size[0], 
          :M.shape[1]-M.shape[1]%blk_size[1]]
    rows = []
    for i in range(0, M.shape[0], blk_size[0]):
        cols = []
        for j in range(0, M.shape[1], blk_size[1]):
            max_ndx = (min(i+blk_size[0], M.shape[0]),
                       min(j+blk_size[1], M.shape[1]))
            cols.append(fun(M[i:max_ndx[0], j:max_ndx[1]]))
        rows.append(np.concatenate(cols, axis=1))
    return np.concatenate(rows, axis=0)


from numpy.lib.stride_tricks import as_strided
def block_view(A, block= (3, 3)):
    """Provide a 2D block view to 2D array. No error checking made.
    Therefore meaningful (as implemented) only for blocks strictly
    compatible with the shape of A."""
    # simple shape and strides computations may seem at first strange
    # unless one is able to recognize the 'tuple additions' involved ;-)
    shape= (A.shape[0]/ block[0], A.shape[1]/ block[1])+ block
    strides= (block[0]* A.strides[0], block[1]* A.strides[1])+ A.strides
    return as_strided(A, shape= shape, strides= strides)

def segmented_stride(M, fun, blk_size=(3,3), overlap=(0,0)):
    # This is some complex function of blk_size and M.shape
    stride = blk_size
    output = np.zeros(M.shape)

    B = block_view(M, block=blk_size)
    O = block_view(output, block=blk_size)

    for b,o in zip(B, O):
        o[:,:] = fun(b);

    return output

def view_process(M, fun=None, blk_size=(16,16), overlap=None):
    # truncate M to a multiple of blk_size
    from itertools import product
    output = np.zeros(M.shape)

    dz = np.asarray(blk_size)
    shape = M.shape - (np.mod(np.asarray(M.shape), 
                          blk_size))
    for indices in product(*[range(0, stop, step) 
                        for stop,step in zip(shape, blk_size)]):
        # Don't overrun the end of the array.
        #max_ndx = np.min((np.asarray(indices) + dz, M.shape), axis=0)
        #slices = [slice(s, s + f, None) for s,f in zip(indices, dz)]
        output[indices[0]:indices[0]+dz[0], 
               indices[1]:indices[1]+dz[1]][:,:] = fun(M[indices[0]:indices[0]+dz[0], 
               indices[1]:indices[1]+dz[1]])

    return output

if __name__ == "__main__":
    R = np.random.rand(128,128)
    squareit = lambda(x):x*2

    from timeit import timeit
    t ={}
    kn = np.array(list(product((8,16,64,128), 
                               (128, 512, 2048, 4096))  ) )

    methods = ("segment_and_concatenate", 
               "view_process", 
               "segmented_stride")    
    t = np.zeros((kn.shape[0], len(methods)))

    for i, (k, N) in enumerate(kn):
        for j, method in enumerate(methods):
            t[i,j] = timeit("""Rprime = %s(R, blk_size=(%d,%d), 
                          overlap = (0,0), 
                          fun = squareit)""" % (method, k, k),
                   setup="""
from segmented_processing import %s
import numpy as np
R = np.random.rand(%d,%d)
squareit = lambda(x):x**2""" % (method, N, N),
number=5
)
        print "k =", k, "N =", N #, "time:", t[i]
        print ("    Speed up (view vs. concat, stride vs. concat): %0.4f, %0.4f" % (
                       t[i][0]/t[i][1], 
                       t[i][0]/t[i][2]))

А вот и результаты:

A comparison of three methods for processing large matrices in a block-wise fashion in numpy Обратите внимание, что метод сегментированного шага выигрывает в 3-4 раза для небольших размеров блока. Только при больших размерах блока (128 x 128) и очень больших матрицах (2048 x 2048 и более) подход к обработке представления выигрывает, и тогда только с небольшим процентом. По результатам конкурса, похоже, что @eat получает галочку! Спасибо вам обоим за хорошие примеры!

3 голосов
/ 09 декабря 2014

Еще позже в игре. Швейцарский пакет обработки изображений под названием Bob доступен по адресу: https://www.idiap.ch/software/bob/docs/releases/last/sphinx/html/index.html Он имеет команду python bob.ip.block, описанную по адресу: https://www.idiap.ch/software/bob/docs/releases/last/sphinx/html/ip/generated/bob.ip.block.html#bob.ip.block который, кажется, делает все, что делает команда Matlab 'blockproc'. Я не проверял это.

Существует также интересная команда bob.ip.DCTFeatures, которая включает в себя возможности «блока» для извлечения или изменения коэффициентов DCT изображения.

3 голосов
/ 17 января 2014

Немного опоздал в игру, но это сделало бы перекрывающиеся блоки. Я не делал этого здесь, но вы могли бы легко приспособиться к размерам шагов для сдвига окна, я думаю:

from numpy.lib.stride_tricks import as_strided
def rolling_block(A, block=(3, 3)):
    shape = (A.shape[0] - block[0] + 1, A.shape[1] - block[1] + 1) + block
    strides = (A.strides[0], A.strides[1]) + A.strides
    return as_strided(A, shape=shape, strides=strides)
0 голосов
/ 05 июля 2016

Я нашел этот учебник - окончательный исходный код обеспечивает именно желаемую функциональность! Это даже должно работать для любой размерности (я не проверял это) http://www.johnvinyard.com/blog/?p=268

Хотя опция "flatten" в самом конце исходного кода кажется немного глючной. Тем не менее, очень хороший кусок программного обеспечения!

...