Shape registration problem have been an active research topic in computational geometry, computer vision, medical image analysis and pattern recognition communities. Also called the shape alignment, it has extensive uses in recognition, indexing, retrieval, generation and other downstream analysis of a set of shapes. There have been a variety of works that approach this problem, with the methods varying mostly in terms of (can be called pillars of registration) the shape representation, transformation and registration criteria that is used. One such method is proposed by Huang et al. in this paper, which uses a novel combination of the three pillars, where an implicit shape representation is used to register an object both globally and locally. For the registration criteria, the proposed method uses Mutual Information based criteria for its global registration phase, while sum-squared differences (SSD) for its local phase. 

The method starts off with defining an implicit, non-parameteric shape representation which is translation, rotation and scale invariant. This makes the first step of the registration pipeline which transforms the input images into a domain where the shape is implicitly defined. The image is first partitioned into three spaces, namely $[\Omega]$ (the image domain), $[R_S]$ (points inside the shape), $[\Omega - R_S]$ (points outside the shape), and $[S]$ (points lying on the shape boundary). Using this partition, a function based upon the Lipschitz function $\phi : \Omega -> \mathbb{R}^+$ is defined as:

\begin{equation}
    \phi_S(x,y)
    \begin{cases} 
         0 & (x,y) \in S \\
         + D((x,y), S)>0 & (x,y) \in [R_s] \\
         - D((x,y), S)<0 & (x,y) \in [\Omega - R_s]
    \end{cases}
\end{equation}
    
Where $D((x,y),S)$ is the distance function which gives the minimum Euclidean distance between point $(x,y)$ and the shape $S$. 

Given the implicit representation, global shape alignment is performed using the Mutual Information (MI) objective function defined between the probability density functions of the pixels in source image and the target image sampled from the domain $\Omega$. 

\begin{equation}

MI(f_{\Omega}, g_{\Omega}^{A}) = \underbrace{\mathcal{H}[p^{f_{\Omega}}(l_1)]}_{\substack{\text{Entropy of the}\\ \text{distribution representing $f_{\Omega}$}}} + \underbrace{\mathcal{H}[p^{g_{\Omega}^{A}}(l_2)]}_{\substack{\text{Entropy of the}\\ \text{distribution representing $g_{\Omega}^{A}$} \\ \text{which is the} \\ \text{transformed source ISR using $A(\theta)$}}} - \underbrace{\mathcal{H}[p^{f_{\Omega}, g_{\Omega}^{A}}(l_1, l_2)]}_{\substack{\text{Entropy of the}\\ \text{joint distribution}\\\text{representing $f_{\Omega}, g_{\Omega}^{A}$}}}

\end{equation}

Following global registration, local registration is performed by embedding a control point grid using the Incremental Free Form Deformation (IFFD) method. The objective function to minimize is used as the sum squared differences (SSD). The local registration is also offset by using a multi-resolution framework, which performs deformations on control points of varying resolution, in order to account for small local deformations in the shape. In case where there is prior information available for feature point correspondence between the two shapes, this prior knowledge can be added as a plugin term in the overall local registration optimization term. 

The method was applied on statistically modeling anatomical structures, 3D face scan and mesh registration. 

Automated segmentation of various anatomical structures of interest from medical images has been a well grounded field of research in medical imaging. One such problem is related to segmenting whole heart region from a sequence of magnetic resonance imaging (MRI), which is currently done manually, and is time consuming and tedious. Although many automated techniques exist for this, the task remains challenging due to the complex nature of the problem, partly because of low contrast between heart and nearby tissue. Moreover many of the methods are unable to incorporate prior information into the process. To this end, Pluempitiwiriyawej et al. proposed a version of active contour energy minimization based method to segment the whole heart region, including the epicardium, and the left and right ventricular endocardia. 

The proposed method follows the framework laid out by Chan and Vese\cite{Chan2001}. However Pluempitiwiriyawej et al. propose a modified energy function, which consists of four energy terms. The energy function is given below, where $C$ is the contour represented as a level set function $\phi(x,y)$:

$J(C) = \lambda_1 J_1(C) + \lambda_2 J_2(C) + \lambda_3 J_3(C) + \lambda_4 J_4(C)$

