Derivative of a TPS warp wrt to its source landmarks
In [ ]:
import numpy as np
import os
import scipy.io as sio
from menpo.shape import PointCloud
import menpo.io as pio
from menpo.transform.tps import TPS
import matplotlib.pyplot as plt
We start by defining the source and target landmarks. Notice that, in this first example source = target!!!
In [ ]:
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...
In [ ]:
tps = TPS(src_landmarks, tgt_landmarks)
np.allclose(tps.apply(src_landmarks).points, tgt_landmarks.points)
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.
In [ ]:
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.
In [ ]:
%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])
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.
In [ ]:
print reshaped[1:5,1:5,0,0]
print reshaped[1:5,1:5,0,1]
The sum of all the derivatives wrt the x coordinates should produce an all 1 image
In [ ]:
summed_x = np.sum(reshaped[:,:,:,0], axis=-1)
np.allclose(np.ones(xx.shape), summed_x)
In [ ]:
plt.imshow(summed_x)
and so should the sum of all derivatives wrt the y coordinates.
In [ ]:
summed_y = np.sum(reshaped[:,:,:,1], axis=-1)
np.allclose(np.ones(xx.shape), summed_y)
In [ ]:
plt.imshow(summed_y)
Finally, the derivatives with respect to the x and y coordinates should be in this case exactly the same!!!
In [ ]:
np.allclose(reshaped[:,:,:,0], reshaped[:,:,:,1])
We now show that when source != target things are a bit different.
In [ ]:
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.
In [ ]:
tps = TPS(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.
In [ ]:
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])
We can be easily see that, in this case, the derivatives wrt the x and y coordinates are not equal!!!