Core Algorithms#

Low-level optical flow computation engine implementing variational optical flow with multi-scale pyramid approach.

Optical Flow#

Optical Flow Implementation#

This module implements the variational optical flow computation using a pyramid-based multi-scale approach with gradient constancy assumption.

This is the default ‘flowreg’ backend implementation that maintains compatibility with the MATLAB Flow-Registration toolbox.

Functions#

get_displacement: Main API function for computing optical flow between two frames
get_motion_tensor_gc: Compute motion tensor components for gradient constancy
get_motion_tensor_gray: Compute motion tensor components for gray-value constancy
get_motion_tensor_cs: Compute motion tensor components for census constancy
imregister_wrapper: Warp an image using computed displacement fields
warpingDepth: Calculate pyramid depth based on image dimensions
add_boundary: Add boundary padding to arrays

Notes

The implementation uses a coarse-to-fine pyramid scheme where flow is computed at progressively finer resolutions, with each level initialized from the previous coarser level.

Main API#

pyflowreg.core.optical_flow.get_displacement(fixed, moving, alpha=(2, 2), update_lag=5, iterations=50, min_level=0, levels=50, eta=0.8, a_smooth=1.0, a_data=0.45, const_assumption='gc', uv=None, weight=None, level_solver_backend=None, gnc_schedule=None, warping_steps=None)[source]#

Compute optical flow displacement field using variational approach.

This function implements the main pyramid-based variational optical flow algorithm using gradient constancy assumption with non-linear diffusion regularization. Flow is computed from coarse to fine scales, with each level initialized from the upsampled result of the previous coarser level.

Parameters:

fixed (np.ndarray) – Reference (fixed) image, shape (H, W) or (H, W, C)
moving (np.ndarray) – Moving image to register, shape (H, W) or (H, W, C)
alpha (tuple of float, default=(2, 2)) – Regularization strength (alpha_x, alpha_y) controlling smoothness. Larger values enforce smoother flow fields.
update_lag (int, default=5) – Number of iterations between updates of non-linearity weights (psi). Smaller values update more frequently (slower, potentially more accurate). Larger values update less frequently (faster convergence).
iterations (int, default=50) – Number of SOR iterations per pyramid level
min_level (int, default=0) – Minimum (finest) pyramid level to compute. 0 = full resolution.
levels (int, default=50) – Maximum number of pyramid levels attempted
eta (float, default=0.8) – Pyramid downsampling factor per level (0 < eta <= 1). Each level is eta times the size of the previous level.
a_smooth (float, default=1.0) – Exponent for generalized Charbonnier penalty on smoothness term.
a_data (float, default=0.45) – Exponent for generalized Charbonnier penalty on data term.
const_assumption (str, default=’gc’) – Constancy assumption. Supported values are ‘gc’/’gradient’, ‘gray’/’brightness’, and ‘cs’/’census’.
uv (np.ndarray, optional) – Initial displacement field (H, W, 2) with [u, v] components.
weight (np.ndarray or list, optional) – Channel weights for multi-channel registration.
level_solver_backend (Callable, optional) – Custom level solver function to use instead of the default CPU solver.
gnc_schedule (sequence of float, optional) – Opt-in stage weights interpolating from quadratic (0.0) to fully robust (1.0). Each stage reruns the pyramid with the previous stage result used as initialization.
warping_steps (int, optional) – Number of warp/relinearize steps per pyramid level in optional GNC mode. If omitted, GNC defaults to 10 steps per level. Ignored when GNC is off.

Warping and Registration#

pyflowreg.core.optical_flow.imregister_wrapper(f2_level, u, v, f1_level, interpolation_method='cubic')[source]#

Backward warp of moving image using displacement field.

Performs backward registration by warping f2_level toward f1_level using displacement field (u, v) with bicubic or bilinear interpolation via cv2.remap. Out-of-bounds pixels are replaced with corresponding pixels from f1_level.

Parameters:

