Как сместить индексы выходной переменной (в ноль) во время присваивания векторизованным способом

Question

Мотивация: предположим, что у меня есть изображение RGB J, и я хочу применить преобразование T (например, поворот) к пикселям J. Я создам новое черное изображение K, с которым связаны его пиксели. к J с помощью K [x, y] = J [T [x, y]]. Теперь проблема в том, что T [x, y] должен быть внутри J, и если я хочу полностью захватить преобразованное изображение J, мне, возможно, придется иметь дело с некоторыми отрицательными значениями x или y или значениями, которые больше, чем размер J. Итак, сначала я должен определить размер K, а затем сместить пиксели K на соответствующий вектор, чтобы избежать отрицательных значений.

Теперь предположим, что я определил соответствующий вектор перевода. Я хочу сделать перевод координат, который отправляет (x, y) в (x + a, y + k).

Цель: Используя циклы for, я хочу сделать следующее:

for i in range(0,J.shape[0]):
    for j in range(0, J.shape[1]):
        K[i+a,j+b] = J[T[i,j]]

Как я могу сделать это векторизованным способом? Любая помощь приветствуется.

---

Edit:

img = face() # dummy RGB data

i,j = np.mgrid[:img.shape[0], :img.shape[1]] # 2d arrays each
i_min, i_max, j_min, j_max = func(*) # assume that these values have been found
i = i + i_min
j = j + j_min
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)
inew,jnew = np.linalg.inv(T) @ [i.ravel(), j.ravel()] # 1d arrays each

inew = np.floor(inew).astype(int)
jnew = np.floor(jnew).astype(int)

out = np.zeros((i_max - i_min, j_max - j_min, 3), dtype=img.dtype)

for i in inew:
    for j in jnew:
        out[i-i_min,j-j_min, :] = img[i,j,:]

Теперь я хочу отменить эффект смещения на i_min и j_min в массиве так же, как код, который я написал с помощью for-loop.

Answer 1

Наивная версия

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

import numpy as np
import matplotlib.pyplot as plt # for plotting the result
from scipy.misc import face # for dummy data
img = face() # dummy RGB data

# transform pixels by 45 degrees
i,j = np.mgrid[:img.shape[0], :img.shape[1]] # 2d arrays each
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)
inew,jnew = T @ [i.ravel(), j.ravel()] # 1d arrays each

# new coordinates now range into negatives, shift back into positives
# and the non-integer pixel indices will be normalized with floor
inew = np.floor(inew - inew.min()).astype(int)
jnew = np.floor(jnew - jnew.min()).astype(int)

# now the new coordinates are all non-negative, this defines the size of the output
out = np.zeros((inew.max() + 1, jnew.max() + 1, 3), dtype=img.dtype)

# fill the necessary indices of out with pixels from img
# reshape the indices to 2d for matching broadcast
inew = inew.reshape(img.shape[:-1])
jnew = jnew.reshape(img.shape[:-1])
out[inew, jnew, :] = img
# OR, alternative with 1d index arrays:
#out[inew, jnew, :] = img.reshape(-1, 3)

# check what we've done
plt.imshow(out)
plt.show()

Суть кода в том, что повернутые координаты пикселя сдвинуты обратно в положительные значения (это соответствует вашему [i+a, j+b] смещению), выделен новый нулевой массив, который будет соответствовать всем новым индексам, и индексирование применяется только на правой стороне ! Это то, что не соответствует вашему коду, но я считаю, что это то, что вы действительно хотите сделать: для каждого пикселя в исходном (неиндексированном) изображении мы устанавливаем его значение RGB в позиции new результирующего массив.

Как видите, на изображении много черных пикселей из-за того, что нецелые преобразованные координаты были округлены до floor. Это нехорошо, поэтому, если мы пойдем по этому пути, мы должны выполнить 2-мерную интерполяцию, чтобы избавиться от этих артефактов. Обратите внимание, что для этого требуется немного памяти и процессорного времени:

import numpy as np
import scipy.interpolate as interp
import matplotlib.pyplot as plt # for plotting the result
from scipy.misc import face # for dummy data
img = face() # dummy RGB data

# transform pixels by 45 degrees
i,j = np.mgrid[:img.shape[0], :img.shape[1]] # 2d arrays each
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)
inew,jnew = T @ [i.ravel(), j.ravel()] # 1d arrays each

