Compare RAFT Optical Flow and Semi-Global Matching for Stereo Reconstruction

Since R2026a

This example shows how to estimate disparity between a pair of stereo images using optical flow from the RAFT deep learning model [1], and compares the results with a reconstruction created using a traditional computer vision technique: semi-global matching (SGM) [2].

Stereo disparity refers to the pixel offset between corresponding points in a pair of stereo images, which comprise the left and right views of the same scene. You can determine a dense depth map using the stereo disparity from a calibrated camera, which is useful in applications such as dense 3-D reconstruction. For more applications of stereo disparity, see Stereo Visual SLAM for UAV Navigation in 3D Simulation and Depth Estimation from Stereo Video.

Traditional computer vision techniques, such as SGM [2] achieve good results on well-textured scenes, but face challenges in scenarios with large uniform regions with no texture, slanted surfaces that are not fronto-parallel to the camera, as well as occlusions and motion-boundaries between the two views.

Accurate optical flow between two images provides dense pixel correspondences and determines the apparent motion between two views [3][4], which provides an estimate of disparity. Optical flow from RAFT is robust even in the presence of textureless regions and motion blur, and can provide more accurate estimates of disparity when traditional techniques cannot provide high-quality results.

Load and Display Stereo Pair

Load the stereo parameters.

load("webcamsSceneReconstruction.mat");

Read in the stereo pair of images.

I1 = imread("sceneReconstructionLeft.jpg");
I2 = imread("sceneReconstructionRight.jpg");

Rectify the images. Rectified images have horizontal epipolar lines, and are row-aligned.

[imLeft, imRight, reprojectionMatrix] = rectifyStereoImages(I1,I2,stereoParams);

Display the two images side-by-side.

figure
montage({imLeft, imRight})
title("Rectified stereo pair image")

Create the stereo anaglyph of the rectified stereo image pair and display it. Observe that the displacement between the two images, after stereo rectification, is entirely due to horizontal camera translation.

A = stereoAnaglyph(imLeft, imRight);
figure
imshow(A)
title("Red-Cyan composite view of the rectified stereo pair image")

Compute Disparity Map Using opticalFlowRAFT

Create an opticalFlowRAFT object.

flowModel = opticalFlowRAFT;

Compute and visualize the optical flow between the stereo image pair.

estimateFlow(flowModel,imLeft);
flow = estimateFlow(flowModel,imRight);
figure
imshow(imLeft)
hold on
plot(flow,DecimationFactor=[40 40],ScaleFactor=0.75,color="g");
hold off
title("Plot of Optical Flow Vectors")

Estimate Disparity Map from Optical Flow

Observe that the magnitude of each optical flow vector is inversely proportional to the distance of the associated point from the camera. Assuming perfectly rectified stereo images, the only displacement of pixels should be along the horizontal axis. The magnitude of the horizontal component of optical flow abs(flow.Vx) thus gives an estimate of the stereo disparity [3][4].

disparityMapRAFT = abs(flow.Vx);
disparityMapRAFT = imresize(disparityMapRAFT, size(imLeft,[1 2]));

Refine Estimated Disparity Map

Discard pixels where the horizontal flow component is smaller than one pixel, or where the flow vectors exceed image bounds, as these are likely to be noisy predictions. This discards areas of extremely low disparity, as well as points on the image boundaries that are not covisible from both the stereo images. Pixels on the image margins are likely to suffer from noisy flow estimates, so discard these as well.

% Modify thresholds based on specific data set
minDisparity = 1;   % in pixels
imageMargin  = 10;  % in pixels

[H,W,~] = size(imLeft);
[X,Y] = meshgrid(1:W,1:H);
X2 = X + imresize(flow.Vx, size(imLeft,[1 2]));
Y2 = Y + imresize(flow.Vy, size(imLeft,[1 2]));

% Determine valid optical flow mask
validFlow = X2>=1 & X2<=W & Y2>=1 & Y2<=H & disparityMapRAFT > minDisparity;

% Replace invalid positions in the disparity map with nan
disparityMapRAFT(~validFlow) = nan;

% Mark image borders as invalid
disparityMapRAFT(:,end-imageMargin:end) = nan;
disparityMapRAFT(end-imageMargin:end,:) = nan;
disparityMapRAFT(1:imageMargin,:) = nan;
disparityMapRAFT(:,1:imageMargin) = nan;

figure
imagesc(disparityMapRAFT)
colorbar
axis image
title("Refined Disparity Map from Optical Flow")

Observe that areas around the image borders and locations of large optical flow magnitudes, that are likely to be noisy, are marked as invalid.

Compute Disparity Map Using disparitySGM

Compute the disparity map using the semi-global matching (SGM) algorithm, which is a classical method for computing the stereo disparity between a pair of rectified images.

imLeftGray = im2gray(imLeft);
imRightGray = im2gray(imRight);
disparityMapSGM = disparitySGM(imLeftGray,imRightGray);

Compare Disparity Maps from opticalFlowRAFT and disparitySGM

Display the disparity maps obtained using opticalFlowRAFT and disparitySGM.

