From 719893af0b83545512e057e9477201553461b88b Mon Sep 17 00:00:00 2001 From: Simon Kiefhaber Date: Fri, 15 Oct 2021 23:42:24 +0200 Subject: [PATCH 1/8] Added implementation of the Horn-Schunck algorithm --- computer_vision/horn_schunck.py | 126 ++++++++++++++++++++++++++++++++ 1 file changed, 126 insertions(+) create mode 100644 computer_vision/horn_schunck.py diff --git a/computer_vision/horn_schunck.py b/computer_vision/horn_schunck.py new file mode 100644 index 000000000000..f6997e98589c --- /dev/null +++ b/computer_vision/horn_schunck.py @@ -0,0 +1,126 @@ +""" + The Horn-Schunck method estimates the optical flow for every single pixel of + a sequence of images. + It works by assuming brightness constancy between two consecutive frames + and smoothness in the optical flow. + + Useful resources: + Wikipedia: https://en.wikipedia.org/wiki/Horn%E2%80%93Schunck_method + Paper: http://image.diku.dk/imagecanon/material/HornSchunckOptical_Flow.pdf +""" + +from typing import Optional + +import numpy as np +from scipy.ndimage.filters import convolve +from typing_extensions import SupportsIndex + + +def warp(image: np.ndarray, u: np.ndarray, v: np.ndarray) -> np.ndarray: + """ + Warps the pixels of an image into a new image using the horizontal and vertical + flows. + Pixels that are warped from an invalid location are set to 0. + + Parameters: + image: Grayscale image + u: Horizontal flow + v: Vertical flow + + Returns: Warped image + """ + flow = np.stack((u, v), 2) + + # Create a grid of all pixel coordinates and subtract the offsets + grid = np.stack( + np.meshgrid(np.arange(0, image.shape[1]), np.arange(0, image.shape[0])), 2 + ) + grid = np.round(grid - flow).astype(np.int) + + # Find the locations outside of the original image + invalid = (grid < 0) | (grid >= np.array([image.shape[1], image.shape[0]])) + grid[invalid] = 0 + + warped = image[grid[:, :, 1], grid[:, :, 0]] + + # Set pixels at invalid locations to 0 + warped[invalid[:, :, 0] | invalid[:, :, 1]] = 0 + + return warped + + +def calc_derivatives( + im0: np.ndarray, im1: np.ndarray +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Calculates the vertical and horizontal image derivatives as well as the + time derivatives used in the Horn-Schunck algorithm + + Parameters: + im0: First image + im1: Warped second image + + Returns: + dX: Horizontal image derivative + dY: Vertical image derivative + dT: Time derivative + """ + # Prepare kernels for the calculation of the derivatives + kernel_x = np.array([[-1, 1], [-1, 1]]) * 0.25 + kernel_y = np.array([[-1, -1], [1, 1]]) * 0.25 + kernel_t = np.array([[1, 1], [1, 1]]) * 0.25 + + # Calculate the derivatives + dX = convolve(im0, kernel_x) + convolve(im1, kernel_x) + dY = convolve(im0, kernel_y) + convolve(im1, kernel_y) + dT = convolve(im0, kernel_t) + convolve(im1, -kernel_t) + + return dX, dY, dT + + +def horn_schunck( + im0: np.ndarray, + im1: np.ndarray, + num_iter: SupportsIndex, + alpha: Optional[float] = None, +) -> tuple[np.ndarray, np.ndarray]: + """ + This function performs the Horn-Schunck algorithm and returns the estimated + optical flow. It is assumed that the input images are grayscale and + normalized to be in [0, 1]. + + Parameters: + im0: First image + im1: Second image + alpha: Regularization constant + num_iter: Number of iterations performed + + Returns: + u: Horizontal flow + v: Vertical flow + """ + if alpha is None: + alpha = 0.1 + + # Initialize flow + u = np.zeros_like(im0) + v = np.zeros_like(im0) + + laplacian_kernel = np.array( + [[1 / 12, 1 / 6, 1 / 12], [1 / 6, 0, 1 / 6], [1 / 12, 1 / 6, 1 / 12]] + ) + + # Iteratively refine the flow + for _ in range(num_iter): + dx, dy, dt = calc_derivatives(warp(im0, u, v), im1) + + avg_u = convolve(u, laplacian_kernel) + avg_v = convolve(v, laplacian_kernel) + + # This updates the flow as proposed in the paper (Step 12) + d = (dx * avg_u + dy * avg_v + dt) / (alpha ** 2 + dx ** 2 + dy ** 2) + + u = avg_u - dx * d + v = avg_v - dy * d + + return u, v From 36b9c2ca25d7e5f89c0667c2db15a0187af783f0 Mon Sep 17 00:00:00 2001 From: Simon Kiefhaber Date: Sat, 16 Oct 2021 00:41:37 +0200 Subject: [PATCH 2/8] Cleaner variable names --- computer_vision/horn_schunck.py | 84 ++++++++++++++------------------- 1 file changed, 35 insertions(+), 49 deletions(-) diff --git a/computer_vision/horn_schunck.py b/computer_vision/horn_schunck.py index f6997e98589c..9233a8ca3ebd 100644 --- a/computer_vision/horn_schunck.py +++ b/computer_vision/horn_schunck.py @@ -16,7 +16,9 @@ from typing_extensions import SupportsIndex -def warp(image: np.ndarray, u: np.ndarray, v: np.ndarray) -> np.ndarray: +def warp( + image: np.ndarray, horizontal_flow: np.ndarray, vertical_flow: np.ndarray +) -> np.ndarray: """ Warps the pixels of an image into a new image using the horizontal and vertical flows. @@ -24,18 +26,19 @@ def warp(image: np.ndarray, u: np.ndarray, v: np.ndarray) -> np.ndarray: Parameters: image: Grayscale image - u: Horizontal flow - v: Vertical flow + horizontal_flow: Horizontal flow + vertical_flow: Vertical flow Returns: Warped image + """ - flow = np.stack((u, v), 2) + flow = np.stack((horizontal_flow, vertical_flow), 2) # Create a grid of all pixel coordinates and subtract the offsets grid = np.stack( np.meshgrid(np.arange(0, image.shape[1]), np.arange(0, image.shape[0])), 2 ) - grid = np.round(grid - flow).astype(np.int) + grid = np.round(grid - flow).astype(np.int32) # Find the locations outside of the original image invalid = (grid < 0) | (grid >= np.array([image.shape[1], image.shape[0]])) @@ -49,38 +52,9 @@ def warp(image: np.ndarray, u: np.ndarray, v: np.ndarray) -> np.ndarray: return warped -def calc_derivatives( - im0: np.ndarray, im1: np.ndarray -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """ - Calculates the vertical and horizontal image derivatives as well as the - time derivatives used in the Horn-Schunck algorithm - - Parameters: - im0: First image - im1: Warped second image - - Returns: - dX: Horizontal image derivative - dY: Vertical image derivative - dT: Time derivative - """ - # Prepare kernels for the calculation of the derivatives - kernel_x = np.array([[-1, 1], [-1, 1]]) * 0.25 - kernel_y = np.array([[-1, -1], [1, 1]]) * 0.25 - kernel_t = np.array([[1, 1], [1, 1]]) * 0.25 - - # Calculate the derivatives - dX = convolve(im0, kernel_x) + convolve(im1, kernel_x) - dY = convolve(im0, kernel_y) + convolve(im1, kernel_y) - dT = convolve(im0, kernel_t) + convolve(im1, -kernel_t) - - return dX, dY, dT - - def horn_schunck( - im0: np.ndarray, - im1: np.ndarray, + image0: np.ndarray, + image1: np.ndarray, num_iter: SupportsIndex, alpha: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: @@ -90,37 +64,49 @@ def horn_schunck( normalized to be in [0, 1]. Parameters: - im0: First image - im1: Second image + image0: First image of the sequence + image1: Second image of the sequence alpha: Regularization constant num_iter: Number of iterations performed Returns: - u: Horizontal flow - v: Vertical flow + horizontal_flow: Horizontal flow + vertical_flow: Vertical flow """ if alpha is None: alpha = 0.1 # Initialize flow - u = np.zeros_like(im0) - v = np.zeros_like(im0) + horizontal_flow = np.zeros_like(image0) + vertical_flow = np.zeros_like(image0) + # Prepare kernels for the calculation of the derivatives and the average velocity + kernel_x = np.array([[-1, 1], [-1, 1]]) * 0.25 + kernel_y = np.array([[-1, -1], [1, 1]]) * 0.25 + kernel_t = np.array([[1, 1], [1, 1]]) * 0.25 laplacian_kernel = np.array( [[1 / 12, 1 / 6, 1 / 12], [1 / 6, 0, 1 / 6], [1 / 12, 1 / 6, 1 / 12]] ) # Iteratively refine the flow for _ in range(num_iter): - dx, dy, dt = calc_derivatives(warp(im0, u, v), im1) + warped_image = warp(image0, horizontal_flow, vertical_flow) + derivative_x = convolve(warped_image, kernel_x) + convolve(image1, kernel_x) + derivative_y = convolve(warped_image, kernel_y) + convolve(image1, kernel_y) + derivative_t = convolve(warped_image, kernel_t) + convolve(image1, -kernel_t) - avg_u = convolve(u, laplacian_kernel) - avg_v = convolve(v, laplacian_kernel) + avg_horizontal_velocity = convolve(horizontal_flow, laplacian_kernel) + avg_vertical_velocity = convolve(vertical_flow, laplacian_kernel) # This updates the flow as proposed in the paper (Step 12) - d = (dx * avg_u + dy * avg_v + dt) / (alpha ** 2 + dx ** 2 + dy ** 2) + update = ( + derivative_x * avg_horizontal_velocity + + derivative_y * avg_vertical_velocity + + derivative_t + ) + update = update / (alpha ** 2 + derivative_x ** 2 + derivative_y ** 2) - u = avg_u - dx * d - v = avg_v - dy * d + horizontal_flow = avg_horizontal_velocity - derivative_x * update + vertical_flow = avg_vertical_velocity - derivative_y * update - return u, v + return horizontal_flow, vertical_flow From 74ec1eb20e41829cc9bfb1a1ac35e093396cb0d5 Mon Sep 17 00:00:00 2001 From: Simon Kiefhaber Date: Sat, 16 Oct 2021 01:02:27 +0200 Subject: [PATCH 3/8] added doctests --- computer_vision/horn_schunck.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/computer_vision/horn_schunck.py b/computer_vision/horn_schunck.py index 9233a8ca3ebd..63b0b0db941d 100644 --- a/computer_vision/horn_schunck.py +++ b/computer_vision/horn_schunck.py @@ -31,10 +31,17 @@ def warp( Returns: Warped image + >>> warp(np.array([[0, 1, 2], [0, 3, 0], [2, 2, 2]]), \ + np.array([[0, 1, -1], [-1, 0, 0], [1, 1, 1]]), \ + np.array([[0, 0, 0], [0, 1, 0], [0, 0, 1]])) + array([[0, 0, 0], + [3, 1, 0], + [0, 2, 3]]) """ flow = np.stack((horizontal_flow, vertical_flow), 2) - # Create a grid of all pixel coordinates and subtract the offsets + # Create a grid of all pixel coordinates and subtract the flow to get the + # target pixels coordinates grid = np.stack( np.meshgrid(np.arange(0, image.shape[1]), np.arange(0, image.shape[0])), 2 ) @@ -69,9 +76,16 @@ def horn_schunck( alpha: Regularization constant num_iter: Number of iterations performed - Returns: - horizontal_flow: Horizontal flow - vertical_flow: Vertical flow + Returns: estimated horizontal & vertical flow + + >>> np.round(horn_schunck(np.array([[0, 0, 2], [0, 0, 2]]), \ + np.array([[0, 2, 0], [0, 2, 0]]), alpha=0.1, num_iter=110)).\ + astype(np.int32) + array([[[ 0, -1, -1], + [ 0, -1, -1]], + + [[ 0, 0, 0], + [ 0, 0, 0]]]) """ if alpha is None: alpha = 0.1 @@ -84,7 +98,7 @@ def horn_schunck( kernel_x = np.array([[-1, 1], [-1, 1]]) * 0.25 kernel_y = np.array([[-1, -1], [1, 1]]) * 0.25 kernel_t = np.array([[1, 1], [1, 1]]) * 0.25 - laplacian_kernel = np.array( + kernel_laplacian = np.array( [[1 / 12, 1 / 6, 1 / 12], [1 / 6, 0, 1 / 6], [1 / 12, 1 / 6, 1 / 12]] ) @@ -95,8 +109,8 @@ def horn_schunck( derivative_y = convolve(warped_image, kernel_y) + convolve(image1, kernel_y) derivative_t = convolve(warped_image, kernel_t) + convolve(image1, -kernel_t) - avg_horizontal_velocity = convolve(horizontal_flow, laplacian_kernel) - avg_vertical_velocity = convolve(vertical_flow, laplacian_kernel) + avg_horizontal_velocity = convolve(horizontal_flow, kernel_laplacian) + avg_vertical_velocity = convolve(vertical_flow, kernel_laplacian) # This updates the flow as proposed in the paper (Step 12) update = ( From 6cff7dfb43f7a33295be229461a2a2bc9d1c0db1 Mon Sep 17 00:00:00 2001 From: Simon Kiefhaber Date: Sat, 16 Oct 2021 13:52:04 +0200 Subject: [PATCH 4/8] Fix doctest --- computer_vision/horn_schunck.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/computer_vision/horn_schunck.py b/computer_vision/horn_schunck.py index 63b0b0db941d..33dc19c1a1f3 100644 --- a/computer_vision/horn_schunck.py +++ b/computer_vision/horn_schunck.py @@ -85,7 +85,7 @@ def horn_schunck( [ 0, -1, -1]], [[ 0, 0, 0], - [ 0, 0, 0]]]) + [ 0, 0, 0]]], dtype=int32) """ if alpha is None: alpha = 0.1 From e705cf2db7c4477e058cae698ed27fc6e48aa6fd Mon Sep 17 00:00:00 2001 From: John Law Date: Mon, 2 May 2022 23:53:25 +0800 Subject: [PATCH 5/8] Update horn_schunck.py --- computer_vision/horn_schunck.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/computer_vision/horn_schunck.py b/computer_vision/horn_schunck.py index 33dc19c1a1f3..fcd258b2b0e7 100644 --- a/computer_vision/horn_schunck.py +++ b/computer_vision/horn_schunck.py @@ -124,3 +124,9 @@ def horn_schunck( vertical_flow = avg_vertical_velocity - derivative_y * update return horizontal_flow, vertical_flow + + +if __name__ == "__main__": + import doctest + + doctest.testmod() From b88c8b6abcdae452234aec1e10c0a98772bd468b Mon Sep 17 00:00:00 2001 From: John Law Date: Mon, 2 May 2022 23:55:49 +0800 Subject: [PATCH 6/8] Update horn_schunck.py --- computer_vision/horn_schunck.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/computer_vision/horn_schunck.py b/computer_vision/horn_schunck.py index fcd258b2b0e7..59519ff1cd2d 100644 --- a/computer_vision/horn_schunck.py +++ b/computer_vision/horn_schunck.py @@ -118,7 +118,7 @@ def horn_schunck( + derivative_y * avg_vertical_velocity + derivative_t ) - update = update / (alpha ** 2 + derivative_x ** 2 + derivative_y ** 2) + update = update / (alpha**2 + derivative_x**2 + derivative_y**2) horizontal_flow = avg_horizontal_velocity - derivative_x * update vertical_flow = avg_vertical_velocity - derivative_y * update From d5da626ca0b5287988dc8d2004a2b16f6a1fdcc1 Mon Sep 17 00:00:00 2001 From: John Law Date: Mon, 2 May 2022 23:59:27 +0800 Subject: [PATCH 7/8] Update horn_schunck.py --- computer_vision/horn_schunck.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/computer_vision/horn_schunck.py b/computer_vision/horn_schunck.py index 59519ff1cd2d..1d2d5fc365ba 100644 --- a/computer_vision/horn_schunck.py +++ b/computer_vision/horn_schunck.py @@ -63,7 +63,7 @@ def horn_schunck( image0: np.ndarray, image1: np.ndarray, num_iter: SupportsIndex, - alpha: Optional[float] = None, + alpha: float | None = None, ) -> tuple[np.ndarray, np.ndarray]: """ This function performs the Horn-Schunck algorithm and returns the estimated From 1d5d7292fc5062ad828ca67182f19a3b553a2de6 Mon Sep 17 00:00:00 2001 From: John Law Date: Tue, 3 May 2022 00:02:34 +0800 Subject: [PATCH 8/8] Update horn_schunck.py --- computer_vision/horn_schunck.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/computer_vision/horn_schunck.py b/computer_vision/horn_schunck.py index 1d2d5fc365ba..1428487d051b 100644 --- a/computer_vision/horn_schunck.py +++ b/computer_vision/horn_schunck.py @@ -9,8 +9,6 @@ Paper: http://image.diku.dk/imagecanon/material/HornSchunckOptical_Flow.pdf """ -from typing import Optional - import numpy as np from scipy.ndimage.filters import convolve from typing_extensions import SupportsIndex