As usual, we first initialise the Python environment and define some helper functions for loading and plotting point sets.
NOTE: If you comment out line 7 by putthin a # in front, you can get 3D plots in an extra window in which you can rotate the view. You may need to restart the IPython kernel in the 'Kernel' menu above.
%matplotlib inline
import numpy as np
from vtk import *
from vtk.util import numpy_support
from matplotlib import pyplot
import pylab
from mpl_toolkits.mplot3d import Axes3D
# Helper method to read point sets from a file
def vtk_point_read(filename):
# load a vtk file as input
reader = vtkPolyDataReader()
reader.SetFileName(filename)
reader.Update()
data = reader.GetOutput()
x = np.zeros(data.GetNumberOfPoints())
y = np.zeros(data.GetNumberOfPoints())
z = np.zeros(data.GetNumberOfPoints())
for i in range(data.GetNumberOfPoints()):
x[i], y[i], z[i] = data.GetPoint(i)
x = np.reshape(x, (1, len(x)))
y = np.reshape(y, (1, len(y)))
z = np.reshape(z, (1, len(z)))
points = np.concatenate((x, y, z), axis=0)
return points
# Helper function for plotting 3D point sets
def myplotin3D(pnt1_i,pnt2_i,title=''):
fig = pylab.figure()
ax = Axes3D(fig)
ax.scatter(pnt1_i[0, 0::5], pnt1_i[1, 0::5], pnt1_i[2, 0::5], c='b')
ax.scatter(pnt2_i[0, 0::5], pnt2_i[1, 0::5], pnt2_i[2, 0::5], c='r')
ax.set_title(title,fontsize=16)
ax.view_init(50,0)
pyplot.show()
Below is a skeleton implementation of the Umeyama method for computing rigid-body transformation between corresponding point sets X and Y. For demonstration purposes, we use a single point set corresponding to a liver surface extracted from a real 3D CT image. The moving point set is generated by transforming the fixed point set and applying some Gaussian noise to it. Thus, the two point sets are not identical.
TASK: Finish the implementation by inserting the missing lines of code. You might want to check out the lecture notes for this.
# Implementation of the Umeyama method for computing rigid-body transformations
# between corresponding point sets X and Y
def umeyama_rigid(X, Y):
# Get dimension and number of points
m, n = X.shape
# Demean the point sets X and Y
X_mean = # INSERT CODE HERE
Y_mean = # INSERT CODE HERE
X_demean = # INSERT CODE HERE
Y_demean = # INSERT CODE HERE
# Computing matrix XY' using demeaned point sets
XY = # INSERT CODE HERE
# Singular value decomposition
U,D,V = np.linalg.svd(XY,full_matrices=True,compute_uv=True)
V=V.T.copy()
# Determine rotation
R = # INSERT CODE HERE
# Determine translation
t = # INSERT CODE HERE
return R,t
# Load point set
filename1 = "./data/liver1_dec.vtk"
pnt1 = vtk_point_read(filename1)
shp1 = pnt1.shape
# Copy the points and apply rigid transform
pnt2=pnt1.copy()
rotation_in = [[-0.3214, -0.5567, 0.7660],[-0.0637, -0.7944, -0.6040],[0.9448, -0.2429, 0.2198]]
trans_in = [10,20,-30]
pnt2 = np.dot(rotation_in,pnt1) + np.tile(trans_in,(shp1[1], 1)).transpose()
# Add Gaussian noise
sigma = max(pnt2.flatten()) * 0.02
noise = np.random.normal(0.0,sigma,pnt2.shape)
pnt2 = pnt2 + noise
# Initial position
myplotin3D(pnt1,pnt2,'initial alignment (moving=blue,target=red)')
# Running Umeyama method
R, t = umeyama_rigid(pnt1,pnt2)
print "Rotation matrix=\n",R,"\nTranslation vector=",t
# After registration
pnt1 = np.dot(R,pnt1) + np.tile(t, (shp1[1], 1)).transpose()
myplotin3D(pnt1,pnt2,'after registration (moving=blue,target=red)')
Rotation matrix= [[-0.32082409 -0.55738439 0.76576403] [-0.06327984 -0.79408172 -0.60450796] [ 0.94502252 -0.24239815 0.21948935]] Translation vector= [ 9.97581107 19.91363448 -30.05124307]
In the lecture, we have only discussed the basic Umeyama method as above. Now, we will implement a more advaned variant that can determine a similarity transformation (i.e. rotations + translation + scaling).
TASK: a) Finish the implementation below by inserting your code from above, b) try to understand what the additional lines of code do. References are given to the equations in the original work by Umeyama which can be found here: http://web.stanford.edu/class/cs273/refs/umeyama.pdf
# Implementation of the Umeyama method for computing similarity transformations
# between corresponding point sets X and Y
def umeyama_similarity(X, Y):
# Get dimension and number of points
m, n = X.shape
# Demean the point sets X and Y
X_mean = # INSERT CODE HERE
Y_mean = # INSERT CODE HERE
X_demean = # INSERT CODE HERE
Y_demean = # INSERT CODE HERE
# Computing matrix XY' using demeaned and NORMALISED point sets (divide by the number of points n)
# See Equation (38) in the paper
XY = # INSERT CODE HERE
# Determine variances of points X and Y, see Equation (36),(37) in the paper
X_var = np.mean(np.sum(X_demean*X_demean, 0))
Y_var = np.mean(np.sum(Y_demean*Y_demean, 0))
# Singular value decomposition
U,D,V = np.linalg.svd(XY,full_matrices=True,compute_uv=True)
V=V.T.copy()
# Determine rotation
R = # INSERT CODE HERE
# Determine the scaling, see Equation (42) in the paper (assume S to be the identity matrix, so ignore)
c = # INSERT CODE HERE
# Determine translation, see Equation (41) in the paper
t = # INSERT CODE HERE
return R,c,t
# Load point set
filename1 = "./data/liver1_dec.vtk"
pnt1 = vtk_point_read(filename1)
shp1 = pnt1.shape
# Copy the points and transform them.
pnt2=pnt1.copy()
rotation_in = [ [-0.3214, -0.5567, 0.7660],
[-0.0637, -0.7944, -0.6040],
[0.9448, -0.2429, 0.2198] ]
scale_in = 2
trans_in = [10,20,-30]
pnt2 = scale_in * np.dot(rotation_in,pnt1) + np.tile(trans_in,(shp1[1], 1)).transpose()
# Add Gaussian noise
sigma = max(pnt2.flatten()) * 0.02
noise = np.random.normal(0.0,sigma,pnt2.shape)
pnt2 = pnt2 + noise
# Initial position
myplotin3D(pnt1,pnt2,'initial alignment (moving=blue,target=red)')
# Running Umeyama method
R, c, t = umeyama_similarity(pnt1,pnt2)
print "Rotation matrix=\n",R,"\nScaling coefficient=",c,"\nTranslation vector=",t
# After registration
pnt1 = c * np.dot(R,pnt1) + np.tile(t, (shp1[1], 1)).transpose()
myplotin3D(pnt1,pnt2,'after registration (moving=blue,target=red)')
Rotation matrix= [[-0.32123678 -0.55677537 0.76603402] [-0.06338237 -0.79445234 -0.60401006] [ 0.94487544 -0.2425833 0.21991757]] Scaling coefficient= 1.99860054983 Translation vector= [ 9.88952992 19.9544177 -30.04901934]