f2_level (np.ndarray) – Moving image to warp, shape (H, W) or (H, W, C)
u (np.ndarray) – Horizontal displacement field, shape (H, W)
v (np.ndarray) – Vertical displacement field, shape (H, W)
f1_level (np.ndarray) – Fixed (reference) image, shape (H, W) or (H, W, C)
interpolation_method (str, default=’cubic’) – Interpolation method: ‘cubic’ (bicubic) or ‘linear’ (bilinear). Defaults to bicubic following Sun et al. best practices.

Returns:

warped (np.ndarray) – Backward-warped image, shape (H, W) or (H, W, C)

Notes

The displacement convention is: warped_pos = original_pos + (u, v) Out-of-bounds regions use values from f1_level to maintain continuity.

Bicubic interpolation is more accurate than bilinear for optical flow estimation and is the recommended default.

References

pyflowreg.core.optical_flow.warpingDepth(eta, levels, m, n)[source]#

Calculate maximum pyramid depth for given dimension and warping factor.

Determines how many pyramid levels can be computed given the downsampling factor eta before the dimension becomes too small (< 10 pixels) for reliable optical flow estimation. At pyramid level i, the dimension size is dim * eta^i, where dim = min(m, n).

Parameters:

eta (float) – Pyramid downsampling factor per level (0 < eta <= 1)
levels (int) – Maximum number of levels to attempt
m (int) – First dimension
n (int) – Second dimension

Returns:

warpingdepth (int) – Maximum pyramid depth satisfying: round(dim * eta^i) >= 10, where dim = min(m, n). Approximately floor(log(10/dim) / log(eta)).

Notes

When called from get_displacement with (m, m) for height and (n, n) for width, this enables independent pyramid depth computation per dimension, allowing narrow ROIs to achieve large displacements along their longer dimension without being limited by the shorter dimension.

Motion Tensor Computation#

pyflowreg.core.optical_flow.get_motion_tensor_gc(f1, f2, hy, hx)[source]#

Compute motion tensor components for gradient constancy assumption.

Calculates the motion tensor components that describe the linearized optical flow constraints under the gradient constancy assumption, where image gradients are assumed to remain constant between frames. The motion tensor encodes motion through temporal derivatives and frame-averaged spatial derivatives between f1 and f2.

Parameters:

f1 (ndarray) – Reference frame (2D array).
f2 (ndarray) – Moving frame (2D array).
hy (float) – Spatial grid spacing in y-direction.
hx (float) – Spatial grid spacing in x-direction.

Returns:

J11, J22, J33, J12, J13, J23 (ndarray) – Motion tensor components used in the optical flow solver. These components encode the linearized gradient constancy constraints with adaptive regularization based on local gradient structure.

Notes

The gradient constancy assumption extends the classic brightness constancy by requiring that spatial gradients also remain constant during motion. This provides additional robustness in textured regions.

The tensor components are computed with adaptive regularization that weights contributions based on local gradient magnitudes, making the estimation robust to noise while preserving motion details.

References

pyflowreg.core.optical_flow.get_motion_tensor_gray(f1, f2, hy, hx)[source]#

Compute motion tensor components for gray-value constancy assumption.

Calculates the motion tensor components that encode the linearized optical flow constraints under the classic brightness (gray-value) constancy assumption. Spatial derivatives are averaged between frames, while the temporal term measures residual brightness changes after warping.

Parameters:

f1 (ndarray) – Reference frame (2D array).
f2 (ndarray) – Moving frame (2D array).
hy (float) – Spatial grid spacing in y-direction.
hx (float) – Spatial grid spacing in x-direction.

Returns:

J11, J22, J33, J12, J13, J23 (ndarray) – Motion tensor components used by the optical flow solver. These components encode the linearized gray-value constancy constraints using averaged spatial gradients and the inter-frame difference.

Notes

Gray-value constancy assumes that pixel intensities remain unchanged between frames. It is sensitive to lighting changes but provides a simple and fast data term that works well when illumination is stable.

Gradients are computed on symmetrically padded images to enforce Neumann boundary conditions, with boundary entries zeroed to avoid wrap-around artifacts.

pyflowreg.core.optical_flow.get_motion_tensor_cs(f1, f2, hy, hx, eps=None)[source]#

Compute motion tensor components for census-based constancy assumption.

