Registration

In this assignment, you will implement the paired point matching (PPM) and iterative closest point (ICP) algorithms for patient-to-image registration as discussed in the “Registration” lecture.

Theory

_images/icp_animation.gif

From Wikipedia, the four steps in ICP are:

  1. For each point in the source point cloud, match the closest point in the reference point cloud (or a selected set).

  2. Estimate the combination of rotation and translation using a root mean square point to point distance metric minimization technique which will best align each source point to its match found in the previous step. This step may also involve weighting points and rejecting outliers prior to alignment.

  3. Transform the source points using the obtained transformation.

  4. Iterate (re-associate the points, and so on).

Paired Point Matching

As discussed in the lecture, ICP is essentially an iterative application of the paired point matching algorithm. Your first task is to implement this algorithm in the file assignments/registration/registration.py. You are given point clouds of corresponding points and must calculate the transformation matrix that maps the source to the target point clouds.

 1def paired_point_matching(source, target):
 2    """
 3    Calculates the transformation T that maps the source to the target point clouds.
 4    :param source: A N x 3 matrix with N 3D points.
 5    :param target: A N x 3 matrix with N 3D points.
 6    :return:
 7        T: 4x4 transformation matrix mapping source to target.
 8        R: 3x3 rotation matrix part of T.
 9        t: 1x3 translation vector part of T.
10    """
11    assert source.shape == target.shape
12    T = np.eye(4)
13    R = np.eye(3)
14    t = np.zeros((1, 3))
15
16    ## TODO: your code goes here
17
18    return T, R, t

Test your implementation by running the script in PyCharm or from the console using python cas/registration/registration.py.

Iterative Closest Point

With the basic building block of ICP in place, implement the iterative part of the algorithm in the same Python file. This script will open a Matplotlib figure to help visualise the results. You can interactively navigate the plot using the commands provided by Matplotlib.

Initial pose

Provide an initial pose by estimating it from the data or selecting it manually.

 1def get_initial_pose(source, target):
 2    """
 3    Calculates an initial rough registration or optionally returns a hand-picked initial pose.
 4    :param source: A N x 3 point cloud.
 5    :param target: A N x 3 point cloud.
 6    :return: An initial 4 x 4 rigid transformation matrix.
 7    """
 8    T = np.eye(4)
 9
10    ## TODO: Your code goes here
11
12    return T

Find the Closest Point

Implement a function to find the closest point in the target data for each point in the source data. It is recommended to use a KD-tree for efficiency (see Resources).

Your function should return a tuple of two arrays: the first containing the distances to the nearest neighbor in the target for each point in the source, and the second containing the indices of these nearest neighbors in the target.

 1def find_nearest_neighbor(source, target):
 2    """
 3    Finds the nearest neighbor in 'target' for every point in 'source'.
 4    :param source: A N x 3 point cloud.
 5    :param target: A N x 3 point cloud.
 6    :return: A tuple containing two arrays: the first array contains the
 7             distances to the nearest neighbor in 'target' for each point
 8             in 'source', and the second array contains the indices of
 9             these nearest neighbors in 'target'.
10    """
11
12    ## TODO: replace this by your code
13    pass

Iterative Matching

Implement the iterative matching process that applies the paired point matching until your error function converges.

 1def icp(source, target, init_pose=None, max_iterations=10, tolerance=0.0001):
 2    """
 3    Iteratively finds the best transformation mapping the source points onto the target.
 4    :param source: A N x 3 point cloud.
 5    :param target: A N x 3 point cloud.
 6    :param init_pose: Initial pose as a 4 x 4 transformation matrix.
 7    :param max_iterations: Maximum iterations.
 8    :param tolerance: Error tolerance.
 9    :return: The optimal 4 x 4 rigid transformation matrix, distances, and registration error.
10    """
11    # Initialisation
12    T = np.eye(4)
13    distances = 0
14    error = np.finfo(float).max
15
16    ## TODO: Your code goes here
17
18    return T, distances, error

Code Submission

Submit a ZIP file named lastname_firstname_assignment3.zip on ILIAS containing:

  1. The modified registration.py as lastname_firstname_assignment3_code.py.

  2. Console output in a text file named lastname_firstname_assignment3_output.txt.

  3. A screenshot of the plots as lastname_firstname_assignment3_screenshot.png.

Online Questions

Complete the “Assignment 3 - Questions” on ILIAS:

  • Answer all questions.

  • Each question has only one correct answer.

  • All questions are equally weighted. Incorrect answers will not result in point deductions.

  • You are allowed only one attempt to complete the test.

Assignment Evaluation

This assignment constitutes 25% of your total assignment grade, split equally between:

  • Code Evaluation (50%): points are awarded as follows:

    • 4 points for a working solution.

    • 3 points for only small errors.

    • 2 points for a substantial effort.

    • 1 point for substantial errors or minimal effort.

    • 0 points for no attempt or plagiarism.

  • Online questions (50%)

Resources

KD-Trees

You don’t neccessarily need to use them, but it will make the search for the closest point faster.

Solving systems of equations