Cardiocyte Shape Representation With Integral Kernels
We wish to retrieve the information embedded in the shapes of a contracting cardiocyte like the ones shown in Figure . In other words, we want to analyze closed planar regions D
2, and their boundaries (finite perimeters), as depicted in Figure . We will describe these regions by binary images, Figure , composed with a suitable class of (continuous and invertible) image domain transformation. The binary images will be represented by a characteristic function
defined for
x ![[set membership]](/corehtml/pmc/pmcents/x2208.gif)
Ω
2, with
D ![[subset or is implied by]](/corehtml/pmc/pmcents/x2282.gif)
Ω, where Ω is the rectangular image domain.
We will define the multi-scale nonlinear features
Rσ, as the convolution of the shape
S with a family of kernels
Kσ, indexed by a scale
σ. More specifically,
σ
+,
K :
2 ×
+ →
+; (
x,
σ)
Kσ (
x). For convenience, we will consider the isotropic Gaussian kernel of the form
We will work with the non-linear features proposed by Hong, Soatto and Vese [
20], which were designed to retain boundary information. They are given by one of the following two expressions,
or the symmetrized version
The shape representation for different scales using the first (simpler) features, Eq. (3), are shown in Figure . The shape representation for different scales using the second (symmetric) features, Eq. (4), are shown in Figure . Both the binary representation (S) and the nonlinear shape features (Rσ) include the original boundary information. However, the nonlinear shape features also encode the local shape information (up to a scale σ), which is not explicitly available in the binary representation.
These shape features have several useful properties (for more details on these properties, please see [
20]): (i) they are very robust under the presence of noise that gets incorporated in the segmentation process; (ii) since their values depend on the local geometry, these features propagate the shape information inside and outside the boundary; (iii) because the value of the features at a point is a local statistic of the shape in a neighborhood of that point, these features capture the context of the particular shape; and (iv) these shape features are very straightforward to compute.
Cardiocyte Contractility Assessment via Shape Matching
In this paper, the contractility analysis is done at the cellular level, thereby only measuring the overall contractions in the cell. This is consistent with the contraction measurements that are being used in our laboratory. We are working on a similar energy conservation and information content approach for assessing contractility in neonatal cardiocytes, where we will measure the fine granular changes in the image. Unlike adult cardiocytes, which are highly organized and quite similar in morphology, the neonatal cardiocyte is in the process of developing its contractile machinery. The neonatal cardiocyte is generally unable to retract its cell boundary during contraction, and noticeable changes occur only within the cell perimeter. For these reasons, it is difficult to perform contractile measurements on this cell type in a manner similar to that of the adult cardiocyte, in which changes in cell boundary are quantified during contraction. (Please see a recent article by Bazan, Torres-Barba, Paolini and Blomgren [
27] for a previously developed computational framework for the quantitative assessment of contractile responses of isolated neonatal cardiac myocytes.)
We are given two shapes of the same topology,
i.e., the shapes of a relaxed and contracted cardiocyte, respectively (Figure ), defined by the two functions
S1,
S2 : Ω → {0, 1}. Our intent is to transform one into the other, and
vice-versa, in a process that resembles that of the contraction/relaxation of the cardiocyte. As proposed in [
20], we will do this by
warping, that is a domain deformation
h : Ω →
2 such that
h (Ω) = Ω and
We are interested in the
distance between the two shapes,
i.e., our proposed measure of contraction. The distance between the shapes will be defined as the energy of the aforementioned warping. Since there are infinitely many warpings that satisfy (5), in order to make this distance unique, it is defined as the one that minimizes the energy in a suitably chosen class [
28]. For instance, as proposed in [
20],