The coefficients $\lambda_{1..4}$ determine the weight of terms $J_{1..4}$. The first term $J_1(C)$ is designed to add stochastic models $\mathcal{M}_1, \mathcal{M}_2$ corresponding to the regions inside and outside of the active contour $C$. The models dictate the probability distribution from which the image intensities making up the inside and outside region of the contour are sampled. The negative log of this term is minimized, which essentially maximizes the probability $p(u | C, \mathcal{M}_1, \mathcal{M}_2)$ given the active contour $C$, and the models $\mathcal{M}_1, \mathcal{M}_2$. 

The second term $J_2(C)$ is designed similar to the classical Snakes\cite{Kass1988} in the sense that it uses edges to guide the contour towards the structure of interest. For this term, a simple edge map is used after convolving with a Gaussian filter which smooths out the noise. The term $J_3(C)$ encodes an shape prior which constraints the contour to follow an elliptical shape, and guides it in conjunction with the region and edge information. 

The final term $J_4(C)$ which encodes the total Euclidean arc length of the contour. This forces the contour to be ``smooth", without rough edges. 

The process of minimizing the energy function follows a three-task approach. The first task is to estimate the stochastic model parameters $\mu_k, \sigma^2_k$, and is performed by fixing the position of initial contour $C$, taking derivatives of $J$ w.r.t stochastic model parameters, and solving by equating to zero. The second task estimates the parameters of the ellipse using least squares method. The third and final task involves the contour using the estimated parameters in task one and two, such that it minimizes the function $J$. The method also performs stochastic relaxation, by dynamically changing the values of parameters $\lambda_1, \lambda_2, \lambda_3, \lambda_4$ as the optimization process proceeds. The intuition is that when the optimization starts, the edge and region terms must guide the contour, and as the process proceeds to it's end, the shape prior and contour length term should carry more weight to regularize the effective shape of the contour. 

The study used 48 MRI studies acquired by imaging rat hearts, and compared the proposed method with two earlier methods, namely Xu and Prince's GVF \cite{ChenyangXu1998}, and Chan and Vese \cite{Chan2001}. The authors also design a new quantitative metric, which is a modification of the Chamfer matching \cite{Barrow} technique. The reported results are observed to be in excellent agreement with the gold standard hand-traced contours. However the similarity values for other methods against human gold-standard were not reported. 

The propagation rule used in this paper is the following:

$$
H^l = \sigma \left(\tilde{D}^{-\frac{1}{2}} \tilde{A} \tilde{D}^{-\frac{1}{2}} H^{l-1} W^l \right)
$$

Where $\tilde{A}$ is the [adjacency matrix][adj] of the undirected graph (with self connections, so has a diagonal of 1s and is symmetric) and $H^l$ are the hidden activations at layer $l$. The $D$ matrices are performing row normalisation, $\tilde{D}^{-\frac{1}{2}} \tilde{A} \tilde{D}^{-\frac{1}{2}}$ is [equivalent to][pygcn] (with $\tilde{A}$ as `mx`):

```
rowsum = np.array(mx.sum(1))            # sum along rows
r_inv = np.power(rowsum, -1).flatten()  # 1./rowsum elementwise
r_mat_inv = sp.diags(r_inv)             # cast to sparse diagonal matrix
mx = r_mat_inv.dot(mx)                  # sparse matrix multiply
```

The symmetric way to express this is part of the [symmetric normalized Laplacian][laplace]. This work, and the [other][hammond] [work][defferrard] on graph convolutional networks, is motivated by convolving a parametric filter over $x$. Convolution becomes easy if we can perform a *graph Fourier transform* (don't worry I don't understand this either), but that requires us having access to eigenvectors of the normalized graph Laplacian (which is expensive).