# new coordinates now range into negatives, shift back into positives
# keep them non-integer for interpolation later
inew -= inew.min()
jnew -= jnew.min()
# (inew, jnew, img) contain the data from which the output should be interpolated


# now the new coordinates are all non-negative, this defines the size of the output
out = np.zeros((int(round(inew.max())) + 1, int(round(jnew.max())) + 1, 3), dtype=img.dtype)
i_interp,j_interp = np.mgrid[:out.shape[0], :out.shape[1]]

# interpolate for each channel
for channel in range(3):
    out[..., channel] = interp.griddata(np.array([inew.ravel(), jnew.ravel()]).T, img[..., channel].ravel(), (i_interp, j_interp), fill_value=0)

# check what we've done
plt.imshow(out)
plt.show()

По крайней мере, результат выглядит намного лучше:

scipy.ndimage: map_coordinates

Подход, который непосредственно соответствует тому, что вы имели в виду, может использовать scipy.ndimage.map_coordinates для выполнения интерполяции с использованием обратного преобразования . Это должно иметь лучшую производительность, чем предыдущая попытка с griddata, поскольку map_coordinates может использовать тот факт, что входные данные определены в сетке. Оказывается, что он действительно использует меньше памяти и гораздо меньше ресурсов процессора:

import numpy as np
import scipy.ndimage as ndi
import matplotlib.pyplot as plt # for plotting the result
from scipy.misc import face # for dummy data

img = face() # dummy RGB data
n,m = img.shape[:-1]

# transform pixels by 45 degrees
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)

# find out the extent of the transformed pixels from the four corners
inew_tmp,jnew_tmp = T @ [[0, 0, n-1, n-1], [0, m-1, 0, m-1]] # 1d arrays each
imin,imax,jmin,jmax = inew_tmp.min(),inew_tmp.max(),jnew_tmp.min(),jnew_tmp.max()
imin,imax,jmin,jmax = (int(round(val)) for val in (imin,imax,jmin,jmax))

# so the pixels of the original map inside [imin, imax] x [jmin, jmax]
# we need an image of size (imax - imin + 1, jmax - jmin + 1) to house this
out = np.zeros((imax - imin + 1, jmax - jmin + 1, 3), dtype=img.dtype)
# indices have to be shifted by [imin, imax]

# compute the corresponding (non-integer) coordinates on the domain for interpolation
inew,jnew = np.mgrid[:out.shape[0], :out.shape[1]]
i_back,j_back = np.linalg.inv(T) @ [inew.ravel() + imin, jnew.ravel() + jmin]

# perform 2d interpolation for each colour channel separately
for channel in range(3):
    out[inew, jnew, channel] = ndi.map_coordinates(img[..., channel], [i_back, j_back]).reshape(inew.shape)

# check what we've done
plt.imshow(out)
plt.show()

Результат все еще хорош:

scipy.ndimage: геометрическое преобразование

Наконец, я понял, что мы можем подняться на один уровень выше и использовать scipy.ndimage.geometric_transform напрямую. В случае с повернутым енотом это кажется медленнее, чем в ручной версии, использующей map_coordinates, но приводит к более чистому коду:

import numpy as np
import scipy.ndimage as ndi
import matplotlib.pyplot as plt # for plotting the result
from scipy.misc import face # for dummy data

img = face() # dummy RGB data
n,m = img.shape[:-1]

# transform pixels by 45 degrees
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)
Tinv = np.linalg.inv(T)

# find out the extent of the transformed pixels from the four corners
inew_tmp,jnew_tmp = T @ [[0, 0, n-1, n-1], [0, m-1, 0, m-1]] # 1d arrays each
imin,imax,jmin,jmax = inew_tmp.min(),inew_tmp.max(),jnew_tmp.min(),jnew_tmp.max()
imin,imax,jmin,jmax = (int(round(val)) for val in (imin,imax,jmin,jmax))

# so the pixels of the original map inside [imin, imax] x [jmin, jmax]
# we need an image of size (imax - imin + 1, jmax - jmin + 1) to house this

def transform_func(output_coords):
    """Inverse transform output coordinates back into input coordinates"""
    inew,jnew,channel = output_coords
    i,j = Tinv @ [inew + imin, jnew + jmin]
    return i,j,channel

out = ndi.geometric_transform(img, transform_func, output_shape = (imax - imin + 1, jmax - jmin + 1, 3))

