In [ ]:
%matplotlib inline
import menpo.io as mio
from menpo.feature import no_op, igo, hog, sparse_hog, double_igo
import numpy as np

# method to load a database
def load_database(path_to_images, crop_percentage, max_images=None):
    images = []
    # load landmarked images
    for i in mio.import_images(path_to_images, max_images=max_images, verbose=True):
        # crop image
        i.crop_to_landmarks_proportion_inplace(crop_percentage)
        
        # convert it to grayscale if needed
        if i.n_channels == 3:
            i = i.as_greyscale(mode='luminosity')
            
        # append it to the list
        images.append(i)
    return images

In [ ]:
training_images = load_database('/mnt/data/nontas/train200/', 0.5)

In [ ]:
from menpofit.builder import normalization_wrt_reference_shape
normalization_diagonal = 100
group = 'PTS'
label = 'all'
reference_shape, normalized_images = normalization_wrt_reference_shape(training_images, group, label, normalization_diagonal, verbose=True)
all_shapes = [i.landmarks[group][label] for i in normalized_images]

In [ ]:
from menpo.shape import Tree

# Star tree
adjacency_array = np.empty((67, 2), dtype=np.int32)
for i in range(68):
    if i < 34:
        adjacency_array[i, 0] = 34
        adjacency_array[i, 1] = i
    elif i > 34:
        adjacency_array[i-1, 0] = 34
        adjacency_array[i-1, 1] = i

root_vertex = 34

tree = Tree(adjacency_array, root_vertex)

In [ ]:
from menpo.shape import Tree

# MST tree
adjacency_array = np.array([[ 0,  1], [ 1,  2], [ 2,  3], [ 3,  4], [ 4,  5], [ 5,  6], [ 6,  7], [ 7,  8], [ 8,  9], 
                            [ 8, 57], [ 9, 10], [57, 58], [57, 56], [57, 66], [10, 11], [58, 59], [56, 55], [66, 67], 
                            [66, 65], [11, 12], [65, 63], [12, 13], [63, 62], [63, 53], [13, 14], [62, 61], [62, 51],
                            [53, 64], [14, 15], [61, 49], [51, 50], [51, 52], [51, 33], [64, 54], [15, 16], [49, 60],
                            [33, 32], [33, 34], [33, 29], [60, 48], [32, 31], [34, 35], [29, 30], [29, 28], [28, 27],
                            [27, 22], [27, 21], [22, 23], [21, 20], [23, 24], [20, 19], [24, 25], [19, 18], [25, 26],
                            [25, 44], [18, 17], [18, 37], [44, 43], [44, 45], [37, 38], [45, 46], [38, 39], [46, 47],
                            [39, 40], [47, 42], [40, 41], [41, 36]])
root_vertex = 0

tree = Tree(adjacency_array, root_vertex)

In [ ]:
# PDM
from menpofit.builder import build_shape_model
shape_model = build_shape_model(all_shapes, 100)

In [ ]:
from cvpr15.builder import _get_relative_locations

relative_locations = _get_relative_locations(all_shapes, tree, '', True)

In [ ]:
def _build_deformation_model(tree, relative_locations, level_str, verbose):
    def_len = 2 * tree.n_vertices
    def_cov = np.zeros((def_len, def_len))
    all_covs = []
    all_means = []
    for e in range(tree.n_edges):
        # get vertices adjacent to edge
        parent = tree.adjacency_array[e, 0]
        child = tree.adjacency_array[e, 1]

        # compute covariance matrix
        edge_cov = np.cov(relative_locations[..., e])
        all_covs.append(edge_cov)
        edge_cov = np.linalg.inv(edge_cov)
        all_means.append(np.mean(relative_locations[..., e], axis=1))
        
        
        for l in range()

        # store its values
        s1 = edge_cov[0, 0]
        s2 = edge_cov[1, 1]
        s3 = edge_cov[0, 1]
        s3_2 = 2 * s3

        # Fill the covariance matrix matrix
        # get indices
        p1 = 2 * parent
        p2 = 2 * parent + 1
        c1 = 2 * child
        c2 = 2 * child + 1

        # up-left block
        def_cov[p1, p1] += s1
        def_cov[p2, p2] += s2
        def_cov[p2, p1] += 2 * s3

        # up-right block
        def_cov[p1, c1] = - s1
        def_cov[p2, c2] = - s2
        def_cov[p1, c2] = - s3
        def_cov[p2, c1] = - s3

        # down-left block
        def_cov[c1, p1] = - s1
        def_cov[c2, p2] = - s2
        def_cov[c1, p2] = - s3
        def_cov[c2, p1] = - s3

        # down-right block
        def_cov[c1, c1] += s1
        def_cov[c2, c2] += s2
        def_cov[c1, c2] += 2 * s3

    return def_cov, all_covs, all_means

In [ ]:
deformation_model, all_covs, all_means = _build_deformation_model(tree, relative_locations, '', True)
print deformation_model.shape
print len(all_covs)
print len(all_means)

In [ ]:
# CHECK IF: E(li-lj) = E(li) - E(lj)
e = 1
parent = tree.adjacency_array[e, 0]
child = tree.adjacency_array[e, 1]

print all_means[e]
print shape_model.mean().points[child, :] - shape_model.mean().points[parent, :]

In [ ]:
# CHECK IF SUM COST AND MATRIX COST ARE EQUAL
from menpo.shape import PointTree

shape = all_shapes[15]

pointtree = PointTree(shape.points, tree.adjacency_array, tree.root_vertex)
cost1 = 0
for e in range(tree.n_edges):
    parent = tree.adjacency_array[e, 0]
    child = tree.adjacency_array[e, 1]
    
    d = pointtree.relative_location_edge(parent, child) - all_means[e]
    
    cost1 += d.T.dot(np.linalg.inv(all_covs[e])).dot(d)

d = shape.as_vector() - shape_model.mean().as_vector()
cost2 = d.T.dot(deformation_model).dot(d)
print cost1, cost2

In [ ]:
# PLOT GAUSSIANS
from cvpr15.utils import plot_gaussian_ellipse
import matplotlib.pyplot as plt

mean_shape = shape_model.mean().points
for e in range(tree.n_edges):
    # find vertices
    parent = tree.adjacency_array[e, 0]
    child = tree.adjacency_array[e, 1]
    
    # relative location mean
    rel_loc_mean = mean_shape[child, :] - mean_shape[parent, :]
    
    # plot ellipse
    plot_gaussian_ellipse(all_covs[e], mean_shape[parent, :] + rel_loc_mean, n_std=2, facecolor='none', edgecolor='r')

# plot mean shape points
shape_model.mean().view_on(plt.gcf().number)

# create and plto edge connections
PointTree(mean_shape, tree.adjacency_array, tree.root_vertex).view_on(plt.gcf().number)

In [ ]: