import numpy as np
import matplotlib.pyplot as plt
def compute_disparity(imgL, imgR, num_disparities, block_size):
"""
Computes the disparity map between a pair of rectified stereo images.
Args:
imgL (np.ndarray): Left grayscale image.
imgR (np.ndarray): Right grayscale image.
num_disparities (int): Maximum disparity range to search.
block_size (int): Size of the square block for block matching.
Returns:
np.ndarray: Disparity map with the same dimensions as the input images.
"""
# Ensure the images are numpy arrays and have float32 type
imgL = np.asarray(imgL, dtype=np.float32)
imgR = np.asarray(imgR, dtype=np.float32)
# Get the dimensions of the images
height, width = imgL.shape
# Initialize the disparity map with zeros
disparity_map = np.zeros((height, width), dtype=np.float32)
# Calculate half of the block size for easier indexing
half_block = block_size // 2
# Loop over each pixel in the left image
for y in range(half_block, height - half_block):
for x in range(half_block, width - half_block):
# Define the block region in the left image
blockL = imgL[y - half_block:y + half_block + 1, x - half_block:x + half_block + 1]
# Initialize variables for the best match
min_ssd = float('inf')
best_disparity = 0
# Search for the best disparity in the range [0, num_disparities)
for d in range(num_disparities):
if x - d - half_block < 0:
continue
# Define the block region in the right image
blockR = imgR[y - half_block:y + half_block + 1, x - d - half_block:x - d + half_block + 1]
# Compute the sum of squared differences (SSD) between blocks
ssd = np.sum((blockL - blockR) ** 2)
# Update the best disparity if SSD is smaller
if ssd < min_ssd:
min_ssd = ssd
best_disparity = d
# Assign the best disparity value to the disparity map
disparity_map[y, x] = best_disparity
return disparity_map
# Load grayscale images as numpy arrays
imgL = plt.imread('./data/tsukuba_l.png')
imgR = plt.imread('./data/tsukuba_r.png')
# Ensure the images are 2D arrays (grayscale)
if imgL.ndim == 3:
imgL = imgL[:, :, 0] # Convert to grayscale if needed
if imgR.ndim == 3:
imgR = imgR[:, :, 0] # Convert to grayscale if needed
# Parameters for disparity computation
num_disparities = 16
block_size = 20
# Compute the disparity map
disparity = compute_disparity(imgL, imgR, num_disparities, block_size)
# Visualization of the results
plt.figure(figsize=(10, 3))
# Original image
plt.subplot(1, 2, 1)
plt.imshow(imgL, cmap='gray')
plt.title('Original Image')
plt.axis("off")
# Disparity map
plt.subplot(1, 2, 2)
plt.imshow(disparity, cmap='inferno')
plt.colorbar()
plt.title('Disparity Map')
plt.axis("off")
# Display the plots
plt.show()
