Жесткое преобразование - Python - ускорение

У меня есть следующий вопрос о более быстром способе вычисления жесткого преобразования (да, я знаю, что могу просто использовать библиотеку, но мне нужно кодировать это самостоятельно).

Мне нужно вычислить x' и y' для каждого x, y в заданном изображении. Мое главное узкое место — скалярное произведение по всем координатам (интерполяция после этого не проблема). На данный момент я реализовал три варианта:

  1. понимание списка

    result = np.array([[np.dot(matrix, np.array([x, y, 1])) for x in xs] for y in ys])

  2. простая двойная for петля

    result = np.zeros((Y, X, 3))
for x in xs:
    for y in ys:
        result[y, x, :] = np.dot(matrix, np.array([x, y, 1]))

  3. np.ndenumerate

    result = np.zeros((Y, X, 3))
for (y, x), value in np.ndenumerate(image):
    result[y, x, :] = np.dot(matrix, np.array([x, y, 1]))

Самый быстрый способ в изображениях 512x512 — это понимание списка (примерно в 1,5 раза быстрее, чем np.ndenumerate, и в 1,1 раза быстрее, чем двойной цикл for).

Есть ли способ сделать это быстрее?


python numpy performance rigid-bodies
person Nefarin    schedule 24.11.2016    source источник

Ответы (2)


arrow_upward
3
arrow_downward

Заимствуя из this post, идею создания массивов 1D вместо сеток 2D или 3D и использования их на лету broadcasted для повышения эффективности использования памяти и, таким образом, для повышения производительности, вот подход -

out = ys[:,None,None]*matrix[:,1] + xs[:,None]*matrix[:,0] + matrix[:,2]

Если вы покрываете все индексы xs и ys для изображений размером 512x512, мы бы создали их с np.arange, например:

ys = np.arange(512)
xs = np.arange(512)

Тест во время выполнения

Определения функций -

def original_listcomp_app(matrix, X, Y): # Original soln with list-compr. 
    ys = np.arange(Y)
    xs = np.arange(X)
    out = np.array([[np.dot(matrix, np.array([x, y, 1])) for x in xs] \
                                                           for y in ys])
    return out    

def indices_app(matrix, X, Y):        # @Eric's soln  
    coords_ext = np.empty((Y, X, 3))
    coords_ext[...,[1,0]] = np.rollaxis(np.indices((Y, X)), 0, start=3)
    coords_ext[...,2] = 1    
    result = np.matmul(coords_ext,matrix.T)
    return result

def broadcasting_app(matrix, X, Y):  # Broadcasting based
    ys = np.arange(Y)
    xs = np.arange(X)
    out = ys[:,None,None]*matrix[:,1] + xs[:,None]*matrix[:,0] + matrix[:,2]
    return out

Сроки и проверка -

In [242]: # Inputs
     ...: matrix = np.random.rand(3,3)
     ...: X,Y = 512, 512
     ...: 

In [243]: out0 = original_listcomp_app(matrix, X, Y)
     ...: out1 = indices_app(matrix, X, Y)
     ...: out2 = broadcasting_app(matrix, X, Y)
     ...: 

In [244]: np.allclose(out0, out1)
Out[244]: True

In [245]: np.allclose(out0, out2)
Out[245]: True

In [253]: %timeit original_listcomp_app(matrix, X, Y)
1 loops, best of 3: 1.32 s per loop

In [254]: %timeit indices_app(matrix, X, Y)
100 loops, best of 3: 16.1 ms per loop

In [255]: %timeit broadcasting_app(matrix, X, Y)
100 loops, best of 3: 9.64 ms per loop
person Divakar    schedule 24.11.2016
comment
Хорошо, что это можно сделать с помощью вещания. Я часто писал здесь xs[None,:,None] для согласованности с ys[:,None,None], даже если в этом нет необходимости. Рад, что мое решение оказалось по крайней мере быстрее, чем составление списка - спасибо за профилирование! - person Eric; 24.11.2016
comment
Спасибо, оба решения работают достаточно быстро — это даже быстрее, чем встроенное преобразование scikit-image (но все же немного медленнее, чем scipy). - person Nefarin; 25.11.2016
comment
@Nefarin Вау, не знал, что у Scipy есть что-то для этой задачи. Ницца! - person Divakar; 25.11.2016

arrow_upward
3
arrow_downward

Вы можете использовать np.indices и np.rollaxis для создания трехмерного массива, где coords[i, j] == [i, j]. Здесь координаты нужно переключить

Затем все, что вам нужно сделать, это добавить 1, которое вы просите, и использовать @

coords_ext = np.empty((Y, X, 3))
coords_ext[...,[1,0]] = np.rollaxis(np.indices((Y, X)), 0, start=3)
coords_ext[...,2] = 1

# convert to column vectors and back for matmul broadcasting
result = (matrix @ coords_ext[...,None])[...,0]

# or alternatively, work with row vectors and do it in the other order
result = coords_ext @ matrix.T
person Eric    schedule 24.11.2016
comment
Пришлось внести несколько изменений в ваш опубликованный код, чтобы сопоставить выходные значения с подходом OP с двойным циклом, необходимым для проверки значений в ходе теста времени выполнения. Надеюсь, они выглядят нормально. - person Divakar; 24.11.2016
comment
@Divakar: Как выглядят мои исправления? - person Eric; 24.11.2016
comment
Потрясающе, спасибо! Обновлены профилирующие тесты. Некоторое улучшение таймингов с этими исправлениями! - person Divakar; 24.11.2016