HomographyCalculator/hgcalc.py

import sys

import numpy as np

def ransac_iters(p, k, z):
    """ Computes the required number of iterations for RANSAC.

    Args:
        p: probability that any given correspondence is valid
        k: number of pairs
        z: total probability of success after all iterations

    Returns:
        minimum number of required iterations
    """

    def log(base, argument):
        """ Computes the logarithm to a given base for a given argument

        :param base: The base of the logarithm.
        :param argument: The argument of the logarithm.
        :return: The logarithm.
        """

        return np.log(argument) / np.log(base)

    """
    Reference: Lecture 8 Slide 23
    Probability that we have no outlier in an iteration: p ** k
    Probability that we have an outlier in an iteration: 1 - p ** k
    Probability that we have an outlier in each iteration: (1 - p ** k) ** number_of_iterations
    (1 - p ** k) ** number_of_iterations should be smaller than or equal to 1 - z
    So number_of_iterations must be at least the logarithm of (1 - z) to the base (1 - p ** k)
    We always round up to the nearest integer using numpy.ceil because we can only execute full loops.
    """
    return int(np.ceil(log(1 - p ** k, 1 - z)))


def ransac(pairs, n_iters, k, threshold):
    """ RANSAC algorithm.

    Args:
        pairs: (l, 4) numpy array containing matched keypoint pairs
        n_iters: number of ransac iterations
        threshold: inlier detection threshold

    Returns:
        H: (3, 3) numpy array, best homography observed during RANSAC
        max_inliers: number of inliers N
        inliers: (N, 4) numpy array containing the coordinates of the inliers
    """

    H = {}
    max_inliers = 0
    inliers = {}

    """
    We need to split pairs in two arrays with two columns so that we can pass it to pick_samples.
    """
    p1, p2 = np.hsplit(pairs, 2)

    """
    We execute the ransac loop n_iters times, so that we have a good chance to have a valid homography.
    """
    for iteration in range(n_iters):
        print(iteration)
        """
        First, we pick a sample of k corresponding point pairs
        """
        sample1, sample2 = pick_samples(p1, p2, k)

        """
        The points in the samples then have to be conditioned, so their values are between -1 and 1.
        We also retrieve the condition matrices from this function.
        """
        sample1_conditioned, T_1 = condition_points(sample1)
        sample2_conditioned, T_2 = condition_points(sample2)

        """
        Using the conditioned coordinates, we calculate the homographies both with respect to the conditioned points
        and the unconditioned points.
        """
        H_current, HC = compute_homography(sample1_conditioned, sample2_conditioned, T_1, T_2)

        """
        We then use the unconditioned points (as clarified in the office hours) to compute the homography distances.
        """
        homography_distances = compute_homography_distance(H_current, p1, p2)

        """
        Using the find_inliers function, we retrieve the inliers for our given homography as well as the number 
        of inliers we currently have.
        """
        number_inliers_current, inliers_current = find_inliers(pairs, homography_distances, threshold)

        """
        If the number of inliers is higher than for a previously computed homography (i.e. our current homography is
        better), we save our current homography as well as the number of inliers and the inliers themselves.
        """
        if number_inliers_current > max_inliers:
            H = H_current
            max_inliers = number_inliers_current
            inliers = inliers_current

    """
    We then return the best homography we found as well as the number of inliers found with it and the inliers 
    themselves.
    """
    return H, max_inliers, inliers


def recompute_homography(inliers):
    """ Recomputes the homography matrix based on all inliers.

    Args:
        inliers: (N, 4) numpy array containing coordinate pairs of the inlier points

    Returns:
        H: (3, 3) numpy array, recomputed homography matrix
    """

    """
    We split the inliers matrix into two arrays with 2 columns each (one array for the points in the left image 
    and one for the points in the right image).
    """
    p1, p2 = np.hsplit(inliers, 2)

    """
    We then condition these points.
    """
    p1_conditioned, T_1 = condition_points(p1)
    p2_conditioned, T_2 = condition_points(p2)

    """
    Using the conditioned points, we recompute the homography. The new homography should give better results than
    the previously computed one, since we only use inliers (correct correspondences) here instead of all 
    correspondences including the wrong ones.
    """
    H, HC = compute_homography(p1_conditioned, p2_conditioned, T_1, T_2)

    return H


