Почему выходные данные содержат только 2 значения, но не смещение для всего изображения?

Question

Я застрял здесь на некоторое время.Я не могу понять, что я делаю неправильно, вычисляя векторы смещения вдоль оси x и оси y, используя метод Lucas Kanade .

Я реализовал это, как указано в приведенной выше ссылке на Википедию.Вот что я сделал:

    import cv2
    import numpy as np


    img_a = cv2.imread("./images/1.png",0)
    img_b = cv2.imread("./images/2.png",0)


    # Calculate gradient along x and y axis
    ix = cv2.Sobel(img_a, cv2.CV_64F, 1, 0, ksize = 3, scale = 1.0/3.0)
    iy = cv2.Sobel(img_a, cv2.CV_64F, 0, 1, ksize = 3, scale = 1.0/3.0)

    # Calculate temporal difference between the 2 images
    it = img_b - img_a


    ix = ix.flatten()
    iy = iy.flatten()
    it = -it.flatten()

    A = np.vstack((ix, iy)).T


    atai = np.linalg.inv(np.dot(A.T,A))
    atb = np.dot(A.T, it)

    v = np.dot(np.dot(np.linalg.inv(np.dot(A.T,A)),A.T),it)

    print(v)

Этот код выполняется без ошибок, но печатает массив из 2 значений!Я ожидал, что матрица v будет того же размера, что и изображение.Почему это происходит?Что я делаю неправильно?

PS: Я знаю, что есть методы, напрямую доступные в OpenCV, но я хочу написать этот простой алгоритм (как также дано в ссылке на Википедию, которой я поделился выше) сам.

Answer 1

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

Это рецепт(обозначение относится к используемому на странице Википедии ):

1. Вычислить градиент изображения ( A ) для первого изображения (ix, iy в ОП), используя любой метод (у Собеля все в порядке, я предпочитаю гауссовские производные).

   ix = cv2.Sobel(img_a, cv2.CV_64F, 1, 0, ksize = 3, scale = 1.0/3.0)
   iy = cv2.Sobel(img_a, cv2.CV_64F, 0, 1, ksize = 3, scale = 1.0/3.0)
   

2. Вычисление тензора структуры ( A ^T WA ): Axx = ix * ix, Axy = ix * iy, Ayy = iy * iy.Каждое из этих трех изображений должно быть сглажено гауссовым фильтром (это окно).Например,

   Axx = cv2.GaussianBlur(ix * ix, (0,0), 5)
   Axy = cv2.GaussianBlur(ix * iy, (0,0), 5)
   Ayy = cv2.GaussianBlur(iy * iy, (0,0), 5)
   

   Эти три изображения вместе образуют тензор структуры, который представляет собой симметричную матрицу 2x2 в каждом пикселе.Для пикселя с (i,j) матрица имеет вид:

   |  Axx(i,j)  Axy(i,j)  |
   |  Axy(i,j)  Ayy(i,j)  |
   

3. Вычислить временной градиент ( b ), вычитая два изображения (it вOP).

   it = img_b - img_a
   

4. Вычислить A ^T Wb : Abx = ix * it, Aby = iy * it и сгладить эти два изображения с помощьютот же фильтр Гаусса, что и выше.

   Abx = cv2.GaussianBlur(ix * it, (0,0), 5)
   Aby = cv2.GaussianBlur(iy * it, (0,0), 5)
   

5. Вычислить инверсию A ^T WA (симметричная положительно-определенная матрица) иумножить на A ^T Wb .Обратите внимание, что эта инверсия имеет матрицу 2x2 в каждом пикселе, а не изображения в целом.Вы можете записать это как набор простых арифметических операций над изображениями Axx, Axy, Ayy, Abx и Aby.

   Обратная матрица A ^T WA определяется как:

   |  Ayy -Axy  |  
   | -Axy  Axx  | / ( Axx*Ayy - Axy*Axy )
   

   , поэтому вы можете записать решение в виде

   norm = Axx*Ayy - Axy*Axy
   vx = ( Ayy * Abx - Axy * Aby ) / norm
   vy = ( Axx * Aby - Axy * Abx ) / norm
   

   Если изображение натуральное, оно будет иметьхоть немного шума, и norm не будет иметь нулей.Но для искусственных изображений norm может иметь нули, что означает, что вы не можете делить на него.Простое добавление к нему небольшого значения позволит избежать деления на ноль ошибок: norm += 1e-6.

Размер фильтра Гаусса выбран как компромисс между точностью и допустимой скоростью движения: большийФильтр будет давать менее точные результаты, но будет работать с большими сдвигами между изображениями.

Обычно vx и vy оцениваются только тогда, когда два собственных значения матрицы A ^T WA достаточно велики (если хотя бы один из них мал, результат является неточным или, возможно, неправильным).

---

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

import PyDIP as dip

img_a = dip.ImageRead("./images/1.png")
img_b = dip.ImageRead("./images/2.png")

A = dip.Gradient(img_a, [1.0])
b = img_b - img_a
ATA = dip.Gauss(A * dip.Transpose(A), [5.0])
ATb = dip.Gauss(A * b, [5.0])
v = dip.Inverse(ATA) * ATb