Hello, I need your help to fix my project. I will be giving you details and the code in quotes and explaining the error that needs fixing. First of all, here is the premise of the project:

"
In this project, we aim to create a solution in Python for merging sub-images by using keypoint description methods (SIFT, SURF, and ORB) and obtain a final panorama image. First of all, we will extract and obtain multiple keypoints from sub-images by using the keypoint description method. Then we will compare and match these key points to merge sub-images into one panorama image. As a dataset, we will use a subset of the HPatches dataset. With this dataset, you get 6 ".png" images and 5 files for ground truth homography.
"

Here are more details about the dataset:

"
There is a reference image (image number 0) and five target images taken under different illuminations and from different viewpoints. For all images, we have the estimated ground truth homography with respect to the reference image.
"

As you can understand, the dataset includes a scene consisting of a reference image and various sub-images with the estimated ground truth homography with respect to the reference image.

The implementation has some restrictions. Here are the implementation details we used to create the project:

"
1. Feature Extraction: We are expected to extract key points in the sub-images by a keypoint extraction method. (SIFT, SURF, and ORB). You can use libraries for this part.

2. Feature Matching: Then we are expected to code a matching function (this can be based on the k-nearest neighbor method) to match extracted keypoints between pairs of sub-images. You can use libraries for this part.

3. Finding Homography: Then you should calculate a Homography Matrix for each pair of sub-images (by using the RANSAC method). For this part, you cannot use OpenCV or any library other than NumPY.

4. Merging by Transformation: Merge sub-images into a single panorama by applying transformation operations to sub-images using the Homography Matrix. For this part, you cannot use OpenCV or any library other than NumPY.
"

With that being said, I hope you understand the main goal of the project. Now the issue I am facing is that the end result seems like a warped image of the last image (sixth image named image5.png). I have setup the tool to test with ground truth homographs however still the result image is not a panorama image. It seems the merge_images function is either completely wrong or has terrible mistakes. Please help me fix the problem. Now I will provide the full code so you can check it and tell me how to fix the project:

"
import numpy as np
import cv2
import os
import glob
import matplotlib.pyplot as plt
import time


def feature_extraction(sub_images, method="SIFT"):
    if method == "SIFT":
        keypoint_extractor = cv2.xfeatures2d.SIFT_create()
    elif method == "SURF":
        keypoint_extractor = cv2.xfeatures2d.SURF_create()
    elif method == "ORB":
        keypoint_extractor = cv2.ORB_create()

    keypoints = []
    descriptors = []

    for sub_image in sub_images:
        keypoint, descriptor = keypoint_extractor.detectAndCompute(sub_image, None)
        keypoints.append(keypoint)
        descriptors.append(descriptor)

    return keypoints, descriptors


def feature_matching(descriptors, matcher_type="BF"):
    if matcher_type == "BF":
        matcher = cv2.BFMatcher()

    matches = []
    for i in range(1, len(descriptors)):
        match = matcher.knnMatch(descriptors[0], descriptors[i], k=2)
        matches.append(match)
        
    return matches

def compute_homography_matrix(src_pts, dst_pts):
    def normalize_points(pts):
        pts_homogeneous = np.hstack((pts, np.ones((pts.shape[0], 1))))
        centroid = np.mean(pts, axis=0)
        scale = np.sqrt(2) / np.mean(np.linalg.norm(pts - centroid, axis=1))
        T = np.array([[scale, 0, -scale * centroid[0]], [0, scale, -scale * centroid[1]], [0, 0, 1]])
        normalized_pts = (T @ pts_homogeneous.T).T
        return normalized_pts[:, :2], T

    src_pts_normalized, T1 = normalize_points(src_pts)
    dst_pts_normalized, T2 = normalize_points(dst_pts)

    A = []
    for p1, p2 in zip(src_pts_normalized, dst_pts_normalized):
        x1, y1 = p1
        x2, y2 = p2
        A.append([0, 0, 0, -x1, -y1, -1, y2 * x1, y2 * y1, y2])
        A.append([x1, y1, 1, 0, 0, 0, -x2 * x1, -x2 * y1, -x2])

    A = np.array(A)
    try:
        _, _, VT = np.linalg.svd(A)
    except np.linalg.LinAlgError:
        return None

    h = VT[-1]

    H_normalized = h.reshape(3, 3)
    H = np.linalg.inv(T2) @ H_normalized @ T1
    return H / H[-1, -1]

