Thermodilution Catheter

10

3D/2D Registration of Mapping Catheter Images for Arrhythmia Interventional Assistance

Abstract

Radiofrequency (RF) catheter ablation has transformed treatment for tachyarrhythmias and has become first-line therapy for some tachycardias. The precise localization of the arrhythmogenic site and the positioning of the RF catheter over that site are problematic: they can impair the efficiency of the procedure and are time consuming (several hours). Electroanatomic mapping technologies are available that enable the display of the cardiac chambers and the relative position of ablation lesions. However, these are expensive and use custom-made catheters. The proposed methodology makes use of standard catheters and inexpensive technology in order to create a 3D volume of the heart chamber affected by the arrhythmia. Further,

We Propose A Novel Method That Uses A Priori 3D

information of the mapping catheter in order to estimate the 3D locations of multiple electrodes across single view C-arm images. The monoplane algorithm is tested for feasibility on computer simulations and initial canine data.

Key words: 3D reconstruction, monoplane imaging, navigation system, catheter ablation, cardiac arrhythmias.

1. Introduction

Severe disorders of the heart rhythm that can cause sudden cardiac death or morbidity, can be treated by radio- frequency (RF) catheter ablation, which consists of inserting a catheter inside the heart, near the area from which originates the abnormal cardiac electrical activity, and delivering RF currents through the catheter tip so as to ablate this arrhythmogenic area. The precise localization of the arrhythmogenic site and positioning of the RF catheter at that site are problematic: they can impair the efficacy of the procedure and the procedure can last many hours, especially for complex arrhythmias. To shorten the duration of RF catheter ablation and increase its efficiency, commercial systems that provide a 3D color display of the

The

arrhythmia have been proposed.

These Systems Incorporate Basket Electrode Arrays

(Constellation, EPT Inc.), catheters with a balloon electrode array (Ensite 3000, Endocardial Solutions Inc.)

(Cartotm, Biosense Webster Inc.). A Complete

navigation and registration framework is also available (CartoMerge, Biosense Webster Inc.). All these systems including purchase of system-specific catheters are costly for hospitals. The first two technologies can map the cardiac activation sequence using data recorded during a single beat whereas the CARTOTM system relies on data recorded point-by-point during numerous beats, which implies that the arrhythmia must remain stable during the procedure.

Other approaches have also been proposed to guide RF ablation therapy, such as the visualization of an optically tracked catheter by making use of magnetic resonance

Catheter , The Combination Of Ultrasound And Pre-

operative computer tomography (CT) , or preoperative

Important

electrophysiological data which allows the interventionist to determine the origin of the arrhythmia. Inverse electrocardiography, an established formulation is the imaging of the activation time map on the entire surface of

The Heart From Ecg Mapping Data. A Few Examples

computed from body surface potential measurements have been incorporated to segmented MRI images by Tilg et al.

And To Biplane Reconstructions Of The Cardiac

geometry by Ghanem et al. . DeBuck et al. have constructed a patient-specific 3D anatomical model from MRI and merged it with fluoroscopic images in an augmented reality environment that enables the transfer of cardiac activation times onto the model. Cristoforetti et al.

have developed a strategy for the spatial registration

Geometrical

reconstruction of the left atria and pulmonary veins retrieved from CT images. Methods based on real-time integration of electroanatomic and tomographic modalities are being developed and offer an intermodal fusion based on semi-automatic registration procedures initialized by

Fiducial Point Pairing . The Above Methods Make

use of preoperative data which do not reflect the real-time movement of the heart at the time of intervention, and introduce a logistical problem as these technologies are not present in the electrophysiological intervention room.

By

obtaining local activation times from a roving catheter whose positions are computed point-by-point from biplane fluoroscopic projections. We also proposed a novel concept by superimposing the isochronal map depicting the cardiac electrical activation sequence directly over the 2D fluoroscopic image of the heart , a method that the

Technology

