Cardiac motion estimation using multi-scale optic flow: A successful application of Mathematica to a life sciences problem

We demonstrate a successful application of the algebraic Mathematica software to a computational life sciences problem. Mathematica was used to analyse cardiac tissue function from magnetic resonance images,and to visualise the motion fields in a vector plot and other important parameters in colour coded raster plots. Disturbances in twisting motion of the cardiac left ventricle during contraction and relaxation have been identified as early signs in cardiac pathologies.Extraction of this type of motion from radiological images requires the use of MR tagging,an imaging protocol that creates stripes in the cardiac muscle by means of spatial modulation of magnetisation,and that move along with the tissue. Because of the deformation of the cardiac muscle,the stripes are stretched and compressed at certain rates throughout the image,with variations in both space and time.The spatial and temporal varying disturbances call for image analysis using a multi-scale paradigm ,i.e. an image is analysed locally at the scale that best suits the local patiotemporal structure. A scale-space representing structure at multiple scales can be constructed by convolution of the image with a bank of filters,in our case Gaussians,with varying standard deviation (s). Motion in image sequences can be estimated using the optic flow constraint equation (OFCE),an intensity conservation law that holds along the flow;in 2D, Lx U + Ly V + Lt = 0 (L is intensity,HU, VLT is the flow vector,and (x, y, t) is the spatiotemporal coordinate;subscripts denote derivation).A combination of both the multi-scale paradigm and the OFCE is used in our application. More specifically,the first order multi-scale polynomial expansion of the OFCE is applied to two independent observations (with perpendicular tag stripes) of the same cardiac motion field. We need two observations (in 2D) to have sufficient equations to solve for all unknowns.In the first order approximation this leads to a system of eight equations with as many unknowns (the system C s = c,with C the coefficient matrix,s the vector HU, VLT,and c the inhomogeneous term), which has to be solved in every pixel tuple of the image sequence. Because of the multi-scale approach, before final estimation of the cardiac motion,optimal-scale selection is required. Our approach uses the condition number of C to find the most suitable scale combination for the local motion estimation. This however requires that the condition number be evaluated,and thus that the system be built,as many times as there are scale combinations (in our case 675,viz. from three temporal and fifteen spatial scales). Finally,at every pixel,the flow is estimated at the selected combination of scales. The above indicates that estimating cardiac motion involves both analytically sophisticated as well as numerically demanding mathematics. For the development of our prototype software we use Mathematica,because of its interactive nature,and rapid development and testing abilities.Because the evaluation of a single slice through the heart takes up to approximately 8 hours of computing,we are currently looking into parallel computation possibilities