# check what we've done
plt.imshow(out)
plt.show()

Результат:

Окончательное исправление: только NumPy

Меня в первую очередь интересовало качество изображения, поэтому все вышеперечисленные решения так или иначе используют интерполяцию. Как вы пояснили в комментариях, это не имеет для вас первостепенного значения. Если это так, мы можем изменить версию, используя map_coordinates, рассчитать приблизительные (округленное целое число) индексы и выполнить векторизованное присваивание:

import numpy as np
import matplotlib.pyplot as plt # for plotting the result
from scipy.misc import face # for dummy data

img = face() # dummy RGB data
n,m = img.shape[:-1]

# transform pixels by 45 degrees
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)

# find out the extent of the transformed pixels from the four corners
inew_tmp,jnew_tmp = T @ [[0, 0, n-1, n-1], [0, m-1, 0, m-1]] # 1d arrays each
imin,imax,jmin,jmax = inew_tmp.min(),inew_tmp.max(),jnew_tmp.min(),jnew_tmp.max()
imin,imax,jmin,jmax = (int(round(val)) for val in (imin,imax,jmin,jmax))

# so the pixels of the original map inside [imin, imax] x [jmin, jmax]
# we need an image of size (imax - imin + 1, jmax - jmin + 1) to house this
out = np.zeros((imax - imin + 1, jmax - jmin + 1, 3), dtype=img.dtype)

# compute the corresponding coordinates on the domain for matching
inew,jnew = np.mgrid[:out.shape[0], :out.shape[1]]
inew = inew.ravel() # 1d array, indices of output array
jnew = jnew.ravel() # 1d array, indices of output array
i_back,j_back = np.linalg.inv(T) @ [inew + imin, jnew + jmin]

# create a mask to grab only those rounded (i_back,j_back) indices which make sense
i_back = i_back.round().astype(int)
j_back = j_back.round().astype(int)
inds = (0 <= i_back) & (i_back < n) & (0 <= j_back) & (j_back < m)
# (i_back[inds], j_back[inds]) maps to (inew[inds], jnew[inds])
# the rest stays black

out[inew[inds], jnew[inds], :] = img[i_back[inds], j_back[inds], :]

# check what we've done
plt.imshow(out)
plt.show()

Результат, хотя и полный однопиксельных неточностей, выглядит достаточно хорошо:

Answer 2

Вы можете использовать функцию карты

for i in range(0,J.shape[0]):
    for j in range(0, J.shape[1]):
        K[i+a,j+b] = J[T[i,j]]

например, вы можете генерировать все кортежи индексов вашей матрицы

indexes = [ (i,j) for i in range(J.shape[0]) for j in range(J.shape[1]) ]

, а затем применить карту с лямбда-функцией

f = lambda coords:  J[T[coords[0],coords[1]]]
resp = list(map(f, indexes))

на данный момент содержит список всех применений f к индексам. Теперь вы должны преобразовать его в хорошую форму. для К

Так что здесь у вас есть две возможности, вы можете сделать список диапазонов размером K, а затем вы можете вернуть ноль, когда это необходимо внутри лямбда-функции

Старый ответ ...

Проблема здесь в том, что вы должны знать размер выходного изображения заранее. Таким образом, есть две возможности: либо вы вычисляете это, либо предполагаете, что оно не будет больше определенной оценки.

Так что, если вы его вычислите, путь зависит от преобразования, которое вы хотите применить. Например, транспонирование означает обмен по осям X и Y. Для поворота размер результата зависит от формы и угла.

So

если вы хотите, чтобы все было очень просто но не обязательно для памяти. Предположим, что ваше преобразование не будет выводить изображение больше, чем в три раза, максимум максимумов X и Y.

С этим вы легко справитесь со своими смещениями

, если ваше изображение NxM с N > M, холст для вашего преобразования будет 3*Nx3*N

теперь допустим, что выходное изображение будет центрировано на этом холсте. В этой ситуации вы должны вычислить смещения a и b, которые вы описали в своем вопросе

Центр преобразованного изображения вдоль вертикальной оси должен совпадать с центром исходного изображения.

if i=N/2 then a+i=3*N/2 это означает, что a=N

то же самое относится к горизонтальной оси и в этом случае

if j=M/2 then b+j=3*N/2 это означает, что b=(3*N - M)/2

Надеюсь понятно