subject to (5), where
h is a diffeomorphism and ||·|| is some chosen norm of integral form. This minimization process can be recast as a variational problem of the form
where H is a suitable function space. The distance between the two shapes, as defined above, is a function of two energy components: Edata (S1, S2 |h), which represents the data fidelity; and Ereg (h), which is a regularization term. Both terms are explained in detail below.
In order to guarantee a meaningful
vis-à-vis correspondence between the two shapes up to a certain scale [
22] (
i.e., a point
x in the interior of the set that defines
S1 is mapped to a point
h (
x) in the interior of the set that defines
S2; and similarly, a point x in the exterior of the set that defines
S1 is mapped to a point
h (
x) in the exterior of the set that defines
S2), we adopt the measure of data fidelity between the two shapes
S1 and
S2 proposed in [
20]:
where Rσ is the shape features (3) or (4).
In our application, the purpose of the regularization term
Ereg (
h) is two-fold. First, it is there to render the problem well posed and is designed to penalize variations of the diffeomorphism function
h in favor of smoothness. Second, it makes the deformations compatible with the deforming material. Several authors have provided insight into the constitutive laws that describe the mechanical responses of resting and contracting cardiac muscle, along with their regional and temporal variations [
23,
29-
31]. The problem is, however, extremely complex and has required of elaborate combinations of multiaxial tissue testing [
32,
33], microstructural morphological modeling [
34,
35], statistical parameter estimation, and validation with measurements [
36,
37]. Very sophisticated numerical methods are also essential for accurate quantitative analysis in all phases of the investigations [
23].
The intact cardiac muscle undergoes finite deformations during the normal cardiac cycle. Thus, the classical linear theory of elasticity is inappropriate for resting myocardial mechanics [
38,
39,
25]. The myocardium is frequently modeled as a finite hyper-elastic material, where the second Piola-Kirchhoff stress tensor components
Pij, are related to the components of the Lagrangian Green's strain tensor
Eij, through the pseudo-strain energy
W, as
where
Cij is the right Cauchy-Green deformation tensor and p is a hydrostatic pressure Lagrange multiplier [
23] (which we assume to be zero in this analysis). Several functional forms have been proposed for
W [
24,
40,
30,
42,
26,
33]. We will adopt the transversely isotropic functional form proposed by Guccione, McCulloch and Waldman [
24] that considers the fibrous structure of the myocardium. The strain energy potential W is an exponential function of the strain components
Eij referred to the fiber coordinates
where
The aforementioned fiber coordinate system has the coordinate directions of the muscle fiber axis, the axis of the myofiber sheets, and the axis normal to the sheets, and is derived by rotating the cardiac coordinate system through the two angles that define the local myofiber-sheet orientation [
24]. In Eq. (9),
E11 is the fiber strain,
E22 is the cross-fiber in-plane strain, and E
12 is the shear strain in the fiber-cross fiber coordinate plane. Omens, MacKenna and McCulloch [
26] have found that the material constants
C = 1.1 kPa,
bf = 9.2 kPa,
bt = 2.0 kPa, and
bfs = 3.7 kPa are appropriate to model the strains measured in the rat midwall. Then, the regularization term,
Ereg (
h), can be written as
The optimal correspondence given by h* is obtained by
The energy minimization is performed in a variational framework using a gradient descent method. The Euler-Lagrange equation corresponding to the energy E = Edata + Ereg yields the gradient direction for h:
where α is a small parameter (Lagrange multiplier). Using the features from Eq. (3),
Rσ (x | S) = S (x) (Kσ * (1 - S (x))), we have
and, in components form, we have
where
With the following short-hand notation
Average Shape of the Relaxed State
In order for us to use the energy of the warping as a measure of the cardiocyte's contractions, we need to determine a baseline that will represent the state when the cell is not contracting (relaxed). We will define this baseline as the average shape of the cardiocyte's relaxed state. In other words, given an ensemble of shapes {S1, S2, ..., Sn}, we are interested in finding the average of the cardiocyte's shapes representing the uncontracted phase.
A video from our lab depicting a contracting cardiocyte will typically comprise about one thousand frames. On average, about ninety percent of these frames will show the cardiocyte in its relaxed state.
Furthermore, the shapes of the cardiocyte in these frames are practically identical, then finding the average shape only guarantees a more unbiased measurement while providing for some regularization. Thus, in the interest of minimizing the overall computing time, we implemented a very simple averaging of contours in lieu of a shape averaging. The algorithm is as follows: 1) Identify the frames from each relaxed phase; 2) Obtain the contours of the shapes; 3) Resample the contours so that they will have the same number of points; 4) Find the average centroid point; 5) Interpolate the points in each contour using splines in polar coordinates; and 6) Interpolate the splines among the contours. Figure shows the effectiveness of this very simple contour averaging approach where the contours of a relaxed shape and a contracted shape were averaged. We show this averaging since the average of two uncontracted contours is practically indistinguishable from the two uncontracted contours.
It is worth noting here that the aforementioned shape averaging could have been accomplished within the same shape deformation approach due to Hong, Soatto and Vese [
20], where they define their shape average
M as the shape that is closest to the ensemble, by simultaneously minimizing their energy functional equivalent to Eq. (6). In other words, they look for
M that minimizes
They perform this optimization via alternating minimization and gradient descent. Implementing the optimization process proposed in [
20] was deemed to be too expensive a proposition for our purposes. Thus, we opted for performing the simpler average of contours which is two orders of magnitudes less expensive to run on a personal computer equipped with MATLAB™.