Introduction

Current density reconstructions (CDRs) of cardiac activation based on non-invasively obtained MCG (magnetocardiogram) and ECG (electrocardiogram) data are a new and promising tool and might support diagnoses in cardiology [6]. A typical result of a CDR is a color-coded activation map representing the magnitude of the current density in a volume or on a surface (e.g. the left ventricular surface, see Fig. 7) or a current dipole distribution in 3-D space. Commonly, these maps are interpreted by the cardiologist and represent the end point of analyses. However, for statistical data analyses (e.g. in group studies) a method is needed which enables a comparison of current density distributions for different time points and within groups of patients or volunteers. One way to achieve this goal is the proper parameterization of current density distributions and the application of statistical analyses to the parameters extracted. Previously, a parameterization based on the visual inspection of up to 8 sub-areas of the heart has been applied [7]. These sub-areas have been manually classified into active or not active ones and statistics have been computed about the number of classified sub-areas. One clear disadvantage of this method is that it yields only a very rough statistical description of the location and extent of the activation maximums and minimums. In this paper, we expand a new technique which has been introduced previously for 2-D planes [8] and full 3-D problems [9] to surface CDRs. This technique is based on the parameterization of current density distributions with the help of equivalent circles. The usefulness of our new technique is demonstrated through patient data.

Methods

A. Equivalent ellipsoid

An equivalent ellipsoid has been defined as a 3-D ellipsoidal object fitted to a current density distribution region in which the magnitude of the currents is above a certain threshold (in the following called supraliminal current density distribution). Reasons for the equivalent ellipsoid technique have been its geometrical simplicity (easy visualization) and a straightforward interpretation of lengths and directions of axes in the current density distributions. Fig. 1 shows an example of CDR results and the fitted ellipsoid. The equivalent ellipsoid is defined by its three orthogonal semi-axes ( a=a1_a, b=b1_b, c=c1_c), where denotes the unit vector connected to the axis . The fitting procedure has been realized in several steps. First, the threshold T_h used for marking the most important region in CDR has been defined by T_h= 100%*( Q_max+Q_mean)/Q_max, where Q_max is the maximum and Q_mean is the mean value of the dipole moments in the current distribution. Then, the center of gravity of the marked region (COG, P_COG) has been calculated. The position of the COG has been used as a new origin of the local coordinate system for the equivalent ellipsoid estimation. The direction of main axis 1_a of the equivalent ellipsoid has been computed on the basis of three different approaches (Fig. 1b). In the first approach, 1_a has been estimated through the normalization of the position of the dipole with the weighted longest distance from the COG (LD). The weighting of the distance has been done by the dipole moment in order to avoid that small and far off dipoles have a large influence on the definition of the main axis. In the second approach, the normalized position of the CDR maximum has been applied (PM). The dominant direction (DD) of the supraliminal region (defined by , where M_i denotes the moment of i-th current dipole and N_T is the number of dipoles in the a supraliminal region) has been employed in the third approach. Please note that DD is based on the sum of the current dipoles (vectors) and therefore represents a completely different approach than the other two approaches. The local coordinate system has been rotated so that the vector 1_a represents the z-axis of the new coordinate system. In the rotated coordinate system, the normalized position of the dipole with the maximal distance from P_COG on the new x-y plane has been used as the direction of the second semi-axis of the equivalent ellipsoid (1_b) (Fig.1c). The direction of the third axis has been determined by (Fig. 1d). The average length of all axes of the equivalent ellipsoid in the new orthogonal coordinate system ( 1_a, 1_b, 1_c) has been calculated iteratively according to

(1)

where .

The factor g is equal to 1 and iteratively increased by 1% until the inside ratio (I_R) is above 70%. The inside ratio  of the equivalent ellipsoid has been introduced in order to estimate how  many current dipoles are located inside the equivalent ellipsoid. Each dipole  has been associated with a voxel representing a small volume around the  dipole. We have defined a voxel as an elementary cube centered around a  current dipole with the edge length equal to the minimum inter distance  between the current dipoles in the distribution. The inside ratio I_R has been defined as the ratio of  volume of dipole voxels located inside the ellipsoid to the volume of all  voxels in the supraliminal distribution.

Figure 1. Construction of the equivalent ellipsoid: current density distribution given as a set of current dipoles ( a), supraliminal distribution and vector 1_a for differentapproaches ( b), vector 1_b for the longest distance approach ( LD) ( c), equivalent ellipsoid estimated on the basis of the longest distance approach ( LD) ( d).