figure

subplot(2,1,1)
imagesc(disparityMapRAFT)
colorbar
axis image
title("Disparity Map using RAFT Optical Flow")

subplot(2,1,2)
imagesc(disparityMapSGM)
colorbar
axis image
title("Disparity Map using Semi-Global Matching")

The disparity map from the disparitySGM method is noisy in several regions, and does not accurately capture the detailed spatial structure of the scene. The disparity map obtained using RAFT optical flow is smoother and more accurate, especially at capturing the finer structures and sharp edges. You can use opticalFlowRAFT to obtain high-quality estimates of stereo disparity in scenarios where disparitySGM does not perform adequately, at the cost of more time and higher compute requirements, such as a GPU.

Reconstruct 3-D Scene from Disparity Maps

Use the estimated stereo disparity maps to reconstruct the scene. The reconstructScene function returns a set of 3-D coordinates for each pixel as an H-by-W-by-3 array, for an image of dimensions H and W.

xyzPointsRAFT = reconstructScene(disparityMapRAFT,reprojectionMatrix);
xyzPointsSGM  = reconstructScene(disparityMapSGM,reprojectionMatrix);

Get depth maps from the reconstructed 3-D points for both methods, and display them. Use a log-scale on the depth map for better visibility of details.

depthMapRAFT = xyzPointsRAFT(:,:,3);
depthMapSGM  = xyzPointsSGM(:,:,3);

figure
subplot(2,1,1)
imagesc(log(1 + depthMapRAFT))
axis image
colorbar
title("Depth Map (log-scale) - RAFT")

subplot(2,1,2)
imagesc(log(1 + depthMapSGM))
axis image
colorbar
title("Depth Map (log-scale) - SGM")

The depth map is inversely proportional to the disparity map: areas closer to the camera, such as the floor in the foreground, have a higher disparity value, and a lower corresponding depth value. The results from opticalFlowRAFT are smoother and have fewer holes, compared to the results from disparitySGM.

Match depth values with RGB pixel colors from the left image.

[H,W,~] = size(imLeft);
[X,Y] = meshgrid(1:W,1:H);
u = X(:);
v = Y(:);
linearIdx = sub2ind([H,W], v, u);
R = imLeft(:,:,1);
G = imLeft(:,:,2);
B = imLeft(:,:,3);
colors = [R(linearIdx), G(linearIdx), B(linearIdx)];

Create a pointCloud object using colors and the 3-D point locations.

ptCloudRAFT = pointCloud(reshape(xyzPointsRAFT, [H*W 3]), Color=colors);
ptCloudSGM  = pointCloud(reshape(xyzPointsSGM, [H*W 3]), Color=colors);

Filter outlier points by setting thresholds along the X, Y and Z axes. Display the thresholded point clouds obtained from stereo reconstruction. The stereo reconstruction is in a metric scale of millimeters.

% Select thresholds using the sliders to remove outliers and retain 
% objects of interest, by visually inspecting the visualized point cloud.
minXThreshold = -3000;
maxXThreshold = 5000;
minYThreshold = -2000;
maxYThreshold = 2000;
maxZThreshold = 6200;
roi = [minXThreshold maxXThreshold minYThreshold  maxYThreshold 0  maxZThreshold];

% Apply selected thresholds on the point clouds
indicesRAFT = findPointsInROI(ptCloudRAFT,roi);
ptCloudRAFTROI = select(ptCloudRAFT,indicesRAFT);
indicesSGM = findPointsInROI(ptCloudSGM,roi);
ptCloudSGMROI = select(ptCloudSGM,indicesSGM);

% Visualize the output point clouds
figure
subplot(2,1,1)
pcshow(ptCloudRAFTROI, VerticalAxis="y", VerticalAxisDir="down");
xlabel("X")
ylabel("Y")
zlabel("Z")
title("Stereo Reconstruction using RAFT")
view(0,-80)

subplot(2,1,2)
pcshow(ptCloudSGMROI, VerticalAxis="y", VerticalAxisDir="down");
xlabel("X")
ylabel("Y")
zlabel("Z")
title("Stereo Reconstruction using SGM")
view(0,-80)

The reconstructed point cloud from RAFT has significantly less noise compared to the point cloud obtained using SGM. Using dense optical flow from the RAFT model provides higher quality metric stereo reconstruction, but at the cost of a longer compute time and greater hardware requirements, such as a GPU.

References

[1] Teed, Z., & Deng, J. "RAFT: Recurrent All-Pairs Field Transforms for Optical Flow." European Conference on Computer Vision, 2020.

[2] Hirschmuller, H. "Accurate and Efficient Stereo Processing by Semi-Global Matching and Mutual Information." IEEE Conference on Computer Vision and Pattern Recognition, 2005.

[3] Longuet-Higgins, C. H., & Prazdny K. "The interpretation of a moving retinal image." Proceedings of the Royal Society of London. Series B. Biological Sciences, 1980.

[4] Bruss, A. R., & Horn B.K.P. "Passive navigation." Computer Vision, Graphics, and Image Processing, 1983.