lacks.

Biplane

fluoroscopy systems are rare in the clinical setting which in turn emphasizes the importance of developing single- view 3D reconstruction algorithms.

The Focus Of This Paper Is A Continuation Of Our

previous work, but with two significant modifications. To begin with, we consider implementing a full perspective camera model instead of the parallel projection concept so as to create a more precise 3D geometry of the heart chamber affected by the arrhythmia. Second, our aim is to answer the following question: is it possible to estimate precisely the 3D locations of the mapping catheter across all C-arm image frames in a cardiac cycle using only a single view? We propose to use a priori 3D information of the mapping ablation catheter positions in order to achieve this. If this is possible, then the interventionist could

Possibly Reduce Procedure Time Accordingly As A

monoplane C-arm fluoroscope is the primary imaging modality used in clinic and can provide real-time data images as well. To our knowledge this is the first work reported on estimating 3D locations using single view sequences in order to assist catheter ablation procedures.

Therefore, we tag this work as a feasibility study by presenting a practical implementation and experimental analysis using computer simulations and initial results on canine data.

2.1 Central Intuition

Clinical accuracy is determined by 3D reconstruction precision in the order of 2 mm or less and we don’t expect to achieve this using a monoplane sequence. As showed in , using only a single 2D C-arm image the recovered catheter depth had an error of about 10 mm. We believe that by incorporating temporal knowledge of the catheter coordinates may help improve depth estimation results.

Initial a priori 3D information of the mapping catheter can be obtained using pre-operative data from CT/MRI or via two-view reconstruction by rotating the monoplane C- arm fluoroscope at two perpendicular angulations and acquiring monoplane sequences. In this paper, we focus on the latter since the principal imaging modality used to

C-Arm

fluoroscope. Like most work in the cardiac field, the 3D a priori coordinates are reconstructed in the diastolic cardiac phase as heart motion is considered to be minimal at this instant. We then proceed to select one of the monoplane sequences acquired and estimate the depth of the mapping catheter at each image frame. The idea is to determine if or when the 3D reconstruction becomes inaccurate. In the end, what will be required is providing a 3D visual aid to the interventionist of the heart chamber and fusing the isochronal times on it in order to help them position

The

arrhythmogenic site.

2.2 C-Arm Fluoroscopy Calibration

Full Perspective Camera Model: Figure 1 shows the full perspective camera model that will be used for the 3D reconstruction problem. If we define a three dimensional point Pworld= [X Y Z 1]T in the world coordinate system, then its 2D projection in an image, m= [u v 1]T, is

(1)

The intrinsic matrix of size [3x3], contains the pixel coordinates of the image center, also known as the principal point (uo, vo), the scaling factor k, which defines

The Number Of Pixels Per Unit Distance In Image

coordinates, and the focal length f of the camera (in meters). The extrinsic matrix of size [3x4] is identified by the transformation needed to align the world coordinate system to the camera coordinate system. This means that a translation vector, t, and a rotation matrix, R, need to be found in order to align the corresponding axis of the two reference frames. Lastly, image resolution (usually mm/pixel) is calculated from the imaging intensifier size divided by the actual size of the image in pixels.

M

Fig 1. The perspective camera model. Any 3D world point can be projected onto a 2D plane and its coordinates would be (u, v)

Ijcsi

Ijcsi

pixels. The camera model is taken from the Epipolar Geometry Toolbox . Image resolution would be equal to C-arm intensifier size divided by image size in pixels.

Parallel & Weak Perspective Camera Models: An

orthographic camera is one that uses parallel projection to generate a two dimensional image of a three dimensional object. The image plane is perpendicular to the viewing direction. Parallel projections are less realistic than full perspective projections, however they have the advantage that parallel lines remain parallel in the projection, and distances are not distorted by perspective foreshortening.

(2)

