In [ ]:
import scipy.io
import numpy as np
from menpo.transform import AffineTransform
from menpo.lucaskanade.image import ImageInverseCompositional
from menpo.image import MaskedImage, BooleanImage
from menpo.lucaskanade.residual import *
import matplotlib.pyplot as plt
from menpo.convolution import log_gabor
from collections import OrderedDict
In [ ]:
def compute_rms_point_error(test_pts, template_affine, M):
iteration_pts = M.apply(template_affine.T)
diff_pts = test_pts - iteration_pts
return np.sqrt(np.mean(diff_pts ** 2))
In [ ]:
def my_test_affine(tdata, pt_offsets, alg_list, n_iters, n_freq_tests, spatial_sigma, warp, verbose):
# tdata has three fields:
# tdata.img1 greyscale target image
# tdata.img2 greyscale template image
# tdata.tmplt [x_start y_start x_end y_end] template rectangle corners
#
# pt_offsets is a N x 6 matrix of (x1, y1, x2, y2, x3, y3) deltas for each
# of the three affine control points. Zero mean, unit variance,
# pt_offsets = np.random.rand(N, 6)
#
# alg_list is an ordered dictionary of similarity measures to use:
# {'String measure name', (MeasureClass, constructor params, constructor keyword arguments)}
# eg.
# {'LeastSquares', (LeastSquares, None, None)}
# Target affine warp control points - triangle on the rectangle
target_affine = np.array([[tdata.tmplt[0], tdata.tmplt[1]],
[tdata.tmplt[2], tdata.tmplt[1]],
[tdata.tmplt[0] + ((tdata.tmplt[2] - tdata.tmplt[0]) / 2) - 0.5, tdata.tmplt[3]]])
# Template image dimensions (formerly had +1)
template_width = tdata.tmplt[2] - tdata.tmplt[0]
template_height = tdata.tmplt[3] - tdata.tmplt[1]
tmplt_shape = (template_height, template_width)
# Template affine warp control points
template_affine = np.array([[0, 0],
[template_width, 0],
[template_width / 2, template_height]])
# Initial warp parameters. Unperturbed translation
p_init = np.zeros([3, 3])
p_init[0, 2] = tdata.tmplt[1] - 1
p_init[1, 2] = tdata.tmplt[0] - 1
# Translate by 0.5 pixels to avoid identity warp. Warping introduces a little
# smoothing and this avoids the case where the first iteration for a forwards
# algorithm is on the "unsmoothed" unwarped image
p_init[0, 2] += 0.5
p_init[1, 2] += 0.5
p_init = p_init[:2, :3].flatten(order='F')
# Pick a total of n_freq_tests point offsets from pt_offsets randomly
ind = np.floor(pt_offsets.shape[0] * np.random.rand(n_freq_tests, 1)).astype(np.uint32)
ind[ind == 0] = 1
# Uncomment this line to randomise the noise added
# pt_offsets1 = pt_offsets[ind, :]
pt_offsets1 = pt_offsets
# Scale point offsets to have required sigma
pt_offsets1 = pt_offsets1 * spatial_sigma
img = MaskedImage(tdata.img1)
results = {}
# Rotate points for new ordering
template_affine = template_affine[:, ::-1]
M1 = np.vstack([template_affine.T, np.ones([1, template_affine.shape[0]])]).T
base_transform = AffineTransform(np.eye(3))
# Run
for offset_idx in xrange(n_freq_tests):
# Test points: apply current point offset to target points
test_pts = target_affine + np.reshape(pt_offsets1[offset_idx, :], [3, 2])
# Rotate points for new ordering
test_pts = test_pts[:, ::-1]
# Solve for affine warp
M2 = np.vstack([test_pts.T, np.ones([1, test_pts.shape[0]])]).T
M = np.linalg.solve(M1, M2).T
# Warp original image to get test "template" image
tmplt = MaskedImage(tdata.img2).warp_to(BooleanImage.blank(tmplt_shape), AffineTransform(M))
if verbose:
print 'Initial RMS: {0}'.format(compute_rms_point_error(test_pts, template_affine.T, AffineTransform.from_vector(p_init)))
# Run each algorithm
for sim_measure in alg_list:
# Allow the passing of arguments in to the instantiated class
clz = alg_list[sim_measure][0]
args = alg_list[sim_measure][1]
kwargs = alg_list[sim_measure][2]
if not args is None and not kwargs is None:
metric = clz(*args, **kwargs)
elif not args is None:
metric = clz(*args)
elif not kwargs is None:
metric = clz(**kwargs)
else:
metric = clz()
final_transform = ImageInverseCompositional(tmplt, metric, base_transform, interpolator='scipy').align(img, p_init, max_iters=n_iters)
rms_pt_error = compute_rms_point_error(test_pts, template_affine.T, final_transform)
if verbose:
print '{0}: {1}'.format(sim_measure, rms_pt_error)
if not sim_measure in results:
results[sim_measure] = []
measure_results = results[sim_measure]
measure_results.append(rms_pt_error)
results[sim_measure] = measure_results
return results
In [ ]:
%matplotlib inline
# Load datasets
np.set_printoptions(suppress=False, linewidth=600)
pt_offset = scipy.io.loadmat('/home/pts08/Dropbox/lk/face_recognition_2011/data/affine_pt_offset.mat')['pt_offset']
yale_cropped = scipy.io.loadmat('/home/pts08/Dropbox/lk/face_recognition_2011/data/myYaleCropped.mat')
example_imgs = yale_cropped['example_imgs']
coords = yale_cropped['coords'][0]
tmplts = yale_cropped['tmplts']
# Pre-compute
tmplt_height = coords[3] - coords[1]
tmplt_width = coords[2] - coords[0]
num_of_subjs = tmplts.shape[2]
num_of_imgs_per_subj = example_imgs.shape[2] / num_of_subjs
# Set up experiment variables
verbose = True
n_iters = 30 # Number of gradient descent iterations
n_freq_tests = 30 # Number of frequency of convergence tests
max_spatial_error = 3.0 # Max location error for deciding convergence
all_spc_sig = np.arange(1, 11) # All spatial sigmas (1,10)
warp = 'scipy'
# Gabor variables
blank_image = np.ones([tmplt_height, tmplt_width])
alg_list = OrderedDict([('LSIntensity', (LSIntensity, None, None)),
('ECC', (ECC, None, None)),
#'IRLS': (IRLS, None, None),
('GaborFourier', (GaborFourier, None, None)),
('GradientImages', (GradientImages, None, None)),
('GradientCorrelation', (GradientCorrelation, None, None))])
results = np.zeros([num_of_subjs, num_of_imgs_per_subj, len(all_spc_sig), len(alg_list)])
In [ ]:
# Run experiment for each subject
for subj in xrange(num_of_subjs):
print 'Subject {0} of {1}'.format(subj, num_of_subjs)
example_imgs_per_subj = example_imgs[:, :, num_of_imgs_per_subj * subj:num_of_imgs_per_subj * (subj + 1)]
tmplt = np.ascontiguousarray(tmplts[:, :, subj])
for subj_img in xrange(num_of_imgs_per_subj):
tdata = lambda x: 0
tdata.img1 = np.ascontiguousarray(example_imgs_per_subj[:, :, subj_img])
tdata.img2 = tmplt
tdata.tmplt = coords
# Matrix S for Gabor-Fourier method
# Thanks to Peter Kovesi's Gabor Filters
# http://www.csse.uwa.edu.au/~pk/
S = log_gabor(blank_image, num_orientations=32, num_scales=32)[2]
alg_list['GaborFourier'] = (GaborFourier, [blank_image.shape], {'filter_bank': S})
# Run tests
for sigma_ind, current_sigma in enumerate(all_spc_sig):
res = my_test_affine(tdata, pt_offset, alg_list, n_iters, n_freq_tests, current_sigma, warp, verbose)
for measure_ind, sim_measure in enumerate(alg_list):
measure_results = res[sim_measure]
# Get whether or not it converges
n_converge = len(filter(lambda error: error < max_spatial_error, measure_results))
results[subj, subj_img, sigma_ind, measure_ind] = n_converge
# Save out results just in case
scipy.io.savemat('results.mat', {'results': results})
In [ ]:
# Plot results
mean_results = np.mean(np.mean(results, 1), 0) / float(n_freq_tests)
line_styles = ['k--D', 'y:^', 'r:*', 'g:^', 'b-s']
lines = []
for i in xrange(mean_results.shape[1]):
lines.append(all_spc_sig)
lines.append(mean_results[:, i])
lines.append(line_styles[i])
p = plt.plot(*lines)
legend_labels = [a for a in alg_list.keys()]
plt.yticks(np.linspace(0, 1, 11))
plt.xticks(all_spc_sig)
plt.ylabel('Frequency of Convergence')
plt.xlabel('Point Standard Deviation')
plt.legend(p, legend_labels)
plt.title('Yale: with Smoothing')