Builds a robust motion tensor using a smoothed census transform that matches the relative ordering of neighboring pixels instead of raw intensity values. The tensor aggregates directional differences over a 3x3 neighborhood to enforce illumination-invariant constraints between f1 and f2.

Parameters:

f1 (ndarray) – Reference frame (2D array).
f2 (ndarray) – Moving frame (2D array).
hy (float) – Spatial grid spacing in y-direction.
hx (float) – Spatial grid spacing in x-direction.
eps (float, optional) – Smoothing width for the smoothed Heaviside function applied to directional differences r = (neighbor - center) / dist. If None, uses 0.1 / 255.0, matching the Hafner/Demetz/Weickert epsilon = 0.1 convention for images scaled from [0, 255] to approximately [0, 1]. When hx or hy are physical units rather than pixel-like units, callers may need to scale eps consistently.

Returns:

J11, J22, J33, J12, J13, J23 (ndarray) – Motion tensor components used in the optical flow solver. These components encode the linearized census constancy constraints, averaged over the eight-connected neighborhood for robustness.

Notes

The hard census transform is invariant to global monotonically increasing grey-value transforms because it depends only on ordering. This implementation uses finite differences, Gaussian-preprocessed inputs, a smoothed Heaviside function, and a linearized motion tensor, so invariance is approximate.

Additive offsets cancel exactly in neighbor-center differences. Positive multiplicative changes are approximately handled only when eps is small relative to the directional-difference scale, or when eps is scaled consistently with image intensity scale. Gamma and other nonlinear monotone transforms preserve hard ordering but not exact smoothed Heaviside values.

References

Boundary Handling#

pyflowreg.core.optical_flow.add_boundary(f)[source]#

Add 1-pixel boundary padding with edge replication.

Pads array with 1 pixel on all sides by replicating edge values, implementing Neumann boundary conditions for the flow solver.

Parameters:: f (np.ndarray) – Input array to pad
Returns:: padded (np.ndarray) – Padded array with shape increased by 2 in each dimension

Notes

Uses np.pad with mode=’edge’, which replicates the values at the edges to provide Neumann (zero-derivative) boundary conditions.

Pyramid Level Solver#

Pyramid Level Solver for Variational Optical Flow#

This module implements the computationally intensive solver for optical flow estimation at each pyramid level. All functions are numba-optimized with JIT compilation for high performance.

The solver uses successive over-relaxation (SOR) to iteratively solve the linearized optical flow equations with non-linear diffusion regularization. The implementation follows the generalized Charbonnier penalty framework for robustness to noise and motion discontinuities.

Functions#

compute_flow: Main iterative solver for flow field computation at a pyramid level
set_boundary_2d: Apply Neumann boundary conditions to 2D arrays
nonlinearity_smoothness_2d: Compute non-linearity weights for smoothness term

Notes

All functions use numba JIT compilation with cache=True for performance. The cache is beneficial since these functions are called repeatedly for each image, but it means the functions will fail if dtypes change between calls.

For optimal performance, input arrays should be contiguous in memory. Numba performs significantly better with contiguous arrays and can be slow with non-contiguous data due to inability to use optimized SIMD instructions.

The solver implements coupled updates for the u and v components of the flow field using SOR with OMEGA=1.95.

Flow Computation#

pyflowreg.core.level_solver.compute_flow(J11, J22, J33, J12, J13, J23, weight, u, v, alpha_x, alpha_y, iterations, update_lag, a_data, a_smooth, hx, hy)[source]#

Iterative solver for optical flow at a single pyramid level.

Solves the linearized optical flow equations using successive over-relaxation (SOR) with non-linear diffusion regularization. The solver minimizes an energy functional combining data fidelity (from motion tensor) and smoothness terms, both with generalized Charbonnier penalties.

Parameters:

J11, J22, J33, J12, J13, J23 (ndarray) – Motion tensor components, shape (m, n, n_channels). Encode the linearized gradient constancy constraints.
weight (ndarray) – Channel weights for multi-channel data, shape (m, n, n_channels).
u, v (ndarray) – Initial flow field components from coarser pyramid level, shape (m, n).
alpha_x, alpha_y (float) – Regularization weights for smoothness term in x and y directions.
iterations (int) – Number of SOR iterations to perform.
update_lag (int) – Update non-linearity weights every update_lag iterations.
a_data (ndarray) – Charbonnier exponents for data term, length n_channels.
a_smooth (float) – Charbonnier exponent for smoothness term.
hx, hy (float) – Spatial grid spacing in x and y directions.

