[Python] Poisson Blending in Image Stitching Results in Blurred, Ghostly Images

I'm working on an image stitching project using Poisson blending to seamlessly blend overlapping images. While the blending quality is generally good, I'm encountering an issue where the resulting stitched image has blurred layers the same size as the masks, leading to a ghostly, blurry result.

Problem:

When stitching multiple images using their corresponding homographies, the final panorama shows ghosting and blurring in the areas where images overlap.

Linear blending result

Poisson blending result

import cv2
import numpy as np
import os
import matplotlib.pyplot as plt
from skimage import img_as_ubyte

def process_images_and_masks(images: list[np.ndarray], homographies: list[np.ndarray]) -> tuple[list[np.ndarray], list[np.ndarray]]:
"""
Process images with their homographies and return the warped images and masks.
"""
def get_transformed_corners(image, H):
h, w = image.shape[:2]
corners = np.float32([[0, 0], [0, h-1], [w-1, h-1], [w-1, 0]]).reshape(-1, 1, 2)
transformed_corners = cv2.perspectiveTransform(corners, H)
return transformed_corners

# Calculate the bounds of the canvas
all_corners = []
for image, H in zip(images, homographies):
transformed_corners = get_transformed_corners(image, H)
all_corners.append(transformed_corners)

all_corners = np.concatenate(all_corners, axis=0).reshape(-1, 2)
min_x, min_y = np.min(all_corners, axis=0).astype(int)
max_x, max_y = np.max(all_corners, axis=0).astype(int)

canvas_width = max_x - min_x
canvas_height = max_y - min_y

# Adjust the homographies to fit all images within the canvas
translation = np.array([[1, 0, -min_x], [0, 1, -min_y], [0, 0, 1]])
adjusted_homographies = [translation @ H for H in homographies]

new_images = []
visible_masks = []

for image, H in zip(images, adjusted_homographies):
warped_image = cv2.warpPerspective(image, H, (canvas_width, canvas_height))

# Create a binary mask for the original image
mask = np.ones(image.shape[:2], dtype=np.uint8) * 255

# Warp the mask using the same homography
warped_mask = cv2.warpPerspective(mask, H, (canvas_width, canvas_height))

# Append the results to the lists
new_images.append(warped_image)
visible_masks.append(warped_mask)

return new_images, visible_masks

def divAB(a, b, fill=np.nan):
""" Element-wise division of a by b, handling division by zero. """
with np.errstate(divide='ignore', invalid='ignore'):
c = np.true_divide(a, b)
if np.isscalar(c):
return c if np.isfinite(c) else fill
else:
c[~np.isfinite(c)] = fill
return c

def computeGrad(img):
# Compute gradient of the image.
out_x = np.zeros(img.shape, np.double)
out_y = np.zeros(img.shape, np.double)
out_x[:, 1:-1, :] = 0.5 * (img[:, 2:, :].astype('double') - img[:, :-2, :].astype('double'))
out_y[1:-1, :, :] = 0.5 * (img[2:, :, :].astype('double') - img[:-2, :, :].astype('double'))
return out_x, out_y

def computeGuidanceField(grads, masks):
# Convert masks to float64 and normalize to [0, 1]
masks = [mask[:, :, np.newaxis].astype(np.float64) / 255.0 for mask in masks]

# Stack the gradients and masks
grad_x_stack = np.stack([grad[0] for grad in grads], axis=0) # Shape: (num_images, H, W, C)
grad_y_stack = np.stack([grad[1] for grad in grads], axis=0)
mask_stack = np.stack(masks, axis=0) # Shape: (num_images, H, W, 1)

# Sum of masks at each pixel
mask_sum = np.sum(mask_stack, axis=0)

# Avoid division by zero
mask_sum[mask_sum == 0] = 1

# Normalize masks to sum to 1 at each pixel
normalized_masks = mask_stack / mask_sum

# Weighted sum of gradients
guide_x = np.sum(grad_x_stack * normalized_masks, axis=0)
guide_y = np.sum(grad_y_stack * normalized_masks, axis=0)

return guide_x, guide_y

def convertDouble2Uint8(img):
img = cv2.normalize(img, None, 0, 255, cv2.NORM_MINMAX)
return img.astype(np.uint8)

def mft(U):
# Define a shortcut for the Fourier transform.
return np.fft.fftshift(np.fft.fft2(U))

def imft(U):
# Define a shortcut for the inverse Fourier transform.
return np.fft.ifft2(np.fft.ifftshift(U)).real

