Acta Polytechnica doi:10.14311/AP.2014.54.0214 Acta Polytechnica 54(3):214–220, 2014 © Czech Technical University in Prague, 2014 available online at http://ojs.cvut.cz/ojs/index.php/ap SEGMENTING TAGGED CARDIAC MRI DATA USING A LOCAL VARIANCE FILTER Tomáš Oberhubera,∗, Slavomír Kučeraa, Jakub Louckýa, Lucie Súkupováb, Jaroslav Tintěrab a Department of Mathematics, Faculty of Nuclear Sciences and Physical Engineering, Czech Technical University in Prague, Trojanova 13, Praha 2, 120 00, Czech Republic b Department of Radiology, Institute for clinical and experimental medicine, Vídeňská 1958/9, Praha 4, 140 21, Czech Republic ∗ corresponding author: tomas.oberhuber@fjfi.cvut.cz Abstract. We present a new method for segmenting tagged MRI data. It is based on the observation that the tagging grid vanishes almost immediately inside the ventricles because of the blood flow. Segmentation can be therefore driven by the intensity of the tagging grid. We estimate the intensity using a local variance filter. The main segmentation is computed by the graph-cuts method, though other segmentation methods can also be used. We compare the results obtained by our method with manual segmentation. Keywords: image segmentation, medical data processing, tagged MRI data, local variance filter, graph cuts. 1. Introduction Cardiac arrhythmia is very common problem in pa- tients. It has a wide spectrum of causes (ischemic disease or various form of cardiomyopathy) and leads to 1–2 % of mortality in developed countries. Exten- sive and expensive efforts are therefore dedicated to the treatment of this disease, e.g implantation of a biventricular cardioverter/defibrillator (ICD). Unfor- tunately, cardiac resynchronization therapy does not always lead to successful treatment. 30% of selected patients do not benefit from the therapy [1, 2] and 40% of them suffer progressive worsening of their heart function after this treatment [2]. Better classification and arrhythmia specification, together with improved navigation of the ICD im- plantation, could significantly help to increase the number of successfully treated patients. A method that has the potential to support the treatment with relevant and important parameters is magnetic reso- nance imaging (MRI) with spatial modulation of the magnetization (SPAMM). In this case, the detection of tagged myocardium on kinetic MRI series during the heart cycle provides important information about the circumferential and radial strain, which is used for assessing of local myocardial motion and wall deforma- tion. Analysis of the tagged myocardium is therefore a method with unique capabilities to describe the level of myocardial dyssynchrony. Myocardial tagging is spatial-selective magnetiza- tion modulation of the myocardium prior to imaging. As a result, a grid moving with the myocardium is measured in the data. An example of tagged cardiac MRI is shown in Figure 1. Tagging is a helpful tool for analysing the deforma- Figure 1. Comparison of common (top) and tagged (bottom) MRI data. tion of myocardium and for detecting dyssynchronies. Especially movement in the direction tangential to the myocardium walls is very difficult to recognize without tagging. Several methods for analysing the myocardium motion have been proposed. These meth- 214 http://dx.doi.org/10.14311/AP.2014.54.0214 http://ojs.cvut.cz/ojs/index.php/ap vol. 54 no. 3/2014 Segmenting Tagged Cardiac MRI Data Using a Local Variance Filter (a) Original image (b) Filtering with n = 2 (c) Filtering with n = 10 (d) Filtering with n = 17 Figure 2. Results of the local variance filter. ods involve movement detection and for postprocess- ing it is usually necessary to segment the ventricles. Widely used segmentation methods based on edge detectors cannot be successful, because the tagging grid creates new spurious edges and it destroys the original ones. One should focus on the grid pattern and not on the intensity itself. Methods proposed so far have mainly been applied to modulation in one direction. In this case, parallel lines appear in the image instead of a rectangular grid. In [3], the authors present a method based on local appearance features like edges, ridges and tagging line breakpoints. Based on these features, the Adaboost algorithm de- tects edges in the tagged data and an active shape model is then used for the segmentation. Grayscale opening and closing as an image segmentation pre- processing phase is used in [4]. Another approach is to remove the tagging lines by the Gabor filter bank technique before the main segmentation as pre- sented in [5]. In [6], a shape learning method is pre- sented. In the new method, we employ the local variance filter (LVF) on the image intensity function. After this preprocessing, common segmentation methods like level-set [7, 8] or graph-cuts [9–11] can be used. We demonstrate our method in combination with the graph-cuts method since it is faster and it gives results comparable with the level-set method. 215 T.Oberhuber, S. Kučera, J. Loucký et al. Acta Polytechnica 2. Contribution of the method We propose a new method involving a local variance filter for segmenting the tagged MRI data. This filter serves as preprocessing for the graph-cuts method. We present results of the method on several time sequences, and also a comparison with manual seg- mentation. The main advantage of our method is that it does not involve any learning phase. Only four parameters need to be set up. One of the parameters can be estimated from the tagging grid period. The method is not very sensitive to the other parameters. 3. Detection of the ventricles from the tagged MRI data In image segmentation, we aim to find one or more regions in the image such that their borders are de- fined by the edges in the image, and the intensity function is homogeneous inside this region. This can- not be fulfilled with tagged MRI data, because the tagging grid itself creates spurious edges which do not define the borders of the ventricles. As can be seen in Figure 1, the tagging grid disappears very quickly inside the ventricles. Blood is tagged in the same way as the cardiac tissue. However, molecules of blood with different magnetic modulation mix together al- most immediately. This observation leads us to the following idea. By making measurements of the local intensity of the tagging grid, we could detect regions inside and outside the ventricles. The grid can be defined as a fast transition of the intensity function from white to dark. This means that there is significant variation of the intensity function. We estimate it by LVF [12]. Let the image consist of pixels (i, j) for i, j ∈ 0 . . . N . Let it be mapped to a domain [0, 1]2, and we assume that it is symmetrically and periodically extrapolated to R2 [13]. For each pixel p of the image with coordi- nates (i, j), we take its squared neighborhood Np Np ≡ { (i′, j′) | i′ ∈ i− n, . . . i+ n, j′ ∈ j − n . . . j + n } , (1) for some n > 0. LVF is defined as Jp = 1 |Np| ∑ q∈Np ( Iq −Np )2 for Np = 1 |Np| ∑ k∈Np Ik . where |Np| denotes the number of pixels in Np, and Iq is the image intensity in pixel q. Figure 2 shows results obtained with different parameter n. If n is small, LVF still recognizes particular meshes of the tagging grid. Spurious edges are still present and segmentation would hardly be possible. As n grows, the grid meshes blend in and we get continuous detection of parts with the tagging grid (white color) and without (dark color). A drawback of filtering with larger n is the loss of sharp edges. A segmentation algorithm must therefore be applied carefully. 4. Segmentation with the graph-cuts Image segmentation by graph-cuts was introduced in [9, 10]. The idea of this method is to build a graph G consisting of so-called non-terminal nodes represented by the pixels of the input image and two terminal nodes represented by a source S and a sink T . The edges of the graph G split into two sets – t- links connecting a non-terminal node with a terminal node and n-links connecting two neighboring non- terminal nodes (pixels). For each interior pixel (i, j), we consider four neighboring pixels (i + 1, j), (i − 1, j), (i, j + 1), (i, j − 1). In a similar way we define neighboring pixels for pixels on the boundaries of the image. To initiate the segmentation process we define seeds, i.e. pixels marked by the user inside the segmented object (which will be denoted as O) and pixels outside the segmented object (which will be denoted as P ). The method of graph-cuts is based on finding either the minimal cut or the maximal flow in a graph with weighted edges. The two problems are dual to each other. Imagine that the edges in the graph corre- spond to pipes and we pump water from the source node. The water passes through the graph and it can leave only at the sink. The maximal flow, that can be reached, depends on the capacities (weights) of the pipes (edges). At a certain point, some set of pipes gets saturated and we cannot pump more water through the graph. This defines the maximal flow. At the same time, the same set of edges is the minimal cut (it is likely that pipes with smaller capacity will be saturated first). With a proper definition of the edge weights in the graph, the minimal-cut/maximal- flow agrees with the segmented object. If an n-link lies near to an edge in the image, we set a smaller capacity (weight) to it, because it is likely a part of the boundary of the segmented object. The t-links coming from the source to the nodes, which are likely a part of the segmented object, get higher capacity. Those connecting nodes outside the segmented object get low capacity. And vice versa for t-links connected with the sink. A summary is presented in Table 1. In this table, functions B(∆I, d), Rs(I) and Rt(I) are defined as B(∆I, d) = exp ( −∆I 2 2σ2 n ) · 1 d , Rs(I) = − lnP (I|O), Rt(I) = − lnP (I|P), ∆I = |Ip − Iq| . (2) which we have adopted from [14]. lnP (I|O) and lnP (I|P ) are conditioned probabilities expressing 216 vol. 54 no. 3/2014 Segmenting Tagged Cardiac MRI Data Using a Local Variance Filter Figure 3. Example of an input non-tagged image with the user defined seeds (sets O and P ) in light blue and yellow color. Link type Edge Capacity n-link (p, q) for p, q ∈ P, d = ‖(p, q)‖ B(∆I, d) t-link (s, p) for p ∈ P \ {O ∪ P} λRs(Ip) for p ∈ O ∞ for p ∈ P 0 (p, t) for p ∈ P \ {O ∪ P} λRt(Ip) for p ∈ O 0 for p ∈ P ∞ Table 1. Summary of edge weights. Ip, Iq denotes: the intensities of pixels p and q respectively (after local variance filtering); d denotes the distance of the pixels. By symbol ∞ we mean any value large enough for the edge never to get saturated while finding the maximal flow. 0 10 20 30 40 0 255 h o (I ) I 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0 255 P (I |O ) I 0 10 20 30 40 50 60 70 0 255 h b (I ) I (a) (b) (c) 0 0.002 0.004 0.006 0.008 0.01 0.012 0 255 P (I |P ) I 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 σn 32 B (∆ I ,1 ) ∆I -40 -20 0 20 40 60 80 100 120 0 255 R (I ) I (d) (e) (f) Figure 4. (a) Histogram of object seeds (O); (b) Approximate probability distribution of object (σt = 10); (c-d) The same for seeds of background (P); Functions used for evaluating edge weights: (e) n-links (σn = 3, d = 1); (f) t-links. The data is obtained from segmentation of the data in Figure 3. 217 T.Oberhuber, S. Kučera, J. Loucký et al. Acta Polytechnica Figure 5. Comparison of the manual segmentation (dotted white lines) and the segmentation by the new method (white lines). The initial seeds are marked as red and blue. 218 vol. 54 no. 3/2014 Segmenting Tagged Cardiac MRI Data Using a Local Variance Filter whether the given pixel belongs to object O and back- ground P , respectively. These probabilities are deter- mined from histograms measured on pixels of object and background respectively. One can show (see [11]) that two t-links (s, p) and (p, t) may be combined into one link (s, t)p with a weight R(I) = λ(Rs(I)−Rt(I)). We demonstrate the computation of the weights on non-tagged data depicted in Figure 3. The histograms computed on the initial seeds are shown in Figure 4 (a) and (c). The histogram is smoothed by a convolution with the Gaussian kernel as follows: P (x) = ∑ i h(i) exp ( − (x−i)2 2σ2 t )√ 2πσ2 t , (3) where h stands for one of the histograms ho or hb (Figure 4 (b) and (d)). After the weights are computed, the minimal- cut/maximal-flow is found by an appropriate optimiza- tion algorithm. At least two approaches are possible. The first is the Ford-Fulkerson algorithm [11, 15, 16], and the second is the Preflow-push algorithm [16, 17]. An advantage of the Preflow-push algorithm is its ca- pability to run in parallel [18]. For our computational study, we chose the Preflow-push algorithm. 5. Results The only parameters that must be set in the new method are n in (1), λ in t-links, σn in (2) and σt in (3). The choice of n has already been discussed and our experiments affirm that the optimal value of n corresponds with the period of the tagging grid. In our data, the period was approximately 35 pixels and that is how we set n. The other parameters were set as λ = 10, σt = 2 and σn = 7 for all experiments. We have segmented data measured on 10 patients. An example is demonstrated in Figure 5. Here we present a comparison of segmentation obtained by our method and by manual segmentation. The images in the same column are for the same spatial slice, and the rows represent the time frames. For each space slice we show the 4th to 7th time frame. The method is able to segment all time frames from the first to the 15th. After this frame, the tagging grid almost disappears. The biggest difference between manual and semi-automatic segmentation is apparent in the first column. In the left bottom part of the segmentation, there is a papilar muscle which was not recognized by the new method. It could easily be fixed by additional red seeds and re-segmentation. In the other space slices (columns 2 and 3), the two segmentations agree well. Acknowledgements This work was partially supported by research project RVO 68407700 (Ministry of Education, Czech Republic), by projects of the Student Grant Agency of the Czech Tech- nical University in Prague No. SGS11/161/OHK4/3T/14 and No. SGS14/206/OHK4/3T/14, by grant project (Min- istry of Health, Czech Republic) for the development of research organization 00023001 (IKEM, Prague, Czech Republic) Institutional Support. References [1] K. Bilchick, V. Dimaano, K. Wu, et al. Cardiac magnetic resonance assessment of dyssynchrony and myocardial scar predicts function class improvement following cardiac resynchronization therapy. Journal of the American College of Cardiology: Cardiovascular imaging 1(5), 2008. doi:10.1016/j.jcmg.2008.04.013. [2] J. White, R. Yee, X. Yuan, et al. Delayed enhancement magnetic resonance imaging predicts response to cardiac resynchronization therapy in patients with intraventricular dyssynchrony. Journal of the American College of Cardiology 48(10), 2006. doi:10.1016/j.jacc.2006.07.046. [3] Z. Qian, D. N. Metaxas, L. Axel. A learning framework for the automatic and accurate segmentation of cardiac tagged MRI images. In Proc. CVIBA Workshop, in Conjunction with ICCV, LNCS 3765, pp. 93–102. 2005. doi:10.1007/11569541_11. [4] A. Montillo, D. Metaxas, L. Axel. Automated segmentation of the left and right ventricles in 4d cardiac spamm images. In Medical Imaging Computing and Computer-Assisted Intervention, pp. 620–633. 2002. doi:10.1007/3-540-45786-0_77. [5] Z. Qian, X. Huang, D. Metaxas, L. Axel. Robust segmentation of 4d cardiac MRI-tagged images via spatio-temporal propagation. In Proceedings of SPIE Medical Imaging: Physiology, Function, and Structure frim Medical Images, vol. 5746, pp. 580–591. 2005. doi:10.1117/12.595862. [6] X. Huang, Z. Li, D. N. Metaxas. Learning couple prior shape and appearance models for segmentation. In MICCAI, vol. 1, pp. 60–69. 2004. [7] M. Beneš, R. Máca. Application of a degenerate diffusion method in 3D medical image processing. In Proceedings of Algoritmy 2012, pp. 282–290. 2012. [8] K. Mikula, A. Sarti. Parametric and geometric deformable models: An application in biomaterials and medical imagery, vol. 2, chap. Parallel co-volume subjective surface method for 3D medical image segmentation, pp. 123–160. 2007. [9] Y. Boykov, V. Kolmogorov. An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision. Pattern Analysis and Machine Intelligence 26(9):1124–1137, 2004. doi:10.1109/TPAMI.2004.60. [10] Y. Boykov, G. Funka-Lea. Graph cuts and efficient N-D image segmentation. International Journal of Computer Vision 70(2):109–131, 2006. doi:10.1007/s11263-006-7934-5. [11] J. Loucký, T. Oberhuber. Graph cuts in segmentation of a left ventricle from mri data. In Y. S. Beneš M., Kimura M. (ed.), Proceedings of Czech-Japanese Seminar in Applied Mathematics 2010. 2010. [12] L. Şendur, S. I. W. Bivariate shrinkage with local variance estimation. IEEE signal processing letters 9(12):438–441, 2002. doi:10.1109/LSP.2002.806054. 219 http://dx.doi.org/10.1016/j.jcmg.2008.04.013 http://dx.doi.org/10.1016/j.jacc.2006.07.046 http://dx.doi.org/10.1007/11569541_11 http://dx.doi.org/10.1007/3-540-45786-0_77 http://dx.doi.org/10.1117/12.595862 http://dx.doi.org/10.1109/TPAMI.2004.60 http://dx.doi.org/10.1007/s11263-006-7934-5 http://dx.doi.org/10.1109/LSP.2002.806054 T.Oberhuber, S. Kučera, J. Loucký et al. Acta Polytechnica [13] G. Aubert, P. Kornprobst. Mathematical problems in image processing. Springer-Verlag, Springer-Verlag New York, Inc., 175 Fifth Avenue, New York, NY 10010, USA, 2002. [14] Y. Boykov, M. P. Jolly. Interactive graph cuts for optimal boundary & region segmentation of objects in N– D images. In Eighth IEEE International Conference on Computer Vision, 2001. ICCV 2001. Proceedings., vol. 1, pp. 105–112. 2001. doi:10.1109/ICCV.2001.937505. [15] L. R. Ford, D. R. Fulkerson. Maximal flow through a network. Canadian Journal of Mathematics 8:399–404, 1956. [16] R. Sedgewick. Algorithms in C++, Part 5 – Graph Algorithms. Addison-Wesley, 2002. [17] A. V. Goldberg, R. E. Tarjan. A new approach to the maximum-flow problem. J ACM 35(4):921–940, 1988. doi:10.1145/48014.61051. [18] V. Vineet, P. J. Narayanan. CUDA cuts: Fast graph cuts on the GPU. In Computer Vision and Pattern Recognition Workshops, 2008. CVPRW ’08. IEEE Computer Society Conference on, pp. 1–8. 2008. doi:10.1109/CVPRW.2008.4563095. 220 http://dx.doi.org/10.1109/ICCV.2001.937505 http://dx.doi.org/10.1145/48014.61051 http://dx.doi.org/10.1109/CVPRW.2008.4563095 Acta Polytechnica 54(3):214–220, 2014 1 Introduction 2 Contribution of the method 3 Detection of the ventricles from the tagged MRI data 4 Segmentation with the graph-cuts 5 Results Acknowledgements References