APMA 1170: HW 4
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
APMA 1170: HW 4
Submit your HW via Gradescope before 8:00pm ET on Friday, October 20, 2023. There are 4 exercises and a total of 50 points
Coding Question: (15 points)
The Kabsch algorithm determines the optimal rotation matrix which most closely aligns a set of N points in space with another set of N points. It was presented by Wolfgang Kabsch in a 1976 paper in the journal Acta Crystallographica - the method was devised in order to compare the atomic arrangement of crystalline solids in different coordinate systems. It is also useful in computer graphics.
The problem setup is as follows. We begin with two sets of data stored in N × 3 matrices P and Q:
Each row of P corresponds to a single point in a data set, where (xi , yi , zi ) are the coordinates in 3-dimensional space. Each row of Q is a single point in a separate data set, with coordinates (ui , vi , wi ).
We define the centroids of the two data sets as:
We search for a rotation matrix R ∈ R3×3 which induces the following transformation for a given data point:
(1)
We seek an “optimal” R in that we seek to minimize the deviation between the points in the data set Q and the points in the data set P after being transformed according to eq. (1).
The algorithm: The algorithm proceeds as follows. (All components are real-valued, so we use transpose instead of adjoint).
1. Calculate the centroids +¯(x) ¯(y) ¯(z),T and+¯(u) ¯(v) ¯(w),T of the two data sets.
2. Define a new matrix P(-) by subtracting the coordinates of the centroid from each data point:
P(-)(i, :) = +xi − ¯(x) yi − ¯(y) zi −¯(z),
3. Similarly, subtract the centroid coordinates (¯(u), ¯(v), ¯(w)) from the data in each row of Q to form a matrix Q(-) .
4. Compute the cross-covariance matrix H = P Q
5. Compute the (full) singular value decomposition H = UΣVT .
6. Calculate a correction factor d = sign.det(VUT )/ which ensures a right-handed coordinate system. (The sign of the determinant - Matlab has built-in functions for computing the determinant and the sign of a number).
(2)
8. The optimal rotation matrix is then given by R = . /−1 .
Implement this algorithm by completing the file “kabsch rotation
.m” . It should take in two N×3 matrices P and Q and return the optimal rotation matrix R ∈ R3×3 . Your code should NOT compute any inverses directly. (No calling Matlab’s “inv” , no backslash operator). The matrix R can be computed using U and V and the identity in eq. (2). You don’t need to actually compute H - you can form R directly with the information given.
Example of optimal rotation for points lying in a two-dimensional plane (i.e., with coordinates zi = wi = 0.
2023-10-21