Проблема с иерархическим методом Лукаса-Канаде для оптического потока ⇐ Python

[Anonymous]

Проблема с иерархическим методом Лукаса-Канаде для оптического потока
Я реализую иерархический метод Лукаса-Канаде в Python на основе
этого руководства. Однако,
применяя этот метод к вращающейся сфере, я получаю неожиданные результаты.
[img]https://i .stack.imgur.com/M0urG.gif[/img]

• Данные можно найти
  здесь.

Пояснение алгоритма
Общая структура алгоритма объясняется ниже. Обратите внимание, что в уравнениях
этого руководства используется
соглашение, согласно которому верхний левый угол равен (0, 0), а нижний правый угол равен

Код: Выделить всё

(width-1, height-1)

, а представленная здесь реализация меняет местами оси x и y. Другими словами, координата верхнего левого угла — (0, 0), а правого нижнего
угла — (высота-1, ширина-1) .

Базовый метод Лукаса-Канаде

С учетом уравнений 19, 20, 23, 29 и 28 базовый метод Лукаса-Канаде
реализовано следующим образом:

Код: Выделить всё

def lucas_kanade(img1, img2):
img1 = np.copy(img1).astype(np.float32)
img2 = np.copy(img2).astype(np.float32)

# Change the window size based on the image's size due to downsampling.
window_size = min(max(3, min(img1.shape[:2]) / 6), 31)
window_size = int(2 * (window_size // 2) + 1)
print("window size: ", window_size)

# Compute image gradients
Ix = np.zeros(img1.shape, dtype=np.float32)
Iy = np.zeros(img1.shape, dtype=np.float32)
Ix[1:-1, 1:-1] = (img1[1:-1, 2:] - img1[1:-1, :-2]) / 2  # pixels on boundry are 0.
Iy[1:-1, 1:-1] = (img1[2:, 1:-1] - img1[:-2, 1:-1]) / 2

# Compute temporal gradient
It = np.zeros(img1.shape, dtype=np.float32)
It = img1 - img2

# Define a (window_size, window_size) kernel for the convolution
kernel = np.ones((window_size, window_size), dtype=np.float32)
# kernel = create_gaussian_kernel(window_size, sigma=1)

# Use convolution to calculate the sum of the window for each pixel
Ix2 = convolve2d(Ix**2, kernel, mode="same", boundary="fill", fillvalue=0)
Iy2 = convolve2d(Iy**2, kernel, mode="same", boundary="fill", fillvalue=0)
Ixy = convolve2d(Ix * Iy, kernel, mode="same", boundary="fill", fillvalue=0)
Ixt = convolve2d(Ix * It, kernel, mode="same", boundary="fill", fillvalue=0)
Iyt = convolve2d(Iy * It, kernel, mode="same", boundary="fill", fillvalue=0)

# Compute optical flow parameters
det = Ix2 * Iy2 - Ixy**2
# Avoid division by zero
u = np.where((det > 1e-6), (Iy2 * Ixt - Ixy * Iyt) / det, 0)
v = np.where((det > 1e-6), (Ix2 * Iyt - Ixy * Ixt) / det, 0)

optical_flow = np.stack((u, v), axis=2)
return optical_flow.astype(np.float32)

Сгенерировать пирамиду Гаусса

Пирамида Гаусса генерируется следующим образом.

Код: Выделить всё

def gen_gaussian_pyramid(im, max_level):
# Return `max_level+1` arrays.
gauss_pyr = [im]
for i in range(max_level):
gauss_pyr.append(cv2.pyrDown(gauss_pyr[-1]))
return gauss_pyr

Повышает дискретизацию потока

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

Код: Выделить всё

def expand(img, dst_size, interpolation=None):
# Increase dimension.
height, width = dst_size[:2]
return cv2.GaussianBlur(
cv2.resize(  # dim: (width, height)
img, (width, height), interpolation=interpolation or cv2.INTER_LINEAR
),
(5, 5),
0,
)

Деформировать изображение потоком на предыдущем уровне

В уравнении 12 нужно сместить правое изображение в зависимости от количества пикселей в
Предыдущий цикл. Я решил использовать функцию opencv.remap для деформации левого изображения, чтобы оно
выровнено по правому изображению.

Код: Выделить всё

def remap(a, flow):
height, width = flow.shape[:2]

# Create a grid of coordinates using np.meshgrid
y, x = np.meshgrid(np.arange(height), np.arange(width), indexing="ij")

# Create flow_map by adding the flow vectors
flow_map = np.column_stack( # NOTE: minus sign on flow
(x.flatten() + -flow[:, :, 0].flatten(), y.flatten() + -flow[:, :, 1].flatten())
)

# Reshape flow_map to match the original image dimensions
flow_map = flow_map.reshape((height, width, 2))

# Ensure flow_map values are within the valid range
flow_map[:, :, 0] = np.clip(flow_map[:, :, 0], 0, width - 1)
flow_map[:, :, 1] = np.clip(flow_map[:, :, 1], 0, height - 1)

# Convert flow_map to float32
flow_map = flow_map.astype(np.float32)

# Use cv2.remap for remapping
warped = cv2.remap(a, flow_map, None, cv2.INTER_LINEAR)

return warped

Собираем все вместе

Определив все основы, мы можем собрать их все вместе. Здесь g_L и d_L — это
переменная в уравнении 7.

Код: Выделить всё

def hierarchical_lucas_kanade(im1, im2, max_level):  # max_level = 4
gauss_pyr_1 = gen_gaussian_pyramid(im1, max_level)  # from finest to roughest
gauss_pyr_2 = gen_gaussian_pyramid(im2, max_level)  # from finest to roughest

g_L = [0 for _ in range(max_level + 1)]  # Every slot will be (h, w, 2) array.
d_L = [0 for _ in range(max_level + 1)]  # Every slot will be (h, w, 2) array.
assert len(g_L) == 5  # 4 + 1 (base)
# Initialzie g_L[0] as (h, w, 2) zeros array
g_L[max_level] = np.zeros(gauss_pyr_1[-1].shape[:2] + (2,)).astype(np.float32)

for level in range(max_level, -1, -1):  # 4, 3, 2, 1, 0
# Warp image 1 by previous flow.
warped = remap(gauss_pyr_1[level], g_L[level])
# Run Lucas-Kanade on warped image and right image.
d_L[level] = lucas_kanade(warped, gauss_pyr_2[level])

# Expand/Upsample the flow so that the dimension can match the finer result.
g_L[level - 1] = 2.0 * expand(
g_L[level] + d_L[level],
gauss_pyr_2[level - 1].shape[:2] + (2,),
interpolation=cv2.INTER_LINEAR,
)
return g_L[0] + d_L[0]

Визуализация
После загрузки данных вы можете запустить их с помощью кода:

Код: Выделить всё

sphere_seq = []
for fname in natsorted(Path("./input/sphere/").rglob("*.ppm")):
sphere_seq.append(cv2.imread(str(fname), cv2.IMREAD_GRAYSCALE))

flows = []
for i in range(len(sphere_seq) - 1):
flows.append(hierarchical_lucas_kanade(sphere_seq[i], sphere_seq[i + 1], max_level=4))
show_flow(sphere_seq[i], flows[i], f"./output/sphere/flow-{i}.png")

Результат выглядит следующим образом:
[img]https://i.stack.imgur .com/1oK0e.gif[/img]

Конкретный вопрос:
В результате возникает несколько проблем:

• • Похоже, что направление X потока правильное, а направление Y — нет. Это
    может быть, мой код визуализации неправильный. Вот код:

    Код: Выделить всё

     def show_flow(img, flow, filename=None):
    x = np.arange(0, img.shape[1], 1)
    y = np.arange(0, img.shape[0], 1)
    x, y = np.meshgrid(x, y)
    plt.figure(figsize=(10, 10))
    fig = plt.imshow(img, cmap="gray", interpolation="bicubic")
    plt.axis("off")
    fig.axes.get_xaxis().set_visible(False)
    fig.axes.get_yaxis().set_visible(False)
    
    num_points_per_axis = 32
    step = int(img.shape[0] / num_points_per_axis)
    plt.quiver(
    x[::step, ::step],
    y[::step, ::step],  # reverse order?
    flow[::step, ::step, 0],
    flow[::step, ::step, 1],  # reverse sign?
    color="r",
    pivot="tail",
    headwidth=2,
    headlength=3,
    )
    if filename is not None:
    plt.savefig(filename, bbox_inches="tight", pad_inches=0)
    

• • На неподвижных пикселях наблюдается ненулевой поток.

Будем очень признательны за любые идеи или предложения по решению этой проблемы. Спасибо
вам!


Подробнее здесь: https://stackoverflow.com/questions/775 ... tical-flow