The weak perspective camera is an approximation of the full perspective camera, with individual depth points Zi replaced by an average depth Zavg. We define the average depth, Zavg as being located at the centroid of the cloud of 3D points in the world coordinate system. The weak

2.3 Mapping Catheter Segmentation

We have recently developed a four-step filter in order to enhance coronary arteries visible in C-arm fluoroscopy images . Optimal filter parameters are presented there and omitted in this paper for brevity. We applied the same filter to enhance the catheters and extract electrode coordinates. Following is a brief description of these filters.

Homomorphic Filtering: A homomorphic filter is first used to denoise the fluoroscopic image. The illumination component of an image is generally characterized by slow spatial variation. The reflectance component of an image tends to vary abruptly. The homomorphic filter tends to decrease the contribution made by the low frequencies and amplify the contribution of high frequencies. The result is simultaneous dynamic range compression and contrast

(4)

with γL < 1 and γH >1. The coefficient c controls the sharpness of the slope at the transition between high and low frequencies, whereas Do is a constant that controls the shape of the filter and D(u,v) is the distance in pixels from the origin of the filter.

Perona-Malik Filtering: The Perona-Malik filter is implemented here in order to reduce and remove both noise and texture from the image, as well as, to preserve and enhance structures. The diffusion equation is given by

(5)

Where I is the input image and c(x, y, t), the diffusion coefficient, will control the degree of smoothing at each pixel point in the image. The diffusion coefficient is a monotonically decreasing function of the image gradient magnitude. It allows for locally adaptive diffusion strengths; edges are selectively smoothed or enhanced based on the evaluation of the diffusion function. Although any monotonically decreasing continuous function of the gradient would suffice as a diffusion function, two

(6)

Κ is referred to as the diffusion constant or the flow constant. The greatest flow is produced when the image gradient magnitude is close to the value of Κ. Therefore, by choosing Κ to correspond to gradient magnitudes produced by noise, the diffusion process can be used to reduce noise in images.

Complex Shock Filtering: The Complex Shock Filter

couples shock and linear diffusion in the discrete domain, showing that the process converges to a trivial constant steady state. To regularize the shock filter, the authors suggest adding a complex diffusion term and using the imaginary value as the controller for the direction of the flow instead of the second derivative. The complex shock

(7)

where a is a parameter that controls the sharpness of the slope, λ= reiθ is a complex scalar, λ̃ is a real scalar, ξ is the direction perpendicular to the gradient and η is the direction of the gradient. As fluoroscopy images have a low signal to noise ratio they tend to be noisy with artifacts present in them. Hence, we believe that applying a complex filter will result in a robust and stable deblurring process for the images as the filter is effective in very noisy environments.

Ijcsi

Ijcsi

Morphological Operation: Morphological filtering was applied as a final image processing step in order to eliminate background elements around the object of interest. The structuring element consists of a pattern specified as the coordinates of a number of discrete points relative to a defined origin. We chose a disk structuring element that has a radius of a few pixels, since the contours of the catheters and electrodes can be modeled as a disk.

Electrode Segmentation & Convex Hull: Using MatLab, the 4-step filter was implemented and applied on both monoplane C-arm datasets, beginning with the diastolic image. The smoothed images were then automatically

Thresholded Using Otsu’S Method. We Then Labeled

connected components in this image using the bwlabel function in MatLab. The centroid of the electrodes in each labeled region was then calculated using the regionprops function so as to obtain the coordinates required in the two

C-Arm Views. The Algorithm Outputs The Centroid

automatically and in case of failure, due to electrodes or catheters overlapping, we manually segment the image to obtain the desired coordinates. The convex hull algorithm is a classical and popular method used to reconstruct a 3D object from an unorganized set of points. The 2D version of this algorithm (which is easier to understand) is described as follows. Let S be a finite set of two- dimensional points in the plane. The convex hull of S is the smallest convex set that contains S. This means that the boundary of the convex hull is a convex polygon, whose vertices are points of S, and whose edges are line segments joining pairs of points of S.

