Заполнение вывода scipy affine_transform для отображения неперекрывающихся областей преобразованных изображений

Если у вас есть два похожих (или одинаковых) изображения и вы хотите их выровнять, вы можете сделать это, используя обе функции rotate и shift :

from scipy.ndimage import rotate, shift

Вам нужно найди сначала разность углов между двумя изображениями angle_to_rotate, при этом вы применить вращение в src:

angle_to_rotate = 25
rotated_src = rotate(src, angle_to_rotate , reshape=True, order=1, mode = "constant")

С помощью reshape=True вы избегаете потери информации из исходной матрицы src и дополняет результат, чтобы изображение можно было перевести вокруг индексов 0,0. Вы можете рассчитать этот перевод как есть (x*cos(angle),y*sin(angle), где x и y — размеры изображения, но, вероятно, это не будет иметь значения.

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

rot_translated_src = shift(rotated_src , [distance_x, distance_y])

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

Но вы можете сделать некоторые дополнения с помощью

np.pad(src, number, mode='constant')

Чтобы рассчитать distance_x и distance_y, вам нужно будет найти точку, которая служит вам ориентиром между rotated_src и пунктом назначения, а затем просто вычислить расстояние по осям x и y.

Резюме

1. Сделайте прокладку в src и dst
2. Найдите угловое расстояние между ними.
3. Повернутьsrc с scipy.ndimage.rotate с использованием reshape=True
4. Найдите горизонтальное и вертикальное расстояниеdistance_x, distance_y между повернутым изображением и dst
5. Перевести ваш 'rotated_src' с scipy.ndimage.shift

Код

from scipy.ndimage import rotate, shift
import matplotlib.pyplot as plt
import numpy as np

Сначала мы делаем целевое изображение:

# make and plot dest
dst = np.ones([40,20])
dst = np.pad(dst,10)
dst[17,[14,24]]=4
dst[27,14:25]=4
dst[26,[14,25]]=4
rotated_dst = rotate(dst, 20, order=1)

plt.imshow(dst) # plot it
plt.imshow(rotated_dst)
plt.show()

Делаем Исходный образ:

# make_src image and plot it
src = np.zeros([40,20])
src = np.pad(src,10)
src[0:20,0:20]=1
src[7,[4,14]]=4
src[17,4:15]=4
src[16,[4,15]]=4
plt.imshow(src)
plt.show()

Затем мы выравниваем src по месту назначения:

rotated_src = rotate(src, 20, order=1) # find the angle 20, reshape true is by default
plt.imshow(rotated_src)
plt.show()
distance_y = 8 # find this distances from rotated_src and dst
distance_x = 12 # use any visual reference or even the corners
translated_src = shift(rotated_src, [distance_y,distance_x])
plt.imshow(translated_src)
plt.show()

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

Анализ сложности

Задача состоит в том, чтобы определить три параметра

Предположим, что у вас есть сетка для углов, смещений по осям x и y, каждое из которых имеет размер O(n), и что ваши изображения имеют размер O(n x n), поэтому вращение, перемещение и сравнение изображений занимают O(n^2), поскольку у вас есть O(n^3) преобразования-кандидаты, чтобы попробовать , вы получаете сложность O(n^5), и, вероятно, поэтому вы задаете вопрос.

Однако часть смещения можно вычислить несколько более эффективно, если вычислить максимальную корреляцию с использованием преобразований Фурье. Преобразования Фурье могут быть выполнены со сложностью O(n log n) по каждой оси, и мы должны выполнить их для двух пространственных измерений, полная корреляционная матрица может быть вычислена в O(n^2 log^2 n), затем мы найдем максимум со сложностью O(n^2), поэтому общая временная сложность определения лучший расклад O(n^2 log^2 n). Однако вы все равно хотите найти лучший ракурс, так как у нас есть O(n) ракурсы-кандидаты, общая сложность этого поиска будет O(n^3 log^2 n). Помните, что мы используем Python, и у нас могут быть значительные накладные расходы, поэтому эта сложность дает нам только представление о том, насколько это будет сложно, и я справлялся с такими проблемами раньше, поэтому я начинаю уверенно.

Подготовка некоторого примера

Я начну с загрузки изображения и применения поворота и центрирования изображения нулями.

def centralized(a, width, height):
    '''
    Image centralized to the given width and height
    by padding with zeros (black)
    '''
    assert width >= a.shape[0] and height >= a.shape[1]
    ap = np.zeros((width, height) + a.shape[2:], a.dtype)
    ccx = (width - a.shape[0])//2
    ccy = (height - a.shape[1])//2
    ap[ccx:ccx+a.shape[0], ccy:ccy+a.shape[1], ...] = a
    return ap
def image_pair(im, width, height, displacement=(0,0), angle=0):
    '''
    this build an a pair of images as numpy arrays
    from the input image.
    Both images will be padded with zeros (black)
    and roughly centralized.
    and will have the specified shape
    
    make sure that the width and height chosen are enough 
    to fit the rotated image
    '''
    a = np.array(im)
    a1 = centralized(a, width, height)
    a2 = centralized(ndimage.rotate(a, angle), width, height)
    a2 = np.roll(a2, displacement, axis=(0,1))
    return a1, a2

def random_transform():
    angle = np.random.rand() * 360
    displacement = np.random.randint(-100, 100, 2)
    return displacement, angle

a1, a2 = image_pair(im, 512, 512, *random_transform())
plt.subplot(121)
plt.imshow(a1)
plt.subplot(122)
plt.imshow(a2)

Поиск смещения

Первое, что нужно сделать, это вычислить корреляцию изображения

def compute_correlation(a1, a2):
    A1 = np.fft.rfftn(a1, axes=(0,1))
    A2 = np.fft.rfftn(a2, axes=(0,1))
    C = np.fft.irfftn(np.sum(A1 * np.conj(A2), axis=2))
    return C

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

displacement, _ = random_transform()
a1, a2 = image_pair(im, 521, 512, displacement, angle=0)
C = compute_correlation(a1, a2)
np.unravel_index(np.argmax(C), C.shape), displacement
a3 = np.roll(a2, np.unravel_index(np.argmax(C), C.shape), axis=(0,1))
assert np.all(a3 == a1)

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

Давайте поместим это в функцию для будущего использования

def get_aligned(a1, a2, angle):
    a1_rotated = ndimage.rotate(a1, angle, reshape=False)
    C = compute_correlation(a2, a1_rotated)
    found_displacement = np.unravel_index(np.argmax(C), C.shape)
    a1_aligned = np.roll(a1_rotated, found_displacement, axis=(0,1))
    return a1_aligned

В поисках угла

Теперь мы можем сделать что-то в два шага,

в одном мы вычисляем корреляцию для каждого угла, затем с углом, который дает максимальную корреляцию, находим выравнивание.

displacement, angle = random_transform()
a1, a2 = image_pair(im, 521, 512, displacement, angle)
C_max = []
C_argmax = []
angle_guesses = np.arange(0, 360, 5)
for angle_guess in angle_guesses:
    a1_rotated = ndimage.rotate(a1, angle_guess, reshape=False)
    C = compute_correlation(a1_rotated, a2)
    i = np.argmax(C)
    v = C.reshape(-1)[i]
    C_max.append(v)
    C_argmax.append(i)

Посмотрим, как будет выглядеть корреляция

plt.plot(angle_guesses, C_max);

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

Давайте применим преобразование к исходному изображению и посмотрим, как оно выглядит.

a1_aligned = get_aligned(a1, a2, angle_guesses[np.argmax(C_max)])
plt.subplot(121)
plt.imshow(a2)
plt.subplot(122)
plt.imshow(a1_aligned)

Отлично, я бы не сделал лучше, чем это вручную.

Я использую изображение подсолнуха из соображений красоты, но процедура одинакова для любого типа изображения. Я использую RGB, показывая, что изображение может иметь одно дополнительное измерение, то есть оно использует вектор признаков, вместо скалярного признака, вы можете изменить форму своих данных на (width, height, 1), если ваш признак является скалярным.

Рабочий код ниже на случай, если кому-то еще понадобятся аффинные преобразования scipy:

def affine_test(angle=0, translate=(0, 0), shape=(200, 100), buffered_shape=(300, 200), nblob=50):
    # Maxiumum translation allowed is half difference between shape and buffered_shape

    np.random.seed(42)

    # Generate a buffered_shape-sized base image
    base = np.zeros(buffered_shape, dtype=np.float32)
    random_locs = np.random.choice(np.arange(2, buffered_shape[0] - 2), nblob * 2, replace=False)
    i = random_locs[:nblob]
    j = random_locs[nblob:]
    for k, (_i, _j) in enumerate(zip(i, j)):
        base[_i - 2 : _i + 2, _j - 2 : _j + 2] = k + 10

    # Impose a rotation and translation on source
    src = rotate(base, angle, reshape=False, order=1, mode = "constant")
    bsc = (np.array(buffered_shape) / 2).astype(int)
    sc = (np.array(shape) / 2).astype(int)
    src = src[
        bsc[0] - sc[0] + translate[0] : bsc[0] + sc[0] + translate[0],
        bsc[1] - sc[1] + translate[1] : bsc[1] + sc[1] + translate[1],
    ]
    # Cut-out destination from the centre of the base image
    dst = base[bsc[0] - sc[0] : bsc[0] + sc[0], bsc[1] - sc[1] : bsc[1] + sc[1]]

    src_y, src_x = src.shape

    def get_matrix_offset(centre, angle, scale):
        """Follows OpenCV.getRotationMatrix2D"""
        angle_rad = angle * np.pi / 180
        alpha = np.round(scale * np.cos(angle_rad), 8)
        beta = np.round(scale * np.sin(angle_rad), 8)
        return (
            np.array([[alpha, beta], [-beta, alpha]]),
            np.array(
                [
                    (1 - alpha) * centre[0] - beta * centre[1],
                    beta * centre[0] + (1 - alpha) * centre[1],
                ]
            ),
        )

    matrix, offset = get_matrix_offset(np.array([((src_y - 1) / 2) - translate[0], ((src_x - 1) / 2) - translate[
    1]]), angle, 1)

    offset += np.array(translate)

    M = np.column_stack((matrix, offset))
    M = np.vstack((M, [0, 0, 1]))
    iM = np.linalg.inv(M)
    imatrix = iM[:2, :2]
    ioffset = iM[:2, 2]

    # Determine the outer bounds of the new image
    lin_pts = np.array([[0, src_y-1, src_y-1, 0], [0, 0, src_x-1, src_x-1]])
    transf_lin_pts = np.dot(matrix, lin_pts) + offset.reshape(2, 1) # - np.array(translate).reshape(2, 1) # both?

    # Find min and max bounds of the transformed image
    min_x = np.floor(np.min(transf_lin_pts[1])).astype(int)
    min_y = np.floor(np.min(transf_lin_pts[0])).astype(int)
    max_x = np.ceil(np.max(transf_lin_pts[1])).astype(int)
    max_y = np.ceil(np.max(transf_lin_pts[0])).astype(int)

    # Add translation to the transformation matrix to shift to positive values
    anchor_x, anchor_y = 0, 0
    if min_x < 0:
        anchor_x = -min_x
    if min_y < 0:
        anchor_y = -min_y

    dot_anchor = np.dot(imatrix, [anchor_y, anchor_x])
    shifted_offset = ioffset - dot_anchor

    # Create padded destination image
    dst_y, dst_x = dst.shape[:2]
    pad_widths = [anchor_y, max(max_y, dst_y) - dst_y, anchor_x, max(max_x, dst_x) - dst_x]
    dst_padded = np.pad(
        dst,
        ((pad_widths[0], pad_widths[1]), (pad_widths[2], pad_widths[3])),
        "constant",
        constant_values=-10,
    )

    dst_pad_y, dst_pad_x = dst_padded.shape
    # Create the aligned and padded source image
    source_aligned = affine_transform(
        src,
        imatrix,
        offset=shifted_offset,
        output_shape=(dst_pad_y, dst_pad_x),
        order=3,
        mode = "constant",
        cval=-10,
    )

Например. Бег:

affine_test(angle=-25, translate=(10, -40))

покажет:

и увеличил:

Извините, код не очень хорошо написан как есть.

Обратите внимание, что, запуская это в дикой природе, я заметил, что он не может обрабатывать какие-либо изменения масштаба изображений, но я не уверен, что это не связано с тем, как я вычисляю преобразование, поэтому предостережение, которое стоит отметить и проверить, если вы выравниваете изображения с разными масштабами.