Background

I’m working on an image stitching project using the Poisson blending method. I found a repository on GitHub that deals with the same problem (https://github.com/erman18/poison-blending-stitching). However, when I run the code with the provided flower data, the resulting image is blurry compared to the original images.

Original image (before stitching):
Original image

Result image (after stitching):
Result image

Problem:

The resulting image after Poisson blending appears blurry. I want the resulting image to be as clear as the original without showing seams.

What I have tried:

• I reviewed the expandImage and poissonSolver functions but couldn’t identify why the blurriness occurs.
• Adjusted the resizing factor but saw no improvement.

Coding:

import re
import os
import cv2
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from typing import List
from skimage import img_as_ubyte


def readFiles(folder_path, filter="mask"):
    # Get all files in the folder.
    files = os.listdir(folder_path)
    
    # Get the file paths.
    file_paths = [os.path.join(folder_path, file) for file in files if re.match(rf"^{filter}.*.png$", file, re.I)]
    
    # Sort the file paths.
    file_paths.sort()
    print(f"Found {len(file_paths)} files in {folder_path}")
    print(file_paths)
    return file_paths
    
def resizeImg(src, factor=0.3):
    #percent by which the image is resized
    if int(factor) == 1:
        return src
    # Compute the new dimensions of the image and resize it
    width = int(src.shape[1] * factor)
    height = int(src.shape[0] * factor)
    dsize = (width, height)
    return cv2.resize(src, dsize)
def convertDouble2Uint8(img):
    
    a = max(abs(img.max()), abs(img.min()))
    print(img.min(), img.max(), "a: ", a)

    _min, _max = -a, a
    minMax = (_max - _min)
    img = 2*(img - _min)
    img = (img / minMax) -1
    cv_image = img_as_ubyte(img)
    return cv_image

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

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

# Expand the image to create symetric 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 symetric to avoid boundary problems
    # My solution: function (g'(-x) = g'(x))
    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 frec
    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 nonzeros 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#/255.0

def computeGrad(img):
    # Compute gradient of the two images.
    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[:,0:-2,:].astype('double'))
    out_y[1:-1,:,:] = 0.5*(img[2:,:,:].astype('double') - img[0:-2,:,:].astype('double'))
    
    return out_x, out_y

def computeGuidanceField(grads, masks):
    # Compute the guidance field.
    guide_x = grads[0][0] * masks[0]
    guide_y = grads[0][1] * masks[0]
    for i in range(1, len(grads)):
        gradc_x = grads[i][0] * masks[i]
        gradc_y = grads[i][1] * masks[i]
        guide_x = np.where(np.abs(guide_x) >= np.abs(gradc_x), guide_x, gradc_x)
        guide_y = np.where(np.abs(guide_y) >= np.abs(gradc_y), guide_y, gradc_y)


    print("Guide X: ", guide_x.shape)
    return guide_x, guide_y

def computeGuidanceFieldAvg(grads, masks):
    # Compute the guidance field.
    guide_x = grads[0][0] * masks[0]
    guide_y = grads[0][1] * masks[0]
    mask_sum = masks[0]
    for i in range(1, len(grads)):
        guide_x += grads[i][0] * masks[i]
        guide_y += grads[i][1] * masks[i]
        mask_sum += masks[i]

    print("Guide X: ", guide_x.shape)
    return divAB(guide_x, mask_sum, fill=0), divAB(guide_y, mask_sum, fill=0)

def main_merge(img_paths: List[str], mask_paths: List[str], scf=0.5) -> None:    
    # Read images and masks.
    imgs = [resizeImg(cv2.imread(img_path), scf) for img_path in img_paths]  
    masks = [resizeImg(cv2.imread(mask_path)//255, scf) for mask_path in mask_paths]
    # masks = [~(resizeImg(cv2.imread(mask_path)//255, scf) == 0) for mask_path in mask_paths]

    # Compute the sum of the masks.
    # Taking a matrix of size 5 as the kernel
    kernel = np.ones((15, 15), np.uint8)
    T = np.zeros_like(imgs[0])
    d_masks = []
    for i, (img, mask) in enumerate(zip(imgs, masks)):
        mask = mask # cv2.dilate(mask, kernel, iterations=1)
        d_masks.append(mask)
        T += img * mask

    # Compute gradient of the images.
    grads = [(computeGrad(img)) for img in imgs]

    # Compute the guidance field.
    guide_x, guide_y = computeGuidanceField(grads, d_masks)
 
    # Compute the Poisson equation.
    img_out = poissonSolver(guide_x, guide_y)
    
    # Convert the double image to uint8.
    img_out = convertDouble2Uint8(img_out)

    return T,guide_x, guide_y, img_out
   
m_path = "poison-blending-stitching/data/flower"
img_paths = readFiles(m_path, "tks")
mask_paths = readFiles(m_path, "mask")

T,guide_x, guide_y, img_out = main_merge(img_paths, mask_paths)

img2 = cv2.cvtColor(img_out, cv2.COLOR_BGR2RGB)
plt.imshow(img2) 
plt.axis("off")

The data: Data link

Question:

Why does the Poisson blending method result in a blurry image, and how can I improve the clarity of the resulting image?

Thank you for your help.