[Hammond's early paper][hammond] suggested getting around this problem by using Chebyshev polynomials for the approximation. This paper takes the approximation even further, using only *first order* Chebyshev polynomial, on the assumption that this will be fine because the modeling capacity can be picked up by the deep neural network.

That's how the propagation rule above is derived, but we don't really need to remember the details. In practice $\tilde{D}^{-\frac{1}{2}} \tilde{A} \tilde{D}^{-\frac{1}{2}} = \hat{A}$ is calculated prior and using a graph convolutional network involves multiplying your activations by this sparse matrix at every hidden layer. If you're thinking in terms of a batch with $N$ examples and $D$ features, this multiplication happens *over the examples*, mixing datapoints together (according to the graph structure). If you want to think of this in an orthodox deep learning way, it's the following:

```
activations = F.linear(H_lm1, W_l)         # (N,D)
activations = activations.permute(1,0)     # (D,N)
activations = F.linear(activations, hat_A) # (D,N)
activations = activations.permute(1,0)     # (N,D)
H_l = F.relu(activations)
```

**Won't this be really slow though, $\hat{A}$ is $(N,N)$!** Yes, if you implemented it that way it would be very slow. But, many deep learning frameworks support sparse matrix operations ([although maybe not the backward pass][sparse]). Using that, a graph convolutional layer can be implemented [quite easily][pygcnlayer].

**Wait a second, these are batches, not minibatches?** Yup, minibatches are future work.

**What are the experimental results, though?** There are experiments showing this works well for semi-supervised experiments on graphs, as advertised. Also, the many approximations to get the propagation rule at the top are justified by experiment.

**This summary is bad.** Fine, smarter people have written their own posts: [the author's][kipf], [Ferenc's][ferenc].

[adj]: https://en.wikipedia.org/wiki/Adjacency_matrix
[pygcn]: https://github.com/tkipf/pygcn/blob/master/pygcn/utils.py#L56-L63
[laplace]: https://en.wikipedia.org/wiki/Laplacian_matrix#Symmetric_normalized_Laplacian
[hammond]: https://arxiv.org/abs/0912.3848
[defferrard]: https://arxiv.org/abs/1606.09375
[sparse]: https://discuss.pytorch.org/t/does-pytorch-support-autograd-on-sparse-matrix/6156/7
[pygcnlayer]: https://github.com/tkipf/pygcn/blob/master/pygcn/layers.py#L35-L68
[kipf]: https://tkipf.github.io/graph-convolutional-networks/
[ferenc]: http://www.inference.vc/how-powerful-are-graph-convolutions-review-of-kipf-welling-2016-2/

The DeepMind guys present an extension to the Neural Turing Machine architecture.

https://i.imgur.com/4TC6yAp.png

They call it a Differentiable Neural Computer (DNC) and it uses 1) an attention mechanism to access information in a matrix that acts as a memory, 2) an attention mechanism to save information to that memory, and 3) a transition matrix that stores information about the order in which rows in the memory are modified, in order to better handle sequential information. They test on the bAbI question answering dataset, a graph inference task, and on solving a puzzle of arranging blocks.

Normally, a Deep Convolutional Network (DCN) has a conditional probabilistic model where we have some parameterised function $f_{\theta}$, and some distribution over the targets to define the loss function (categorical for classification or Gaussian for regression). But, that conditional function is treated as a black box (probabilistically speaking, it's just fit by maximum likelihood). This paper breaks the entire network up into a number of latent factors.

The latent factors are designed to represent familiar parts of DCNs; for example, max-pooling selects from a set of activations and we can model the uncertainty in this selection with a categorical random variable. To recreate every pixel you can imagine selecting a set of paths backwards through the network to reproduce that pixel activation. That's not how their model is parameterised, all I can do is point you to the paper for the real DRMM model definition.

# Inference

The big trick of this paper is that they have designed the probabilistic model so that max-sum-product message passing (also introduced in this paper) inference equates to the forward prop in a DCN. What does that mean? Well, since the network structure is now a hierarchical probabilistic model we can hope to throw better learning algorithms at it, which is what they do.

# Learning

Using this probabilistic formulation, you can define a loss function that includes a reconstruction loss generating the image from responsibilities you can estimate during the forward pass. Then, you can have reconstruction gradients _and_ gradients wrt your target loss, so you can train semi-supervised or unsupervised but still work with what is practically a DCN.

In addition, it's possible to derive a full EM algorithm in this case, and the M-step corresponds to just solving a generalised Least Squares problem. Which means gradient-free and better principled training of neural networks.

# Experimental Results

Not all of the theory presented in this paper is covered in the experiments. They do demonstrate their method working well on MNIST and CIFAR-10 semi-supervised (SOTA with additional variational tweaks). But, there are not yet any experiments using the full EM algorithm they describe (only reporting that results _appear promising_); all experiments use gradient optimisation. They report that their network will train 2-3x faster, but only demonstrate this on MNIST (and what about variation in the chosen optimizer?).

Also, the model can be sampled from, but we don't have any visualisations of the kind of images we would get.