Returns:

flow (ndarray) – Computed flow field, shape (m, n, 2) where flow[:,:,0] is u-component and flow[:,:,1] is v-component.

Notes

Uses SOR with relaxation parameter OMEGA=1.95 for convergence acceleration. The solver updates flow increments (du, dv) relative to the input (u, v) using a coupled Gauss-Seidel scheme with immediate u-updates affecting v.

Boundary Conditions#

pyflowreg.core.level_solver.set_boundary_2d(f)[source]#

Apply Neumann boundary conditions to 2D array in-place.

Replicates edge values to provide zero-derivative (Neumann) boundary conditions for the flow solver. This enables one-sided finite differences at boundaries (forward differences at left/top edges, backward differences at right/bottom edges) without going out of bounds.

Parameters:: f (ndarray) – 2D array to apply boundary conditions to (modified in-place).

Nonlinearity Functions#

pyflowreg.core.level_solver.nonlinearity_smoothness_2d(psi_smooth, u, du, v, dv, m, n, a, hx, hy)[source]#

Compute derivative of smoothness penalty (non-linearity weights).

Calculates ψ(s²) = a(s² + ε)^(a-1), the derivative of the Charbonnier smoothness penalty ρ(s²) = (s² + ε)^a, where s² = ||∇u||² + ||∇v||². This non-linearity arises from the variational derivative of the energy functional.

Parameters:

psi_smooth (ndarray) – Output array for smoothness weights (m, n), modified in-place.
u (ndarray) – Current u-component of flow field (m, n).
du (ndarray) – Incremental update to u-component (m, n).
v (ndarray) – Current v-component of flow field (m, n).
dv (ndarray) – Incremental update to v-component (m, n).
m (int) – Height of arrays.
n (int) – Width of arrays.
a (float) – Charbonnier penalty exponent for smoothness term.
hx (float) – Spatial grid spacing in x-direction.
hy (float) – Spatial grid spacing in y-direction.

DISO Backend#

Dense Inverse Search Optical Flow (DIS) backend using OpenCV. Provides an alternative optical flow computation method to the variational approach.

class pyflowreg.core.diso_optical_flow.DisoOF(preset=2, finest_scale=2, gradient_descent_iterations=12, patch_size=8, patch_stride=4, use_mean_normalization=True, use_spatial_propagation=True)[source]#

Bases: object

Dense Inverse Search Optical Flow (DIS) implementation using OpenCV.

This class provides a callable interface compatible with get_displacement for computing optical flow between fixed and moving images.

Uses lazy initialization to ensure pickle compatibility for multiprocessing.

Parameters:

preset (int) – DIS preset mode: - cv2.DISOPTICAL_FLOW_PRESET_ULTRAFAST: Ultrafast preset - cv2.DISOPTICAL_FLOW_PRESET_FAST: Fast preset - cv2.DISOPTICAL_FLOW_PRESET_MEDIUM: Medium preset (default)
finest_scale (int) – Finest scale for the image pyramid (0 = original scale)
gradient_descent_iterations (int) – Number of gradient descent iterations at each pyramid level
patch_size (int) – Size of the patch for matching (default: 8)
patch_stride (int) – Stride between neighbor patches (default: 4)
use_mean_normalization (bool) – Whether to use mean normalization (default: True)
use_spatial_propagation (bool) – Whether to use spatial propagation (default: True)

set_preset(preset)[source]#

Update the DIS preset configuration.

Parameters:: preset (int) – One of cv2.DISOPTICAL_FLOW_PRESET_*

get_params()[source]#

Get current DIS parameters.

Return type:: Dict[str, Any]
Returns:: dict – Dictionary containing current parameter values

set_params(**params)[source]#

Update DIS parameters.

Parameters:: **params (dict) – Parameters to update