Poisson Blending in Image Stitching Results in Ghostly Images

[khoatran02]

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

[LiteBrite82]

Thanks for posting in the GitHub Community, @khoatran02 !

We’ve moved your post to our Programming Help 🧑‍💻 category, which is more appropriate for this type of discussion.

Please review our guidelines about the Programming Help category for more information.

[khoatran02]

Thank you @LiteBrite82