def pick_samples(p1, p2, k):
    """ Randomly select k corresponding point pairs.

    Args:
        p1: (n, 2) numpy array, given points in first image
        p2: (m, 2) numpy array, given points in second image
        k:  number of pairs to select

    Returns:
        sample1: (k, 2) numpy array, selected k pairs in left image
        sample2: (k, 2) numpy array, selected k pairs in right image
    """

    """
    We assume that m = n (as confirmed in the office hours)
    """

    """
    Reference: https://numpy.org/doc/stable/reference/random/index.html
    We calculate the number of points n (= m) and generate k unique random integers in [0; n) (random_numbers).
    We then use the random numbers as indices to select points from p1 and the corresponding points from p2.
    """
    n = p1.shape[0]
    generator = np.random.default_rng()
    random_numbers = generator.choice(n, size=k, replace=False)
    # random_numbers = np.array([96, 93, 63, 118])
    sample1 = p1[random_numbers]
    sample2 = p2[random_numbers]

    return sample1, sample2


def condition_points(points):
    """ Conditioning: Normalization of coordinates for numeric stability
    by substracting the mean and dividing by half of the component-wise
    maximum absolute value.
    Further, turns coordinates into homogeneous coordinates.
    Args:
        points: (l, 2) numpy array containing unnormailzed cartesian coordinates.

    Returns:
        ps: (l, 3) numpy array containing normalized points in homogeneous coordinates.
        T: (3, 3) numpy array, transformation matrix for conditioning
    """

    """
    First, we calculate half of the component-wise maximum absolute value for all points, as discussed in the 
    office hours.
    """
    s = np.max(np.abs(points), axis=0) / 2
    sx = s[0]
    sy = s[1]

    """
    We then compute the component-wise mean of all points.
    """
    t = np.mean(points, axis=0)
    tx = t[0]
    ty = t[1]

    """
    We then calculate the transformation matrix T like on slide 18 in Lecture 8.
    We use the x components of s and t in the first row and the y components of s and t in the second row.
    """
    T = np.array(([1 / sx, 0, - tx / sx],
                  [0, 1 / sy, - ty / sy],
                  [0, 0, 1]))

    """
    We now iterate over each point and transform it into homogeneous coordinates.
    Then we transform the point using T and append the result to the list ps.
    Finally, the list ps is converted to a numpy array and returned together with the transformation matrix T.
    """
    ps = []
    for point in points:
        homogeneous_point = np.append(point, 1)
        homogeneous_transformed_point = np.matmul(T, homogeneous_point)
        ps.append(homogeneous_transformed_point)

    return np.array(ps), T


def compute_homography(p1, p2, T1, T2):
    """ Estimate homography matrix from point correspondences of conditioned coordinates.
    Both returned matrices shoul be normalized so that the bottom right value equals 1.
    You may use np.linalg.svd for this function.

    Args:
        p1: (l, 3) numpy array, the conditioned homogeneous coordinates of interest points in img1
        p2: (l, 3) numpy array, the conditioned homogeneous coordinates of interest points in img2
        T1: (3,3) numpy array, conditioning matrix for p1
        T2: (3,3) numpy array, conditioning matrix for p2

    Returns:
        H: (3, 3) numpy array, homography matrix with respect to unconditioned coordinates
        HC: (3, 3) numpy array, homography matrix with respect to the conditioned coordinates
    """

    """
    Reference: Lecture 8 Slide 9, 13 ff.
    We did something similar in Assignment 1
    """

    """
    We calculate the number of rows in p1 (which must be the same as number of rows in p2).
    """
    N = p1.shape[0]
    # Create a (2*N, 9) matrix filled with zeros
    A = np.zeros((2 * N, 9))

    # We iterate from 0 to N (N is the number of correspondences)
    for i in range(0, N):
        # The point from the first image (divided by the last component so that x and y have the non-homogeneous
        # value. We could cut off the last component (which is now 1), but that step is not necessary here.)
        point1 = p1[i] / p1[i][2]

        # The corresponding point from the second image (also divided by the last component to get the
        # non-homogeneous values for x and y)
        point2 = p2[i] / p2[i][2]

        # We always calculate two rows of the matrix A in each loop iteration
        # Here, we calculate the first row according to the lecture slides Lecture 8 slide 14
        first_row = np.concatenate((np.zeros(3), point1, - point2[1] * point1))

        # Calculate the second row as in Lecture 8 slide 14
        second_row = np.concatenate((- point1, np.zeros(3), point2[0] * point1))

        # Calculate the row indices at matrix A, where to store the two calculated rows
        # When we iterate from 0 to N, we get indices 0, 2, 4, 6... etc. for the first row index
        first_row_index = i * 2

        # The second row index is simply the first row index plus one, so the rows are below each other
        second_row_index = first_row_index + 1

        """
        Finally, we save the values in the corresponding rows in matrix A
        """
        A[first_row_index] = first_row
        A[second_row_index] = second_row

    """
    Now we perform a single value decomposition on A and retrieve the last right singular vector
    """
    u, s, vt = np.linalg.svd(A)
    h = vt[-1]

    """
    Since we used conditioned points, we retrieve the homography matrix HC regarding the conditioned coordinates 
    by reshaping the last right singular vector .
    """
    HC = np.reshape(h, (3, 3))

    """
    The matrix is then normalized so that the bottom right value equals 1.
    """
    HC = HC / HC[2, 2]

    """
    We calculate the homography matrix with respect to the unconditioned coordinates by using the formula from 
    lecture 8 on slide 18.
    """
    H = np.matmul(np.linalg.inv(T2), np.matmul(HC, T1))

    """
    The homography matrix is also normalized.
    """
    H = H / H[2, 2]

    return H, HC