2.4 A Priori 3D Monoplane Algorithm

We suppose that we have at our disposal a set of n 3D ablation catheter electrode coordinates (X0n, Y0n, Z0n) at time t = 0 obtained from two-view fluoroscopic data.

These a priori coordinates are expressed in the camera

Reference Frame In Order To Have A Z-Direction

corresponding to catheter electrode depth. Secondly, we

Have At Our Disposal The C-Arm Fluoroscope Gantry

parameters which can be extracted from the image header DICOM files. These parameters enable us to construct a projection matrix for a specific viewing angle. By selecting one of the two acquired monoplane datasets, we can solve for the 3D displacements (dxi,i+1, dyi,i+1, dzi,i+1) between consecutive C-arm image frames beginning with the first image i=1. Expanding equation (1) and using an additional C-arm image i =2, we obtain our first two

(8)

Both equations describe the pixel coordinates in the second image (u2, v2) and the twelve coefficients mi=1:12 are the values of the projection matrix. By adding an additional C-arm image frame i = 3 we obtain two new

(9)

These four equations take into account the spatial positions of a projected 3D world point on the acquired C-arm images. As we have four equations with six unknowns we can extract two additional equations based on the fact that the Euclidean distances in pixels, d, between catheter electrode points in two consecutive images are known. It is to note that the distance between 3D points is not the same as the distance between their projected image points. Thus, we consider orthogonal projection estimations in this case and we deem that this approximation is suitable enough for the proposed analysis. We arrive at the following two

(10)

We can now solve for the three dimensional displacements.

A Levenberg-Marquardt Optimization Scheme Can Be

used here in order to solve for the unknown displacements. For the optimization scheme, initial approximations are a requirement to initialize the process. Hence, a suitable approximation for the displacements dx and dy can be obtained if we consider a parallel back projection of the 2D image points into the world coordinate system. As for the displacements dz, if we assume that the average depth of the catheter electrodes remains relatively constant in consecutive time frames (i.e. weak perspective camera model), then we can calculate the average depth of the 3D

Points (X0, Y0, Z0). This Average Depth Should Be

relatively the same at future time instants t = 2, 3, etc., signifying that depths dz will be equal to zero for the optimization scheme. However, for the sake of a more exhaustive analysis, we also consider depth displacements dz ∈ millimeters as well. To justify the range of approximations, we note that the acquisition frame rate of

The Fluoroscopy Images Was Such That The Movement

between catheter electrodes in two adjacent images

Ijcsi

Ijcsi

resulted in displacements less than ten pixels. Finally, we note that the proposed equations do not incorporate nonrigid constraints implying that the monoplane 3D reconstruction of rigid moving objects (i.e. catheter electrodes) might be feasible.

2.5 Clinical Data Acquisition

A mongrel dog was anesthetized and laid on its right side on a fluoroscopy table (Integris Allura, Philips Inc.). A reference catheter and a pacing catheter were inserted into the right ventricle, close to the septal wall. The role of the reference catheter was to define an origin for our 3D coordinate system. This was important as motion artifacts are ever present during the experiments (heart beat, respiration, etc.), hence we deemed it appropriate to position it near a rigid landmark so that it experiences less movement due to artifacts. The role of the pacing catheter was to produce a simple electrical activation sequence so as to validate the isochronal maps. Finally, a standard RF ablation catheter was inserted from the femoral vein into the left ventricle (LV) of the dog. During the course of the

Experiment, This Mapping Catheter Was Moved To 20

different sites (i.e. point-by-point as the CARTOTM system) within the ventricle in order to obtain electrical and geometrical data from sufficient sites to map the activation sequence. The 20 landmarks were selected to reflect as closely as possible the entire volume of the ventricle.

The

CardioMap software system (Research Center at Sacré Coeur Hospital, Montreal, Canada). Specifications for the acquisition software were: 1000 samples/second, 0.05 Hz

