Discontinuity Preserving Registration using Truncated L1 Regularization and Minimum Spanning Tree based Motion Clustering

Dongxiao Li^{1,2}, Juerong Wu^{1}, Kofi M. Deh^{2}, Thanh D. Nguyen^{2}, Martin R. Prince^{2}, Yi Wang^{2,3}, and Pascal Spincemaille^{2}

The motion field was
modeled using multi-level 3D B-spline meshes. At each level, the registration
problem was formulated as a MRF labelling problem: $${\bf u}=\substack{{\tt\large
arg\,min}\\{{\bf u}_p,\,p\in\mathcal P}}\left(\sum_{p\in\mathcal
P}\|I\left({\bf x}_p\right)-J\left({\bf x}_p+{\bf
u}_p\right)\|_1+\alpha\sum_{\left(p,q\right)\in\mathcal
N}\min\left(\frac{\|{\bf u}_p-{\bf u}_q\|_1}{\|{\bf x}_p-{\bf x}_q\|_2},\tau\right)\right),$$
where $$$p,\,q$$$ were two neighboring control nodes in the mesh, $$${\bf x}_p,\,{\bf x}_q$$$ were their
coordinates, $$${\bf u}_p,\,{\bf u}_q$$$ were their motion vectors (MVs), $$$I\left({\bf
x}_p\right)$$$ was the patch partitioned to node $$$p$$$ in the target image $$$I$$$,
$$$J\left({\bf x}_p+{\bf u}_p\right)$$$ was the corresponding patch in the
source image $$$J$$$, $$$\tau$$$ was a threshold, and $$$\alpha$$$ was a
weighting parameter. In this MRF model, each vertex represented a control node,
and its label was a quantized MV. The source image was upsampled to support
sub-voxel MV search. The first term in this cost function represented the
matching error between the target and the transformed source, and the second
pair-wise term enforced smoothness of the MVs of neighboring control nodes. To
allow for motion discontinuity, we used a truncated L1 regularization. Each vertex
was initially connected with its 6 immediate neighbors. The optimization was
initialized using the L2-norm of differences between neighboring patches in the
target image as edge costs and a single MST was constructed using Prim’s^{5}
algorithm. After that, we used the L2-norm of finite MV differences between
neighboring vertices as edge costs to construct a MST and used the Maximum
Standard Deviation Reduction (MSDR)^{6} algorithm to partition the vertices
into clusters. Next, a new MRF cost function was constructed using only edges completely
within these clusters and solved by fast message passing^{7}. A
multi-resolution approach was used to speed up convergence. Cluster guided
interpolation was used to generate the initial motion field for the next resolution
level and the label space was successively reduced using incremental MV search.
At each level, motion clustering and regularization were alternated for several
iterations.

Data were acquired in one volunteer using a 1.5T
scanner, 8-channel cardiac coil and LAVA-Flex sequence with IRB approval.
Imaging parameters were: voxel size 0.7813x0.7813x2 mm and matrix size 512x512x88.
Five breath-holds M0,…,M4 at different exhalation levels were obtained,
approximately evenly distributed between maximum exhalation at residual volume
(M0) to maximum inhalation at total lung capacity (M4). Registration of M1,…,M4 to
M0 was performed using the proposed method as well as two previously published
nonrigid registration algorithms Elastix^{8} and deedsMST^{2} for comparison.

Fig. 1. Example (a) axial, (b) coronal and (c) sagittal
section of image M0 acquired at the maximum exhalation state. The arrows show
the motion field estimated in registration of M2 with M0. The contours in red
show boundaries of the segmentation generated by the proposed motion clustering
method.

Fig. 2. Comparison of spatial
matching error between M0 and M1~M4 without registration and after using Elastix,
deedsMST and the proposed method: (a) Mean Square Error (MSE) and (b) Normalized
Mutual Information (NMI). In MSE calculation, voxel values were linearly scaled
by mapping the mean of M0 to 0.5.

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)

0781