def compute_homography_distance(H, p1, p2):
    """ Computes the pairwise symmetric homography distance.

    Args:
        H: (3, 3) numpy array, homography matrix
        p1: (l, 2) numpy array, interest points in img1
        p2: (l, 2) numpy array, interest points in img2

    Returns:
        dist: (l, ) numpy array containing the distances
    """

    """
    Determine the number of points in p1 (should be the same as in p2).
    """
    l = p1.shape[0]

    """
    Calculate the inverse of the homography H.
    """
    H_inverse = np.linalg.inv(H)

    """
    Create a zero-filled array which we later fill with the distances.
    """
    distances = np.zeros(l)

    """
    We calculate the transformed points for the given homography.
    """
    p1_transformed = transform_pts(np.array(p1), H)
    p2_transformed = transform_pts(np.array(p2), H_inverse)

    """
    We calculate the symmetric squared distances using the formula in the assignment sheet for every point pair.
    """
    for index in range(0, l):
        x1 = p1[index]
        x2 = p2[index]
        x1_transformed = p1_transformed[index]
        x2_transformed = p2_transformed[index]
        distances[index] = np.linalg.norm(x1_transformed - x2) ** 2 + np.linalg.norm(x1 - x2_transformed) ** 2

    return distances


def find_inliers(pairs, dist, threshold):
    """ Return and count inliers based on the homography distance.

    Args:
        pairs: (l, 4) numpy array containing keypoint pairs
        dist: (l, ) numpy array, homography distances for k points
        threshold: inlier detection threshold

    Returns:
        N: number of inliers
        inliers: (N, 4)
    """

    inliers_list = []

    """
    We iterate over each pair and use its index in the array and its value.
    We then use the index to retrieve the corresponding distance for the pair in dist. If the distance is smaller
    than or equal to the given threshold, we add the pair to the inliers_list.  
    """
    for pair_index, pair_value in enumerate(pairs):
        if dist[pair_index] <= threshold:
            inliers_list.append(pair_value)

    """
    We then get the number of inliers as the length of the list. 
    Furthermore, we convert the inliers_list to a numpy array.
    Both values are then returned.
    """
    N = len(inliers_list)
    inliers = np.array(inliers_list)

    return N, inliers


def transform_pts(p, H):
    """ Transform p through the homography matrix H.

    Args:
        p: (l, 2) numpy array, interest points
        H: (3, 3) numpy array, homography matrix

    Returns:
        points: (l, 2) numpy array, transformed points
    """

    """
    Reference: Lecture 8 Slide 8
    """

    points_list = []

    """
    Iterate over every point in p.
    For each point, we first append a 1 to transform the point into homogeneous coordinates.
    Then, we perform a matrix multiplication of the point with the homography H, as shown on slide 8 in lecture 8.
    The transformed point is then converted back into non-homogeneous coordinates by dividing it by the last 
    value and then cutting off the last value (which is then 1 again).
    The transformed points are returned as a numpy array.
    """
    for point_index, point_value in enumerate(p):
        homogeneous_p = np.append(point_value, 1)
        homogeneous_transformed_p = np.matmul(H, homogeneous_p)
        normalized_homogeneous_transformed_p = homogeneous_transformed_p / homogeneous_transformed_p[2]
        transformed_p = normalized_homogeneous_transformed_p[0:2]
        points_list.append(transformed_p)

    return np.array(points_list)


def main():
    if len(sys.argv) != 3:
        print('Usage: hgcalc <pairs_file_path> <homography_file_path>')
    else:
        pairs_file_path = sys.argv[1]
        homography_file_path = sys.argv[2]
        # RANSAC Parameters
        ransac_threshold = 0.02  # inlier threshold
        p = 0.35  # probability that any given correspondence is valid
        k = 4  # number of samples drawn per iteration
        z = 0.99  # total probability of success after all iterations

        pairs = np.loadtxt(pairs_file_path)

        n_iters = ransac_iters(p, k, z)
        print('Interations to be done: ', n_iters)
        H, num_inliers, inliers = ransac(pairs, n_iters, k, ransac_threshold)
        print('Number of inliers:', num_inliers)

        # recompute homography matrix based on inliers
        H = recompute_homography(inliers)
        print(H)
        np.savetxt(homography_file_path, H)

if __name__ == "__main__":
    main()