# Expand the image to create symmetry to avoid boundary problems
def expandImage(original_image, d="x"):
# Get the dimensions of the original image
rows, cols, ch = np.shape(original_image)

# Create an empty array to hold the final image
final_image = np.empty((rows * 2, cols * 2, ch), dtype=original_image.dtype)

if d == "x":
# Place the original and mirror images in the final image
final_image[0:rows, 0:cols, :] = original_image
final_image[0:rows, cols:cols*2, :] = -np.flip(original_image, axis=1)
final_image[rows:rows*2, 0:cols, :] = np.flip(original_image, axis=0)
final_image[rows:rows*2, cols:cols*2, :] = -np.flip(original_image, axis=(0,1))
elif d == "y":
# Place the original and mirror images in the final image
final_image[0:rows, 0:cols, :] = original_image
final_image[0:rows, cols:cols*2, :] = np.flip(original_image, axis=1)
final_image[rows:rows*2, 0:cols, :] = -np.flip(original_image, axis=0)
final_image[rows:rows*2, cols:cols*2, :] = -np.flip(original_image, axis=(0,1))

return final_image

def poissonSolver(gx, gy):
print("poissonSolver")
# Initialization of the output image
I = np.zeros(gx.shape)

# Expand the image to create symmetry to avoid boundary problems
gx = expandImage(gx, "x")
gy = expandImage(gy, "y")

H,W,C = gx.shape

# Define frequency domain
wx, wy = np.meshgrid(np.arange(1, W+1, 1), np.arange(1, H+1, 1))

wx0 = np.floor(W/2)+1
wy0 = np.floor(H/2)+1 # Zero frequency
wx = wx - wx0
wy = wy - wy0

cx = ((1j*2*np.pi)/W)*wx
cy = ((1j*2*np.pi)/H)*wy
d = (cx)**2 + (cy)**2
print("---", gx.shape)

print(f"{'-->':>4} Print zeros : {np.argwhere(np.abs(d) == 0)}, Center: ({int(wy0)}, {int(wx0)})")

del wx, wy

for c in range(0, C):
Gx = gx[:,:,c]
Gy = gy[:,:,c]

Vx = (cx) * mft(Gx)
Vy = (cy) * mft(Gy)

# FT_I = ( Vx + Vy ) / ( d )
FT_I = np.zeros_like(Vx)
np.divide(( Vx + Vy ), d, out=FT_I, where=d!=0) # Only divide non-zeros else 1

FT_I[int(wy0-1), int(wx0-1)] = 10 # Set DC value (undefined in the previous div.)

Aux = imft(FT_I)
I[:,:,c] = Aux[0:int(H/2), 0:int(W/2)] # Keep the original portion of the space.

cv2.normalize(I, I, 0, 1, cv2.NORM_MINMAX)
return I

def poisson_blending(images_list, homographies_list) -> np.ndarray:
new_images, visible_masks = process_images_and_masks(images_list, homographies_list)
grads = [computeGrad(img) for img in new_images]
guide_x, guide_y = computeGuidanceField(grads, visible_masks)

img_out = poissonSolver(guide_x, guide_y)
img_out = convertDouble2Uint8(img_out)

return guide_x, guide_y, img_out

# Directory where images and H matrices are stored
output_dir = 'samples/output_images'

# List to store images and homographies
images = []
homographies = []

# Load the images
image_filenames = sorted([f for f in os.listdir(output_dir) if f.endswith('.png')])
for image_filename in image_filenames:
image_path = os.path.join(output_dir, image_filename)
image = cv2.imread(image_path)
images.append(image)

# Load the H matrices (homographies)
H_matrices_filename = os.path.join(output_dir, 'H_matrices.npy')
homographies = np.load(H_matrices_filename, allow_pickle=True)

guide_x, guide_y, img_out = poisson_blending(images, homographies)

# Display the resulting image
plt.axis("off") # Hide axis
plt.imshow(cv2.cvtColor(img_out, cv2.COLOR_BGR2RGB)) # Convert from BGR to RGB for correct display
plt.show()



Additional Information: For each image, there is a corresponding homography that positions the image in 2D space. The images and homographies are loaded from the samples/output_images directory. The masks are created as binary masks (all ones) and warped using the homographies.

Data link: output_images - Google Drive

How can I adjust the Poisson blending process to eliminate the blurring and achieve a seamless panorama?

Continue reading...