def filter_matches(matches, ratio_thres=0.7):
    filtered_matches = []
    for match in matches:
        good_match = []
        for m, n in match:
            if m.distance < ratio_thres * n.distance:
                good_match.append(m)
        filtered_matches.append(good_match)
        
    return filtered_matches


def find_homography(keypoints, filtered_matches):
    homographies = []
    skipped_indices = []  # Keep track of skipped images and their indices
    
    for i, matches in enumerate(filtered_matches):
        src_pts = np.float32([keypoints[0][m.queryIdx].pt for m in matches]).reshape(-1, 1, 2)
        dst_pts = np.float32([keypoints[i + 1][m.trainIdx].pt for m in matches]).reshape(-1, 1, 2)
        
        H = ransac_homography(src_pts, dst_pts)
        
        if H is not None:
            H = H.astype(np.float32)
            homographies.append(H)
        else:
            print(f"Warning: Homography computation failed for image pair (0, {i + 1}). Skipping.")
            skipped_indices.append(i + 1) # Add indices of skipped images to the list
            continue
        
    return homographies, skipped_indices

def ransac_homography(src_pts, dst_pts, iterations=2000, threshold=3):
    best_inlier_count = 0
    best_homography = None

    for _ in range(iterations):
        indices = np.random.choice(len(src_pts), 4, replace=True)
        src_subset = src_pts[indices].reshape(-1, 2)
        dst_subset = dst_pts[indices].reshape(-1, 2)

        homography = compute_homography_matrix(src_subset, dst_subset)
        if homography is None:
            continue

        inliers = 0
        for i in range(len(src_pts)):
            projected_point = np.dot(homography, np.append(src_pts[i], 1))
            projected_point = projected_point / projected_point[-1]
            distance = np.linalg.norm(projected_point - np.append(dst_pts[i], 1))

            if distance < threshold:
                inliers += 1

        if inliers > best_inlier_count:
            best_inlier_count = inliers
            best_homography = homography

    return best_homography

def read_ground_truth_homographies(dataset_path):
    H_files = sorted(glob.glob(os.path.join(dataset_path, "H_*")))
    ground_truth_homographies = []
    
    for filename in H_files:
        H = np.loadtxt(filename)
        ground_truth_homographies.append(H)
    
    return ground_truth_homographies

def warp_perspective(img, H, target_shape):
    h, w = target_shape
    target_y, target_x = np.meshgrid(np.arange(h), np.arange(w))
    target_coordinates = np.stack([target_x.ravel(), target_y.ravel(), np.ones(target_x.size)])

    source_coordinates = np.dot(np.linalg.inv(H), target_coordinates)
    source_coordinates /= source_coordinates[2, :]

    valid = np.logical_and(np.logical_and(0 <= source_coordinates[0, :], source_coordinates[0, :] < img.shape[1] - 1),
                            np.logical_and(0 <= source_coordinates[1, :], source_coordinates[1, :] < img.shape[0] - 1))
    
    # Change the dtype to int instead of float64 to save memory
    valid_source_coordinates = np.round(source_coordinates[:, valid].astype(np.float32)[:2]).astype(int)
    valid_target_coordinates = target_coordinates[:, valid].astype(int)[:2]

    valid_source_coordinates[0] = np.clip(valid_source_coordinates[0], 0, img.shape[1] - 1)
    valid_source_coordinates[1] = np.clip(valid_source_coordinates[1], 0, img.shape[0] - 1)
    valid_target_coordinates[0] = np.clip(valid_target_coordinates[0], 0, w - 1)
    valid_target_coordinates[1] = np.clip(valid_target_coordinates[1], 0, h - 1)

    warped_image = np.zeros((h, w, 3), dtype=np.uint8)
    for i in range(3):
        warped_image[..., i][valid_target_coordinates[1], valid_target_coordinates[0]] = img[..., i][valid_source_coordinates[1], valid_source_coordinates[0]]

    return warped_image

