- Email: [email protected]
Joint edge detection and motion estimation of cardiac MR image sequence by a phase field method
ARTICLE IN PRESS Computers in Biology and Medicine 40 (2010) 21–28
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Joint edge detection and motion estimation of cardiac MR image sequence by a phase field method Abouzar Eslami a,, Mehran Jahed a, Tobias Preusser b a b
Electrical Engineering Department, Sharif University of Technology, Azadi Ave., Tehran, Iran Jacobs University Bremen and Fraunhofer MEVIS, Bremen, Germany
a r t i c l e in fo
abstract
Article history: Received 24 January 2009 Accepted 9 October 2009
In this paper a variational framework for joint segmentation and motion estimation is employed for inspecting heart in Cine MRI sequences. A functional including Mumford–Shah segmentation and optical flow based dense motion estimation is approximated using the phase-field technique. The minimizer of the functional provides an optimum motion field and edge set by considering both spatial and temporal discontinuities. Exploiting calculus of variation principles, multiple partial differential equations associated with the Euler–Lagrange equations of the functional are extracted, first. Next, the finite element method is used to discretize the resulting PDEs for numerical solution. Several simulation runs are used to test the convergence and the parameter sensitivity of the method. It is further applied to a comprehensive set of clinical data in order to compare with conventional cascade methods. Developmental constraints are identified as memory usage and computational complexities, which may be resolved utilizing sparse matrix manipulations and similar techniques. Based on the results of this study, joint segmentation and motion estimation outperforms previously reported cascade approaches especially in segmentation. Experimental results substantiated that the proposed method extracts the motion field and the edge set more precisely in comparison with conventional cascade approaches. This superior result is the consequence of simultaneously considering the discontinuity in both motion field and image space and including consequent frames (usually five) in our joint process functional. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Cardiac Cine MRI Segmentation Motion estimation Optical flow Mumford-Shah’s functional
1. Introduction The breath-held Cine MRI technique provides a detailed viewing of the beating heart consistent with the standard echocardiography scenes. It is used for accurate estimation of the ejection fraction and the stroke volume, investigation of regional wall motion abnormalities, as well as scar and infarct imaging. ECG gated and fast imaging technique of Cine MRI along with its ability in representing soft tissues has made it an important tool for inspection of cardiovascular system [1–4]. Cardiac motion estimation and segmentation are two intensively addressed topics in the biomedical literature. Inspection of heart deformation and strain requires the motion of all pixels (i.e. dense motion estimation), that can be estimated by using global motion estimation methods such as optical flow techniques. A popular optical flow algorithm is the Lucas–Kanade method, which estimates the motion locally, assuming that the velocity field is constant within a window [5]. The introduction of a global variational approach for the optical flow problem by Horn and Schunck [6] was followed by developments which successfully
Corresponding author.
E-mail addresses: [email protected], [email protected] (A. Eslami). 0010-4825/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2009.10.004
resolved issues such as smoothing of discontinuities and computational cost [7–9]. Mailloux et al. extended the optical flow algorithm by adding a linearity constraint to the motion field [10]. Later, Zini et al. added an additional incompressibility constraint [11]. Besides dense motion estimation, fine segmentation of cardiac chambers is of much interest in cardiac diagnosis [12–14]. Various methods such as active appearance and deformable models [12,15–18], fuzzy connectedness [19], and hybrid methods combining edge, region and shape information [20,21] are applied for segmentation and tracking of heart boundaries. Generally, employing priori information in the statistical approaches restricts the performance due to the dependence on sample size and training set comprehensiveness. Hence, knowledge free methods are considered in many applications and studies like this paper. Two-phase approaches with cascade processing, where either segmentation or motion estimation is performed first, are applied to the cardiac deformation analysis through various modalities in [22–25]. Theoretically, segmentation associates with discontinuity extraction in either the motion field or the image field. Moving objects with smooth boundary and distinctive objects without motion cause each one of the above two phase approaches to fail. Consequently, all these approaches suffer from an unreliable
ARTICLE IN PRESS 22
A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
decision on extracting discontinuities either on the motion or the image fields. Alternatively, joint processing of these tasks eliminates the necessity of such a decision. The idea of joint processing is proposed to cope with the ambiguity of the tasks’ priority in applications. It has been applied to mathematical image and motion processing such as joint segmentation and registration or motion estimation. Simultaneous segmentation and registration are recently studied by Droske and Ring in [26,27]. Kornprobst et al. considered piecewise smooth motion patterns as well as piecewise smooth objects on image sequences [28–30]. Their work focuses on the functions of bounded variation (BV) and proposes effective numerical implementation. Basically, two different measurements are employed in motion estimation. One is the mean square error (MSE) distortion between two frames after motion compensation (i.e. intensity matching). This is based on intensity constancy assumption, which depicts motion as the only source for intensity variation. A functional based on this measurement implements the warping theory also used in registration. Brox et al. minimized a functional of constancy measure (intensity matching) and segmentation approximated by level set method for joint segmentation and motion estimation [31]. Alternatively, in optical flow theory, by using the Taylor’s series expansion, the constancy measure is replaced with its first order approximation. As a consequence of this approximation, a minimization of the functional leads to linear partial differential equations that can be solved with finite element methods and numerical linear algebra, which is much faster than gradient descent methods required for solving PDEs associated with warping based functionals. Nesi related Mumford–Shah segmentation with the optical flow motion estimation in [32] and Keeling and Ring investigated the optical flow extraction with respect to optimization in [33]. Papenberg et al. recently proposed a total variation regularization of the optical flow involving higher order gradients [34]. In this paper an accurate and reliable knowledge-free method is developed in the spirit of [7] and applied to Cine MRI sequences. The functional for joint motion estimation and segmentation is developed based on Mumford–Shah segmentation and optical flow motion estimation. The minimization of the functional is accomplished by employing a phase-field approximation by deriving partial differential equations through the mathematical calculus of variation [35,36]. The minimizer corresponds to the optimal motion field and edge-set, considering discontinuities in both image and motion field. The phase-field approximation to the joint functional is inherited from Ambrosio and Tortorelli’s approach to the Mumford–Shah functional [37]. By exploiting a finite element method for solving the resulting partial differential equations, the whole process converges in just a few iterations. The paper is organized as follows: the next section provides a brief insight into the Mumford–Shah functional for image segmentation, as well as the mathematical principles for the joint segmentation and motion estimation approach. While the formulas are derived for 2D+T case, they can be easily generalized to 3D+ T. Section 4 discusses experimental results of applying the joint tasks method to both synthetic and clinical data. It also reports process features such as convergence rate and parameter selection considerations. The paper will be concluded in Section 5, while the appendix provides a concise definition of the matrices for an implementation of the method with finite elements.
2. Joint segmentation and motion estimation functional Image and motion segmentation by global optimization approaches comprises two major issues: mathematical descrip-
tion of the problem as an energy functional or probabilistic framework (Bayes or minimum description length) and its numerical solution. Typical probabilistic criteria involve conditional probability of segmentation that usually require powerful optimization algorithms such as simulated annealing and graduated non-convexity [38–40]. Shape prior for variational segmentation is also proposed based on a probabilistic definition of shapes [41]. This technique is usually employed for region segmentation (volumes in 3D and surfaces in 2D) and may be used for segmentation of myocardium and cardiac cavities. In this paper, we focus on extracting cardiac boundaries (epicardium and endocardium) based on a discontinuity extracting functional. The well referenced image edge segmentation functional is introduced by Mumford and Shah [42] as Z Z l m ðI I0 Þ2 dx þ jrIj2 dx þuHd1 ðGÞ ð1Þ EðI; GÞ ¼ 2 O 2 O\G The Mumford–Shah functional partitions the input image (I0 defined on an image domain O Rd ) into a piecewise smooth representation I and an edge set G. While I approximates I0 in a least mean square error sense, it shall be smooth distant from discontinuities (G). In addition, G should be smooth and small with respect to the (d 1)-dimensional Hausdorff measure ðHd1 Þ. A treatmeant for this problem has been proposed in the space of functions of bounded variation BV and SBV [43]. The level set and the phase field method are two mathematically proven tools for interface problems, which have been applied to the Mumford– Shah functional [37]. In the level set method, the edge set is represented by the zero-level set of a smooth function. Alternatively the phase field function ðzÞ is defined over the image such that it is small on edges and almost equal to 1 elsewhere. Based on the phase field representation of the edge set, the functional (1) can be written as Z l m 2 n ðI I0 Þ2 þ ðz þ ke ÞjrIj2 þ nejrzj2 þ ð1 zÞ2 dx EðI; zÞ ¼ 2 4e O 2 ð2Þ Here the parameter e controls the width of edges in the phasefield: the larger e, the coarser and more blurred the edge set is represented by the phase field. The small positive regularizing parameter ke ¼ OðeÞ 5 1, mathematically ensures strict coercivity 2 with respect to I. Since z vanishes over edges, the second term in Eq. (2) measures the smoothness of I distant from discontinuities ðO\GÞ. The third term measures the area of the edge set and the 2 fourth term limits the phase field to z 1 away from discontinuities. It has been proven that Eq. (2) converges to Eq. (1) for small e. Now suppose that @O is Lipschitz and I(x,t) represents an image sequence in a finite interval t A ½0; T. In this condition, the I and G parts of image sequence are the regular and singular parts with respect to space-time, respectively. So that for any frame at time t0 Gð; t0 Þ represents the edge set in that frame. In order to build a joint segmentation and motion estimation functional, we follow [7] by adding a motion estimation cost function to the Mumford–Shah functional. The optical flow equation relates the intensity variation of the object with its texture and velocity. Assuming that w ¼ ½1; m is the augmented motion vector of a pixel, from the brightness constancy assumption of the optical flow (ðIðt þ e; x þ emÞ ¼ const:Þ), we have
rI w ¼ 0
ð3Þ
The first order condition (Eq. (3)) is traditionally used in optical flow estimation functionals and further utilized in this paper. Higher order conditions associating with other terms in Taylor’s series approximation of the brightness constancy condition can also be included in the optical flow functional [34].
ARTICLE IN PRESS A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
Combining the Mumford–Shah functional with the optical flow condition and motion field regularity, the functional for joint segmentation and motion estimation can be defined as Z Z Z l m a EðI; G; wÞ ¼ ðI I0 Þ2 dw þ jrIj2 dw þ jw rIj2 dw 2 2 2 O
þ
b
Z
2
O\G
jrwj2 dw þ u
Z
T
O\G
Hd1 ðGðtÞÞdt
Since the phase-field represents the singularity of image sequence, the last term in the functional computes the motion field regularity excluding the edge set. It should be mentioned that in this functional, the discontinuity of the image and the motion field are represented simultaneously by one phase-field. Alternatively, it is possible to decouple motion and image singularities by involving a much finer scale parameter e^ 5 e for representing motion singularities. However, practically, it is hard to discretize the resulting functional since e is itself in the range of pixel size [7].
3. Numerical solution Let /df E; gS denote the variation of the energy functional E in direction g with respect to f. Assuming sufficiently smooth I, z, and w the variations of our functional (5) with respect to the unknowns are Z Z /dI E; J S ¼ lðI I0 ÞJ dw þ aðw rIÞðw O O Z 2 rJ Þ dw þ mðz þ ke ÞrI rJ dw O
/dz E; xS ¼
Z O
mzjrIj2 x dw þ
Z
bzjrwj2 x dw O
/du E; uS ¼
u ðz 1Þx dw 2e
2uerz rx þ
Z
aIx ðw rIÞu dw þ
Z
O
Z
bðz2 þ k A Þru ru dw
O
aIy ðw rIÞv dw þ O
ð4Þ
The integration in Eq. (4) is over space and time domain and dw is the spatio-temporal element [7]. The first two terms and the last one were previously defined in the Mumford–Shah functional (1) for segmentation. Additional motion dependent terms pertain to the optical flow (third term) and the motion field smoothness (fourth term). All the constant weights are supposed to be positive. The minimizer is the set of a piecewise smooth image sequence (I), a piecewise smooth motion vector field (w) and a discontinuity set (G). However, the variational minimization of the functional is difficult because of its dependency on the geometry of the discontinuity set. Some approaches are proposed in the literature for minimizing such a functional [44,45]. Brox et al. proposed a level set approximation that does not encode the motion concentrated on edges [31]. In this paper, to gain more flexibility and incorporate a simple scalable treatment, the phase-field approximation is employed based on the spirit of the works presented earlier [37,46,7]. Similar to the basic Mumford–Shah problem, the edge set is represented by an auxiliary variable (the phase-field variable) that vanishes on the edges and is almost equal to 1 distant from it. Defining e as the scale parameter, the functional yielding the smoothed image sequence (I), the phase field variable ðzÞ, and the motion field (w) is Z Z Z l a m 2 EðI; z; wÞ ¼ ðI I0 Þ2 dw þ jw rIj2 dw þ ðz þ ke ÞjrIj2 dw O2 O2 O2 Z Z u b 2 ðz þ ke Þjrwj2 dw þ uejrzj2 þ ð1 zÞ2 dw þ 4e O O2 ð5Þ
Z O
/dv E; vS ¼
0
O\G
þ
23
Z
bðz2 þ k A Þrv rv dw O
All J , x, u and v are scalar test functions and the functional (5) is rewritten expanding w ¼ ½1; u; vT in 2D space. So, the Euler– Lagrange equations associated with I and z are m a 2 I div ðz þ ke ÞrI þ wðw rIÞ ¼ I0
l
l
z
m 1 1 jrIj2 þ 2 þ bjrwj2 Dz ¼ 2 2ev 4e 4e
Functional (5) is strictly convex and even quadratic with respect to these two variables for each feasible and fixed motion field (w). Hence, a fast and effective numerical solution to the set of linear PDEs can be achieved through the finite element method [47]. Discretizing the spatio-temporal domain T O by a regular hexahedral grid leads to the system of equations associated with the Euler–Lagrange equations: ðAðw; zÞ þ CÞI ¼ CI 0
ð6Þ
ðCðI; wÞ þDÞZ ¼ E
ð7Þ
The finite element method is also employed for the phase-field approximation of the Mumford–Shah functional in [48,49] which is an extension to the approach first presented by Chambolle and Dal Maso [50]. Finally, gradient descent of EðI; z; oÞ with respect to the motion field variables ðu and vÞ is given by
b du ¼ ðw rIÞ@x1 I þ divððz2 þke ÞruÞ
ð8Þ
b dv ¼ ðw rIÞ@x2 I þ divððz2 þ ke ÞrvÞ
ð9Þ
a
a
where @xj I ¼ @I=@xj . For a fixed e the minimization strategy can be summarized as Step 0: Initialize I ¼ I0 ; z 1; u 0 and v 0. Step 1: Solve Eq. (6) for fixed u; v and z. Step 2: Solve Eq. (7) for fixed u, v, I. Step 3: Compute descent of E with respect to u from Eq. (8) when z; v; I are fixed. Step 4: Compute descent direction of E with respect to v from Eq. (8) when z; u; I are fixed. Step 5: Return to step 1. In contrast to the level set approach, the initialization of the phasefield is simpler as there is no need for a precise initialization close to the total edge set. Besides, the phase-field approximation comes along with a natural scale parameter. The smoothness of the energy depends on the scale parameter e and larger values of e yield coarse and smooth approximations of the images, the phase field and the deformation. In a multiscale strategy, the whole minimization process (steps 0–5) will be computed starting from a large value of e. This results in a coarse approximation as the stationary point in the simplified energy landscape. The process will be followed iteratively, reducing the scale parameter by taking the solution of the previous scale as the new initial guess on the next finer scale. Finally, e has to be of the same order as the grid size for the finest scale. It should be noted that due to the special definition of the functional, the motion field regularity is applied to the smooth part of image. Consequently, it is critical to
ARTICLE IN PRESS 24
A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
finalize the multiscale approach with small e in order to reach a regular motion field.
4. Experimental results The proposed joint segmentation and motion estimation model is implemented by using an Intels Core2TM Duo processor (T7300) at 2.0 GHz with a MATLAB platform and a Windows VistaTM operating system. The implemented solution is based on just one scale parameter instead of the multiscale strategy mentioned before, due to the small size of images in the clinical data. In the small frames of clinical Cine MRI, the anatomical objects are so close that they may be merged by choosing large scale parameters. Hence, in our application it is impractical to exploit a multiscale approach starting from large scales. Defining h as the grid width, experimentally preferred values of the parameters are b ¼ l ¼ 1, a ¼ 105 h2 , m ¼ h2 , and e ¼ h=4 [7]. The grid width of the FEM in this case can be computed from h ¼ minfhx ; hy ; ht g where hx ¼ 1=ðNx 1Þ, hy ¼ 1=ðNy 1Þ and ht ¼ 1=ðNt 1Þ are the inverses of number of pixels in each direction and the inverse of frames number. Ahead of the validation tests, it is necessary to discuss implementation considerations such as memory usage and convergence speed. Basically, the proposed joint method of this paper requires plenty of memory because the motion field, the phase field and all finite element matrices should be available simultaneously. It is strongly recommended to exploit sparse matrices and their operations for the finite element solution. But it still needs much more memory in comparison with tracking methods to contain the motion and phase-field variables. The most effective scheme is to limit the number of frames. It should be considered that for each pixel there are three variables namely, one phase-field and two motion field variables. So, omitting the finite element matrices, by increasing the number of frames by one, the number of unknown variables grows by three times of image pixel number. Fig. 1 depicts the convergence rate of the method from the initial condition of I= I0, z 1; u 0 and v 0. The evolution of
different terms in the functional is depicted in Fig. 1 for a normal Cine MRI of five frames. While the total energy and the optical R flow term ð a2jw rIj2 dwÞ are decreasing, the motion field Rb 2 regularity ð 2ðz þke Þjrwj2 dwÞ increases during evolution of motion field from its fully regular initial state (u 0 and v 0) to its minimizer. The decay of the total energy implies convergence speed of proposed method. In order to check the performance of the joint method, it is applied to a synthesized image sequence in which a disk with smooth gradient moves north west with a motion vector of m ¼ ½1; 1T . The sequence is degraded by inhomogeneous light to roughly simulate inappropriate imaging. The light decreases from bottom-right to upper-left where the object boundary completely fades. Fig. 2 shows the extracted edge set of the proposed joint process in comparison with a one task segmentation process [51]. In theory, all image based processes fail to extract the whole disk because of a lack of the intensity gradient in dark sites of image. Alternatively, the joint process can extract the disk completely relying on motion discontinuities in the absence of image discontinuities (Fig. 2(b)). As a consequence of considering both spatial and temporal discontinuities simultaneously, the extracted edge set successfully copes with the light condition and results in a more reliable edge set especially where the image gradient degrades. In the optical flow literature, the motion field regularity term is cited as desirable and it is shown to even reciprocate flickering [7]. But in the case of cardiac MR with rapid heart deformation and low frame rate, this term results in fictitious motion vectors especially during the systole phase. Instant contraction of the left ventricle at the beginning of the systole (much like all the other phase transitions) can be missed due to the motion smoothness term. On the other hand, the smoothness term is necessary in order to cope with the functional’s ill-posedness and noise which has destructive effects in the optical flow approach. In the optical flow equation (Eq. (3)), the variation of the pixels’ intensity is interpreted as their motion. By omitting the motion smoothness term, the estimated motion of adjacent pixels can be completely disparate. In order to inspect the latency and ability of the proposed method in tracking erratic motions, a simulated
Fig. 1. Evolution of different energy terms for a normal Cine MRI with five frames. (a) Optical flow term Total energy (Eq. (5)).
Ra Rb 2 2 2 2 jw; r Ij dw, (b) Motion field regularity 2 ðz ; k A Þjrwj dw, (c)
ARTICLE IN PRESS A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
25
Fig. 2. Extracted edge set from a synthetic image sequence with light inhomogeneity: (a) moving disk with light degeneration, (b) color coded motion field extracted by joint process, (c) segmentation result of [51] and (d) the phase-field function (z) of joint process.
Fig. 3. The effect of motion field regularity in tracking of rapid and sudden motion variations. Solid: ideal, dashed: b = 1, dotted: b =0.6.
pulse-like motion is synthesized in an images sequence of 10 frames. After three frames a textured object begins to move north east with m ¼ ½2; 2T and keeps its motion until the end of the sequence. Fig. 3 shows the latency caused by the motion smoothness term with b ¼ 1. While the solid line in Fig. 3 shows the norm of the real motion, the other lines show the object’s estimated motion norm with different b. Reducing b increases the proposed method’s ability in tracking of objects with rapid and sudden movements. However, there is a trade off as very small values of b reduce the robustness and can cause undesirable and irregular motion fields. For validation, the proposed method is applied to a set of clinical cardiac Cine MR image sequences. The dataset consists of
17 cases (seven females and 10 males) and a total of five healthy hearts. Each image sequence has 20 frames of 125 190 images from one heart cycle in long axes view. Our parameter setting for the clinical data was l = 1, b = 0, 5, a =105 h 2, m = 10 1 h 2 and e = h/10. The values of b and e were reduced to extract a fine edge set and to follow rapid movements. Processing each image sequence (size of 125 190 20) took 160s in average (with a maximum of 173s) using the implementation specifications stated before. Fig. 4 shows a color coded motion field extracted by the joint method associated with a normal heart. In order to provide benchmark, a specialist in interpreting the cardiac MRI delineated all the cardiac boundaries in all frames. Fig. 5 provides a comparison between manual and automatic
ARTICLE IN PRESS 26
A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
contours at four different events: end of diastole (ED), beginning of systole (BS), mid-systole (MS) and end of systole (ES). While the phase-field function is depicted by gray, white regions represent its difference from manual boundaries. Exploiting contrast agent, blood flow can be seen as well as wall motion in Cine MRI and can be extracted by dense motion estimation strategy in Eq. (5). Unfortunately there is no known reference for quantitative test of estimated motion by using Cine MRI, solely. But, it is possible to inspect the accuracy of the endocardium motion estimation by computing the difference
between the acquired edge set and the manual border delineated by an expert. Since the atrio-ventricular valves are not discriminative in Cine MR images, the heart is partitioned just in left and right sides, and Tanimoto’s index is used instead of ejection fraction [19]. Assuming that Oe and Om are extracted and manually delineated regions, respectively, Tanimoto’s index is defined as T1 ¼
JOe \ Om J JOe [ Om J
This true positive criterion determines the common area of automatically and manually extracted objects and implies both segmentation and motion estimation accuracies [52]. Tables 1 and 2 include the mean and standard deviation of T1 at four different states of beating heart in comparison with two other tracking based approaches [23,25]. The acquired results of Table 1 Mean of TI at ED, ES and MS phases.
Joint method Cascade method of [23] Level set and tracking of [25]
ED
BS
MS
ES
95.77 91.13 92.52
90.56 90.21 91.02
96.36 92.52 93.43
96.87 92.65 94.78
Table 2 Standard deviation of TI at ED, ES and MS phases.
Fig. 4. (a) Extracted motion associating with a normal heart at mid-systole besides applied velocity color map.
Joint method Cascade method of [23] Level set and tracking of [25]
ED
BS
MS
ES
1.03 1.56 1.24
3.10 3.62 2.36
2.12 2.39 2.43
1.19 2.68 1.32
Fig. 5. Extracted edge set in comparison with manually delineated boundary at four events of heart cycle: (a) end of diastole, (b) beginning of systole, (c) mid-systole and (d) end of systole.
ARTICLE IN PRESS A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
the joint process are superior to the results of cascade processes as a consequence of considering more than one frame for segmentation. But, no privilege is gained at the beginning of systole when a capricious contraction happens. This has been previously investigated as the effect of motion regularity term that may degrade the results in two or three frames, utmost (cf. Fig. 3). Although it is recommended to reduce the effect of the motion regularity term when the frame rate is low, it should not be suppressed to reach in appropriate smooth motion fields.
Ci;j ¼
Di;j ¼
Ei ¼
Z
27
m 1 jrIj2 þ 2 þ bjrwj2 ci ci dw 2ev 4e
Z
1 4e2
rci rcj dw Z
ci dw
It should be mentioned that w is spatio-temporal variable.
5. Conclusion
References
In this paper, a mathematical method for the simultaneous segmentation and motion estimation of cardiac MR image sequences is proposed. The method is based on the calculus of variation and minimizes a functional defined by using a combination of an optical flow functional and the Mumford–Shah segmentation functional. It is shown that finite element discretization and phase-field approximation provide an appropriate solution for the joint segmentation and motion estimation problem due to the fast convergence and the simple initialization. The edge set and the dense motion field associated with the minimizer of the functional are optimum and the proposed method outperforms conventional approaches. Consideration of spatial and temporal discontinuities simultaneously is the most important feature of proposed method. The tests performed on simulated image sequences imply the ability of the joint process in extracting even weak edges by using motion discontinuities. We have discussed the selection of the method’s parameters and especially the effect of the motion field regularity term on the method’s ability in tracking erratic movements. It can be inferred from the clinical experiments that the proposed method outperforms conventional cascade approaches without requiring any pre- or post-processing and even fine initialization. Since it is completely automatic and free from a priori knowledge, it can be employed for various heart shapes and functionality.
[1] D.J. Pennell, U.P. Sechtem, C.B. Higgins, W.J. Manning, G.M. Pohost, F.E. Rademakers, A.C. Van Rossum, L.J. Shaw, E.K. Yucel, Clinical indications for cardiovascular magnetic resonance (CMR): consensus panel report, Elsevier Eur. Heart J. 25 (2004) 1940–1965. [2] /http://www.scmr.org/S. [3] /http://atlas.scmr.org/cineplayer.htmlS. [4] /http://en.wikipedia.org/wiki/Cardiovascular_magnetic_resonanceS. [5] B. Lucas, T. Kanade, An iterative image restoration technique with an application to stereo vision, in: Proceedings of Darpa Uworkshop, 1981, pp. 121–130. [6] B. Horn, B. Schunk, Determining optical flow, Artif. Intell. 17 (1981) 185–204. [7] T. Preusser, M. Droske, C.S. Garbe, A. Telea, M. Rumpf, A phase field method for joint denoising, edge detection, and motion estimation in image sequence processing, SIAM J. Appl. Math. 68 (3) (2007) 599–618. [8] H.H. Nagel, W. Enkelmann, An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences, IEEE Trans. Pattern Anal. Mach. Intell. 8 (1986) 565–593. [9] A. Bruhn, J. Weickert, A confidence measure for variational optic flow methods, in: R. Klette, R. Kozera, L. Noakes, J. Weickert (Eds.), Geometric Properties from Incomplete Data, Springer, Dordrecht, The Netherlands, 2006, pp. 283–297. [10] G.E. Mailloux, F. Langlois, P.Y. Simard, M. Bertrand, Restoration of the velocity field of the heart from two-dimensional echocardiograms, IEEE Trans. Med. Imag. 8 (2) (1989) 143–153. [11] G. Zini, A. Sarti, C. Lamberti, Application of continuum theory and multi-grid methods to motion evaluation from 3D echocardiography, IEEE Trans. Ultrason. Ferroelect. Freq. Control 44 (3) (1997) 297–308. [12] J.G. Bosch, S.C. Mitchell, B.P.F. Lelieveldt, F. Nijland, O. Kamp, M. Sonka, J.H.C. Reiber, Automatic segmentation of echocardiographic sequences by active appearance motion models, IEEE Trans. Med. Imag. 21 (11) (2002) 1374–1383. [13] C. Pluempitiwiriyawej, J.M.F. Moura, Y. Lin Wu, C. Ho, Stacs, New active contour scheme for cardiac MR image segmentation, IEEE Trans. Med. Imag. 24 (5) (2005) 593–603. [14] Y. Wang, Y. Jia, Segmentation of the left ventricle from MR images via snake models incorporating shape similarities, IEEE Int. Conf. Image Process. 45 (2006) 213–216. [15] H. Van Assen, M. Danilouchkine, F. Behloul, H. Lamb, R. Van Der Geest, J. Reiber, B. Lelieveldt, Cardiac LV segmentation using a 3d active shape model driven by fuzzy inference, in: Proceedings of Medical Image Computing And Computer-Assisted Intervention (Miccai’03), Montreal, QC, Canada, November 16–18, 2003, pp. 533–540. [16] S. Mitchell, B. Lelieveldt, R. Van Der Geest, H. Bosch, J. Reiber, M. Sonka, Multistage hybrid active appearance model matching: segmentation of left and right ventricles in cardiac MR images, IEEE Trans. Med. Imag. 20 (5) (2001) 415–423. [17] M. Kaus, J. Von Berg, W. Niessen, V. Pekar, Automated segmentation of the left ventricle in cardiac MRI, Med. Image Anal. 8 (2004) 245–254. [18] J. Lotjonen, S. Kivisto, J. Koikkalainen, D. Smutek, K. Laurema, Statistical shape model of atria, ventricles and epicardium from short- and long-axis MR images, Med. Image Anal. 8 (2004) 371–386. [19] A. Pednekar, U. Kurkure, R. Muthupillai, S. Flamm, I.A. Kakadiaris, Automated left ventricular segmentation in cardiac MRI, IEEE Trans. Biomed. Eng. 53 (7) (2006) 1425–1428. [20] N. Paragios, Shape-based segmentation and tracking in cardiac image analysis, IEEE Trans. Med. Imag. 22 (6) (2003) 773–776. [21] A.S. Pednekar, I.A. Kakadiaris, U. Kurkure, R. Muthupillai, S.D. Flamm, Intensity and morphology-based energy minimization for the automatic segmentation of the myocardium, in: Proceedings of IEEE Workshop On Variational, Geometric and Level Set Methods in Computer Vision (Vlsm’03), Nice, France, October 2003, pp. 185–192. [22] D. Metaxas, L. Axel, A. Montillo, K. Park, Automated Segmentation, Motion Estimation of LV/RV Motion From MRI, in: EMBS/BMES Conference, 2002. Proceedings of the Second Joint, vol. 2, 23–26 October 2002, pp. 1099–1100. [23] J. Bosch, S. Mitchell, B. Lelieveldt, F. Nijland, O. Kamp, M. Sonka, J. Reiber, Automatic segmentation of echocardiographic sequences by active appearance motion models, IEEE Trans. Med. Imag. 21 (11) (2002) 1374–1383. [24] L. Legrand, C. Bordier, A. Lalande, P. Walker, F. Brunotte, C. Quantin, Magnetic resonance image segmentation and heart motion tracking with an active mesh based system, Comput. Cardiol. 22–25 (2002) 177–180.
Conflict of interest statement None declared.
Acknowledgments The authors would like to thank Mr. Alinaghizadeh and the Nour clinic for providing Cine MRI sequences and Dr. Mohebbi for delineation and tracking of heart chambers. The authors would also in particular like to thank Dr. Razvan for giving mathematical insight into variational theory.
Appendix A The system of equations from the finite element method contains Z I0 ci dw I 0 ¼ ½I0i ¼
C ¼ ½Ci;j ¼
Z
ci cj dw
Z m 2 A ¼ ½Aðw; zÞi;j ¼ ðz þke Þrci l Z a rcj dw þ ðrci wÞðrcj wÞ dw
l
ARTICLE IN PRESS 28
A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
¨ [25] M. Suhling, M. Arigovindan, C. Jansen, P. Hunziker, M. Unser, Myocardial motion analysis from B-mode echocardiograms, IEEE Trans. Image Process. 14 (4) (2005) 525–536. [26] M. Droske, W. Ring, A Mumford–Shah level-set approach for geometric image registration, SIAM J. Appl. Math. 66 (2006) 2127–2148. [27] M. Droske, W. Ring, M. Rumpf, Mumford–Shah based registration: a comparison of a level set and a phase field approach, Comput. Vis. Sci. 12 (3) (2008) 101–114. [28] P. Kornprobst, R. Deriche, G. Aubert, Image sequence analysis via partial differential equations, J. Math. Imag. Vision 11 (1999) 5–26. [29] G. Aubert, P. Kornprobst, A mathematical study of the relaxed optical flow problem in the space Bv(O), SIAM J. Math. Anal. 30 (1999) 1282–1308. [30] G. Aubert, R. Deriche, P. Kornprobst, Computing optical flow via variational techniques, SIAM J. Appl. Math. 60 (1999) 156–182. [31] T. Brox, A. Bruhn, J. Weickert, Variational motion segmentation with level sets, in: H. Bischof, A. Leonardis, A. Pinz (Eds.), Computer Vision–Eccv 2006, Lecture Notes in Computer Science, vol. 3951, Springer, Berlin, 2006, pp. 471–483. [32] P. Nesi, Variational approach to optical flow estimation managing discontinuities, Image Vision Comput. 11 (1993) 419–439. [33] S.L. Keeling, W. Ring, Medical image registration and interpolation by optical flow with maximal rigidity, J. Math. Imag. Vision 23 (2005) 47–65. [34] N. Papenberg, A. Bruhn, T. Brox, S. Didas, J. Weickert, Highly accurate optic flow computation with theoretically justified warping, Int. J. Comput. Vision 67 (2006) 141–158. ¨ [35] B.V. Brunt, The Calculus of Variations, Birkhauser, 2004. [36] H. Sagan, Introduction to the Calculus of Variations, Courier Dover Publications, 1992. [37] L. Ambrosio, V.M. Tortorelli, On the approximation of free discontinuity problems, Boll. Un. Mat. Ital. B 7 (6) (1992) 105–123. [38] N. Vasconcelos, A. Lippman, Empirical Bayesian motion segmentation, IEEE Trans. Pattern Anal. Mach. Intell. 23 (2) (2001) 217–221. [39] S.C. Zhu, A. Yuille, Region competition: unifying snakes, region growing, and Bayes/MDLfor multiband image segmentation, IEEE Trans. Pattern Anal. Mach. Intell. 18 (9) (1996) 884–900.
[40] Y. Weiss, E.H. Adelson, A unified mixture framework for motion segmentation: incorporating spatial coherence and estimating the number of models, Proc. IEEE Conf. Comput. Vision Pattern Recognition 45 (1996) 321–326. [41] D. Cremers, F.R. Schmidt, F. Barthel, Shape priors in variational image segmentation: convexity, lipschitz continuity and globally optimal solutions, Proc. IEEE Conf. Comput. Vision Pattern Recognition 45 (2008) 1–6. [42] D. Mumford, J. Shah, Optimal approximation by piecewise smooth functions and associated variational problems, Comm. Pure Appl. Math. 42 (1989) 577–685. [43] L. Ambrosio, N. Fusco, D. Pallara, Functions of Bounded Variation and Free Discontinuity Problems, Oxford University Press, New York, 2000. [44] D. Cremers, S. Soatto, Motion competition: a variational approach to piecewise parametric motion segmentation, Int. J. Comput. Vision 62 (2005) 249–265. [45] D. Cremers, T. Kohlberger, C. SchnOrr, Nonlinear shape statistics in Mumford– Shah based segmentation, in: Proceedings of the Seventh European Conference on Computer Vision, Copenhagen, Lecture Notes in Computer Science, vol. 2351, Springer, Berlin, 2002, pp. 93–108. [46] G. Aubert, R. Deriche, P. Kornprobst, Computing optical flow via variational techniques, SIAM J. Appl. Math. 60 (1999) 156–182. [47] P. Solin, Partial Differential Equations and the Finite Element Method, WileyInterscience, 2005. [48] B. Bourdin, Image segmentation with a finite element method, M2AN Math. Model. Numer. Anal. 33 (1999) 229–244. [49] B. Bourdin, A. Chambolle, Implementation of an adaptive finite-element approximation of the Mumford–Shah functional, Numer. Math. 85 (2000) 609–646. [50] A. Chambolle, G. Dal Maso, Discrete approximation of the Mumford–Shah functional in dimension two, M2AN Math. Model. Numer. Anal. 33 (1999) 651–672. [51] M. Mora, C. Tauber, H. Batatia, Robust level set for heart cavities detection in ultrasound images, Comput. Cardiol. 45 (2005) 235–238. [52] S. Theodoridis, K. Koutroumbas, Pattern Recognition, Academic Press, USA, 1999 p. 366.
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Joint edge detection and motion estimation of cardiac MR image sequence by a phase field method Abouzar Eslami a,, Mehran Jahed a, Tobias Preusser b a b
Electrical Engineering Department, Sharif University of Technology, Azadi Ave., Tehran, Iran Jacobs University Bremen and Fraunhofer MEVIS, Bremen, Germany
a r t i c l e in fo
abstract
Article history: Received 24 January 2009 Accepted 9 October 2009
In this paper a variational framework for joint segmentation and motion estimation is employed for inspecting heart in Cine MRI sequences. A functional including Mumford–Shah segmentation and optical flow based dense motion estimation is approximated using the phase-field technique. The minimizer of the functional provides an optimum motion field and edge set by considering both spatial and temporal discontinuities. Exploiting calculus of variation principles, multiple partial differential equations associated with the Euler–Lagrange equations of the functional are extracted, first. Next, the finite element method is used to discretize the resulting PDEs for numerical solution. Several simulation runs are used to test the convergence and the parameter sensitivity of the method. It is further applied to a comprehensive set of clinical data in order to compare with conventional cascade methods. Developmental constraints are identified as memory usage and computational complexities, which may be resolved utilizing sparse matrix manipulations and similar techniques. Based on the results of this study, joint segmentation and motion estimation outperforms previously reported cascade approaches especially in segmentation. Experimental results substantiated that the proposed method extracts the motion field and the edge set more precisely in comparison with conventional cascade approaches. This superior result is the consequence of simultaneously considering the discontinuity in both motion field and image space and including consequent frames (usually five) in our joint process functional. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Cardiac Cine MRI Segmentation Motion estimation Optical flow Mumford-Shah’s functional
1. Introduction The breath-held Cine MRI technique provides a detailed viewing of the beating heart consistent with the standard echocardiography scenes. It is used for accurate estimation of the ejection fraction and the stroke volume, investigation of regional wall motion abnormalities, as well as scar and infarct imaging. ECG gated and fast imaging technique of Cine MRI along with its ability in representing soft tissues has made it an important tool for inspection of cardiovascular system [1–4]. Cardiac motion estimation and segmentation are two intensively addressed topics in the biomedical literature. Inspection of heart deformation and strain requires the motion of all pixels (i.e. dense motion estimation), that can be estimated by using global motion estimation methods such as optical flow techniques. A popular optical flow algorithm is the Lucas–Kanade method, which estimates the motion locally, assuming that the velocity field is constant within a window [5]. The introduction of a global variational approach for the optical flow problem by Horn and Schunck [6] was followed by developments which successfully
Corresponding author.
E-mail addresses: [email protected], [email protected] (A. Eslami). 0010-4825/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2009.10.004
resolved issues such as smoothing of discontinuities and computational cost [7–9]. Mailloux et al. extended the optical flow algorithm by adding a linearity constraint to the motion field [10]. Later, Zini et al. added an additional incompressibility constraint [11]. Besides dense motion estimation, fine segmentation of cardiac chambers is of much interest in cardiac diagnosis [12–14]. Various methods such as active appearance and deformable models [12,15–18], fuzzy connectedness [19], and hybrid methods combining edge, region and shape information [20,21] are applied for segmentation and tracking of heart boundaries. Generally, employing priori information in the statistical approaches restricts the performance due to the dependence on sample size and training set comprehensiveness. Hence, knowledge free methods are considered in many applications and studies like this paper. Two-phase approaches with cascade processing, where either segmentation or motion estimation is performed first, are applied to the cardiac deformation analysis through various modalities in [22–25]. Theoretically, segmentation associates with discontinuity extraction in either the motion field or the image field. Moving objects with smooth boundary and distinctive objects without motion cause each one of the above two phase approaches to fail. Consequently, all these approaches suffer from an unreliable
ARTICLE IN PRESS 22
A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
decision on extracting discontinuities either on the motion or the image fields. Alternatively, joint processing of these tasks eliminates the necessity of such a decision. The idea of joint processing is proposed to cope with the ambiguity of the tasks’ priority in applications. It has been applied to mathematical image and motion processing such as joint segmentation and registration or motion estimation. Simultaneous segmentation and registration are recently studied by Droske and Ring in [26,27]. Kornprobst et al. considered piecewise smooth motion patterns as well as piecewise smooth objects on image sequences [28–30]. Their work focuses on the functions of bounded variation (BV) and proposes effective numerical implementation. Basically, two different measurements are employed in motion estimation. One is the mean square error (MSE) distortion between two frames after motion compensation (i.e. intensity matching). This is based on intensity constancy assumption, which depicts motion as the only source for intensity variation. A functional based on this measurement implements the warping theory also used in registration. Brox et al. minimized a functional of constancy measure (intensity matching) and segmentation approximated by level set method for joint segmentation and motion estimation [31]. Alternatively, in optical flow theory, by using the Taylor’s series expansion, the constancy measure is replaced with its first order approximation. As a consequence of this approximation, a minimization of the functional leads to linear partial differential equations that can be solved with finite element methods and numerical linear algebra, which is much faster than gradient descent methods required for solving PDEs associated with warping based functionals. Nesi related Mumford–Shah segmentation with the optical flow motion estimation in [32] and Keeling and Ring investigated the optical flow extraction with respect to optimization in [33]. Papenberg et al. recently proposed a total variation regularization of the optical flow involving higher order gradients [34]. In this paper an accurate and reliable knowledge-free method is developed in the spirit of [7] and applied to Cine MRI sequences. The functional for joint motion estimation and segmentation is developed based on Mumford–Shah segmentation and optical flow motion estimation. The minimization of the functional is accomplished by employing a phase-field approximation by deriving partial differential equations through the mathematical calculus of variation [35,36]. The minimizer corresponds to the optimal motion field and edge-set, considering discontinuities in both image and motion field. The phase-field approximation to the joint functional is inherited from Ambrosio and Tortorelli’s approach to the Mumford–Shah functional [37]. By exploiting a finite element method for solving the resulting partial differential equations, the whole process converges in just a few iterations. The paper is organized as follows: the next section provides a brief insight into the Mumford–Shah functional for image segmentation, as well as the mathematical principles for the joint segmentation and motion estimation approach. While the formulas are derived for 2D+T case, they can be easily generalized to 3D+ T. Section 4 discusses experimental results of applying the joint tasks method to both synthetic and clinical data. It also reports process features such as convergence rate and parameter selection considerations. The paper will be concluded in Section 5, while the appendix provides a concise definition of the matrices for an implementation of the method with finite elements.
2. Joint segmentation and motion estimation functional Image and motion segmentation by global optimization approaches comprises two major issues: mathematical descrip-
tion of the problem as an energy functional or probabilistic framework (Bayes or minimum description length) and its numerical solution. Typical probabilistic criteria involve conditional probability of segmentation that usually require powerful optimization algorithms such as simulated annealing and graduated non-convexity [38–40]. Shape prior for variational segmentation is also proposed based on a probabilistic definition of shapes [41]. This technique is usually employed for region segmentation (volumes in 3D and surfaces in 2D) and may be used for segmentation of myocardium and cardiac cavities. In this paper, we focus on extracting cardiac boundaries (epicardium and endocardium) based on a discontinuity extracting functional. The well referenced image edge segmentation functional is introduced by Mumford and Shah [42] as Z Z l m ðI I0 Þ2 dx þ jrIj2 dx þuHd1 ðGÞ ð1Þ EðI; GÞ ¼ 2 O 2 O\G The Mumford–Shah functional partitions the input image (I0 defined on an image domain O Rd ) into a piecewise smooth representation I and an edge set G. While I approximates I0 in a least mean square error sense, it shall be smooth distant from discontinuities (G). In addition, G should be smooth and small with respect to the (d 1)-dimensional Hausdorff measure ðHd1 Þ. A treatmeant for this problem has been proposed in the space of functions of bounded variation BV and SBV [43]. The level set and the phase field method are two mathematically proven tools for interface problems, which have been applied to the Mumford– Shah functional [37]. In the level set method, the edge set is represented by the zero-level set of a smooth function. Alternatively the phase field function ðzÞ is defined over the image such that it is small on edges and almost equal to 1 elsewhere. Based on the phase field representation of the edge set, the functional (1) can be written as Z l m 2 n ðI I0 Þ2 þ ðz þ ke ÞjrIj2 þ nejrzj2 þ ð1 zÞ2 dx EðI; zÞ ¼ 2 4e O 2 ð2Þ Here the parameter e controls the width of edges in the phasefield: the larger e, the coarser and more blurred the edge set is represented by the phase field. The small positive regularizing parameter ke ¼ OðeÞ 5 1, mathematically ensures strict coercivity 2 with respect to I. Since z vanishes over edges, the second term in Eq. (2) measures the smoothness of I distant from discontinuities ðO\GÞ. The third term measures the area of the edge set and the 2 fourth term limits the phase field to z 1 away from discontinuities. It has been proven that Eq. (2) converges to Eq. (1) for small e. Now suppose that @O is Lipschitz and I(x,t) represents an image sequence in a finite interval t A ½0; T. In this condition, the I and G parts of image sequence are the regular and singular parts with respect to space-time, respectively. So that for any frame at time t0 Gð; t0 Þ represents the edge set in that frame. In order to build a joint segmentation and motion estimation functional, we follow [7] by adding a motion estimation cost function to the Mumford–Shah functional. The optical flow equation relates the intensity variation of the object with its texture and velocity. Assuming that w ¼ ½1; m is the augmented motion vector of a pixel, from the brightness constancy assumption of the optical flow (ðIðt þ e; x þ emÞ ¼ const:Þ), we have
rI w ¼ 0
ð3Þ
The first order condition (Eq. (3)) is traditionally used in optical flow estimation functionals and further utilized in this paper. Higher order conditions associating with other terms in Taylor’s series approximation of the brightness constancy condition can also be included in the optical flow functional [34].
ARTICLE IN PRESS A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
Combining the Mumford–Shah functional with the optical flow condition and motion field regularity, the functional for joint segmentation and motion estimation can be defined as Z Z Z l m a EðI; G; wÞ ¼ ðI I0 Þ2 dw þ jrIj2 dw þ jw rIj2 dw 2 2 2 O
þ
b
Z
2
O\G
jrwj2 dw þ u
Z
T
O\G
Hd1 ðGðtÞÞdt
Since the phase-field represents the singularity of image sequence, the last term in the functional computes the motion field regularity excluding the edge set. It should be mentioned that in this functional, the discontinuity of the image and the motion field are represented simultaneously by one phase-field. Alternatively, it is possible to decouple motion and image singularities by involving a much finer scale parameter e^ 5 e for representing motion singularities. However, practically, it is hard to discretize the resulting functional since e is itself in the range of pixel size [7].
3. Numerical solution Let /df E; gS denote the variation of the energy functional E in direction g with respect to f. Assuming sufficiently smooth I, z, and w the variations of our functional (5) with respect to the unknowns are Z Z /dI E; J S ¼ lðI I0 ÞJ dw þ aðw rIÞðw O O Z 2 rJ Þ dw þ mðz þ ke ÞrI rJ dw O
/dz E; xS ¼
Z O
mzjrIj2 x dw þ
Z
bzjrwj2 x dw O
/du E; uS ¼
u ðz 1Þx dw 2e
2uerz rx þ
Z
aIx ðw rIÞu dw þ
Z
O
Z
bðz2 þ k A Þru ru dw
O
aIy ðw rIÞv dw þ O
ð4Þ
The integration in Eq. (4) is over space and time domain and dw is the spatio-temporal element [7]. The first two terms and the last one were previously defined in the Mumford–Shah functional (1) for segmentation. Additional motion dependent terms pertain to the optical flow (third term) and the motion field smoothness (fourth term). All the constant weights are supposed to be positive. The minimizer is the set of a piecewise smooth image sequence (I), a piecewise smooth motion vector field (w) and a discontinuity set (G). However, the variational minimization of the functional is difficult because of its dependency on the geometry of the discontinuity set. Some approaches are proposed in the literature for minimizing such a functional [44,45]. Brox et al. proposed a level set approximation that does not encode the motion concentrated on edges [31]. In this paper, to gain more flexibility and incorporate a simple scalable treatment, the phase-field approximation is employed based on the spirit of the works presented earlier [37,46,7]. Similar to the basic Mumford–Shah problem, the edge set is represented by an auxiliary variable (the phase-field variable) that vanishes on the edges and is almost equal to 1 distant from it. Defining e as the scale parameter, the functional yielding the smoothed image sequence (I), the phase field variable ðzÞ, and the motion field (w) is Z Z Z l a m 2 EðI; z; wÞ ¼ ðI I0 Þ2 dw þ jw rIj2 dw þ ðz þ ke ÞjrIj2 dw O2 O2 O2 Z Z u b 2 ðz þ ke Þjrwj2 dw þ uejrzj2 þ ð1 zÞ2 dw þ 4e O O2 ð5Þ
Z O
/dv E; vS ¼
0
O\G
þ
23
Z
bðz2 þ k A Þrv rv dw O
All J , x, u and v are scalar test functions and the functional (5) is rewritten expanding w ¼ ½1; u; vT in 2D space. So, the Euler– Lagrange equations associated with I and z are m a 2 I div ðz þ ke ÞrI þ wðw rIÞ ¼ I0
l
l
z
m 1 1 jrIj2 þ 2 þ bjrwj2 Dz ¼ 2 2ev 4e 4e
Functional (5) is strictly convex and even quadratic with respect to these two variables for each feasible and fixed motion field (w). Hence, a fast and effective numerical solution to the set of linear PDEs can be achieved through the finite element method [47]. Discretizing the spatio-temporal domain T O by a regular hexahedral grid leads to the system of equations associated with the Euler–Lagrange equations: ðAðw; zÞ þ CÞI ¼ CI 0
ð6Þ
ðCðI; wÞ þDÞZ ¼ E
ð7Þ
The finite element method is also employed for the phase-field approximation of the Mumford–Shah functional in [48,49] which is an extension to the approach first presented by Chambolle and Dal Maso [50]. Finally, gradient descent of EðI; z; oÞ with respect to the motion field variables ðu and vÞ is given by
b du ¼ ðw rIÞ@x1 I þ divððz2 þke ÞruÞ
ð8Þ
b dv ¼ ðw rIÞ@x2 I þ divððz2 þ ke ÞrvÞ
ð9Þ
a
a
where @xj I ¼ @I=@xj . For a fixed e the minimization strategy can be summarized as Step 0: Initialize I ¼ I0 ; z 1; u 0 and v 0. Step 1: Solve Eq. (6) for fixed u; v and z. Step 2: Solve Eq. (7) for fixed u, v, I. Step 3: Compute descent of E with respect to u from Eq. (8) when z; v; I are fixed. Step 4: Compute descent direction of E with respect to v from Eq. (8) when z; u; I are fixed. Step 5: Return to step 1. In contrast to the level set approach, the initialization of the phasefield is simpler as there is no need for a precise initialization close to the total edge set. Besides, the phase-field approximation comes along with a natural scale parameter. The smoothness of the energy depends on the scale parameter e and larger values of e yield coarse and smooth approximations of the images, the phase field and the deformation. In a multiscale strategy, the whole minimization process (steps 0–5) will be computed starting from a large value of e. This results in a coarse approximation as the stationary point in the simplified energy landscape. The process will be followed iteratively, reducing the scale parameter by taking the solution of the previous scale as the new initial guess on the next finer scale. Finally, e has to be of the same order as the grid size for the finest scale. It should be noted that due to the special definition of the functional, the motion field regularity is applied to the smooth part of image. Consequently, it is critical to
ARTICLE IN PRESS 24
A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
finalize the multiscale approach with small e in order to reach a regular motion field.
4. Experimental results The proposed joint segmentation and motion estimation model is implemented by using an Intels Core2TM Duo processor (T7300) at 2.0 GHz with a MATLAB platform and a Windows VistaTM operating system. The implemented solution is based on just one scale parameter instead of the multiscale strategy mentioned before, due to the small size of images in the clinical data. In the small frames of clinical Cine MRI, the anatomical objects are so close that they may be merged by choosing large scale parameters. Hence, in our application it is impractical to exploit a multiscale approach starting from large scales. Defining h as the grid width, experimentally preferred values of the parameters are b ¼ l ¼ 1, a ¼ 105 h2 , m ¼ h2 , and e ¼ h=4 [7]. The grid width of the FEM in this case can be computed from h ¼ minfhx ; hy ; ht g where hx ¼ 1=ðNx 1Þ, hy ¼ 1=ðNy 1Þ and ht ¼ 1=ðNt 1Þ are the inverses of number of pixels in each direction and the inverse of frames number. Ahead of the validation tests, it is necessary to discuss implementation considerations such as memory usage and convergence speed. Basically, the proposed joint method of this paper requires plenty of memory because the motion field, the phase field and all finite element matrices should be available simultaneously. It is strongly recommended to exploit sparse matrices and their operations for the finite element solution. But it still needs much more memory in comparison with tracking methods to contain the motion and phase-field variables. The most effective scheme is to limit the number of frames. It should be considered that for each pixel there are three variables namely, one phase-field and two motion field variables. So, omitting the finite element matrices, by increasing the number of frames by one, the number of unknown variables grows by three times of image pixel number. Fig. 1 depicts the convergence rate of the method from the initial condition of I= I0, z 1; u 0 and v 0. The evolution of
different terms in the functional is depicted in Fig. 1 for a normal Cine MRI of five frames. While the total energy and the optical R flow term ð a2jw rIj2 dwÞ are decreasing, the motion field Rb 2 regularity ð 2ðz þke Þjrwj2 dwÞ increases during evolution of motion field from its fully regular initial state (u 0 and v 0) to its minimizer. The decay of the total energy implies convergence speed of proposed method. In order to check the performance of the joint method, it is applied to a synthesized image sequence in which a disk with smooth gradient moves north west with a motion vector of m ¼ ½1; 1T . The sequence is degraded by inhomogeneous light to roughly simulate inappropriate imaging. The light decreases from bottom-right to upper-left where the object boundary completely fades. Fig. 2 shows the extracted edge set of the proposed joint process in comparison with a one task segmentation process [51]. In theory, all image based processes fail to extract the whole disk because of a lack of the intensity gradient in dark sites of image. Alternatively, the joint process can extract the disk completely relying on motion discontinuities in the absence of image discontinuities (Fig. 2(b)). As a consequence of considering both spatial and temporal discontinuities simultaneously, the extracted edge set successfully copes with the light condition and results in a more reliable edge set especially where the image gradient degrades. In the optical flow literature, the motion field regularity term is cited as desirable and it is shown to even reciprocate flickering [7]. But in the case of cardiac MR with rapid heart deformation and low frame rate, this term results in fictitious motion vectors especially during the systole phase. Instant contraction of the left ventricle at the beginning of the systole (much like all the other phase transitions) can be missed due to the motion smoothness term. On the other hand, the smoothness term is necessary in order to cope with the functional’s ill-posedness and noise which has destructive effects in the optical flow approach. In the optical flow equation (Eq. (3)), the variation of the pixels’ intensity is interpreted as their motion. By omitting the motion smoothness term, the estimated motion of adjacent pixels can be completely disparate. In order to inspect the latency and ability of the proposed method in tracking erratic motions, a simulated
Fig. 1. Evolution of different energy terms for a normal Cine MRI with five frames. (a) Optical flow term Total energy (Eq. (5)).
Ra Rb 2 2 2 2 jw; r Ij dw, (b) Motion field regularity 2 ðz ; k A Þjrwj dw, (c)
ARTICLE IN PRESS A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
25
Fig. 2. Extracted edge set from a synthetic image sequence with light inhomogeneity: (a) moving disk with light degeneration, (b) color coded motion field extracted by joint process, (c) segmentation result of [51] and (d) the phase-field function (z) of joint process.
Fig. 3. The effect of motion field regularity in tracking of rapid and sudden motion variations. Solid: ideal, dashed: b = 1, dotted: b =0.6.
pulse-like motion is synthesized in an images sequence of 10 frames. After three frames a textured object begins to move north east with m ¼ ½2; 2T and keeps its motion until the end of the sequence. Fig. 3 shows the latency caused by the motion smoothness term with b ¼ 1. While the solid line in Fig. 3 shows the norm of the real motion, the other lines show the object’s estimated motion norm with different b. Reducing b increases the proposed method’s ability in tracking of objects with rapid and sudden movements. However, there is a trade off as very small values of b reduce the robustness and can cause undesirable and irregular motion fields. For validation, the proposed method is applied to a set of clinical cardiac Cine MR image sequences. The dataset consists of
17 cases (seven females and 10 males) and a total of five healthy hearts. Each image sequence has 20 frames of 125 190 images from one heart cycle in long axes view. Our parameter setting for the clinical data was l = 1, b = 0, 5, a =105 h 2, m = 10 1 h 2 and e = h/10. The values of b and e were reduced to extract a fine edge set and to follow rapid movements. Processing each image sequence (size of 125 190 20) took 160s in average (with a maximum of 173s) using the implementation specifications stated before. Fig. 4 shows a color coded motion field extracted by the joint method associated with a normal heart. In order to provide benchmark, a specialist in interpreting the cardiac MRI delineated all the cardiac boundaries in all frames. Fig. 5 provides a comparison between manual and automatic
ARTICLE IN PRESS 26
A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
contours at four different events: end of diastole (ED), beginning of systole (BS), mid-systole (MS) and end of systole (ES). While the phase-field function is depicted by gray, white regions represent its difference from manual boundaries. Exploiting contrast agent, blood flow can be seen as well as wall motion in Cine MRI and can be extracted by dense motion estimation strategy in Eq. (5). Unfortunately there is no known reference for quantitative test of estimated motion by using Cine MRI, solely. But, it is possible to inspect the accuracy of the endocardium motion estimation by computing the difference
between the acquired edge set and the manual border delineated by an expert. Since the atrio-ventricular valves are not discriminative in Cine MR images, the heart is partitioned just in left and right sides, and Tanimoto’s index is used instead of ejection fraction [19]. Assuming that Oe and Om are extracted and manually delineated regions, respectively, Tanimoto’s index is defined as T1 ¼
JOe \ Om J JOe [ Om J
This true positive criterion determines the common area of automatically and manually extracted objects and implies both segmentation and motion estimation accuracies [52]. Tables 1 and 2 include the mean and standard deviation of T1 at four different states of beating heart in comparison with two other tracking based approaches [23,25]. The acquired results of Table 1 Mean of TI at ED, ES and MS phases.
Joint method Cascade method of [23] Level set and tracking of [25]
ED
BS
MS
ES
95.77 91.13 92.52
90.56 90.21 91.02
96.36 92.52 93.43
96.87 92.65 94.78
Table 2 Standard deviation of TI at ED, ES and MS phases.
Fig. 4. (a) Extracted motion associating with a normal heart at mid-systole besides applied velocity color map.
Joint method Cascade method of [23] Level set and tracking of [25]
ED
BS
MS
ES
1.03 1.56 1.24
3.10 3.62 2.36
2.12 2.39 2.43
1.19 2.68 1.32
Fig. 5. Extracted edge set in comparison with manually delineated boundary at four events of heart cycle: (a) end of diastole, (b) beginning of systole, (c) mid-systole and (d) end of systole.
ARTICLE IN PRESS A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
the joint process are superior to the results of cascade processes as a consequence of considering more than one frame for segmentation. But, no privilege is gained at the beginning of systole when a capricious contraction happens. This has been previously investigated as the effect of motion regularity term that may degrade the results in two or three frames, utmost (cf. Fig. 3). Although it is recommended to reduce the effect of the motion regularity term when the frame rate is low, it should not be suppressed to reach in appropriate smooth motion fields.
Ci;j ¼
Di;j ¼
Ei ¼
Z
27
m 1 jrIj2 þ 2 þ bjrwj2 ci ci dw 2ev 4e
Z
1 4e2
rci rcj dw Z
ci dw
It should be mentioned that w is spatio-temporal variable.
5. Conclusion
References
In this paper, a mathematical method for the simultaneous segmentation and motion estimation of cardiac MR image sequences is proposed. The method is based on the calculus of variation and minimizes a functional defined by using a combination of an optical flow functional and the Mumford–Shah segmentation functional. It is shown that finite element discretization and phase-field approximation provide an appropriate solution for the joint segmentation and motion estimation problem due to the fast convergence and the simple initialization. The edge set and the dense motion field associated with the minimizer of the functional are optimum and the proposed method outperforms conventional approaches. Consideration of spatial and temporal discontinuities simultaneously is the most important feature of proposed method. The tests performed on simulated image sequences imply the ability of the joint process in extracting even weak edges by using motion discontinuities. We have discussed the selection of the method’s parameters and especially the effect of the motion field regularity term on the method’s ability in tracking erratic movements. It can be inferred from the clinical experiments that the proposed method outperforms conventional cascade approaches without requiring any pre- or post-processing and even fine initialization. Since it is completely automatic and free from a priori knowledge, it can be employed for various heart shapes and functionality.
[1] D.J. Pennell, U.P. Sechtem, C.B. Higgins, W.J. Manning, G.M. Pohost, F.E. Rademakers, A.C. Van Rossum, L.J. Shaw, E.K. Yucel, Clinical indications for cardiovascular magnetic resonance (CMR): consensus panel report, Elsevier Eur. Heart J. 25 (2004) 1940–1965. [2] /http://www.scmr.org/S. [3] /http://atlas.scmr.org/cineplayer.htmlS. [4] /http://en.wikipedia.org/wiki/Cardiovascular_magnetic_resonanceS. [5] B. Lucas, T. Kanade, An iterative image restoration technique with an application to stereo vision, in: Proceedings of Darpa Uworkshop, 1981, pp. 121–130. [6] B. Horn, B. Schunk, Determining optical flow, Artif. Intell. 17 (1981) 185–204. [7] T. Preusser, M. Droske, C.S. Garbe, A. Telea, M. Rumpf, A phase field method for joint denoising, edge detection, and motion estimation in image sequence processing, SIAM J. Appl. Math. 68 (3) (2007) 599–618. [8] H.H. Nagel, W. Enkelmann, An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences, IEEE Trans. Pattern Anal. Mach. Intell. 8 (1986) 565–593. [9] A. Bruhn, J. Weickert, A confidence measure for variational optic flow methods, in: R. Klette, R. Kozera, L. Noakes, J. Weickert (Eds.), Geometric Properties from Incomplete Data, Springer, Dordrecht, The Netherlands, 2006, pp. 283–297. [10] G.E. Mailloux, F. Langlois, P.Y. Simard, M. Bertrand, Restoration of the velocity field of the heart from two-dimensional echocardiograms, IEEE Trans. Med. Imag. 8 (2) (1989) 143–153. [11] G. Zini, A. Sarti, C. Lamberti, Application of continuum theory and multi-grid methods to motion evaluation from 3D echocardiography, IEEE Trans. Ultrason. Ferroelect. Freq. Control 44 (3) (1997) 297–308. [12] J.G. Bosch, S.C. Mitchell, B.P.F. Lelieveldt, F. Nijland, O. Kamp, M. Sonka, J.H.C. Reiber, Automatic segmentation of echocardiographic sequences by active appearance motion models, IEEE Trans. Med. Imag. 21 (11) (2002) 1374–1383. [13] C. Pluempitiwiriyawej, J.M.F. Moura, Y. Lin Wu, C. Ho, Stacs, New active contour scheme for cardiac MR image segmentation, IEEE Trans. Med. Imag. 24 (5) (2005) 593–603. [14] Y. Wang, Y. Jia, Segmentation of the left ventricle from MR images via snake models incorporating shape similarities, IEEE Int. Conf. Image Process. 45 (2006) 213–216. [15] H. Van Assen, M. Danilouchkine, F. Behloul, H. Lamb, R. Van Der Geest, J. Reiber, B. Lelieveldt, Cardiac LV segmentation using a 3d active shape model driven by fuzzy inference, in: Proceedings of Medical Image Computing And Computer-Assisted Intervention (Miccai’03), Montreal, QC, Canada, November 16–18, 2003, pp. 533–540. [16] S. Mitchell, B. Lelieveldt, R. Van Der Geest, H. Bosch, J. Reiber, M. Sonka, Multistage hybrid active appearance model matching: segmentation of left and right ventricles in cardiac MR images, IEEE Trans. Med. Imag. 20 (5) (2001) 415–423. [17] M. Kaus, J. Von Berg, W. Niessen, V. Pekar, Automated segmentation of the left ventricle in cardiac MRI, Med. Image Anal. 8 (2004) 245–254. [18] J. Lotjonen, S. Kivisto, J. Koikkalainen, D. Smutek, K. Laurema, Statistical shape model of atria, ventricles and epicardium from short- and long-axis MR images, Med. Image Anal. 8 (2004) 371–386. [19] A. Pednekar, U. Kurkure, R. Muthupillai, S. Flamm, I.A. Kakadiaris, Automated left ventricular segmentation in cardiac MRI, IEEE Trans. Biomed. Eng. 53 (7) (2006) 1425–1428. [20] N. Paragios, Shape-based segmentation and tracking in cardiac image analysis, IEEE Trans. Med. Imag. 22 (6) (2003) 773–776. [21] A.S. Pednekar, I.A. Kakadiaris, U. Kurkure, R. Muthupillai, S.D. Flamm, Intensity and morphology-based energy minimization for the automatic segmentation of the myocardium, in: Proceedings of IEEE Workshop On Variational, Geometric and Level Set Methods in Computer Vision (Vlsm’03), Nice, France, October 2003, pp. 185–192. [22] D. Metaxas, L. Axel, A. Montillo, K. Park, Automated Segmentation, Motion Estimation of LV/RV Motion From MRI, in: EMBS/BMES Conference, 2002. Proceedings of the Second Joint, vol. 2, 23–26 October 2002, pp. 1099–1100. [23] J. Bosch, S. Mitchell, B. Lelieveldt, F. Nijland, O. Kamp, M. Sonka, J. Reiber, Automatic segmentation of echocardiographic sequences by active appearance motion models, IEEE Trans. Med. Imag. 21 (11) (2002) 1374–1383. [24] L. Legrand, C. Bordier, A. Lalande, P. Walker, F. Brunotte, C. Quantin, Magnetic resonance image segmentation and heart motion tracking with an active mesh based system, Comput. Cardiol. 22–25 (2002) 177–180.
Conflict of interest statement None declared.
Acknowledgments The authors would like to thank Mr. Alinaghizadeh and the Nour clinic for providing Cine MRI sequences and Dr. Mohebbi for delineation and tracking of heart chambers. The authors would also in particular like to thank Dr. Razvan for giving mathematical insight into variational theory.
Appendix A The system of equations from the finite element method contains Z I0 ci dw I 0 ¼ ½I0i ¼
C ¼ ½Ci;j ¼
Z
ci cj dw
Z m 2 A ¼ ½Aðw; zÞi;j ¼ ðz þke Þrci l Z a rcj dw þ ðrci wÞðrcj wÞ dw
l
ARTICLE IN PRESS 28
A. Eslami et al. / Computers in Biology and Medicine 40 (2010) 21–28
¨ [25] M. Suhling, M. Arigovindan, C. Jansen, P. Hunziker, M. Unser, Myocardial motion analysis from B-mode echocardiograms, IEEE Trans. Image Process. 14 (4) (2005) 525–536. [26] M. Droske, W. Ring, A Mumford–Shah level-set approach for geometric image registration, SIAM J. Appl. Math. 66 (2006) 2127–2148. [27] M. Droske, W. Ring, M. Rumpf, Mumford–Shah based registration: a comparison of a level set and a phase field approach, Comput. Vis. Sci. 12 (3) (2008) 101–114. [28] P. Kornprobst, R. Deriche, G. Aubert, Image sequence analysis via partial differential equations, J. Math. Imag. Vision 11 (1999) 5–26. [29] G. Aubert, P. Kornprobst, A mathematical study of the relaxed optical flow problem in the space Bv(O), SIAM J. Math. Anal. 30 (1999) 1282–1308. [30] G. Aubert, R. Deriche, P. Kornprobst, Computing optical flow via variational techniques, SIAM J. Appl. Math. 60 (1999) 156–182. [31] T. Brox, A. Bruhn, J. Weickert, Variational motion segmentation with level sets, in: H. Bischof, A. Leonardis, A. Pinz (Eds.), Computer Vision–Eccv 2006, Lecture Notes in Computer Science, vol. 3951, Springer, Berlin, 2006, pp. 471–483. [32] P. Nesi, Variational approach to optical flow estimation managing discontinuities, Image Vision Comput. 11 (1993) 419–439. [33] S.L. Keeling, W. Ring, Medical image registration and interpolation by optical flow with maximal rigidity, J. Math. Imag. Vision 23 (2005) 47–65. [34] N. Papenberg, A. Bruhn, T. Brox, S. Didas, J. Weickert, Highly accurate optic flow computation with theoretically justified warping, Int. J. Comput. Vision 67 (2006) 141–158. ¨ [35] B.V. Brunt, The Calculus of Variations, Birkhauser, 2004. [36] H. Sagan, Introduction to the Calculus of Variations, Courier Dover Publications, 1992. [37] L. Ambrosio, V.M. Tortorelli, On the approximation of free discontinuity problems, Boll. Un. Mat. Ital. B 7 (6) (1992) 105–123. [38] N. Vasconcelos, A. Lippman, Empirical Bayesian motion segmentation, IEEE Trans. Pattern Anal. Mach. Intell. 23 (2) (2001) 217–221. [39] S.C. Zhu, A. Yuille, Region competition: unifying snakes, region growing, and Bayes/MDLfor multiband image segmentation, IEEE Trans. Pattern Anal. Mach. Intell. 18 (9) (1996) 884–900.
[40] Y. Weiss, E.H. Adelson, A unified mixture framework for motion segmentation: incorporating spatial coherence and estimating the number of models, Proc. IEEE Conf. Comput. Vision Pattern Recognition 45 (1996) 321–326. [41] D. Cremers, F.R. Schmidt, F. Barthel, Shape priors in variational image segmentation: convexity, lipschitz continuity and globally optimal solutions, Proc. IEEE Conf. Comput. Vision Pattern Recognition 45 (2008) 1–6. [42] D. Mumford, J. Shah, Optimal approximation by piecewise smooth functions and associated variational problems, Comm. Pure Appl. Math. 42 (1989) 577–685. [43] L. Ambrosio, N. Fusco, D. Pallara, Functions of Bounded Variation and Free Discontinuity Problems, Oxford University Press, New York, 2000. [44] D. Cremers, S. Soatto, Motion competition: a variational approach to piecewise parametric motion segmentation, Int. J. Comput. Vision 62 (2005) 249–265. [45] D. Cremers, T. Kohlberger, C. SchnOrr, Nonlinear shape statistics in Mumford– Shah based segmentation, in: Proceedings of the Seventh European Conference on Computer Vision, Copenhagen, Lecture Notes in Computer Science, vol. 2351, Springer, Berlin, 2002, pp. 93–108. [46] G. Aubert, R. Deriche, P. Kornprobst, Computing optical flow via variational techniques, SIAM J. Appl. Math. 60 (1999) 156–182. [47] P. Solin, Partial Differential Equations and the Finite Element Method, WileyInterscience, 2005. [48] B. Bourdin, Image segmentation with a finite element method, M2AN Math. Model. Numer. Anal. 33 (1999) 229–244. [49] B. Bourdin, A. Chambolle, Implementation of an adaptive finite-element approximation of the Mumford–Shah functional, Numer. Math. 85 (2000) 609–646. [50] A. Chambolle, G. Dal Maso, Discrete approximation of the Mumford–Shah functional in dimension two, M2AN Math. Model. Numer. Anal. 33 (1999) 651–672. [51] M. Mora, C. Tauber, H. Batatia, Robust level set for heart cavities detection in ultrasound images, Comput. Cardiol. 45 (2005) 235–238. [52] S. Theodoridis, K. Koutroumbas, Pattern Recognition, Academic Press, USA, 1999 p. 366.