WIP: calculate homography with python (hardcoded sample)

.idea/HomographyCalculator.iml

@@ -2,7 +2,7 @@
 <module type="PYTHON_MODULE" version="4">
   <component name="NewModuleRootManager">
     <content url="file://$MODULE_DIR$" />
-    <orderEntry type="inheritedJdk" />
+    <orderEntry type="jdk" jdkName="Python 3.10 (HomographyCalculator)" jdkType="Python SDK" />
     <orderEntry type="sourceFolder" forTests="false" />
   </component>
   <component name="PyDocumentationSettings">

.idea/misc.xml

@@ -1,4 +1,4 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <project version="4">
-  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.9 (HomographyCalculator)" project-jdk-type="Python SDK" />
+  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.10 (HomographyCalculator)" project-jdk-type="Python SDK" />
 </project>

main.py

@@ -1,15 +1,438 @@
 import numpy as np
-# This is a sample Python script.
 
-# Press Umschalt+F10 to execute it or replace it with your code.
-# Press Double Shift to search everywhere for classes, files, tool windows, actions, and settings.
+def print_hi(name):
+    print(f'Hi, {name}')
 
+    # RANSAC Parameters
+    ransac_threshold = 5.0  # 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
 
-def print_hi(name):
-    # Use a breakpoint in the code line below to debug your script.
-    print(f'Hi, {name}')  # Press Strg+F8 to toggle the breakpoint.
-    matrix = np.loadtxt(c)
+    pairs = np.loadtxt('C:\\Git\\git.tk.informatik.tu-darmstadt.de\\StreetLight\\Assets\\StreamingAssets\\pairs.csv')
+
+    n_iters = ransac_iters(p, k, z)
+    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)
+
+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):
+        """
+        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)
 
 # Press the green button in the gutter to run the script.
 if __name__ == '__main__':