In order to assess the quality of the equivalent ellipsoid we have introduced a goodness factor G_o, which has been defined as the volume ratio of the dipole voxels located inside the ellipsoid and the equivalent ellipsoid. The equivalent ellipsoid technique has been tested through the use of the two problems depicted in the following sections.

B. Equivalent circle (Dali's object)

The idea of equivalent circle has been influenced by the Salvador Dali's picture "The persistence of memory", which shows three clocks fitted to the curved surfaces (Fig. 2).

Figure 2. The persistence of memory, Salvador Dali (1936) (Salvador Dali Art Gallery, http://www.dali-gallery.com).

In order to apply the equivalent circle technique the reconstructed current density distribution has to be defined as a set of current dipoles on a surface. The information about the neighborhood of every current dipole or about the surface itself has to be given too. This is not very strong limitation because usually, the surface on which the CD is reconstructed has to be defined first. Standard output files from Curry ( .cdr files) contain information about the neighborhood. To realize the equivalent circle technique a following algorithm has been formulated:

• A supraliminal current density distribution is defined according to the description in section A.
• From the supraliminal current density distribution, a mesh of triangular elements is constructed similar. The outer line of this mesh is smoothed according to a suitable angle criterion. The constructed mesh can also be used for plotting the magnitude of the CDR in the region of interest (Fig. 3 and 4)

Figure 3. Supraliminal current density distribution (left) and triangular mesh with outline points (right)

.

Figure 4. Triangular mesh with smoothed outline (left) and the magnitude distribution of reconstructed CD over non-smoothed mesh (right).

• The center of gravity (COG) is calculated for the supraliminal mesh. Next, an average normal to the surface is computed as the average sum of normals found for every node in the supraliminal mesh. The COG is projected on the supraliminal surface using the average normal direction as a projection vector. The projected COG is used as the center of the equivalent circle.
• The set of 32 planes crossing the projected and original COG is set up with the constant angle between the planes (11.25^o). These planes define meridian lines (traces) on the supraliminal (Fig. 5)

Figure 5. Meridian traces with the pole located in projected COG.

• An average concentration radius of the CDR is calculated on the basis of the average lengths of the traces and defines the radius of the equivalent circle.
• Using the meridian traces and starting from the projected COG, five equidistant curves lying evenly with a parallel of latitude are set up. In this process a geodesic measure has been applied producing a set of points shown in Fig.6 (left). If the distance from the pole is larger than the radius of the equivalent circle the position of point lying outside the distribution is determined on the basis of the last part of the meridian trace.

Figure 6. Set of equidistant points lying evenly with a parallel of latitude (left) and equivalent circle ( Dali's object) (right).

• The set of equidistant points is used to create an equivalent circle fitted to the 3-D surface on which CD is reconstructed. The constructed equivalent circle (in a form of triangular mesh) is now named as Dali's object.

Several parameters accompanying a Dali's object can be used in statistical comparisons e.g. the position of the projected COG, the average concentration radius, the area of the supraliminal region, maximum and minimum length of the meridian traces determined on the basis of geodesic or Euclidian measure.

C. Patient data

The measurements have been taken in a magnetically shielded room (AK3b, Vacuumschmelze, Hanau, Germany) at the Biomagnetic Center in Jena, Germany. The magnetic field has been recorded with a twin dewar biomagnetometer system (2x31 channels) with first order axial gradiometers (Philips, Hamburg, Germany) [1]. Using system described above, we have measured the magnetocardiagram of a patient aged 70 years who had non- sustained ventricular tachycardia that developed after anterior left ventricular myocardial infarction and apical aneurysm at the Biomagnetic Center Jena. The subject has been lying in a supine position and the two dewars have been positioned above the thorax so that they covered the field maximums. We have recorded 600 s of signals at a sampling rate of 1000 Hz. To stabilize the baseline, an analog high pass filter (first order, 0.036 Hz) has been applied to the analog signals. We have performed a two step averaging procedure as described in [5] in order to improve the signal-to-noise ratio. A noise level of 50 fT has been estimated in both dewars. The last 40 ms of the bi-directional 30 Hz highpass filtered depolarization signal (late potentials, LP) have been used for inverse computations. A 3-D MRI data set of the chest of the patient has been obtained. A BEM model consisting of the left and right lungs as well as the outer torso surface has been applied for the magnetic field computations (forward model). We have used a triangular mesh with 2990 nodes and a linear potential approach for each triangle [3]. The surface of the lungs has been eroded by 3 mm in order to avoid numerical problems arising from a too short distance between the currents on the left ventricle (lv) and the surface of the lungs [4]. The conductivity ratio of torso to lungs was 5 to 1. The surface of the lv has been segmented from the MRI data set and subsequently used for the restriction of the source space. Two source configurations have been investigated. The first configuration has consisted of dipoles distributed on the surface of the left ventricle (1022 dipoles with an average dipole spacing of 4.7 mm). In the second one, a regular grid of dipolar sources (5 mm grid spacing in x, y, and z) has been placed into the lv resulting in 1715 current dipoles. The source parameters (dipole strength and orientation) have been determined through a minimum norm least squares algorithm (L2 norm) [2] for all time points. The equivalent ellipsoids have been estimated on the basis of the distribution containing the dipoles with the maximum strength over the last 40 ms time interval (LP).

Results

Figure 7 shows a diaphragmal view on the segmented 3-D surfaces and the maximum magnitude of the reconstructed dipole distribution in the late potential interval (LP image). The area with high activation corresponds to the infarcted area [6]. Figure 8 shows the equivalent ellipsoids found for the CDRs (minimum L₂ norm) on the surface of the lv and on the regular 3-D grid of points located inside the left ventricle.

Figure 7. Diaphragmal view on the 3-D magnetic resonance torso image (a), magnetic late potentials (LP) heart signals and analyzed time interval (gray) (b), zoomed diaphragmal surface of the left ventricle: LP image (c). The scale indicates the maximum dipole magnitude at every location during the LP interval.

 

Figure 8. Equivalent ellipsoids fitted to the reconstructed current dipole distributions located on the surface of the left ventricle (lower row) and on the regular 3-D grid points inside the left ventricle (upperrow) using different approaches. Diaphragmal view (left column), apical view (middle column), and magnified apical view for DD approach (right column).

The PM approach has yielded larger ellipsoids for the surface-based reconstructions than the other two approaches (Fig. 8). The results in Table I demonstrate that LD and DD yield similar volumes and semi-axes lengths, while PM clearly differs. The PM approach also exhibits the lowest goodness factor for the surface-based reconstruction in Table I, while the DD approach has the best goodness factor. The difference between the surface and 3-D grid- based COG for the supraliminal region has been equal to 4.3 mm.

TABLE I. Equivalent ellipsoids parameters for maximum LB

For the surface CDR the Dali's object has been constructed (Fig. 9). The area of supraliminal region found for the threshold T_h= 69.7% is equal to 1601.1 mm². The length of the minimum and maximum meridian traces is equal 17.0 mm and 32.0 mm respectively. The comparison with the minimum and maximum distances of the outline points to the projected COG (which are equal 15.9 mm and 29.6 mm) gives an imagination about the relatively small local curvature of the supraliminal region. The average convergence radius around the projected COG located in [26.9, -26.2, -24.6] mm has been found as 23.0 mm. The area of Dali's object calculated as an area of the circle using the convergence radius is equal 1655.5 mm² and based on a real object 1524.1 mm². The percentage ratio of the real Dali's object area to the supraliminal region area can be treated as the fitting quality factor in the equivalent circle technique (95.1% in the presented case).

Figure 9. Dali's object fitted to the reconstructed current dipole distribution located on the surface of the left ventricle (d), meridian "legs" of Dali's object (c), supraliminal region found for T_h=69,7% (a) and current density magnitude distribution (b).

Conclusions

In this paper we have introduced new techniques which enable the postprocessing of 3D CDRs. The techniques have been found to be useful for the analysis of inverse cardiac problems. Three different approaches of the algorithm have been tested. For 3-D grid-based CDRs all approaches yielded similar results. However, for a surface-based CDR the DD approach has been found to be the most appropriate one. One major advantage of the equivalent ellipsoid/circle techniques proposed is that they extract parameters from CDRs which can be used for statistical analyses (volume of the equivalent ellipsoid, length and direction of the semi-axes, COG, convergence radius, area of supraliminal surface). This allows both comparisons of different algorithms and comparisons of groups of patients or volunteers. Furthermore, the both procedures provide an easy and straightforward visualization of the foci of even very complex 3-D current density distributions. There are some limitations of the presented methods. The algorithms depicted have been applied only to CDRs with a single focus. In order to deal with multi-focal distributions an iterative strategy should be used. One possible strategy is to find local supraliminal regions and to employ the algorithms separately to each focus.

References