High-Pass And 450 Hz Low-Pass Frequencies. Local

activation time was measured (ms) as the difference in the times of the fastest negative deflections (dV/dt) seen in the two electrograms recorded with the reference catheter and the mapping catheter. The fluoroscopic image acquisition rate was set to 60 fps so as to minimize motion artifacts as interframe 2D images are thus closer to each other and yield a smaller displacement between the objects in

Consecutive Images. Images Were Recorded During

approximately 2 seconds at the end of the expiration. The monoplane fluoroscopic C-arm was rotated by 90˚ to acquire two images for each mapping site: a left lateral view (the C-arm in a vertical position) and a posterior view (the C-arm in a horizontal position). Images were recorded with a 512 x 512 pixel resolution. We selected the diastolic frame for each mapping site so as to minimize motion blur by identifying the frame having the smallest root-mean-square difference with the preceding frame.

2.6 Evaluation

Computer Simulations: In order to validate out method, we first tested our proposed procedure on synthetic experimentation. We created a 3D helix containing 30 coordinate points so as to model the shape of a catheter.

Then we created two C-arm gantry setups that represented the posterior/anterior (PA) and left lateral (LAT) views of the heart. The focal length of the fluoroscopic X-ray system was equal to 1 meter and the helix location was set to 50 cm along the focal axis. The primary angles were equal to (90˚, 0˚) respectively for the PA and LAT views, whereas the secondary angles were equal to (0˚, 0˚) for both views respectively. The image sizes were set to [512x512] pixels and the intensifier size was chosen to be [178x178] millimeters. This allowed for a resolution 0.347mm/pixel. We could now calculate two projection matrices and we projected the 3D helix points on a first set of biplane images. These two-view images represent a first time instant at t=0. Subsequent biplane images are calculated by applying rigid movement on the 3D helix coordinates using the following rigid motion equation:

(11)

The 3D angles were set to [θ, Φ, φ] ∈ [0.5˚, 0.5˚, 0.5˚] and the 3D translations were set to [Tx, Ty, Tz] ∈ [1, 1, 1] millimeters each. This accounted for 2D interframe displacements in the images of (-6.82, -6.03) pixels in the

Left Lateral View And (4.97, 7.12) Pixels For The

posterior/anterior views. These values are similar to what was perceived using the 60 fps acquisition rate in the clinical images. A new set of biplane images and 3D helix points are obtained at t=1. In a similar fashion, equation (11) is reapplied to produce subsequent biplane data.

Excluding the a priori biplane set at t=0, a total of five biplane datasets are generated. For good measure we added errors of up to 2 mm to the coordinates.

Clinical Experiment: Regarding the clinical data, as we had a total of 20 biplane datasets, we calculated the 3D catheter points using a conventional biplane reconstruction algorithm (triangulation method). Details of the biplane algorithm are omitted here as it can be deduced that by using the two projection matrices for the X-ray gantry settings and the 2D electrode coordinates, we can easily estimate the 3D world points by solving a set of 4 projection equations with 3 unknowns (X, Y, Z) when considering a single point. A complete description is available in . For a specific dataset, we first extracted the diastolic images and used our 4-step filter to extract the electrodes. Then we performed two-view reconstruction.

This 3D a priori model was then used to initialize our monoplane equations. Since we are supposing that the C- arm fluoroscope gantry parameters remain fixed in a monoplane analysis, we can then calculate projection matrices that will allow the 3D points to be projected in the images temporally. We optimized the equations using a total of 12 visible electrodes including the actual mapping tip-electrode; the one that is in contact with the ventricle wall and directs heat to the arrhythmogenic site. With

Related Journal Articles & DOI Links

Selected peer-reviewed publications relevant to 12 Lead ECG Acquisition. Click the DOI to access the full paper (may require institutional access).

Why Choose Us?

Bangalore guidance for robotics, Spectre and autonomous systems projects.