Page 166 - Advances in Biomechanics and Tissue Regeneration

162 8. TOWARDS THE REAL-TIME MODELING OF THE HEART

Next, the GMM centroids are reformulated in terms of a velocity function, v, which is used to update the position of
Y as:

TðY,vÞ¼ Y + vðYÞ: (8.39)
To find an estimate of v, a negative log-likelihood function is minimized while imposing independent and identically
distributed data assumptions and adding a regularization term, ϕ(v):

N λ CPD
2 X
f ðv,σ Þ¼ logðpðxÞÞ + ϕðvÞ, (8.40)
2
n¼1
2
with λ CPD being the constant regulating the contribution of the regularization term. In order to find v and σ , the expec-
tation maximization algorithm is employed. The regularization term, ϕ(v), is chosen to take the form of:
Z

, (8.41)
j vðsÞj
ϕðvÞ¼
 D GðsÞ

where s is the so-called frequency domain variable, v, the Fourier transformation of v, and G is taken to be a Gaussian

kernel. A Gaussian kernel is chosen because it is positive definite symmetric and provides a mean to regulate spatial
smoothness. More importantly, it also allows the regularization term to be equivalent to the one presented in the
motion coherence theory, which involves defining a coherent velocity field across a set of points with no prior data
of their motion [89].
The expectation maximization algorithm usually consists of two steps, the E step and the M step. According to
CPD
Myronenko et al. [47], the M step leads in finding W m from the following expression:
1
ðG + λ CPD 2 1 CPD ¼ diagðP1Þ PX Y: (8.42)
σ ðdiagðPÞ1Þ ÞW
W CPD , of size M D, is the list of weight constants. G is known as the Gram matrix with the size of M M and whose
components are given by

2

1
y i y j
(8.43)
2 β
G ij ¼ e
CPD
,
where β CPD is a constant and k k is the norm. With Eq. (8.39), Y is updated through Y update ¼T ðY,W CPD Þ¼ Y + GW CPD
while the matrix of posterior probabilities, P, can be evaluated during the E step of the expectation maximization algo-
rithm through the following equation derived from Eq. (8.40):
1
exp 2σ 2 kx n ðy m + Gðm, ÞW CPD Þk 2

,
P mn ¼
M CPD 2 D=2 (8.44)
X 1 CPD w M
exp 2σ 2 kx n ðy k + Gðm, ÞW Þk + ð2πσ Þ

1 w CPD N
k¼1
N M
P P old 2
where N P ¼ n¼1 m¼1 P ðmjx n Þ. Finally, σ is obtained with the following expression [47]:
1
2 T T T T
σ ¼ ðtrðX diagðP Þ1XÞ 2trððPXÞ TÞ +trðT diagðPÞ1TÞÞ: (8.45)
N P D
To demonstrate CPD registration, consider an idealized 3D human LV created from a half-cut ellipsoid using dimen-
sional data from the literature, as listed in Fig. 8.20.
The LV is then discretized using tetrahedral finite elements to obtain a spatial distribution of points. To represent the
template and the data geometry, two different mesh configurations are utilized. The template mesh consists of 974
nodes and 4060 elements while the data mesh has 771 nodes and 3059 elements. Besides the mesh densities, also
the heart anatomy of the data mesh is altered in terms of translation, stretch, and rotation. The template mesh and
the data mesh are given in Figs. 8.21 and 8.22, respectively.
The CPD algorithm used in this research is based on a MATLAB code written by Myronenko, one of the authors of
Ref. [48]. The code can be publicly accessed from: https://sites.google.com/site/myronenko/research/cpd. To start
the registration procedure, the MATLAB’s code default parameters of w CPD ¼ 0.7, λ CPD ¼ 2, and β CPD ¼ 3 are assigned.
Following the CPD calculation, the results are then analyzed by comparing the deformed data mesh with the original
one. In this case, the undeformed data mesh configuration is considered the exact one as its shape is the same as the
template mesh. A visual representation of the obtained registered mesh as compared with the template one is shown in
Fig. 8.23.

I. BIOMECHANICS