def merge_images(sub_images, homographies, skipped_indices):
    ref_img = sub_images[0]

    for i, H_i in enumerate(homographies):
        if i + 1 in skipped_indices:
            print(f"Image {i + 1} was skipped due to homography computation failure.")
            continue

        img_i = sub_images[i + 1]

        # Get corners of the second (transferred) image
        h, w, _ = img_i.shape
        corners = np.array([[0, 0, 1], [w-1, 0, 1], [w-1, h-1, 1], [0, h-1, 1]], dtype=np.float32)
        corners_transformed = H_i @ corners.T
        corners_transformed = corners_transformed / corners_transformed[2]
        corners_transformed = corners_transformed.T[:, :2]

        # Calculate size of the stitched image
        x_min = int(min(corners_transformed[:, 0].min(), 0))
        x_max = int(max(corners_transformed[:, 0].max(), ref_img.shape[1]))
        y_min = int(min(corners_transformed[:, 1].min(), 0))
        y_max = int(max(corners_transformed[:, 1].max(), ref_img.shape[0]))

        # Get the transformation to shift the origin to min_x, min_y
        shift_mtx = np.array([[1, 0, -x_min],
                              [0, 1, -y_min],
                              [0, 0, 1]], dtype=np.float32)

        # Apply the transformation on both images
        img_i_transformed = warp_perspective(img_i, shift_mtx @ H_i, (x_max-x_min, y_max-y_min))
        ref_img_transformed = warp_perspective(ref_img, shift_mtx, (x_max-x_min, y_max-y_min))

        # Calculate blend masks
        mask_img_i = (cv2.cvtColor(img_i_transformed, cv2.COLOR_BGR2GRAY) > 0).astype(np.float32)
        mask_ref_img = (cv2.cvtColor(ref_img_transformed, cv2.COLOR_BGR2GRAY) > 0).astype(np.float32)
        mask_overlap = np.logical_and(mask_img_i, mask_ref_img).astype(np.float32)

        # Normalize weights and combine images
        total_weight = mask_img_i + mask_ref_img + 2 * mask_overlap
        ref_img = ((mask_img_i[..., None] * img_i_transformed) + (mask_ref_img[..., None] * ref_img_transformed) + (2 * mask_overlap[..., None] * img_i_transformed * ref_img_transformed)) / total_weight[..., None]

        # Crop black regions
        gray = cv2.cvtColor(ref_img.astype(np.uint8), cv2.COLOR_BGR2GRAY)
        _, thresh = cv2.threshold(gray, 1, 255, cv2.THRESH_BINARY)
        contours, _ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)[-2:]
        cnt = max(contours, key=cv2.contourArea)
        x, y, w, h = cv2.boundingRect(cnt)
        ref_img = ref_img[y:y+h, x:x+w]

    ref_img = ref_img.astype(np.uint8)
    return ref_img


def main(dataset_path):
    filenames = sorted(glob.glob(os.path.join(dataset_path, "*.png")))
    sub_images = []

    for filename in filenames:
        img = cv2.imread(filename, cv2.IMREAD_COLOR) # Load images as color
        sub_images.append(img)
        
    ground_truth_homographies = read_ground_truth_homographies(dataset_path)
    skipped_indices = []
    panorama = merge_images(sub_images, ground_truth_homographies, skipped_indices)
    plt.figure()
    plt.imshow(panorama, cmap="gray", aspect="auto")
    plt.title(f"Panorama - Test")
    return

    methods = ["SIFT", "SURF", "ORB"]
    for method in methods:
        start_time = time.time()
        
        keypoints, descriptors = feature_extraction(sub_images, method=method)
        matches = feature_matching(descriptors)
        filtered_matches = filter_matches(matches)
        homographies, skipped_indices = find_homography(keypoints, filtered_matches)
        
        panorama = merge_images(sub_images, homographies, skipped_indices)
        
        end_time = time.time()
        runtime = end_time - start_time

        print(f"Method: {method} - Runtime: {runtime:.2f} seconds")

        for idx, (image, kp) in enumerate(zip(sub_images, keypoints)):
            feature_plot = cv2.drawKeypoints(image, kp, None)
            plt.figure()
            plt.imshow(feature_plot, cmap="gray")
            plt.title(f"Feature Points - {method} - Image {idx}")

        for i, match in enumerate(filtered_matches):
            matching_plot = cv2.drawMatches(sub_images[0], keypoints[0], sub_images[i + 1], keypoints[i + 1], match, None, flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)
            plt.figure()
            plt.imshow(matching_plot, cmap="gray")
            plt.title(f"Feature Point Matching - {method} - Image 0 - Image {i + 1}")

        plt.figure()
        plt.imshow(panorama, cmap="gray", aspect="auto")
        plt.title(f"Panorama - {method}")

        print("\nGround truth homographies")
        for i, H_gt in enumerate(ground_truth_homographies):
            print(f"Image 0 to {i+1}:")
            print(H_gt)

        print("\nComputed homographies")
        for i, H_est in enumerate(homographies):
            print(f"Image {i} to {i+1}:")
            print(H_est)

    plt.show()


if __name__ == "__main__":
    dataset_path = "dataset/v_bird"
    main(dataset_path)
"