Derivative of a TPS warp wrt to its source landmarks
import os
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt
from menpo.shape import PointCloud
import menpo.io as mio
from menpo.transform import ThinPlateSplines
We start by defining the source and target landmarks. Notice that, in this first example source = target!!!
src_landmarks = PointCloud(np.array([[-1, -1],
[-1, 1],
[ 1, -1],
[ 1, 1]]))
tgt_landmarks = PointCloud(np.array([[-1, -1],
[-1, 1],
[ 1, -1],
[ 1, 1]]))
The warp can be effectively computed, although the rendering will not appear to be correct...
tps = ThinPlateSplines(src_landmarks, tgt_landmarks)
np.allclose(tps.apply(src_landmarks).points, tgt_landmarks.points)
True
The next step is to define the set of points at which the derivative of the previous TPS warp must be evaluated. In this case, we use the function meshgrid to generate points inside the convex hull defined by the source landmarks.
x = np.arange(-1, 1, 0.01)
y = np.arange(-1, 1, 0.01)
xx, yy = np.meshgrid(x, y)
points = np.array([xx.flatten(1), yy.flatten(1)]).T
We evaluate the derivative, reshape the output, and visualize the result.
%matplotlib inline
dW_dxy = tps.jacobian_source(points)
reshaped = dW_dxy.reshape(xx.shape + (4,2))
#dW_dx
plt.subplot(240)
plt.imshow(reshaped[:,:,0,0])
plt.subplot(241)
plt.imshow(reshaped[:,:,1,0])
plt.subplot(242)
plt.imshow(reshaped[:,:,2,0])
plt.subplot(243)
plt.imshow(reshaped[:,:,3,0])
#dW_dy
plt.subplot(244)
plt.imshow(reshaped[:,:,0,1])
plt.subplot(245)
plt.imshow(reshaped[:,:,1,1])
plt.subplot(246)
plt.imshow(reshaped[:,:,2,1])
plt.subplot(247)
plt.imshow(reshaped[:,:,3,1])
<matplotlib.image.AxesImage at 0x171f5748>
If everything goes as expected, the upper corner of the images defining the derivative of the warp wrt the x and y coordinates of the first of the source landmarks should both contain values close to 1.
print(reshaped[1:5,1:5,0,0])
print(reshaped[1:5,1:5,0,1])
[[ 0.99240803 0.98853625 0.98459474 0.98059178] [ 0.98853625 0.98468246 0.98075813 0.9767711 ] [ 0.98459474 0.98075813 0.97685217 0.97288316] [ 0.98059178 0.9767711 0.97288316 0.96893276]] [[ 0.99240803 0.98853625 0.98459474 0.98059178] [ 0.98853625 0.98468246 0.98075813 0.9767711 ] [ 0.98459474 0.98075813 0.97685217 0.97288316] [ 0.98059178 0.9767711 0.97288316 0.96893276]]
The sum of all the derivatives wrt the x coordinates should produce an all 1 image
summed_x = np.sum(reshaped[:,:,:,0], axis=-1)
np.allclose(np.ones(xx.shape), summed_x)
True
plt.imshow(summed_x)
<matplotlib.image.AxesImage at 0x179a2320>
and so should the sum of all derivatives wrt the y coordinates.
summed_y = np.sum(reshaped[:,:,:,1], axis=-1)
np.allclose(np.ones(xx.shape), summed_y)
True
plt.imshow(summed_y)
<matplotlib.image.AxesImage at 0x17d80d30>
Finally, the derivatives with respect to the x and y coordinates should be in this case exactly the same!!!
np.allclose(reshaped[:,:,:,0], reshaped[:,:,:,1])
True
We now show that when source != target things are a bit different.
tgt_landmarks = PointCloud(np.array([[-0.6, -1.3],
[-0.8, 1.2],
[ 0.7, -0.8],
[ 1.3, 0.5]]))
The warp can be computed and visualized.
tps = ThinPlateSplines(src_landmarks, tgt_landmarks)
np.allclose(tps.apply(src_landmarks).points, tgt_landmarks.points)
tps.view()
and so can its derivative, which we evaluate using the the same set of points used in the first example.
dW_dxy = tps.jacobian_source(points)
reshaped = dW_dxy.reshape(xx.shape + (4,2))
#dW_dx
plt.subplot(240)
plt.imshow(reshaped[:,:,0,0])
plt.subplot(241)
plt.imshow(reshaped[:,:,1,0])
plt.subplot(242)
plt.imshow(reshaped[:,:,2,0])
plt.subplot(243)
plt.imshow(reshaped[:,:,3,0])
#dW_dy
plt.subplot(244)
plt.imshow(reshaped[:,:,0,1])
plt.subplot(245)
plt.imshow(reshaped[:,:,1,1])
plt.subplot(246)
plt.imshow(reshaped[:,:,2,1])
plt.subplot(247)
plt.imshow(reshaped[:,:,3,1])
<matplotlib.image.AxesImage at 0x189d2b38>
We can be easily see that, in this case, the derivatives wrt the x and y coordinates are not equal!!!