MATLAB Function for Pyramid-Based LK Optical Flow

This article presents a complete MATLAB implementation of the Lucas-Kanade (LK) optical flow algorithm using a pyramidal (coarse-to-fine) approach. The function estimates the apparent motion vector field between two successive grayscale image frames.

Main Algorithm Function

function [flowX, flowY] = computePyramidLK(frameA, frameB, numLevels, winSz)
% computePyramidLK Pyramid-based Lucas-Kanade optical flow estimator.
% Inputs:
%   frameA, frameB - Consecutive grayscale image frames (2D matrices).
%   numLevels      - Number of pyramid levels.
%   winSz          - Size of the local analysis window (odd integer).
% Outputs:
%   flowX, flowY   - Horizontal and vertical flow component matrices.

% Convert to double precision for numerical stability.
frameA = im2double(frameA);
frameB = im2double(frameB);

% Build Gaussian image pyramids.
pyramidA = cell(1, numLevels);
pyramidB = cell(1, numLevels);
pyramidA{1} = frameA;
pyramidB{1} = frameB;

for idx = 2:numLevels
    pyramidA{idx} = impyramid(pyramidA{idx-1}, 'reduce');
    pyramidB{idx} = impyramid(pyramidB{idx-1}, 'reduce');
end

% Initialize flow at the coarsest pyramid level.
flowX = zeros(size(pyramidA{numLevels}));
flowY = zeros(size(pyramidA{numLevels}));

% Coarse-to-fine iteration.
for currLevel = numLevels:-1:1
    I1 = pyramidA{currLevel};
    I2 = pyramidB{currLevel};

    % Upscale flow from a coarser level, if not at the top.
    if currLevel < numLevels
        flowX = 2 * imresize(flowX, size(I1), 'bilinear');
        flowY = 2 * imresize(flowY, size(I1), 'bilinear');
    end

    % Warp the second image using the current flow estimate.
    I2_warped = applyWarp(I2, flowX, flowY);

    % Compute incremental flow correction at this level.
    [dX, dY] = estimateSingleLevelFlow(I1, I2_warped, winSz);

    % Update the complete flow estimate.
    flowX = flowX + dX;
    flowY = flowY + dY;
end

% Optional: Visualize the computed flow field.
visualizeFlowField(flowX, flowY, frameA);
end

Supporting Function: Image Warping

function imgOut = applyWarp(imgIn, warpU, warpV)
% applyWarp Warps an input image according to flow vectors (u,v).
[Xgrid, Ygrid] = meshgrid(1:size(imgIn,2), 1:size(imgIn,1));
imgOut = interp2(Xgrid, Ygrid, imgIn, ...
                 Xgrid + warpU, Ygrid + warpV, 'linear', 0);
end

Supporting Function: Single-Level LK Flow Estimation

function [u, v] = estimateSingleLevelFlow(imgRef, imgTarget, winRadius)
% estimateSingleLevelFlow Computes optical flow for one pyramid level.
% Computes spatial and temporal derivatives.
[gradX, gradY] = gradient(imgRef);
gradT = imgTarget - imgRef; % Temporal derivative.

u = zeros(size(imgRef));
v = zeros(size(imgRef));
halfWin = floor(winRadius / 2);

% Iterate over central pixels with a full neighborhood.
for row = halfWin+1:size(imgRef,1)-halfWin
    for col = halfWin+1:size(imgRef,2)-halfWin
        % Extract local window patches.
        patchGx = gradX(row-halfWin:row+halfWin, col-halfWin:col+halfWin);
        patchGy = gradY(row-halfWin:row+halfWin, col-halfWin:col+halfWin);
        patchGt = gradT(row-halfWin:row+halfWin, col-halfWin:col+halfWin);

        % Form the linear system A * d = b.
        A = [patchGx(:), patchGy(:)]; % Design matrix.
        b = -patchGt(:);              % Target vector.

        % Solve via normal equations (A'*A) * d = A' * b.
        if rank(A'*A) >= 2
            delta = (A'*A) \ (A'*b);
            u(row, col) = delta(1);
            v(row, col) = delta(2);
        end
    end
end

% Post-processing: median filtering to reduce outliers.
filterSize = [5 5];
u = medfilt2(u, filterSize);
v = medfilt2(v, filterSize);
end

Supporting Function: Flow Visualization

function visualizeFlowField(U, V, backgroundImg)
% visualizeFlowField Overlays flow vectors on an image.
figure;
imshow(backgroundImg);
hold on;

% Subsample for a clearer visualization.
sampleStep = 10;
[gridX, gridY] = meshgrid(1:sampleStep:size(backgroundImg,2), ...
                          1:sampleStep:size(backgroundImg,1));
U_sampled = U(1:sampleStep:end, 1:sampleStep:end);
V_sampled = V(1:sampleStep:end, 1:sampleStep:end);

% Normalize arrow lengths for consistent display.
maxMag = max(sqrt(U_sampled(:).^2 + V_sampled(:).^2));
if maxMag > 0
    normFactor = 3 / maxMag;
else
    normFactor = 1;
end

quiver(gridX(:), gridY(:), ...
       U_sampled(:)*normFactor, V_sampled(:)*normFactor, ...
       0, 'y', 'LineWidth', 1.5);
title('Optical Flow Vector Field');
hold off;
end

Example Usage

% Load two consecutive video frames.
frame1 = imread('video_frame_001.png');
frame2 = imread('video_frame_002.png');

% Ensure images are grayscale.
if size(frame1, 3) == 3
    frame1 = rgb2gray(frame1);
    frame2 = rgb2gray(frame2);
end

% Set algorithm parameters and compute flow.
pyramidDepth = 4;   % Number of pyramid levels.
windowLength = 17; % Local window size (odd).
[horizontalFlow, verticalFlow] = computePyramidLK(frame1, frame2, ...
                                                  pyramidDepth, windowLength);

Algorithm Breakdown and Notes

1. Pyramidal Framework: A Gaussian pyramid is constructed using impyramid. Flow is computed at the coarsest level and then upsampled to initialize the computation at the next finer level. This allows the algorithm to handle larger displacements.

2. Core LK Method: For each level, the standard LK equations are solved within a local window. The system models brightness constancy, leading to the linear equation shown in the original material.

3. Iterative Warping: At each pyramid level, the second image is warped towards the first using the current flow estimate. The incremental flow (dX, dY) is then computed between the first image and the warped second image. This process helps refine the estimate.

4. Parameter Guidance:

• Pyramid Levels (numLevels): 3 to 5 levels are typical. More levels can help capture very large motions.
• Window Size (winSz): A size like 15x15 or 21x21 is common. Larger windows provide more spatial support but can oversmooth motion boundaries.

5. Potantial Enhancements: The basic implementation can be extended by adding iterative refinement within each pyramid level, employing robust cost functions, performing sub-pixel interpolation, or incorporating global smoothness constraints.