Я реализую иерархический метод Лукаса-Канаде в Python на основе
этого руководства. Однако,
применяя этот метод к вращающейся сфере, я получаю неожиданные результаты.
[img]https://i .stack.imgur.com/M0urG.gif[/img]
- Данные можно найти
здесь.
Общая структура алгоритма объясняется ниже. Обратите внимание, что в уравнениях
этого руководства используется
соглашение, согласно которому верхний левый угол равен (0, 0), а нижний правый угол равен
Код: Выделить всё
(width-1, height-1)угла — (высота-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)
- Похоже, что направление X потока правильное, а направление Y — нет. Это
-
- На неподвижных пикселях наблюдается ненулевой поток.
вам!
Подробнее здесь: https://stackoverflow.com/questions/775 ... tical-flow