Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Dynamic Coronary Roadmapping via Catheter Tip Tracking in X-ray Fluoroscopy with Deep Learning Based Bayesian Filtering

Dynamic Coronary Roadmapping via Catheter Tip Tracking in X-ray Fluoroscopy with Deep Learning... Percutaneous coronary intervention (PCI) is typically performed with image guidance using X-ray an- giograms in which coronary arteries are opacified with X-ray opaque contrast agents. Interventional cardi- ologists typically navigate instruments using non-contrast-enhanced fluoroscopic images, since higher use of contrast agents increases the risk of kidney failure. When using fluoroscopic images, the interventional cardiologist needs to rely on a mental anatomical reconstruction. This paper reports on the development of a novel dynamic coronary roadmapping approach for improving visual feedback and reducing contrast use during PCI. The approach compensates cardiac and respiratory induced vessel motion by ECG alignment and catheter tip tracking in X-ray fluoroscopy, respectively. In particular, for accurate and robust tracking of the catheter tip, we proposed a new deep learning based Bayesian filtering method that integrates the detec- tion outcome of a convolutional neural network and the motion estimation between frames using a particle filtering framework. The proposed roadmapping and tracking approaches were validated on clinical X-ray images, achieving accurate performance on both catheter tip tracking and dynamic coronary roadmapping ex- periments. In addition, our approach runs in real-time on a computer with a single GPU and has the potential to be integrated into the clinical workflow of PCI procedures, providing cardiologists with visual guidance during interventions without the need of extra use of contrast agent. Keywords: dynamic coronary roadmapping, X-ray fluoroscopy, catheter tip tracking, deep learning, Bayesian filtering, particle filter. 1. Introduction A guiding catheter is firstly positioned into the os- tium of the coronary artery. Through the guiding 1.1. Clinical Background catheter, a balloon catheter carrying a stent is intro- Percutaneous coronary intervention (PCI) is a duced over a guidewire to the stenosed location. The minimally invasive procedure for treating patients balloon is then inflated and the stent is deployed to with coronary artery disease. During these proce- prevent the vessel from collapsing and restenosing. dures, medical instruments inserted through a guid- ing catheter are advanced to treat coronary stenoses. PCI is typically performed with image-guidance using X-ray angiography (XA). Coronary arteries are This manuscript has been accepted for publication in Med- visualized with X-ray opaque contrast agent. During ical Image Analysis. For the published journal article, please refer to DOI: https://doi.org/10.1016/j.media.2020.101634. the procedure, interventional cardiologists may re- yy c 2020. This manuscript version is made available un- peatedly inject contrast agent to visualize the vessels, der the CC-BY-NC-ND 4.0 license. For license details, please as the opacification of coronary arteries only lasts for check via http://creativecommons.org/licenses/by-nc-nd/4.0/. a short period. The amount of periprocedural con- Corresponding author Email address: huama.research@gmail.com (Hua Ma) trast use has been correlated to operator experience, Preprint submitted to Medical Image Analysis January 14, 2020 arXiv:2001.03801v1 [eess.IV] 11 Jan 2020 procedural complexity, renal function and imaging with the live fluoroscopic sequence by aligning the setup ([32]). Furthermore, the risk for contrast in- R waves of their corresponding ECG signals. This duced nephropathy has been associated to contrast system compensates the cardiac motion of vessels, volume ([41]). Manoeuvring guidewires and mate- yet does not correct the respiratory motion during in- rial, however, typically occurs without continuous terventions. Two later studies by [50] and [26] in- contrast injections. In these situations, the navigation troduced image-based respiratory motion compensa- of devices is guided with ”vessel-free” fluoroscopic tion methods for DCR. Their methods assumed an images. Cardiologists have to mentally reconstruct ane respiratory motion model in ECG-gated fluo- the position of vessels and stenosis based on previ- roscopic frames and recovered the respiratory mo- ous angiograms. tion from soft tissues with special handling of static structures. The limitation of these approaches is that 1.2. Dynamic Coronary Roadmapping they require relevant tissue to be suciently visi- ble in the field of view for reliable motion compen- Dynamic coronary roadmapping (DCR) is a sation which is not always the case. In addition, promising solution towards improving visual feed- they require to be run on cardiac-gated frames. In back and reducing usage of contrast medium during a more recent work by [20], binary vessel masks PCI ([12, 50, 26, 20]). This approach dynamically were created as the roadmaps from at least one car- superimposes images or models of coronary arter- diac cycle of angiographic images. Temporal align- ies onto live X-ray fluoroscopic sequences. The dy- ment of the roadmaps and the fluoroscopic sequence, namic overlay serves as a roadmap that provides im- which compensated the cardiac motion of vessels, mediate feedback to cardiologists during the inter- was performed by registering ECG signals using vention, so as to assist in navigating a guidewire into cross-correlation. Additionally, the respiratory mo- the appropriate coronary branch and proper place- tion was corrected by aligning the guidewire center- ment of a stent at the stenotic site with reduced ap- line in the fluoroscopy to the contour of vessels in the plication of contrast agent. Studies with a phantom angiogram where the roadmaps were created. The setup using research software ([20]) or on first cases system has been shown useful in a phantom-based of clinical interventions using commercially avail- study, nevertheless no accuracy evaluation of the spa- able systems ([9, 48, 40]) have investigated the use- tiotemporal alignment was presented. Furthermore, fulness of DCR in assisting PCI, reporting that DCR the spatial registration relies on robust extraction of helps to reduce procedure time, radiation dose and vessels and guidewires which is often challenging for contrast volume. X-ray images. To develop a DCR system, it is important but yet a challenge to accurately overlay a roadmap of coro- Unlike direct roadmapping, the model-based ap- nary arteries onto an X-ray fluoroscopic image, as proaches build a model to predict motion in fluoro- limited information of vessels is present in the target scopic frames. The motion models are often func- fluoroscopic image for inferring the compensation of tions that relate the motion of roadmaps to surro- the vessel motion resulting from patient respiration gate signals derived from images or ECG, so that and heartbeat. The methods that have been proposed once the surrogates for fluoroscopic frames are ob- on motion compensation for DCR can be generally tained, the motion can be computed by the model. grouped into two categories: direct roadmapping and For cardiac interventions including PCI, the organ model-based approaches. motion is mainly a ected by respiratory and car- Direct roadmapping methods use information diac motion. Many previous works often built a mo- from X-ray images and ECG signals to directly cor- tion model parameterized by a cardiac signal derived rect the motion caused by respiration and heartbeat. from ECG and a respiratory signal obtained from di- The first DCR system ([12]) used digital subtrac- aphragm tracking ([37, 42, 13]) or automatic PCA- tion of a contrast sequence and a mask sequence based surrogate ([15]). Some other works model to create a full cardiac cycle of coronary roadmaps. only the respiratory motion in cardiac-gated images The roadmaps were stored and later synchronized ([35, 21, 30]). For a complete review on respira- 2 tory motion modeling, we refer readers to the sur- fitting the B-spline to measurements based on gray vey article by [27]. One limitation of the model- intensity and vesselness image. [5] proposed a con- based approaches is that the motion models are often volutional neural network (CNN) to detect catheter patient-specific, which requires training the model electrodes in X-ray images, which treated catheter every time for a new subject. Additionally, once detection as a segmentation problem. The method the surrogate values during inference are out of the used a weighted cross-entropy loss to cope with the surrogate range for building the model, e.g. for ab- class imbalancing problem due to the small size of normal motion, extrapolation is needed, which may the target. [22] and [11] tracked surgical instruments hamper accurate motion compensation. using a deep network having an encoder-decoder ar- chitecture. Their approaches combined instrument 1.3. Interventional / Surgical Tool Tracking segmentation and detection in a multi-task learning problem to make the tool detection in a cluttered Tracking interventional tools is relevant for mo- background more robust. tion compensation ([35, 7, 25, 4, 2]). In particular for Di erent from tracking by detection, which relies PCI, the guiding catheter tip typically remains within solely on the current image, temporal tracking also the coronary ostium which is in the field of view uses information from previous frames. The tempo- during the largest part of the intervention, making ral information can reduce the search space for de- it a suitable landmark for tracking. [4] have shown tection, or put additional constraints in the model, that catheter tip motion during PCI can be modeled making the tracking more robust. as a combination of cardiac and respiratory motion. As using catheter tip displacement can only correct Temporal information has been used in various translational motion, [4] further showed that, com- ways. Some methods mainly relied on a detection pared to a rigid motion model for the respiratory mo- model, but incorporate temporal information in the tion, modeling only the translational part of the res- preprocessing ([7]) or post-processing ([16]) step or piratory motion deteriorated the accuracy marginally, in the input ([33, 1]). Approaches based on back- which confirms the observations by [36] that the rota- ground estimation have been used for catheter ([49]) tional part of respiratory motion is small. These find- or guidewire ([31]) tracking. In these approaches, the ings motivate motion compensation for DCR through background was recursively updated for every frame, tracking the catheter tip in X-ray fluoroscopic se- and was used for enhancing the foreground contain- quences. ing instruments. Apart from those, many works Many works have been proposed to address the adopted a Bayesian framework for tracking instru- problem of tracking interventional or surgical tools ments via a maximum a posteriori (MAP) formu- in medical images for various applications. The lation. Representations based on key points ([47]), tracking methods from these works can be generally B-splines ([45, 29, 19, 18]), or segment-like fea- categorized into two kinds of approaches: tracking tures ([43]) have been used to model catheters or by detection, and temporal tracking. guidewires. Markov random field (MRF) was used The tracking by detection approaches treat track- to model the position or deformation of the con- ing as a detection problem, which rely on features trol points in the B-spline ([29, 19, 18, 47]). In only from the current image without using informa- the work by [43], temporal information was incor- tion from previous frames. For example in electro- porated in the prior term using Kalman filter. Partic- physiology procedure, as the catheters present spe- ularly, learning-based approaches were used in sev- cific features in shape or intensity, ad hoc meth- eral works to obtain the likelihood for a more ro- ods were proposed based on, e.g. blob detection, bust measurement using probabilistic boosting tree shape constrained searching and model-/template- ([45, 46]) or support vector regression ([29]). In ad- based detection ([25, 24]). [8] modeled the catheter dition, temporal tracking models based on Bayesian tracking problem by optimizing the posterior in a filtering were also a popular approach for instrument Bayesian framework, in which the catheter was rep- tracking. [2] used a hidden Markov model (HMM) resented as a B-spline tube model and was tracked by to track catheter tip in a 3D vessel tree, for which 3 the likelihood was obtained based on the 3D-2D reg- 3. We evaluate the proposed approach visually istration outcome. [38] used particle filters to track and quantitatively on clinical X-ray sequences, surgical tools in medical images. They used a likeli- achieving low errors on both tracking and hood based on the segmentation of instruments, and roadmapping tasks. a dynamic model that incorporates samples from two 4. The proposed DCR method runs in real-time previous time steps. In a later work, [39] used a with a modern GPU, thus can potentially be multi-object particle filter to track multiple instru- used during PCI in real clinical settings. ment regions simultaneously, in which a particle is the concatenation of the states of several objects. 2. Scenario Setup and Method Overview Despite of many existing works on inverventional The proposed method assumes that the scenario of or surgical tool tracking in medical images, an au- performing dynamic coronary roadmapping to guide tomatic approach for tracking the tip of guiding a PCI procedure consists of an oine phase and an catheter in X-ray fluoroscopy for PCI has not been online phase. An overview of the method is shown investigated yet. The challenges of this task are: in Figure 1. (1) di erent from the catheters for EP that can be viewed as blobs or a circle, the guiding catheter for 2.1. Oine Phase PCI presents a dark tubular appearance which shows This phase (Step 0 in Figure 1) is performed o - no prominent features; (2) the shape of the guiding line before the actual roadmapping is conducted. In catheter tip segment varies depending on the orienta- this stage, roadmaps of coronary arteries contain- tion of the C-arm, making feature-/model- based de- ing multiple cardiac phases are created from an X- tection challenging; (3) the background may contain ray angiography sequence acquired with injection of structures that have similar appearance to a catheter contrast agent. A roadmap can be a vessel model tip, such as vertebral structures or residual contrast in the form of centerlines, contours, masks, etc. It agent, which makes robust tracking dicult. may also contain information of clinical interest, e.g. stenosis. Since the main focus of this paper is on 1.4. Contributions accurate overlay of a roadmap, we do not investigate We propose and evaluate a novel approach for dy- how to create the most suitable roadmaps, but use the namic coronary roadmapping. The approach com- images containing only vessels and catheters that are pensates changes in vessel shapes and cardiac motion created using the layer separation method by [23] as by selecting the roadmap of the same cardiac phase the roadmaps to show the concept of dynamic coro- through ECG alignment, and corrects the respiratory nary roadmapping. Along with the XA sequence, induced motion via tracking the tip of the guiding ECG signals are also acquired and stored for later catheter. Our contributions are: selecting a roadmap that has similar cardiac phase to 1. We develop a new way to perform dynamic a given X-ray fluoroscopic frame in the online phase coronary roadmapping on free breathing, non- (see details in Section 3). Once the image sequence cardiac-gated X-ray fluoroscopic sequences. and ECG signals are acquired, the catheter tip loca- Particularly, the respiratory-induced vessel mo- tion in every frame is obtained to serve as a refer- tion is robustly compensated via the displace- ence point for roadmap transformation. In this work ment of catheter tip. we manually annotated the catheter tip in the oine 2. We proposed a deep learning based method XA sequence. In real clinical scenarios, the annota- within a Bayesian filtering framework for online tion work can be done by the clinician or a person detection and tracking of guiding catheter tip in who assists the intervention, such as a technician or X-ray fluoroscopic images. The method models a nurse. the likelihood term of Bayesian filtering with a 2.2. Online Phase convolutional neural network, and integrates it with particle filtering in a comprehensive man- This is when the dynamic roadmapping is actu- ner, leading to more robust tracking. ally performed. In this phase, non-contrast X-ray 4 Offline Phase Online Phase A fluoroscopic frame ECG image Manual annotation Angiographic sequence Tips Tips Track tip images image Tip Tip image Roadmaps Roadmaps ECG Deep learning based Translation Bayesian filtering Correlation Layer separation score Tip Tip Transform ECG ECG Roadmap Roadmap Roadmap STEP 0 STEP 1: Select Roadmap STEP 2: Track Catheter Tip STEP 3: Transform Roadmap Section 2.1 Section 3 Section 4 Section 2.2 Figure 1: The overview of the proposed dynamic coronary roadmapping method. The colored blocks with a dash line border denote objects acquired in the online phase; the colored blocks with a solid line border are objects originated from the oine phase. fluoroscopic images with the same view angles as oroscopic image, which compensates the di erence the roadmaps created during the oine phase are ac- of vessel shape and pose induced by cardiac motoin. quired sequentially. At the same time, ECG signals An approach similar to the ECG matching method by along with the roadmapping frames are also obtained [20] is used to accomplish this task. and are compared with the stored ECG to select the To select roadmaps images based on ECG, a tem- most matched roadmap (Step 1 in Figure 1; see de- poral mapping between X-ray images and ECG sig- tails in Section 3). This is to compensate the change nal points needs to be built first. We assume that of vessel shape and position between frames due ECG signals and X-ray images are well synchronized to cardiac motion. Simultaneously, the catheter tip during acqusition. In the oine phase, the begin- location in the acquired X-ray fluoroscopic images ning and the end of the image sequence are aligned is tracked online using the proposed deep learning with the start and end ECG signal points; the XA based Bayesian filtering method in Section 4 (Step frames in between are then evenly distributed on the 2 in Figure 1). The displacement of catheter tip be- timeline of ECG. This way, a mapping between the tween the current image and the selected roadmap stored sequence images and its ECG signal can be image is then obtained and are applied to transform set up: for each image, the closest ECG signal point the roadmap. Finally, the transformed roadmap is to the location of the image on the timeline can be overlaid on the current non-contrast frame to guide found; for each ECG point, an image that is closest the procedure (Step 3 in Figure 1). to this point on the timeline can be similarly located. Once the mapping is available, all images with good vessel contrast filling and the ECG points that are as- 3. ECG matching for Roadmap Selection sociated to these images are selected from the XA Roadmap selection in this work is achieved by sequence for the pool of roadmaps. In this process, comparing the ECG signal associated with the flu- at least one heartbeat of frames should be acquired, oroscopic image and the ECG of the angiographic which is generally the case in our data. In the on- sequence, such that the most suitable candidate line phase, similar to the approach of [20], for ac- roadmap is selected where the best match of the ECG quisition of each image, a block of N latest ECG ECG signals is found. The selected roadmap has the same signal points is constantly stored and updated in the (or very similar) cardiac phase with the X-ray flu- history bu er. This is considered as the ECG signal Select roadmap Overlay corresponding to the fluoroscopic frame. the likelihood p(z jx ) from which x and z can be k k k k obtained by sampling. To compare the ECG signals associated with the angiographic sequence and the online fluoroscopic With these definitions and p(x ), the inital be- lief of x , Bayesian filtering seeks an estimation of image, a temporal registration of the two signals us- ing cross-correlation is applied ([20]). The two ECG x (k  1) based on the set of all available observa- tions z = fz ; i = 0; :::; kg up to time k via recur- signals are first cross-correlated for every possible 0:k i sively computing the posterior probability p(x jz ) position on the signals, resulting in a 1D vector of k 0:k as Eq.(1) ([3]): correlation scores. The candidate frame for dynamic overlay is then selected as the one associated with the point on the ECG of the angiographic sequence that p(x jz ) / p(z jx ) p(x jx )p(x jz )dx : k 0:k k k k k1 k1 0:k1 k1 is corresponding to the highest correlation score. | {z } p(x jz ) k 0:k1 (1) 4. Bayesian Filtering for Catheter Tip Tracking Assuming the initial probability p(x jz ) = p(x ) is 0 0 0 known, based on Eq.(1), Bayesian filtering runs in Bayesian filtering is a state-space approach aiming cycles of two steps: prediction and update. In the at estimating the true state of a system that changes prediction step, the prior probability p(x jz ), the k 0:k1 over time from a sequence of noisy measurement initial belief of x given previous observations, is made on the system ([3]). One popular application estimated by computing the integral in Eq.(1). In area of this approach is tracking objects in a series of the update step, the prior probability is corrected by images. the current likelihood p(z jx ) to obtain the posterior k k p(x jz ). k 0:k 4.1. Theory of Bayesian Filtering In Section 4.2, we will firstly introduce how to Bayesian filtering typically includes the following model the likelihood. Then in Section 4.3, a way components: hidden system states, a state transition to represent and eciently approximate the posterior model, observations and a observation model. Let will be discussed. Finally in Section 4.4, a summary x 2 R (k = f0; 1; 2; :::g) denote the state, the loca- of the complete catheter tip tracking method will be tion of guiding catheter tip in the k-th frame, a 2D given. vector representing the coordinates in the X-ray im- 4.2. A Deep Learning based Likelihood age space. The transition of the system from one state to the next state is given by the state transition Directly modeling the likelihood p(z jx ) is chal- k k model x = f (x ; v ), where v 2 R is an in- lenging due to (1) the complexity of the generation k k k1 k1 k1 dependent and identically distributed (i.i.d.) process process h and (2) the computational complexity of 2 2 2 2 noise, f : R  R ! R is a possibly nonlinear p(z jx ) for every value x 2 R . In this work, we k k k k function that maps the previous state x to the cur- simplify the problem by only computing the likeli- k1 rent state x with noise v . The observation z in hood p(z jx ) in the image pixel space, i.e. the in- k k1 k k k this work is defined as the k-th X-ray image of a se- teger pixel coordinate. For a subpixel x , the value wh quence, so that z 2 R , where w and h are the of p(z jx ) can possibly be approximated by inter- k k k width and height of an X-ray image. We further de- polation. To this end, we propose to use a deep fine the observation model as z = h (x ; n ), where neural network D to approximate p(z jx ) for inte- k k k k k k n 2 R is an i.i.d measurement noise (n is the di- ger pixel locations. The network takes an image z k k k 2 n wh mension of n ), h : R R ! R is a highly non- as input and outputs a probability of observing the k k linear function that generates the observation z from input z for every pixel location x . Therefore, the k k k the state x with noise n . The state transition model approximated likelihood is a function of x , denoted k k k f and the observation model h , respectively, can as D (x ). Since x is defined within the scope of k k z k k also be equivalently represented using probabilistic the image pixel space, D (x ) is essentially a proba- z k forms, i.e. the state transition prior p(x jx ) and bility map having the same dimension and size with k k1 6 [28], residual blocks ([17]) are adopted at each reso- lution level in the encoder and decoder to ease gradi- ent propagation in a deep network. The encoder con- sists of 4 down blocks in which a residual block fol- lowed by a stride-2 convolution is used for extraction and down-scaling of feature maps. The number of (a) (b) (c) feature maps is doubled in each downsampling step. The decoder has 4 up blocks where a transposed con- Figure 2: Input and ground truth labels for the deep neural net- volution of stride-2 is used for upsampling of the in- work: (a) an input X-ray fluoroscopic image, (b) the binary catheter mask of (a) for catheter segmentation, (c) a 2D Gaus- put feature maps. Dropout is used in the residual unit sian PDF ( = 4 px) for likelihood estimation for (a). of the up block for regularization of the network. Be- tween the encoder and the decoder, another residual block is used to process the feature maps extracted the input image z , in which the entry at each loca- by the encoder. The detailed network architecture is tion x ( j = 1; 2; : : : ; wh) in the map represents the shown in Figure 3. probability of observing z given x . It is worth men- Due to similar appearance between a guiding tioning that the deep neural network is used for ap- catheter tip and corners of a background structure, proximation of p(z jx ), which should be clearly dis- k k such as vertebral bones, lung tissue, stitches or tinguished from the generation model h that maps guidewires, ambiguity may exist when the network an x to z . The existence of h is merely for the con- k k k is expected to output only one blob in the probability venience of definition, its explicit form, however, is map. To alleviate the issue, we adopt a similar strat- not required in the context of this work. egy as [22], using a catheter mask (Figure 2b) as an To obtain the training labels, we assume that there additional label to jointly train the network to output exists a mapping h , such that the training label can both the catheter segmentation heatmap and the like- be defined as a distance-based probability map, i.e. lihood probability map. The segmentation heatmap the farther away x is from the ground truth tip loca- is obtained by applying a 1  1 convolution with tion in the image z , the less possible it is to observe ReLU activation on the feature maps of the last up z given x through the process h . This definition block. To compute the likelihood probability map, a k k k matches the intuition that from a location x that is residual block is firstly applied on the feature maps far from the ground truth tip location, the probability of the last up block. The output feature maps are of observing a z with the catheter tip being located then concatenated with the segmentation heatmap as at the ground truth position should be low. For sim- one additional channel, followed by a 1  1 convo- plicity, a 2D Gaussian probability density function lution. Finally, to ensure the network detection out- 0 2 (PDF) N (x ; x ;  I) centered at the ground truth tip put fits the definition of a probability map on image 0 2 location x with variance  I in the image space is locations, following the 1  1 convolution, a spatial used as the label to train the network (Figure 2c). softmax layer is computed as Eq.(2): Note that this training label makes the estimation of k;l p(z jx ) equivalent to a catheter tip detection problem k k D = P ; (2) k;l i; j such that the deep neural network learns features of i; j catheter tip and outputs high probability at locations where A is the output feature map of the 1  1 con- where the features are present. Due to this reason, volution, A denotes the value of A at location (i; j), i; j we also call p(z jx ) “detection output” or “detec- k k D is the final output of the detection network, a 2D tion probability” and call the estimation of p(z jx ) k k probability map representing p(z jx ). The details k k “catheter tip detection” in the context of this paper. are shown in Figure 3. The network that we use follows a encoder- The training loss is defined as a combination of the decoder architecture with skip connections similar to segmentation loss and the detection loss. The seg- U-net ([34]). Additionally, similar to the work by mentation loss L in this work is a Dice loss defined 7 1 c c c 1 1x1 Conv + Relu down up S 2c 2c 3x3 Conv + c 1x1 Conv + 1 Bn + Relu Spatial down up Softmax 4c 4c residual D Concat down up 8c 8c Add down up 16c 3x3 Conv + 3x3 Conv Relu Bn + Relu + Bn residual residual Add Concat Ch / 2 Add 3x3 Conv + 3x3 Conv + 3x3 Conv + Relu 3x3 Conv + 3x3 Conv Relu Bn + Relu Bn + Dp + Bn + Dp 3x3 stride-2 Bn + Relu + Bn Relu Conv + Bn 3x3 stride-2 + Relu transposed Conv (Ch x 2) + Bn (Ch /2) down up Figure 3: A joint segmentation and detection network for catheter tip detection. This figure shows an example network with 4 levels of depth (the number of down or up blocks). Meaning of abbreviations: Conv, 2D convolution; Bn, batch normalization; Relu, ReLU activation; Dp, dropout; Concat, concatenation; Ch, number of channels; S, segmentation output; D, detection output. The number above an image or feature maps indicates the number of channels; the number of channels in the residual network in a block is shown above the block; c is the basic number of channels, the channel number in the first down block. The number next to a rectangle denotes the size of the image or feature maps. Red arrows indicate a change of number of channels. by Eq.(3): 4.3. Approximation of the Posterior with Particle Filter 2 M S i; j i; j i; j Once the deep neural network in Section 4.2 is L = 1 (3) P P 2 2 M + S trained, its weights are fixed during inference for i; j i; j i; j i; j computing the posterior p(x jz ) for new data. Ide- k 0:k where M denotes the ground truth binary catheter aly, the network detection output p(z jx ) should be a k k masks, S is the segmentation heatmap. The loss Gaussian PDF during inference, as it is trained with function for detection L is mean square error (MSE) labels of Gaussian PDFs. However, due to simi- given by Eq.(4): lar appearance of background structures or contrast residual, the detection output is unlikely to be a per- fect Gaussian (possibly non-Gaussian or having mul- L = jT D j (4) d i; j i; j w h tiple modes), which prevents the posterior p(x jz ) k 0:k iw; jh in Eq.(1) being solved with an analytical method. In where T denotes the ground truth PDF, w and h are practice, the posterior can be approximated using a the width and height of an image. The total training particle filter method ([3]). loss L is defined as Eq.(5): Particle filter methods approximate the posterior PDF by a set of N random samples with associated i i s weights fx ; w g ([3]). As N becomes very large, L = L + L (5) k k i=1 s d this discrete representation approaches the true pos- where  is a weight to balance L and L . terior. According to [3], the approximation of the s d 256 posterior p(x jz ) is given by Eq.(6): weights to be 1=N , so the number of e ective sam- k 0:k s ples which have actual contribution to approximate p(x jz ) is maximized ([3]). If the resampling is ap- k 0:k i i p(x jz )  w (x x ) (6) k 0:k k k k plied at every time step, the particle filter becomes a i=1 sampling importance resampling (SIR) filter, and the where () is the Dirac delta function. The weight weight update rule follows Eq.(9). w can be computed in a recursive manner as Eq.(7) i i once w is known ([3]): w / p(z jx ) (9) k1 k k k i i i The final decision on catheter tip location in frame k p(z jx )p(x jx ) k k k1 i i w / w (7) can then be computed as the expectation of x , x ˆ = k k1 i i k k q(x jx ; z ) k k1 x p(x jz )dx , which is in this case, the weighted k k 0:k k where q(x jx ; z ) is an importance density from sum of all samples: k k k1 which it should be possible to sample x easily. For simplicity, a good and convenient choice of i i x = w x : (10) i k k k the importance density is the prior p(x jx ) ([3]), k1 i=1 so that the weight update rule (7) becomes w / i i w p(z jx ). k 4.4. Summary k1 k A sample can be drawn from p(x jx ) in the fol- k1 The overall catheter tip tracking using a deep lowing way. First, a process noise sample v is k1 learning based Bayesian filtering method is summa- sampled from p (v ), the PDF of v ; then x is v k1 k1 rized in Algorithm 1. generated from x via the state transition model k1 i i i x = f (x ; v ). In this work, p (v ) is set to be a k v k1 k k1 k1 2 5. Experimental Setup GaussianN (0;  I). The choice of motion model for f is important for an accurate representation of the 5.1. Data true state transition prior p(x jx ). A random mo- k k1 Anonymized clinical imaging data were used for tion cannot characterize well the motion of catheter our experiments. The data were acquired with stan- tip in XA frames. In this work, we estimated the dard clinical protocol using Siemens AXIOM-Artis motion from adjacent frames using an optical flow system, and are from 55 patients who underwent a method, as this approach 1) takes into account of the PCI procedure at the Department of Cardiology at observation z , which results in a better guess of the Erasmus MC in Rotterdam, Netherlands. Out of catheter tip motion, and 2) enables estimation of a these data, we selected data from 37 patients which dense motion field where the motion of a sample x were acquired since the year 2014 to develop our can be eciently obtained. Therefore, f is defined method, and used the data from the other 18 patients as Eq.(8): acquired before the year 2013 for evaluation. The detailed information about the data is listed in Table x = x + u (x ) + v (8) k k1 k1 k1 k1 where u () is the motion from frame k1 to frame In order to evaluate the proposed roadmapping k1 k estimated with optical flow using the method by method, for which an o -line angiographic sequence [14], u (x ) is the motion from state x . is required for roadmap preparation and an on- k1 k1 k1 Once samples are drawn and their weights are line fluoroscopic sequence taken from the same C- updated, the so-called “resampling” of the samples arm position is needed for performing the actual should be performed to prevent the degenaracy prob- roadmapping (see Section 2), we selected the con- lem, where all but one sample will have negligi- trast frames from a real clinical sequence to simu- ble weight after a few iterations ([3]). The resam- late the o -line sequence, and chose the non-contrast pling step resamples the existing samples according frames from the same clinical sequence to simulate to their updated weights and then resets all sample the online sequence. The selected contrast sequence 9 Algorithm 1 Deep learning based Bayesian filtering for online tracking of catheter tip in X-ray fluoroscopy Require: fz ; : : : ; z g (sequentially observed frames), D (A trained network from Section 4.2), p(x ) (the 0 T 0 initial PDF),  (the variance of v ; k = 1; : : : ; T ), T (number of frames for tracking), N (number of k1 s samples) i i 1: Draw x  p(x ), set w = 1=N , 8 i = 1; : : : ; N 0 s s 0 0 2: for k = 1 to T do 3: Compute u from z to z using the optical flow method of [14] k1 k1 k 4: for i = 1 to N do i 2 5: Draw v  N (0;  I) k1 i i i 6: Compute the motion of x : u = u (x ) k1 k1 k1 k1 i i i i i i 7: Draw x  p(x jx ): x = x + u + v k k1 k k1 k1 k1 i i i 8: Update weight w = p(z jx ) = D (x ) k z k k k 9: end for i i s i 10: Normalize w w = w , 8 i = 1; : : : ; N k k i=1 k s i i 11: Prediciton in frame k: x ˆ = w x i=1 k k i i s i 12: Resample fx ; w g using the method of [3] (so all w are set to 1=N again) k k i=1 k 13: end for were ensured suciently long to cover at least one complete cardiac cycle. 5.2. Data Split for Catheter Tip Detection and Table 1: Basic information of the acquired X-ray image data Tracking for our experiments. The number in the parenthesis next to the To develop the catheter tip tracking method, 1086 pixel size indicates the possible image size. X-ray fluoroscopic images selected from 260 non- contrast sequences of 25 patients from the develop- Data Development Evaluation ment set were used for training the network from Fig- No. patients 37 18 ure 3; 404 images from 94 non-contrast sequences of No. sequences 354 34 another 12 patients from the development set were Frame rate (fps) 15 15 used as validation set for the network model and hy- Image size (px) 512  512 512  512 perparameter selection. In the training and validation 600  600 600  600 sets, 4-5 frames were randomly selected from each 776  776 776  776 sequence, which are not necessarily continuous. To 960  960 1024  1024 tune the parameters for tracking, 1583 images from 1024  1024 88 sequences out of the 94 from the same 12 patients Pixel size (mm) 0.108 (1024) 0.139 (1024) of the validation set were used (6 sequences were 0.139 (1024) 0.184 (600) not selected for this task due to very short sequence 0.184 (600) 0.184 (776) length not more than 5 frames). Finally, to evaluate 0.184 (776) 0.184 (1024) catheter tip tracking accuracy, 1355 images from 34 0.184 (960) 0.216 (512) non-contrast sequences of 18 patients from the eval- 0.184 (1024) 0.279 (512) uation set were used for testing. The frames selected 0.216 (512) for tracking from each sequence must be continuous; the number of selected frames for tracking might vary, depending on the number of the non-contrast frames in the sequences. Details of the datasets for training, validation and test are listed in Table 2. 10 Table 2: Dataset of training, validation and test for detection and tracking of catheter tip in X-ray fluoroscopic frames. Training Validation Validation Test (detection) (detection) (tracking) (tracking) No. patients 25 12 12 18 No. sequences 260 94 88 34 No. frames 1086 404 1583 1355 Continous frames? No No Yes Yes 5.3. Experimental Settings for Training the Deep e.g. coronary arteries in our case, is not directly vis- Network ible in the target image. One possible choice intro- duced by [50] is to use the guidewire as a surrogate This section describes the basic experimental set- of the target vessel centerline in non-contrast images, tings for training the deep neural network. Details of as guidewire is always inside vessels and commonly the training setup can be found in Appendix A. present in image sequences during interventions. In this work, we follow a similar strategy to evaluate the 5.3.1. Preprocessing accuracy of dynamic coronary roadmapping. As the image data have di erent size ranging from The first step is to select frames for roadmapping 512  512 to 1024  1024, all images were resam- evaluation. From each non-contrast sequence in the pled to a grid of 256 256 before being processed by test set for tracking in Section 5.2, we uniformly se- the neural network. In addition, the image intensities lect 8-20 frames to annotate guidewire. The num- were rescaled to the range from 0 to 1. ber of the selected frames from each sequence de- pends on the sequence length, the frame interval size 5.3.2. Training label and guidewire visibility. For some rare cases in our The standard deviation  of the Gaussian PDF for data where no guidewire is present in the image, we the training label of the detection network was set to discarded that non-contrast frame, and chose those 4 pixels in the resampled image space (256  256). frames with little vessel contrast from the same se- This choice corresponds to the estimation of the quence and annotated the vessel centerline. The se- maximal possible catheter tip radius. An example lection results in 409 frames from 34 sequences in to- of the Gaussian PDF is shown in Figure 2c. tal. Once the target non-contrast frames for evaluat- ing roadmapping are chosen, their corresponding an- 5.3.3. Evaluation Metric giographic frames were found using the ECG match- To select hyperparameters and model weights in ing method in Section 3. We then annotated the cen- training, an evaluation metric is required. As the terline of the vessel corresponding to the guidewire deep network is essentially a catheter tip detector, ac- in the non-contrast frames. curate detection of the tip location is desired. There- The next step is performing the transformation of fore, we chose the location with the highest value the labelled vessel centerline from the angiographic in the detection output, and computed the Euclidean frame to its corresponding target non-contrast frame distance between the chosen location and the ground via displacement of catheter tip in the two frames. truth tip coordinate as the evaluation metric to tune This step simulates the roadmapping transformation the deep network. in the last step in Figure 1. Finally, the distance between the guidewire anno- 5.4. Setup for Evaluating Dynamic Coronary tation in the target frame and the transformed vessel Roadmapping centerline is reported as the roadmapping accuracy. It is in general a challenge to evaluate the In order to compute the distance between two point roadmapping accuracy, as the structure of interest, sets of annotations (e.g. Figure 4a), point-point cor- 11 respondence between the two sets is required (Fig- 6. Experiments and Results ure 4b). The point sets were firstly resampled with the point interval being 1 mm. We then followed The following experiments are performed to as- the approach of [44] to find such correspondences sess the methods. First, In Section 6.1, the training which minimizes the sum of the Euclidean distance of the deep neural network is described. Then in Sec- of all valid point-point correspondence paths. This tion 6.2, the accuracy of catheter tip tracking using way guarantees no cross-over connection and each the optimized trained network and the tuned particle point in one set is connected to at least one point in filter is presented. Section 6.3 describes the accuracy the other set. As the annotated point sets may have evaluation of dynamic coronary roadmapping via the di erent size, the point correspondences to endpoints proposed catheter tip tracking method. Finally, in are excluded such that we only focused on the dis- Section 6.4, we measure the processing time of the tance between corresponding sections, not the entire proposed DCR approach. centerlines (Figure 4c). Once the point-point corre- spondence is available, the distance between the two 6.1. Training the Deep Neural Network points in a pair can be used for evaluating the accu- racy of DCR. The purpose of this experiment is to train the deep neural network to output reasonable likelihood prob- ability map. The network hyperparameters were tuned to optimize the detection performance. The training and validation data for detection men- tioned in Section 5.2 were used for training the deep neural network. The evaluation metric mentioned in Section 5.3, the mean Euclidean distance between the ground truth and the predicted tip location av- (a) (b) (c) eraged over all validation frames, was used as the validation criteria for selecting the optimal train- Figure 4: Correspondence between the labelled guidewire (green) and the transformed vessel centerline (red). The yellow ing epoch and the network hyperparameters. When lines connecting the two point sets illustrate the correspondence we evaluated hyperparameter settings, we firstly se- between red and green points. lected the training epoch with the lowest mean vali- dation error for each setting, then the settings were compared using the model weights (trainable net- work parameters) of their chosen epochs. 5.5. Implementation The network hyperparameters we investigated in the experiments include (1) the basic channel num- ber, i.e. the number of channels or feature maps in The proposed method was developed in Python. the first down block, (2) the network depth level, the The framework used for developing the deep learn- number of down or up blocks, and (3) the dropout ing approach for likelihood approximation is Py- probability. Torch. The major experiments of dynamic coronary roadmapping were performed on a computer with an The validation errors for di erent hyperparameter Intel Xeon E5-2620 v3 2.40 GHz CPU and 16 GB settings using the experimental settings in Section RAM running Ubuntu 16.04. The deep neural net- 5.3 are shown in Table 3. The table shows that the work and the optical flow method were running on hyperparameter setting with the lowest mean error, an NVIDIA GeForce GTX 1080 GPU. The approach which has 4 level in depth and 64 channels in the for evaluating dynamic coronary roadmapping was first down block, achieves a validation error of about developed and running in MeVisLab on a computer 2 mm. The table also shows other good choices of with an Intel Core i7-4800MQ 2.70 GHz CPU and network architecture that have a small validation er- 16 GB RAM running Windows 7. ror (shown in red in Table 3): 32 channels in the first 12 in Section 5.2 (see Appendix B for details). Then 1.2 training training 25 in Section 6.2.1, we evaluated the tracking accuracy 1.0 validation validation with the tuned optimal parameter setting (see Table 0.8 0.6 4) on the test dataset, and compared the proposed 0.4 5 tracking method with alternative approaches using 0.2 only the detection network in Section 4.2 or using 0 20 40 60 80 100 0 20 40 60 80 100 only optical flow. Finally, in Section 6.2.2, we inves- (a) Total loss (b) Detection error (mm) 0.5 tigated tracking accuracy with di erent ways of tip training 0.07 training 0.4 validation 0.06 initialization in the first frame. validation 0.05 0.3 0.04 0.2 0.03 6.2.1. Tracking Methods Evaluation 0.02 0.1 0.01 In this experiment, the proposed tracking method 0.0 0.00 0 20 40 60 80 100 0 20 40 60 80 100 in Algorithm 1 uses the ground truth tip probabil- (c) Segmentation loss (d) Detection loss ity map of the first frame as the initial PDF p(x ) to Figure 5: Learning curves for the chosen hyperparameter set- draw samples. This method is referred to as “Track- ting. ing”. In addition, we compared the proposed method with three alternatives. The first one tracks catheter down block with 4 or 5 levels in depth, or 64 chan- tip using only the detection network in Section 4.2 nels with 3 or 4 depth levels. The dropout regulariza- with the chosen network architecture and trained pa- tion improves the accuracy of the model, compared rameters in Section 6.1, therefore, no temporal infor- to the ones without dropout. mation is used. This method is referred to as “De- The learning curves of the training process with tection (Net)”. The other two methods in this ex- the chosen hyperparameter setting are illustrated in periment use only optical flow to track catheter tip Figure 5. The curves indicate that both segmentation starting from the ground truth tip position in the first and detection reach convergence after training 100 frame. The motion field towards the current frame, estimated by the two methods, was based on the de- epochs. We did not investigate a model with more than 64 formation from the previous frame or the first frame channels or 5 depth levels, because (1) it will fur- in the sequence, respectively. The same implemen- ther increase the processing time which makes on- tation setting as in Appendix B.1 was used for these line applications less feasible; (2) the results in Table two methods. They are called “Optical Flow (pre- 3 show that such a setting (64 channels, 5 depth lev- vious)” and “Optical Flow (first)”, or in short form, els) starts increasing the validation error compared to “OF (pre)” and “OF (1st)”. Additionally, we refer those less complex models. the interested readers to Appendix C.1 where the in- The subsequent experiments will be based on the fluence of catheter segmentation on the detection and tracking approaches is reported. network trained with the chosen hyperparameter set- ting indicated in Table 3 (64 channels, 4 depth levels, The tracking accuracies of all methods reported in dropout 0.2, also see Table 4). this section were obtained on the test set from Table 2. The mean, the median and the maximal tracking 6.2. Catheter Tip Tracking error between the predicted and the ground truth tip The purpose of this experiment is to assess the position of all test images are reported in Table 5. accuracy of catheter tip tracking with the proposed In addition, as the sequences in the test set have dif- method in Section 4. Guiding catheter tip is tracked ferent lengths, we also computed the mean and the in X-ray fluoroscopy using Algorithm 1 based on a median error per sequence, and report the the aver- trained network with the optimal hyperparameter set- age of the sequence mean and median errors, so that ting from Section 6.1. First, the parameters of the each sequence contributes equally in these metrics. optical flow method used in Algorithm 1 and particle Table 5 shows that the results from the detection net- filter were tuned on the validation data for tracking work have large average errors which are caused by 13 Table 3: Validation errors (mm) for di erent hyperparameter settings. Red cells show the settings with the 10 smallest validation errors. bold number indicates the setting with the lowest error. Basic Number Depth Dropout of Channels Level none 0.1 0.2 0.3 0.4 0.5 8 3 5.43 4.99 5.02 5.37 4.38 4.24 4 4.17 4.45 4.25 5.04 4.75 4.36 5 3 4.14 3.53 4.28 3.95 4.11 16 3 3.74 4.29 3.57 4.11 3.74 3.4 4 3.36 3.11 3.63 3.33 3.36 3.78 5 3.38 2.89 3.16 2.52 2.71 2.74 32 3 2.99 3.02 3.26 2.82 3.26 2.56 4 2.87 2.34 2.46 2.6 2.65 2.27 5 3.04 2.51 2.21 2.29 2.3 2.25 64 3 2.19 2.54 2.34 2.27 2.26 2.49 4 2.55 2.31 2.04 2.44 2.22 2.27 5 2.42 2.29 2.73 2.77 2.61 2.85 Table 4: The chosen (hyper-)parameters for di erent building blocks of the catheter tip tracking algorithm. The parameters of the optical flow method can be found in Appendix B.1. Building block (hyper-)parameters value Deep learning Basic channel number 64 Depth 4 Dropout 0.2 OF (pre) OF (1st) Detection Tracking OF (pre) OF (1st) Detection Tracking (a) Overall view of tracking er- (b) A zoom-in view of (a) Particle filter  (px) 5 rors N 1000 Figure 6: Tracking errors for the 4 methods on all test images. some completely failed cases. The proposed tracking racy when the target is on the track (row 4), they method has median errors of about 1 mm and mean errors of about 1.3 mm. It achieves the lowest errors present periodic error patterns in two sequences due to large cardiac motion. The detection method shows compared to the other 3 methods on all listed evalu- ation criteria. peaks of large errors, this is because temporal rela- tion between frames is not modeled by the approach, Figure 6 illustrates the boxplots of tracking er- thus the detection on di erent frames is independent rors made by the 4 methods on all test images. It of each other. The proposed tracking method over- shows that the proposed tracking approach outper- comes the problems that other methods have and forms the detection method by avoiding making ex- presents accurate detection on these 4 sequences. tremely large errors (Figure 6a); meanwhile, it main- The tracking results of the 4 methods on example tains as accurate as the detection method for cases frames from the 4 sequences are illustrated in Figure with small errors, and is more accurate than the methods based solely on optical flow (Figure 6b). Figure 7 shows longitudinal views of tracking er- Figure 9 illustrates how the proposed tracking rors of the 4 methods on 4 example sequences. Al- method works on the same 4 frames in Figure 8. though the optical flow methods show high accu- It shows that the prior hypotheses (samples) assists Ground truth / Prediction Distance (mm) Ground truth / Prediction Distance (mm) Table 5: Catheter tip tracking errors (mm) of the 4 methods on the test (tracking) dataset. y indicates that the di erence between that method and the “Tracking” method are statistically highly significant with the two-sided Wilcoxon signed-rank test (p < 0:001). Optical Flowy Optical Flowy Detection Nety Tracking Evaluation Metrics (previous) (first) (Section 4.2) Maximal error of all images 29.16 20.83 108.20 17.72 Median error of all images 1.78 1.22 0.96 0.96 Mean error of all images 3.74  4.93 3.05  4.05 5.62  15.91 1.29  1.76 Average of sequence median error 2.35  2.52 2.64  3.52 6.26  17.11 1.03  0.49 Average of sequence mean error 2.59  2.69 3.31  2.81 6.83  13.88 1.29  0.94 to focus on the correct target location and results views (Figure 7), the outliers for the tracking with au- in reliable posterior estimation, especially when the tomatic initialization are more consistent over time. detection produces ambiguity in cases of multiple Figure 11 shows the temporal change of tracking er- catheters or contrast residual presented in images. rors for the 6 sequences with outliers using the track- ing with automatic initialization. For the 3 sequences 6.2.2. Catheter Tip Initialization on the top row, the tracking with automatic initial- ization makes large errors at the beginning, but be- In this experiment, the initial PDF p(x ) from comes accurate very fast in a few frames; for the 3 which samples are drawn in the proposed tracking sequences on the bottom row, however, the tracking is investigated (Algorithm 1). In particular, we ex- errors remain large till the end of the sequences. plored and evaluated the tracking accuracy with an Figure 12 shows example frames to give an in- automatic initialization using the probability map ob- sight of the tracking with automatic initialization on tained from the trained detection network in Section the 6 sequences in Figure 11. For the 3 sequences 4.2 with the chosen setting in Section 6.1. on the top row (Figure 12a), although the initializa- Figure 10 shows the boxplot of tracking errors on tion on the first frame (frame 0) is overall not cor- all test images with automatic initialization (Auto) rect, the true tip positions are still covered by some and manual initialization (Manual) for which the samples; once the detection in subsequent frames is ground truth tip probability map of the first frame correct, the tracker can still converge to the right tar- was used. The tracking with automatic initialization get. For the 3 sequence on the bottom row (Figure presents an accuracy similar to the one with manual 12b), the initializations of samples are ambiguous in initialization for small tracking errors, but has more frame 0; the detection in subsequent frames focuses large tracking errors which influence the mean error on a wrong area also given by the initial samples due over all test images (Table 6). We, therefore, de- to residual of contrast agent or multiple catheters, the fined the tracking errors on the right side of the gap tracker then tends to find the wrong target. in the boxplot (> 40 mm) as outliers, and explored the statistics without those outliers. 6.3. Dynamic Coronary Roadmapping Table 6 indicates that, the mean and median er- ror of the tracking with automatic initialization ex- In this experiment, the accuracy of dynamic coro- cluding the outliers are only slightly higher than the nary roadmapping using the proposed method with tracking with manual initialization and the detection manual tip initialization was evaluated. For roadmap method. While the tracking with automatic initial- selection with ECG matching (Section 3), the num- ization has 100 outliers in total from 6 sequences, ber of online ECG signal points N was manu- ECG the detection method that has 10 sequences contain- ally determined so that the ECG signal stored in the ing 106 outliers. bu er corresponding to 12 X-ray frames (0.8 second Unlike the detection method for which the outliers in acquisition time). Following the setup in Section are mainly presented as the peaks in the longitudinal 5.4, we used the distance between the two points in 15 Optical Flow (previous) Optical Flow (first) Detection Tracking 80 80 80 60 60 60 40 40 40 20 20 0 0 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 80 80 80 80 60 60 60 60 40 40 40 20 20 20 20 0 0 0 0 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 80 80 80 60 60 60 60 40 40 40 20 20 20 20 0 0 0 0 0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25 80 80 60 60 60 40 40 40 40 20 20 20 20 0 0 0 0 0 20 40 60 80 100 120 0 20 40 60 80 100 120 0 20 40 60 80 100 120 0 20 40 60 80 100 120 Figure 7: Longitudinal view of tracking errors made by the 4 methods on 4 test sequences (one sequence per row). The x-axis denotes the time steps of a sequence, the y-axis is the tracking error (mm). each point pair as the evaluation metric for DCR (the and the frame mean point distances of all 409 evalua- length of a yellow line segment in Figure 4). As each tion frames are illustrated in Figure 13. The compar- frame may have di erent numbers of point pairs, de- ison of the three DCR approaches from Table 7 and pending on the length of the target guidewire, the av- Figure 13 indicates that the accuracy of the proposed erage point pair distance per frame was also com- DCR method has shown improvement over the DCR puted for evaluation. These distances were evaluated without tracking, and is only slightly less than the on 409 selected frames with manual annotation of DCR with manual tip tracking (although the di er- guidewires and vessel centerlines (Section 5.4). ence is statistically significant). Additionally, inter- ested readers are referred to Appendix C.2 where the In the experiment, we compared the DCR with the influence of catheter segmentation on the accuracy of proposed tracking method to those with manual tip DCR is investigated. tracking and without tracking. All three approaches were based on the same ECG matching method (Sec- Table 8 shows how the frame mean point distances tion 3) for selecting roadmaps. The accuracy of the of the 409 evaluation frames are distributed. The DCR without tracking in Table 7 shows that the mean DCR with the proposed method has similar error dis- distances are reduced to less than 3 mm by com- tribution as the one with manual tip tracking: they pensating only cardiac motion via roadmap selection both have about 1/3 of the distances less than 1 mm with ECG matching. Table 7 also shows that the and 1/3 of the distances between 1 and 2 mm. The DCR with the proposed method achieves median dis- proposed method has slightly more distances larger tances of about 1.4 mm and mean distances of about than 5 mm than manual tip tracking. Both methods 2 mm. The boxplots of the distances of all point pairs are more accurate than the DCR without tracking on 16 Sequence 1 Sequence 2 Sequence 3 Sequence 4 Optical Flow Optical Flow Detection Input Tracking Ground Truth (previous) (first) Net Figure 8: Tracking results on example frames from the same 4 sequences in Figure 7. The blue point indicates the predicted catheter tip location; the red point shows the ground truth location. (Best viewed in color) intervals of small errors (< 2 mm). network and the optical flow component of the track- ing method were running on the GPU. Figure 14 shows overlays of selected roadmaps on example frames of 4 sequences with the three DCR In the experiments, the runtimes for roadmap se- approaches. The DCR without tracking presents mis- lection (step 1) and roadmap transformation (step 3) match of catheters, guidewires or residual of contrast in Figure 1 were negligible (< 1 ms / frame). The agent in the images, whereas the other methods im- runtime of the proposed catheter tip tracking method prove the alignment and show good match between is shown in Table 9 and Figure 15. The average the structures in the original X-ray image and the time to compute the likelihood with the deep learning roadmaps. Compared to the DCR with manual tip setup (DL) is 31.5 ms / frame. The particle filtering tracking, the proposed method show similar visual (PF) step, which consists of the optical flow estima- alignment of the roadmaps to the original X-ray im- tion, sample propagation, sample weight update and ages. For a dynamic view of a roadmapping exam- normalization, prediction and resampling, takes on ple, we refer readers to the supplemental material. average 23 ms / frame. Therefore, the average track- ing time in total is 54.5 ms / frame. The total av- 6.4. Processing Time erage time of the proposed DCR including roadmap The processing time of all steps in the proposed selection, catheter tip tracking and roadmap transfor- DCR method was measured with the hardware and mation is still less than the acquisition time of our software setup in Section 5.5. The ECG matching data (66.7 ms / frame, 15 fps), indicating that the method for roadmap selection was running in Python proposed DCR method would run in real-time with on the CPU of the Linux machine; the deep neural our setup. 17 Sequence 1 Sequence 2 Sequence 3 Sequence 4 Manual Auto Segmentation Detection Particles Particles Prediction Input Ground Truth (heatmap) (likelihood) (prior) (posterior) (expectation) Figure 9: Workflow of the proposed tracking method on the same 4 frames in Figure 8. The high probability is shown with bright color in the detection map. Samples or particles are presented as green dots. The blue point indicates the predicted catheter tip location; the red point shows the ground truth location. (Best viewed in color) 50 60 40 60 0 0 0 0 20 40 60 80 100 0 5 10 15 20 25 0 10 20 30 40 50 60 70 60 80 60 60 40 30 40 20 20 0 20 40 60 80 100 0 0 0 Tracking error (mm) 0 3 6 9 12 15 18 0 10 20 30 40 50 60 70 0 1 2 3 4 5 6 7 Figure 10: Catheter tip tracking errors (mm) with manual and Figure 11: Longitudinal views of tracking errors (mm) for the automatic initialization. 6 sequences with outliers using automatic initialization. 7. Discussion using a proposed deep learning based Bayesian fil- We have presented a new approach to perform on- tering. The proposed tracking method represents and line dynamic coronary roadmapping on X-ray flu- tracks the posterior of catheter tip via a particle fil- oroscopic sequences for PCI procedures. The ap- ter, for which a likelihood probability map is com- proach compensates the cardiac-induced vessel mo- puted for updating the particle weights using a con- tion via selecting oine-stored roadmaps with ap- volutional neural network. In the experiments, the propriate cardiac phase using ECG matching, and proposed DCR approach has been trained and eval- corrects the respiratory motion of vessels by online uated on clinical X-ray sequences for both tracking tracking of guiding catheter tip in X-ray fluoroscopy and roadmapping tasks. 18 Table 6: Catheter tip tracking errors (mm) of detection and tracking with manual and automatic initialization Tracking Detection Manual init. Automatic init. Maximal error 108.20 17.23 98.58 Median error 0.96 0.96 0.96 Mean error 5.62  15.91 1.29  1.76 5.16  13.91 No. of outliers (> 40 mm) 106 0 100 No. of sequences with outliers 10 0 6 Maximal error of inliers 31.06 17.23 28.28 Median error of inliers 0.96 0.96 0.96 Mean error of inliers 1.17  1.78 1.29  1.76 1.34  2.15 Table 7: The statistics of DCR accuracy (mm) with three di erent tracking scenarios. With the two-sided Wilcoxon signed-rank test: y denotes that the di erence between the DCR without tracking and that with the proposed tracking method is statistically highly significant (p < 0:001); * indicates a statistically significantly di erence between the DCR using manual tip tracking and that with the proposed tracking approach (p < 0:05). Without Trackingy Proposed Tracking Method Manual Tip Tracking* All point pairs Maximal distance 27.19 20.24 13.12 Median distance 1.97 1.43 1.35 Mean distance 2.94  2.83 2.07  2.08 1.85  1.72 Frame mean distance Median distance 2.11 1.42 1.38 Average distance 2.76  2.08 1.91  1.52 1.75  1.30 One prerequisite of accurate tracking with the pro- formance does not increase much compared to those posed approach is to obtain a reasonably good like- with simpler settings, implying that the model starts lihood estimation, which requires to train the deep overfitting on our dataset. neural network to detect catheter tip well. In this In addition to the deep neural network, the other work, we have investigated the influence of three important component of the proposed tracking ap- network hyperparameters on the performance of the proach is the sampling in the particle filter that yields detection network (Section 6.1): the basic channel the samples for representing the prior and the poste- number and network depth level are model capac- rior of catheter tip position. First, a sucient number ity parameters, the dropout adds regularization to the of samples in the whole sample space are required model. The experiment showed that the detection ac- to well characterize the probability distributions (see curacy improves when the basic channel number and Appendix B.2). Second, the sample dynamics plays the network depth level increase (Table 3). This ob- an important role in tracking, in particular, as indi- servation matches the expectation that a more com- cated by Eq.(8), the process noise and the sample plex model has higher capacity to model the variation motion. The process noise influences the tracking in the data, hence results in better accuracy. When accuracy, according to Table 10 in Appendix B.2. the complexity reaches a certain level, e.g. 64 ba- Additionally, sample motion is another key aspect of sic channels and 5 level of depth, the network per- sample dynamics. Motion estimation has previouly 19 Table 8: Distribution of frame mean point distances of the 409 evaluation frames. Error Intervals (mm) Tracking Methods of DCR < 1 1-2 2-3 3-4 4-5  5 Without tracking 81 115 69 47 31 66 Proposed tracking method 131 145 61 32 17 23 Manual tip tracking 139 144 61 35 20 10 Table 9: Statistics of the runtime of catheter tip tracking (ms / frame) on the test (tracking) dataset. Deep Learning Particle Filtering Total Tracking Time Mean 31.5  10.3 23.0  8.7 54.5  12.3 Median 35.1 22.8 57.7 been incorporated in a motion-based particle filter, pecially when the CNN detector fails to find the cor- such as adaptive block matching ([6]). In our work, rect target area. The association of knowledge from optical flow was chosen for motion estimation, as two sources together improves the tracking accuracy its non-parametric nature allows to characterize the compared to each single source. complexity of motion in X-ray fluoroscopy well. In The initial state is a also key component of track- addition, the advantage of such approach from a the- ing approaches. In the context of Bayesian filtering, oretical point of view is that it takes into account of the initial state provides the prior knowledge of the the current observation, leading to a more optimal tracking target. Most tracking algorithms assume a importance density ([3]) compared to random mo- known initial state from which the target is tracked, tion. e.g. our proposed method with manual initialization The tracking results in Section 6.2.1 show that in Section 6.2.2. In this case, the prior knowledge is the proposed tracking approach is able to track the provided by human. In Section 6.2.2, we also investi- catheter tip in X-ray fluoroscopy accurately with an gated a scenario where the initial state is given by the average tracking error of about 1.3 mm. It also shows detection network, so that the complete tracking pro- advantages over methods based only on optical flow cess is fully automated. The results indicate that, the or the detection network. The OF (pre) method relies proposed tracking method with automatic initializa- heavily on tracking in the previous frame, hence the tion works reasonably well on most sequences even error could accumulate. The OF (first) method may when the initialization is sometimes incorrect (Fig- su er from large motion from the first frame to the ure 12a). This is because (1) the true position is cov- current frame. The detection method uses informa- ered by a few samples, and (2) the correct detection tion only from the current frame, no temporal rela- in later frames corrects the initial mistake in the first tion between frames is utilized; therefore, it results in frame. The automatic initialization fails when (1) spikes in the longitudinal view, as shown in Figure 7. a wrong position is covered by a few samples and The proposed tracking method has a CNN to provide (2) the wrong detection in subsequent frames con- an accurate observation on the current frame which firms the mistake in the initial frame (Figure 12b). improves the accuracy of optical flow tracking within This happens when there is contrast agent remaining the framework of Bayesian filtering. In the mean- in the image or there are multiple catheters, which time, the optical flow based particle filter maintains are the major sources causing ambiguity in detec- and propagates the prior knowledge from the initial tion. In practice, the automatic initialization would tip position to provide a constraint on searching for work well when contrast agent is washed out and the potentially correct positions, which is useful es- only one catheter is present in the field of view, oth- 20 Detection Initial Particles Detection Posterior Frame 0 Frame 10 Frame 0 Frame 10 Frame 0 Frame 0 Frame 5 Frame 5 Figure 13: Accuracy (mm) of DCR with three di erent tracking scenarios. Frame 20 Frame 0 Frame 20 Frame 0 (a) Sequence 1-3 on the top row in Figure 11 Detection Initial Particles Detection Posterior Original X-ray Without Catheter With The Proposed With Manual Tip Image Tip Tracking Method Tracking Frame 0 Frame 0 Frame 4 Frame 4 Frame 0 Frame 15 Frame 0 Frame 15 Frame 7 Frame 0 Frame 0 Frame 7 (b) Sequence 4-6 on the bottom row in Figure 11 Figure 12: Examples frames from the 6 sequences in Figure 11. The high probability in the detection heatmap is highlighted as bright color. Particles are presented as green dots. The red dots in the last column indicate the ground truth tip location. (Best viewed in color) Figure 14: Examples of superimposition of selected roadmaps (red) on X-ray fluoroscopic frames. (Best viewed in color) 21 Image acquisition data. Additionally, variation exists between di erent DL time 66.7 ms / frame cardiac cycles ([27]), therefore, choosing a roadmap from another cycle may cause inaccuracy for cardiac PF motion compensation. Finally, the way of DCR eval- Total uation in Section 5.4 might also introduce inaccura- 0 20 40 60 80 100 cies in the error measurement. Since guidewires of- Processing Time (ms) ten attach to the inner curves of a vessel to take the Figure 15: Runtime of catheter tip tracking (ms / frame) on all shortest path, the small di erence between the an- test frames. notation of guidewire and vessel centerlines was ig- nored in the evaluation. erwise manual initialization would be needed which In addition to accuracy, processing speed is also requires only one click to initiate tracking. critical for intraoperative applications. The results in Section 6.4 indicate that the total processing time of Dynamic coronary roadmapping is the direct ap- the proposed DCR approach is less than the image plication of the catheter tip tracking results. In our acquisition time, meaning that it runs in real-time on experiments, the DCR was performed with manual our setup. To build a real-time system for PCI in tip initialization to show the potential of the proposed practice, the overall latency of the complete system tracking method, and was compared with the DCR needs to be considered. It is also worth noticing that without tracking and with manual tracking. The the DL and PF steps of the proposed tracking method results indicate that using catheter tip tracking can are independent from each other. In practice, in case improve DCR accuracy, as the respiratory-induced more than one GPU are available, the proposed DCR vessel motion is corrected by the displacement of approach can be further accelerated by paralleling catheter tip in addtion to cardiac motion correction. the DL and PF steps, making them running on dif- The results also show that the proposed DCR reaches ferent GPUs. a good accuracy (mean error is about 2 mm) and per- Compared to the previous works on DCR, the pro- forms only slightly worse than its best case, the DCR with manual tip tracking which is not applicable for posed approach in this paper shows advancement in intraoperative use. Additionally, according to a pre- several aspects. First, our systems works on non- vious study by [10], the average lumen diameters of cardiac-gated sequences which does not require ad- human coronary arteries are between 1.9 mm (dis- ditional setups for cardiac motion gating that were tal left anterior descending artery) and 4.5 mm (left needed for some methods ([50, 26]). Second, our main artery). This means that the accuracy achieved approach compensates both respiratory- and cardiac- with the proposed approach is comparable with the induced vessel motion, which is more accurate than size of coronary arteries. systems that correct only cardiac motion ([12]). In addition, the proposed DCR approach follows a data- Apart from catheter tip tracking, several other pos- driven paradigm that learns target feature from se- sible factors in di erent steps of the experiments may quences acquired from di erent patients and various influence the final DCR accuracy. First, in the oine view angles, making it more robust than the method phase, the signal of contrast agent may become too that relies on traditional vesselness filtering ([20]) strong and completely cover the catheter tip, com- or methods that require specific tissue being present plicating the tip visibility in some cases. In this sit- ([50, 26]). These are the major advantages of the uation, the uncertainty in the manual tip annotation proposed DCR over the existing direct roadmapping may result in errors in roadmap transformation. Sec- systems. Compared to model-based motion compen- ond, in the roadmap selection step, the oine-stored roadmaps are only discrete samples of complete car- sation, our approach does not require the extraction of motion surrogate signals and train a motion model diac cycles which might not fully characterize every for each new patient, but can be directly run with a possible change in the cardiac motion. This problem trained model. could possibly be addressed in the future by inter- polating frames between the existing frames in the The proposed deep learning based Bayesian filter- 22 ing method has several advantages over the existing In the future, it may be worth investigating the instrument tracking methods. First, the deep learning following directions related to this work. As the component enables a more general framework to de- data used in this study was acquired from one hos- tect instruments in medical images than methods tai- pital using a machine from a single vendor, it would lored for specific tools ([25, 24]). Compared to the be interesting to evaluate the proposed approach on existing detection methods based on deep learning multi-center data acquired with machines from dif- ([5, 22, 11]), our approach takes into account of the ferent vendors. Next, since the ECG signals of our information between frames; the Bayesian filtering data appear to be regular, it may be necessary in a framework allows interaction between temporal in- future study to acquire data with irregular ECG that formation and the detection of a convolutional neural could be obtained in practice, and validate the pro- network, making the tracking more robust. Bayesian posed approach on those data. Besides, it would frameworks have been used in many previous tem- be also interesting to validate our approach during poral instrument tracking methods. Particularly, the PCI procedures in an environment simulating the real likelihood term in some works was designed based clinical settings. Additionally, from a methodolog- on registration or segmentation outcomes ([2, 38]) or ical point of view, although the proposed tracking traditional machine learning approaches with hand- method is invariant under di erent view angles, the crafted features ([45, 46, 29]). In our method, we whole DCR approach works only when the oine approximated the likelihood with a deep neural net- and online phase have the same view angle, i.e. it work learned from the clinical data which exempts is a 2D roadmapping system. Therefore, one future the need of feature engineering but yet possesses direction would be to develop a 3D DCR system that more discriminative power; the network directly out- would work with various view angles in the online puts the probability map, making it more straightfor- phase. ward to use. Finally, compared to the existing instru- ment tracking approaches based on Bayesian filter- 8. Conclusion ing ([2, 38, 39]), the state transition in our method was based on the motion estimated from two adja- We have developed and validated a novel approach cent frames, which is more reliable than totally ran- to perform dynamic coronary roadmapping for PCI dom motion or artificially-designed state transition image guidance. The approach compensates car- models. diac motion through ECG alignment and respiratory motion by guiding catheter tip tracking during fluo- From a practical point of view, the proposed DCR roscopy with a deep learning based Bayesian filter- approach could potentially fit into the clinical work- ing method. The proposed tracking and roadmap- flow of PCI. The oine phase of the method can ping approaches were trained and evaluated on clini- be done eciently by a person who assists the pro- cal X-ray image datasets and were proved to perform cedures: selecting and creating roadmaps from an accurately on both catheter tip tracking and dynamic angiography acquisition, annotating the catheter tip coronary roadmapping tasks. Our approach also runs (one point) in the images. This phase is typically in real-time on a setup with a modern GPU and thus done before a fluoroscopy acquisition during which has the potential to be integrated into routine PCI the guidewire advancement and stent placement are procedures, assisting the operator with real-time vi- performed. In the online phase, when a fluoroscopic sual image guidance. image is acquired, the proposed system selects the most suitable roadmap, tracks the catheter tip and transforms the roadmap to prospectively show a ves- Acknowledgment sel overlay on the fluorosocpic image. The online updated coronary roadmap can provide real-time vi- This work was supported by NWO-domain TTW sual guidance to cardiologists to manipulate inter- (The Division of Applied and Engineering Sciences ventional tools during the procedure without the need in The Netherlands Organisation for Scientific Re- of administering extra contrast agent. search), IMAGIC project under the iMIT program 23 (grant number 12703). Ries and Simon van Walsum [12] Elion, J.L., 1989. Dynamic coronary roadmapping. US Patent 4,878,115. are acknowledged for their contribution in the man- [13] Faranesh, A.Z., Kellman, P., Ratnayaka, K., Lederman, ual annotations. R.J., 2013. Integration of cardiac and respiratory motion into MRI roadmaps fused with X-ray. Medical physics References [14] Farneback, G., 2003. Two-frame motion estimation based on polynomial expansion, in: Scandinavian conference on [1] Ambrosini, P., Ruijters, D., Niessen, W.J., Moelker, A., Image analysis, Springer. pp. 363–370. van Walsum, T., 2017a. Fully automatic and real-time [15] Fischer, P., Faranesh, A., Pohl, T., Maier, A., Rogers, T., catheter segmentation in X-ray fluoroscopy, in: Inter- Ratnayaka, K., Lederman, R., Hornegger, J., 2018. An national Conference on Medical Image Computing and MR-based model for cardio-respiratory motion compen- Computer-Assisted Intervention, Springer. pp. 577–585. sation of overlays in X-ray fluoroscopy. IEEE transactions [2] Ambrosini, P., Smal, I., Ruijters, D., Niessen, W.J., on medical imaging 37, 47–60. Moelker, A., Van Walsum, T., 2017b. A hidden markov [16] Garc´ ıa-Peraza-Herrera, L.C., Li, W., Gruijthuijsen, C., model for 3D catheter tip tracking with 2D X-ray catheter- Devreker, A., Attilakos, G., Deprest, J., Vander Poorten, ization sequence and 3D rotational angiography. IEEE E., Stoyanov, D., Vercauteren, T., Ourselin, S., 2016. Trans. Med. Imaging 36, 757–768. Real-time segmentation of non-rigid surgical tools based [3] Arulampalam, M.S., Maskell, S., Gordon, N., Clapp, on deep learning and tracking, in: International Workshop T., 2002. A tutorial on particle filters for online on Computer-Assisted and Robotic Endoscopy, Springer. nonlinear/non-Gaussian Bayesian tracking. IEEE Trans- pp. 84–95. actions on signal processing 50, 174–188. [17] He, K., Zhang, X., Ren, S., Sun, J., 2016. Deep resid- [4] Baka, N., Lelieveldt, B., Schultz, C., Niessen, W., van ual learning for image recognition, in: Proceedings of the Walsum, T., 2015. Respiratory motion estimation in X- IEEE conference on computer vision and pattern recogni- ray angiography for improved guidance during coronary tion, pp. 770–778. interventions. Physics in Medicine & Biology 60, 3617. [18] Heibel, H., Glocker, B., Groher, M., Pfister, M., Navab, [5] Baur, C., Albarqouni, S., Demirci, S., Navab, N., N., 2013. Interventional tool tracking using discrete opti- Fallavollita, P., 2016. Cathnets: detection and single- mization. IEEE transactions on medical imaging 32, 544– view depth prediction of catheter electrodes, in: Interna- tional Conference on Medical Imaging and Virtual Real- [19] Honnorat, N., Vaillant, R., Paragios, N., 2011. Graph- ity, Springer. pp. 38–49. based geometric-iconic guide-wire tracking, in: Inter- [6] Bouaynaya, N., Schonfeld, D., 2009. On the optimality national Conference on Medical Image Computing and of motion-based particle filtering. IEEE transactions on Computer-Assisted Intervention, Springer. pp. 9–16. circuits and systems for video technology 19, 1068–1072. [20] Kim, D., Park, S., Jeong, M.H., Ryu, J., 2018. Regis- [7] Brost, A., Liao, R., Strobel, N., Hornegger, J., 2010. Res- tration of angiographic image on real-time fluoroscopic piratory motion compensation by model-based catheter image for image-guided percutaneous coronary interven- tracking during EP procedures. Medical Image Analysis tion. International journal of computer assisted radiology 14, 695–706. and surgery 13, 203–213. [8] Chang, P.L., Rolls, A., De Praetere, H., Vander Poorten, [21] King, A.P., Boubertakh, R., Rhode, K.S., Ma, Y., Chin- E., Riga, C.V., Bicknell, C.D., Stoyanov, D., 2016. Ro- chapatnam, P., Gao, G., Tangcharoen, T., Ginks, M., bust catheter and guidewire tracking using B-spline tube Cooklin, M., Gill, J.S., et al., 2009. A subject-specific model and pixel-wise posteriors. IEEE Robotics and Au- technique for respiratory motion correction in image- tomation Letters 1, 303–308. guided cardiac catheterisation procedures. Medical Image [9] Dannenberg, L., Polzin, A., Bullens, R., Kelm, M., Zeus, Analysis 13, 419–431. T., 2016. On the road: First-in-man bifurcation percu- [22] Laina, I., Rieke, N., Rupprecht, C., Vizca´ ıno, J.P., Eslami, taneous coronary intervention with the use of a dynamic A., Tombari, F., Navab, N., 2017. Concurrent segmenta- coronary road map and stentboost live imaging system. tion and localization for tracking of surgical instruments, International journal of cardiology 215, 7–8. in: International conference on medical image computing [10] Dodge, J.T., Brown, B.G., Bolson, E.L., Dodge, H.T., and computer-assisted intervention, Springer. pp. 664– 1992. Lumen diameter of normal human coronary arter- ies. influence of age, sex, anatomic variation, and left ven- [23] Ma, H., Dibildox, G., Banerjee, J., Niessen, W., Schultz, tricular hypertrophy or dilation. Circulation 86, 232–246. C., Regar, E., van Walsum, T., 2015. Layer separation for [11] Du, X., Kurmann, T., Chang, P.L., Allan, M., Ourselin, S., vessel enhancement in interventional X-ray angiograms Sznitman, R., Kelly, J.D., Stoyanov, D., 2018. Articulated using morphological filtering and robust PCA, in: Work- multi-instrument 2D pose estimation using fully convolu- shop on Augmented Environments for Computer-Assisted tional networks. IEEE transactions on medical imaging Interventions, Springer. pp. 104–113. 24 [24] Ma, Y., Gogin, N., Cathier, P., Housden, R.J., Gijsbers, image-guided cardiac interventions, in: 2010 IEEE Com- G., Cooklin, M., O’Neill, M., Gill, J., Rinaldi, C.A., puter Society Conference on Computer Vision and Pattern Razavi, R., et al., 2013. Real-time X-ray fluoroscopy- Recognition, IEEE. pp. 2948–2954. based catheter detection and tracking for cardiac electro- [36] Shechter, G., Ozturk, C., Resar, J.R., McVeigh, E.R., physiology interventions. Medical physics 40. 2004. Respiratory motion of the heart from free breath- [25] Ma, Y., King, A.P., Gogin, N., Gijsbers, G., Rinaldi, C.A., ing coronary angiograms. IEEE transactions on medical Gill, J., Razavi, R., Rhode, K.S., 2012. Clinical evalua- imaging 23, 1046. tion of respiratory motion compensation for anatomical [37] Shechter, G., Shechter, B., Resar, J.R., Beyar, R., 2005. roadmap guided cardiac electrophysiology procedures. Prospective motion correction of X-ray images for coro- IEEE Transactions on Biomedical Engineering 59, 122– nary interventions. IEEE Transactions on Medical Imag- 131. ing 24, 441–450. [26] Manhart, M., Zhu, Y., Vitanovski, D., 2011. Self- [38] Speidel, S., Delles, M., Gutt, C., Dillmann, R., 2006. assessing image-based respiratory motion compensation Tracking of instruments in minimally invasive surgery for fluoroscopic coronary roadmapping, in: Biomedical for surgical skill analysis, in: International Workshop on Imaging: From Nano to Macro, 2011 IEEE International Medical Imaging and Virtual Reality, Springer. pp. 148– Symposium on, IEEE. pp. 1065–1069. 155. [27] McClelland, J.R., Hawkes, D.J., Schae ter, T., King, [39] Speidel, S., Kuhn, E., Bodenstedt, S., Rohl, ¨ S., Ken- A.P., 2013. Respiratory motion models: a review. Medi- ngott, H., Muller ¨ -Stich, B., Dillmann, R., 2014. Visual cal image analysis 17, 19–42. tracking of da vinci instruments for laparoscopic surgery, [28] Milletari, F., Navab, N., Ahmadi, S.A., 2016. V-net: Fully in: Medical Imaging 2014: Image-Guided Procedures, convolutional neural networks for volumetric medical im- Robotic Interventions, and Modeling, International Soci- age segmentation, in: 3D Vision (3DV), 2016 Fourth In- ety for Optics and Photonics. p. 903608. ternational Conference on, IEEE. pp. 565–571. [40] Takimura, H., Muramatsu, T., Tsukahara, R., Nakano, M., [29] Pauly, O., Heibel, H., Navab, N., 2010. A machine learn- Nishio, S., Takimura, Y., Yabe, T., Kawano, M., Hada, T., ing approach for deformable guide-wire tracking in fluo- 2018. Usefulness of novel roadmap guided percutaneous roscopic sequences, in: International Conference on Med- coronary intervention for bifurcation lesions (the world’s ical Image Computing and Computer-Assisted Interven- first cases). Journal of the American College of Cardiol- tion, Springer. pp. 343–350. ogy 71, A1365. [30] Peressutti, D., Penney, G.P., Housden, R.J., Kolbitsch, C., [41] Tehrani, S., Laing, C., Yellon, D.M., Hausenloy, D.J., Gomez, A., Rijkhorst, E.J., Barratt, D.C., Rhode, K.S., 2013. Contrast-induced acute kidney injury following pci. King, A.P., 2013. A novel Bayesian respiratory mo- European journal of clinical investigation 43, 483–490. tion model to estimate and resolve uncertainty in image- [42] Timinger, H., Krueger, S., Dietmayer, K., Borgert, J., guided cardiac interventions. Medical image analysis 17, 2005. Motion compensated coronary interventional navi- 488–502. gation by means of diaphragm tracking and elastic motion [31] Petkovic, ´ T., Loncari ˇ c, ´ S., 2010. Guidewire tracking with models. Physics in Medicine & Biology 50, 491. projected thickness estimation, in: Biomedical Imaging: [43] Vandini, A., Glocker, B., Hamady, M., Yang, G.Z., 2017. From Nano to Macro, 2010 IEEE International Sympo- Robust guidewire tracking under large deformations com- sium on, IEEE. pp. 1253–1256. bining segment-like features (seglets). Medical image [32] Piayda, K., Kleinebrecht, L., Afzal, S., Bullens, R., ter analysis 38, 150–164. Horst, I., Polzin, A., Veulemans, V., Dannenberg, L., [44] van Walsum, T., Schaap, M., Metz, C.T., van der Giessen, Wimmer, A.C., Jung, C., et al., 2018. Dynamic coronary A.G., Niessen, W.J., 2008. Averaging centerlines: mean roadmapping during percutaneous coronary intervention: shift on paths, in: International Conference on Medical a feasibility study. European journal of medical research Image Computing and Computer-Assisted Intervention, 23, 36. Springer. pp. 900–907. [33] Rieke, N., Tan, D.J., di San Filippo, C.A., Tombari, F., [45] Wang, P., Chen, T., Zhu, Y., Zhang, W., Zhou, S.K., Co- Alsheakhali, M., Belagiannis, V., Eslami, A., Navab, N., maniciu, D., 2009. Robust guidewire tracking in fluo- 2016. Real-time localization of articulated surgical instru- roscopy, in: 2009 IEEE Conference on Computer Vision ments in retinal microsurgery. Medical image analysis 34, and Pattern Recognition, IEEE. pp. 691–698. 82–100. [46] Wu, W., Chen, T., Strobel, N., Comaniciu, D., 2012. Fast [34] Ronneberger, O., Fischer, P., Brox, T., 2015. U-net: tracking of catheters in 2D fluoroscopic images using an Convolutional networks for biomedical image segmenta- integrated CPU-GPU framework, in: Biomedical Imag- tion, in: International Conference on Medical image com- ing (ISBI), 2012 9th IEEE International Symposium on, puting and computer-assisted intervention, Springer. pp. IEEE. pp. 1184–1187. 234–241. [47] Wu, X., Housden, J., Ma, Y., Razavi, B., Rhode, K., [35] Schneider, M., Sundar, H., Liao, R., Hornegger, J., Xu, C., Rueckert, D., 2015. Fast catheter segmentation from 2010. Model-based respiratory motion compensation for echocardiographic sequences based on segmentation from 25 corresponding X-ray fluoroscopy for cardiac catheteriza- B. Parameter Tuning for Catheter Tip Tracking tion interventions. IEEE transactions on medical imaging This section describes the details of tuning the pa- 34, 861–876. rameters of optical flow and particle filter for catheter [48] Yabe, T., Muramatsu, T., Tsukahara, R., Nakano, M., tip tracking. Takimura, Y., Takimura, H., Kawano, M., Hada, T., 2018. The impact of percutaneous coronary intervention using the novel dynamic coronary roadmap system. Journal of B.1. Tuning Optical Flow Parameters the American College of Cardiology 71, A1103. The approach of [14] was used as the optical flow [49] Yatziv, L., Chartouni, M., Datta, S., Sapiro, G., 2012. To- implementation in Algorithm 1. A grid search to find ward multiple catheters detection in fluoroscopic image the optimal parameter setting was done on the fol- guided interventions. IEEE Transactions on Information Technology in Biomedicine 16, 770–781. lowing parameters of the method: (1) the image scale [50] Zhu, Y., Tsin, Y., Sundar, H., Sauer, F., 2010. Image- to build the pyramids, (2) the number pyramid levels, based respiratory motion compensation for fluoroscopic (3) the averaging window size, (4) the number of it- coronary roadmapping, in: International Conference on erations, (5) the size of the pixel neighborhood used Medical Image Computing and Computer-Assisted Inter- to find polynomial expansion in each pixel, and fi- vention, Springer. pp. 287–294. nally (6) the standard deviation of the Gaussian that is used to smooth derivatives used as a basis for the Appendix polynomial expansion. The above parameters were tuned independently A. Details of the Training Setup of the deep neural network, as optical flow di- A.1. Data Augmentation rectly estimates the catheter tip motion between two To increase the number of training samples and frames. To tune the parameters, we tracked the their diversity, data augmentation was used. The aug- catheter tip in X-ray fluoroscopy starting from the mentation includes geometric transformation such as ground truth tip position in the first frame using the flipping (left-right, up-down), rotation of multiple of motion field between two adjacent frames estimated 90 degrees, random ane transformation (translation with optical flow. The average and median distance -10 to 10 px, scaling 0.9 to 1.1, rotation -5 to 5 de- between the tracked tip position and the ground truth grees, shear -5 to 5 px), random elastic deformation were used as the evaluation criteria for the tuning. (deformation range -4 to 4 px, grid size of control The method of [14] was implemented by using the points 64 px). A training sample has 0.5 probability OpenCV function calcOpticalFlowFarneback. of being processed with one of the transformations. With consideration of the suggested parameter val- The probability for applying each transformation is: ues from the documentation, the parameter setting flipping left-right (1/24), flipping up-down (1/24), chosen for optical flow from the grid search is rotation of multiple of 90 degrees (1/12), ane trans- pyr scale = 0:5, levels = 3, winsize = 10, formation (1/6), elastic deformation (1/6), no trans- iterations = 30, poly n = 5, poly sigma = 1:1. formation (1/2). To make the trained model robust to Details of the parameters can be found on the func- noise, in addition to the geometric transformations, tion documentation page . we also augmented data by adding Gaussian noise to the pixel value with a zero mean and a standard de- B.2. Tuning Particle Filter Parameters viation between 0.01 and 0.03. The probability of The parameters to tune for the particle filter are adding the noise is 0.5. the number of samples N and the variance of pro- cess noise  . When tuning them, we fixed the pa- A.2. Training Settings v rameters of the trained network and the optical flow The  value in the training loss Eq. (5) was set method, and used their optimal parameter settings to 10 to make the scale of the two terms similar. during this experiment. Following Algorithm 1, we Adam optimizer was used to minimize the loss func- tion with a learning rate 0.0001. The number of train- ing samples in a batch is 4. The network was trained https://docs.opencv.org/2.4/modules/video/ with 100 epochs to ensure convergence. doc/motion_analysis_and_object_tracking.html? 26 tracked the catheter tip from the ground truth posi- Detection 6 Detection 120 Tracking Tracking tion (probability map) in the first frame, and used the mean and median distance between the tracked and the true position as the validation metric. The tracking results on the validation (tracking) 20 1 set are shown in Table 10. The table shows that 100 With Segmentation Without Segmentation With Segmentation Without Segmentation samples are suboptimal, while 1000 samples seem (a) Overall view of tracking er- (b) A zoom-in view of (a) sucient, as 10000 samples result in tracking accu- rors racies similar to 1000 samples. It also shows that the optimal choices of the standard deviation of the Figure 16: Tracking errors on all test images with and without process noise are 4 or 5 px for the downsampled im- catheter segmentation. ages. One possible reason for such choices may be that they are similar to the size of guiding catheters. C.1. Catheter tip detection and tracking In general, good choices for N are 1000 and 10000, The same metrics in Table 5 are used to report for  are 4 and 5. By considering the mean, the the accuracy of catheter tip detection and tracking standard deviation and the median of tracking errors, without catheter segmentation. Table 11 and Fig- the parameter setting  = 5, N = 1000 was chosen. v s ure 16 both manifest that the segmentation sub-task improves the accuracy of catheter tip detection and tracking. Although the improvement on the track- C. Detection, tracking and roadmapping without ing task is marginal and not statistically significant catheter segmentation (p = 0:06), the segmentation helps to reduce the magnitude and amount of outliers (large errors). Training of the network in Figure 3 requires catheter labels for detection and segmentation. As C.2. Dynamic coronary roadmapping the segmentation labels are often more expensive to In this experiment, the same setup in Section acquire than the detection label in practice, we also 6.3 was used to assess the accuracy of DCR using investigated the performance of catheter tip detec- catheter tip tracking without segmenting the catheter. tion, tracking and dynamic coronary roadmapping Table 12 indicate that tracking the catheter tip with without segmenting the catheter. To this end, we catheter segmentation improves the DCR accuracy used a similar encoder-decoder network architecture compared to that without catheter segmentation. Al- as Figure 3 which computes only the detection output though the improvement is not statistically signifi- directly after the last up block of the decoder with a cant (p = 0:43), the approach with segmentation is 1 1 convolution followed by a spatial softmax layer, more robust by making less large roadmapping er- instead of having two outputs. We then followed rors (Figure 17). the same way as the approach using the network with segmentation to search for (hyper-)parameters for the approach without segmentation. The follow- ing parameter setting was chosen for the experiments in this section: for deep learning, the basic channel number is 64, the depth is 5, the dropout rate is zero; for particle filtering,  = 3, N = 10000. With this v s setup, we compared the performance of the approach without catheter segmentation to the proposed ap- proach with segmentation on catheter tip detection and tracking (Appendix C.1) and dynamic coronary roadmapping (Appendix C.2) on the test set from Ta- ble 2. Ground truth / Prediction Distance (mm) Ground truth / Prediction Distance (mm) Table 10: Catheter tip tracking errors (mm) on the validation (tracking) dataset of di erent parameter settings for particle filter. The tracked tip point was rounded to the pixel center. The error of all images (mean  std / median) are presented. Red cells show the good choices of parameters; bold number indicates the chosen setting. v s (px) 100 1000 10000 3 1.52  2.19 / 0.79 1.49  2.18 / 0.79 1.48  2.18 / 0.79 4 1.50  2.17 / 0.79 1.46  2.17 / 0.79 1.47  2.18 / 0.79 5 1.52  2.21 / 0.79 1.47  2.17 / 0.74 1.47  2.19 / 0.74 6 1.53  2.39 / 0.79 1.49  2.33 / 0.79 1.48  2.29 / 0.74 7 1.56  2.42 / 0.79 1.50  2.29 / 0.74 1.50  2.39 / 0.74 8 1.58  2.41 / 0.79 1.51  2.40 / 0.74 1.51  2.42 / 0.74 9 1.56  2.22 / 0.79 1.53  2.43 / 0.79 1.52  2.45 / 0.61 10 2.25  6.18 / 0.79 1.54  2.46 / 0.79 1.53  2.47 / 0.61 Table 11: Catheter tip tracking errors (mm) with and without catheter segmentation on the test (tracking) dataset. y indicates that the di erence between the detection with and without segmentation is statistically highly significant with the two-sided Wilcoxon signed-rank test (p < 0:001). No statistically significantly di erence is observed between the tracking with and without segmenta- tion using the two-sided Wilcoxon signed-rank test (p = 0:06). With Segmentation Without Segmentation Evaluation Metrics Detectiony Tracking Detection Tracking Maximal error of all images 108.20 17.72 133.94 23.2 Median error of all images 0.96 0.96 0.96 0.96 Mean error of all images 5.62  15.91 1.29  1.76 9.32  19.73 1.75  3.17 Average of sequence median error 6.26  17.11 1.03  0.49 9.34  18.64 1.42  2.14 Average of sequence mean error 6.83  13.88 1.29  0.94 10.41  15.94 1.69  1.97 Figure 17: Accuracy (mm) of DCR via catheter tip tracking with and without catheter segmentation. 28 Table 12: The statistics of DCR accuracy (mm) via catheter tip tracking with and without catheter segmentation. With the two-sided Wilcoxon signed-rank test, no statistically significantly di erence is observed between the DCR with and without segmentation (p = 0:43). With Segmentation Without Segmentation All point pairs Maximal distance 20.24 25.20 Median distance 1.43 1.43 Mean distance 2.07  2.08 2.44  3.15 Frame mean distance Median distance 1.42 1.40 Average distance 1.91  1.52 2.23  2.59 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Electrical Engineering and Systems Science arXiv (Cornell University)

Dynamic Coronary Roadmapping via Catheter Tip Tracking in X-ray Fluoroscopy with Deep Learning Based Bayesian Filtering

Loading next page...
 
/lp/arxiv-cornell-university/dynamic-coronary-roadmapping-via-catheter-tip-tracking-in-x-ray-rMBuDmqOsh

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

ISSN
1361-8415
eISSN
ARCH-3348
DOI
10.1016/j.media.2020.101634
Publisher site
See Article on Publisher Site

Abstract

Percutaneous coronary intervention (PCI) is typically performed with image guidance using X-ray an- giograms in which coronary arteries are opacified with X-ray opaque contrast agents. Interventional cardi- ologists typically navigate instruments using non-contrast-enhanced fluoroscopic images, since higher use of contrast agents increases the risk of kidney failure. When using fluoroscopic images, the interventional cardiologist needs to rely on a mental anatomical reconstruction. This paper reports on the development of a novel dynamic coronary roadmapping approach for improving visual feedback and reducing contrast use during PCI. The approach compensates cardiac and respiratory induced vessel motion by ECG alignment and catheter tip tracking in X-ray fluoroscopy, respectively. In particular, for accurate and robust tracking of the catheter tip, we proposed a new deep learning based Bayesian filtering method that integrates the detec- tion outcome of a convolutional neural network and the motion estimation between frames using a particle filtering framework. The proposed roadmapping and tracking approaches were validated on clinical X-ray images, achieving accurate performance on both catheter tip tracking and dynamic coronary roadmapping ex- periments. In addition, our approach runs in real-time on a computer with a single GPU and has the potential to be integrated into the clinical workflow of PCI procedures, providing cardiologists with visual guidance during interventions without the need of extra use of contrast agent. Keywords: dynamic coronary roadmapping, X-ray fluoroscopy, catheter tip tracking, deep learning, Bayesian filtering, particle filter. 1. Introduction A guiding catheter is firstly positioned into the os- tium of the coronary artery. Through the guiding 1.1. Clinical Background catheter, a balloon catheter carrying a stent is intro- Percutaneous coronary intervention (PCI) is a duced over a guidewire to the stenosed location. The minimally invasive procedure for treating patients balloon is then inflated and the stent is deployed to with coronary artery disease. During these proce- prevent the vessel from collapsing and restenosing. dures, medical instruments inserted through a guid- ing catheter are advanced to treat coronary stenoses. PCI is typically performed with image-guidance using X-ray angiography (XA). Coronary arteries are This manuscript has been accepted for publication in Med- visualized with X-ray opaque contrast agent. During ical Image Analysis. For the published journal article, please refer to DOI: https://doi.org/10.1016/j.media.2020.101634. the procedure, interventional cardiologists may re- yy c 2020. This manuscript version is made available un- peatedly inject contrast agent to visualize the vessels, der the CC-BY-NC-ND 4.0 license. For license details, please as the opacification of coronary arteries only lasts for check via http://creativecommons.org/licenses/by-nc-nd/4.0/. a short period. The amount of periprocedural con- Corresponding author Email address: huama.research@gmail.com (Hua Ma) trast use has been correlated to operator experience, Preprint submitted to Medical Image Analysis January 14, 2020 arXiv:2001.03801v1 [eess.IV] 11 Jan 2020 procedural complexity, renal function and imaging with the live fluoroscopic sequence by aligning the setup ([32]). Furthermore, the risk for contrast in- R waves of their corresponding ECG signals. This duced nephropathy has been associated to contrast system compensates the cardiac motion of vessels, volume ([41]). Manoeuvring guidewires and mate- yet does not correct the respiratory motion during in- rial, however, typically occurs without continuous terventions. Two later studies by [50] and [26] in- contrast injections. In these situations, the navigation troduced image-based respiratory motion compensa- of devices is guided with ”vessel-free” fluoroscopic tion methods for DCR. Their methods assumed an images. Cardiologists have to mentally reconstruct ane respiratory motion model in ECG-gated fluo- the position of vessels and stenosis based on previ- roscopic frames and recovered the respiratory mo- ous angiograms. tion from soft tissues with special handling of static structures. The limitation of these approaches is that 1.2. Dynamic Coronary Roadmapping they require relevant tissue to be suciently visi- ble in the field of view for reliable motion compen- Dynamic coronary roadmapping (DCR) is a sation which is not always the case. In addition, promising solution towards improving visual feed- they require to be run on cardiac-gated frames. In back and reducing usage of contrast medium during a more recent work by [20], binary vessel masks PCI ([12, 50, 26, 20]). This approach dynamically were created as the roadmaps from at least one car- superimposes images or models of coronary arter- diac cycle of angiographic images. Temporal align- ies onto live X-ray fluoroscopic sequences. The dy- ment of the roadmaps and the fluoroscopic sequence, namic overlay serves as a roadmap that provides im- which compensated the cardiac motion of vessels, mediate feedback to cardiologists during the inter- was performed by registering ECG signals using vention, so as to assist in navigating a guidewire into cross-correlation. Additionally, the respiratory mo- the appropriate coronary branch and proper place- tion was corrected by aligning the guidewire center- ment of a stent at the stenotic site with reduced ap- line in the fluoroscopy to the contour of vessels in the plication of contrast agent. Studies with a phantom angiogram where the roadmaps were created. The setup using research software ([20]) or on first cases system has been shown useful in a phantom-based of clinical interventions using commercially avail- study, nevertheless no accuracy evaluation of the spa- able systems ([9, 48, 40]) have investigated the use- tiotemporal alignment was presented. Furthermore, fulness of DCR in assisting PCI, reporting that DCR the spatial registration relies on robust extraction of helps to reduce procedure time, radiation dose and vessels and guidewires which is often challenging for contrast volume. X-ray images. To develop a DCR system, it is important but yet a challenge to accurately overlay a roadmap of coro- Unlike direct roadmapping, the model-based ap- nary arteries onto an X-ray fluoroscopic image, as proaches build a model to predict motion in fluoro- limited information of vessels is present in the target scopic frames. The motion models are often func- fluoroscopic image for inferring the compensation of tions that relate the motion of roadmaps to surro- the vessel motion resulting from patient respiration gate signals derived from images or ECG, so that and heartbeat. The methods that have been proposed once the surrogates for fluoroscopic frames are ob- on motion compensation for DCR can be generally tained, the motion can be computed by the model. grouped into two categories: direct roadmapping and For cardiac interventions including PCI, the organ model-based approaches. motion is mainly a ected by respiratory and car- Direct roadmapping methods use information diac motion. Many previous works often built a mo- from X-ray images and ECG signals to directly cor- tion model parameterized by a cardiac signal derived rect the motion caused by respiration and heartbeat. from ECG and a respiratory signal obtained from di- The first DCR system ([12]) used digital subtrac- aphragm tracking ([37, 42, 13]) or automatic PCA- tion of a contrast sequence and a mask sequence based surrogate ([15]). Some other works model to create a full cardiac cycle of coronary roadmaps. only the respiratory motion in cardiac-gated images The roadmaps were stored and later synchronized ([35, 21, 30]). For a complete review on respira- 2 tory motion modeling, we refer readers to the sur- fitting the B-spline to measurements based on gray vey article by [27]. One limitation of the model- intensity and vesselness image. [5] proposed a con- based approaches is that the motion models are often volutional neural network (CNN) to detect catheter patient-specific, which requires training the model electrodes in X-ray images, which treated catheter every time for a new subject. Additionally, once detection as a segmentation problem. The method the surrogate values during inference are out of the used a weighted cross-entropy loss to cope with the surrogate range for building the model, e.g. for ab- class imbalancing problem due to the small size of normal motion, extrapolation is needed, which may the target. [22] and [11] tracked surgical instruments hamper accurate motion compensation. using a deep network having an encoder-decoder ar- chitecture. Their approaches combined instrument 1.3. Interventional / Surgical Tool Tracking segmentation and detection in a multi-task learning problem to make the tool detection in a cluttered Tracking interventional tools is relevant for mo- background more robust. tion compensation ([35, 7, 25, 4, 2]). In particular for Di erent from tracking by detection, which relies PCI, the guiding catheter tip typically remains within solely on the current image, temporal tracking also the coronary ostium which is in the field of view uses information from previous frames. The tempo- during the largest part of the intervention, making ral information can reduce the search space for de- it a suitable landmark for tracking. [4] have shown tection, or put additional constraints in the model, that catheter tip motion during PCI can be modeled making the tracking more robust. as a combination of cardiac and respiratory motion. As using catheter tip displacement can only correct Temporal information has been used in various translational motion, [4] further showed that, com- ways. Some methods mainly relied on a detection pared to a rigid motion model for the respiratory mo- model, but incorporate temporal information in the tion, modeling only the translational part of the res- preprocessing ([7]) or post-processing ([16]) step or piratory motion deteriorated the accuracy marginally, in the input ([33, 1]). Approaches based on back- which confirms the observations by [36] that the rota- ground estimation have been used for catheter ([49]) tional part of respiratory motion is small. These find- or guidewire ([31]) tracking. In these approaches, the ings motivate motion compensation for DCR through background was recursively updated for every frame, tracking the catheter tip in X-ray fluoroscopic se- and was used for enhancing the foreground contain- quences. ing instruments. Apart from those, many works Many works have been proposed to address the adopted a Bayesian framework for tracking instru- problem of tracking interventional or surgical tools ments via a maximum a posteriori (MAP) formu- in medical images for various applications. The lation. Representations based on key points ([47]), tracking methods from these works can be generally B-splines ([45, 29, 19, 18]), or segment-like fea- categorized into two kinds of approaches: tracking tures ([43]) have been used to model catheters or by detection, and temporal tracking. guidewires. Markov random field (MRF) was used The tracking by detection approaches treat track- to model the position or deformation of the con- ing as a detection problem, which rely on features trol points in the B-spline ([29, 19, 18, 47]). In only from the current image without using informa- the work by [43], temporal information was incor- tion from previous frames. For example in electro- porated in the prior term using Kalman filter. Partic- physiology procedure, as the catheters present spe- ularly, learning-based approaches were used in sev- cific features in shape or intensity, ad hoc meth- eral works to obtain the likelihood for a more ro- ods were proposed based on, e.g. blob detection, bust measurement using probabilistic boosting tree shape constrained searching and model-/template- ([45, 46]) or support vector regression ([29]). In ad- based detection ([25, 24]). [8] modeled the catheter dition, temporal tracking models based on Bayesian tracking problem by optimizing the posterior in a filtering were also a popular approach for instrument Bayesian framework, in which the catheter was rep- tracking. [2] used a hidden Markov model (HMM) resented as a B-spline tube model and was tracked by to track catheter tip in a 3D vessel tree, for which 3 the likelihood was obtained based on the 3D-2D reg- 3. We evaluate the proposed approach visually istration outcome. [38] used particle filters to track and quantitatively on clinical X-ray sequences, surgical tools in medical images. They used a likeli- achieving low errors on both tracking and hood based on the segmentation of instruments, and roadmapping tasks. a dynamic model that incorporates samples from two 4. The proposed DCR method runs in real-time previous time steps. In a later work, [39] used a with a modern GPU, thus can potentially be multi-object particle filter to track multiple instru- used during PCI in real clinical settings. ment regions simultaneously, in which a particle is the concatenation of the states of several objects. 2. Scenario Setup and Method Overview Despite of many existing works on inverventional The proposed method assumes that the scenario of or surgical tool tracking in medical images, an au- performing dynamic coronary roadmapping to guide tomatic approach for tracking the tip of guiding a PCI procedure consists of an oine phase and an catheter in X-ray fluoroscopy for PCI has not been online phase. An overview of the method is shown investigated yet. The challenges of this task are: in Figure 1. (1) di erent from the catheters for EP that can be viewed as blobs or a circle, the guiding catheter for 2.1. Oine Phase PCI presents a dark tubular appearance which shows This phase (Step 0 in Figure 1) is performed o - no prominent features; (2) the shape of the guiding line before the actual roadmapping is conducted. In catheter tip segment varies depending on the orienta- this stage, roadmaps of coronary arteries contain- tion of the C-arm, making feature-/model- based de- ing multiple cardiac phases are created from an X- tection challenging; (3) the background may contain ray angiography sequence acquired with injection of structures that have similar appearance to a catheter contrast agent. A roadmap can be a vessel model tip, such as vertebral structures or residual contrast in the form of centerlines, contours, masks, etc. It agent, which makes robust tracking dicult. may also contain information of clinical interest, e.g. stenosis. Since the main focus of this paper is on 1.4. Contributions accurate overlay of a roadmap, we do not investigate We propose and evaluate a novel approach for dy- how to create the most suitable roadmaps, but use the namic coronary roadmapping. The approach com- images containing only vessels and catheters that are pensates changes in vessel shapes and cardiac motion created using the layer separation method by [23] as by selecting the roadmap of the same cardiac phase the roadmaps to show the concept of dynamic coro- through ECG alignment, and corrects the respiratory nary roadmapping. Along with the XA sequence, induced motion via tracking the tip of the guiding ECG signals are also acquired and stored for later catheter. Our contributions are: selecting a roadmap that has similar cardiac phase to 1. We develop a new way to perform dynamic a given X-ray fluoroscopic frame in the online phase coronary roadmapping on free breathing, non- (see details in Section 3). Once the image sequence cardiac-gated X-ray fluoroscopic sequences. and ECG signals are acquired, the catheter tip loca- Particularly, the respiratory-induced vessel mo- tion in every frame is obtained to serve as a refer- tion is robustly compensated via the displace- ence point for roadmap transformation. In this work ment of catheter tip. we manually annotated the catheter tip in the oine 2. We proposed a deep learning based method XA sequence. In real clinical scenarios, the annota- within a Bayesian filtering framework for online tion work can be done by the clinician or a person detection and tracking of guiding catheter tip in who assists the intervention, such as a technician or X-ray fluoroscopic images. The method models a nurse. the likelihood term of Bayesian filtering with a 2.2. Online Phase convolutional neural network, and integrates it with particle filtering in a comprehensive man- This is when the dynamic roadmapping is actu- ner, leading to more robust tracking. ally performed. In this phase, non-contrast X-ray 4 Offline Phase Online Phase A fluoroscopic frame ECG image Manual annotation Angiographic sequence Tips Tips Track tip images image Tip Tip image Roadmaps Roadmaps ECG Deep learning based Translation Bayesian filtering Correlation Layer separation score Tip Tip Transform ECG ECG Roadmap Roadmap Roadmap STEP 0 STEP 1: Select Roadmap STEP 2: Track Catheter Tip STEP 3: Transform Roadmap Section 2.1 Section 3 Section 4 Section 2.2 Figure 1: The overview of the proposed dynamic coronary roadmapping method. The colored blocks with a dash line border denote objects acquired in the online phase; the colored blocks with a solid line border are objects originated from the oine phase. fluoroscopic images with the same view angles as oroscopic image, which compensates the di erence the roadmaps created during the oine phase are ac- of vessel shape and pose induced by cardiac motoin. quired sequentially. At the same time, ECG signals An approach similar to the ECG matching method by along with the roadmapping frames are also obtained [20] is used to accomplish this task. and are compared with the stored ECG to select the To select roadmaps images based on ECG, a tem- most matched roadmap (Step 1 in Figure 1; see de- poral mapping between X-ray images and ECG sig- tails in Section 3). This is to compensate the change nal points needs to be built first. We assume that of vessel shape and position between frames due ECG signals and X-ray images are well synchronized to cardiac motion. Simultaneously, the catheter tip during acqusition. In the oine phase, the begin- location in the acquired X-ray fluoroscopic images ning and the end of the image sequence are aligned is tracked online using the proposed deep learning with the start and end ECG signal points; the XA based Bayesian filtering method in Section 4 (Step frames in between are then evenly distributed on the 2 in Figure 1). The displacement of catheter tip be- timeline of ECG. This way, a mapping between the tween the current image and the selected roadmap stored sequence images and its ECG signal can be image is then obtained and are applied to transform set up: for each image, the closest ECG signal point the roadmap. Finally, the transformed roadmap is to the location of the image on the timeline can be overlaid on the current non-contrast frame to guide found; for each ECG point, an image that is closest the procedure (Step 3 in Figure 1). to this point on the timeline can be similarly located. Once the mapping is available, all images with good vessel contrast filling and the ECG points that are as- 3. ECG matching for Roadmap Selection sociated to these images are selected from the XA Roadmap selection in this work is achieved by sequence for the pool of roadmaps. In this process, comparing the ECG signal associated with the flu- at least one heartbeat of frames should be acquired, oroscopic image and the ECG of the angiographic which is generally the case in our data. In the on- sequence, such that the most suitable candidate line phase, similar to the approach of [20], for ac- roadmap is selected where the best match of the ECG quisition of each image, a block of N latest ECG ECG signals is found. The selected roadmap has the same signal points is constantly stored and updated in the (or very similar) cardiac phase with the X-ray flu- history bu er. This is considered as the ECG signal Select roadmap Overlay corresponding to the fluoroscopic frame. the likelihood p(z jx ) from which x and z can be k k k k obtained by sampling. To compare the ECG signals associated with the angiographic sequence and the online fluoroscopic With these definitions and p(x ), the inital be- lief of x , Bayesian filtering seeks an estimation of image, a temporal registration of the two signals us- ing cross-correlation is applied ([20]). The two ECG x (k  1) based on the set of all available observa- tions z = fz ; i = 0; :::; kg up to time k via recur- signals are first cross-correlated for every possible 0:k i sively computing the posterior probability p(x jz ) position on the signals, resulting in a 1D vector of k 0:k as Eq.(1) ([3]): correlation scores. The candidate frame for dynamic overlay is then selected as the one associated with the point on the ECG of the angiographic sequence that p(x jz ) / p(z jx ) p(x jx )p(x jz )dx : k 0:k k k k k1 k1 0:k1 k1 is corresponding to the highest correlation score. | {z } p(x jz ) k 0:k1 (1) 4. Bayesian Filtering for Catheter Tip Tracking Assuming the initial probability p(x jz ) = p(x ) is 0 0 0 known, based on Eq.(1), Bayesian filtering runs in Bayesian filtering is a state-space approach aiming cycles of two steps: prediction and update. In the at estimating the true state of a system that changes prediction step, the prior probability p(x jz ), the k 0:k1 over time from a sequence of noisy measurement initial belief of x given previous observations, is made on the system ([3]). One popular application estimated by computing the integral in Eq.(1). In area of this approach is tracking objects in a series of the update step, the prior probability is corrected by images. the current likelihood p(z jx ) to obtain the posterior k k p(x jz ). k 0:k 4.1. Theory of Bayesian Filtering In Section 4.2, we will firstly introduce how to Bayesian filtering typically includes the following model the likelihood. Then in Section 4.3, a way components: hidden system states, a state transition to represent and eciently approximate the posterior model, observations and a observation model. Let will be discussed. Finally in Section 4.4, a summary x 2 R (k = f0; 1; 2; :::g) denote the state, the loca- of the complete catheter tip tracking method will be tion of guiding catheter tip in the k-th frame, a 2D given. vector representing the coordinates in the X-ray im- 4.2. A Deep Learning based Likelihood age space. The transition of the system from one state to the next state is given by the state transition Directly modeling the likelihood p(z jx ) is chal- k k model x = f (x ; v ), where v 2 R is an in- lenging due to (1) the complexity of the generation k k k1 k1 k1 dependent and identically distributed (i.i.d.) process process h and (2) the computational complexity of 2 2 2 2 noise, f : R  R ! R is a possibly nonlinear p(z jx ) for every value x 2 R . In this work, we k k k k function that maps the previous state x to the cur- simplify the problem by only computing the likeli- k1 rent state x with noise v . The observation z in hood p(z jx ) in the image pixel space, i.e. the in- k k1 k k k this work is defined as the k-th X-ray image of a se- teger pixel coordinate. For a subpixel x , the value wh quence, so that z 2 R , where w and h are the of p(z jx ) can possibly be approximated by inter- k k k width and height of an X-ray image. We further de- polation. To this end, we propose to use a deep fine the observation model as z = h (x ; n ), where neural network D to approximate p(z jx ) for inte- k k k k k k n 2 R is an i.i.d measurement noise (n is the di- ger pixel locations. The network takes an image z k k k 2 n wh mension of n ), h : R R ! R is a highly non- as input and outputs a probability of observing the k k linear function that generates the observation z from input z for every pixel location x . Therefore, the k k k the state x with noise n . The state transition model approximated likelihood is a function of x , denoted k k k f and the observation model h , respectively, can as D (x ). Since x is defined within the scope of k k z k k also be equivalently represented using probabilistic the image pixel space, D (x ) is essentially a proba- z k forms, i.e. the state transition prior p(x jx ) and bility map having the same dimension and size with k k1 6 [28], residual blocks ([17]) are adopted at each reso- lution level in the encoder and decoder to ease gradi- ent propagation in a deep network. The encoder con- sists of 4 down blocks in which a residual block fol- lowed by a stride-2 convolution is used for extraction and down-scaling of feature maps. The number of (a) (b) (c) feature maps is doubled in each downsampling step. The decoder has 4 up blocks where a transposed con- Figure 2: Input and ground truth labels for the deep neural net- volution of stride-2 is used for upsampling of the in- work: (a) an input X-ray fluoroscopic image, (b) the binary catheter mask of (a) for catheter segmentation, (c) a 2D Gaus- put feature maps. Dropout is used in the residual unit sian PDF ( = 4 px) for likelihood estimation for (a). of the up block for regularization of the network. Be- tween the encoder and the decoder, another residual block is used to process the feature maps extracted the input image z , in which the entry at each loca- by the encoder. The detailed network architecture is tion x ( j = 1; 2; : : : ; wh) in the map represents the shown in Figure 3. probability of observing z given x . It is worth men- Due to similar appearance between a guiding tioning that the deep neural network is used for ap- catheter tip and corners of a background structure, proximation of p(z jx ), which should be clearly dis- k k such as vertebral bones, lung tissue, stitches or tinguished from the generation model h that maps guidewires, ambiguity may exist when the network an x to z . The existence of h is merely for the con- k k k is expected to output only one blob in the probability venience of definition, its explicit form, however, is map. To alleviate the issue, we adopt a similar strat- not required in the context of this work. egy as [22], using a catheter mask (Figure 2b) as an To obtain the training labels, we assume that there additional label to jointly train the network to output exists a mapping h , such that the training label can both the catheter segmentation heatmap and the like- be defined as a distance-based probability map, i.e. lihood probability map. The segmentation heatmap the farther away x is from the ground truth tip loca- is obtained by applying a 1  1 convolution with tion in the image z , the less possible it is to observe ReLU activation on the feature maps of the last up z given x through the process h . This definition block. To compute the likelihood probability map, a k k k matches the intuition that from a location x that is residual block is firstly applied on the feature maps far from the ground truth tip location, the probability of the last up block. The output feature maps are of observing a z with the catheter tip being located then concatenated with the segmentation heatmap as at the ground truth position should be low. For sim- one additional channel, followed by a 1  1 convo- plicity, a 2D Gaussian probability density function lution. Finally, to ensure the network detection out- 0 2 (PDF) N (x ; x ;  I) centered at the ground truth tip put fits the definition of a probability map on image 0 2 location x with variance  I in the image space is locations, following the 1  1 convolution, a spatial used as the label to train the network (Figure 2c). softmax layer is computed as Eq.(2): Note that this training label makes the estimation of k;l p(z jx ) equivalent to a catheter tip detection problem k k D = P ; (2) k;l i; j such that the deep neural network learns features of i; j catheter tip and outputs high probability at locations where A is the output feature map of the 1  1 con- where the features are present. Due to this reason, volution, A denotes the value of A at location (i; j), i; j we also call p(z jx ) “detection output” or “detec- k k D is the final output of the detection network, a 2D tion probability” and call the estimation of p(z jx ) k k probability map representing p(z jx ). The details k k “catheter tip detection” in the context of this paper. are shown in Figure 3. The network that we use follows a encoder- The training loss is defined as a combination of the decoder architecture with skip connections similar to segmentation loss and the detection loss. The seg- U-net ([34]). Additionally, similar to the work by mentation loss L in this work is a Dice loss defined 7 1 c c c 1 1x1 Conv + Relu down up S 2c 2c 3x3 Conv + c 1x1 Conv + 1 Bn + Relu Spatial down up Softmax 4c 4c residual D Concat down up 8c 8c Add down up 16c 3x3 Conv + 3x3 Conv Relu Bn + Relu + Bn residual residual Add Concat Ch / 2 Add 3x3 Conv + 3x3 Conv + 3x3 Conv + Relu 3x3 Conv + 3x3 Conv Relu Bn + Relu Bn + Dp + Bn + Dp 3x3 stride-2 Bn + Relu + Bn Relu Conv + Bn 3x3 stride-2 + Relu transposed Conv (Ch x 2) + Bn (Ch /2) down up Figure 3: A joint segmentation and detection network for catheter tip detection. This figure shows an example network with 4 levels of depth (the number of down or up blocks). Meaning of abbreviations: Conv, 2D convolution; Bn, batch normalization; Relu, ReLU activation; Dp, dropout; Concat, concatenation; Ch, number of channels; S, segmentation output; D, detection output. The number above an image or feature maps indicates the number of channels; the number of channels in the residual network in a block is shown above the block; c is the basic number of channels, the channel number in the first down block. The number next to a rectangle denotes the size of the image or feature maps. Red arrows indicate a change of number of channels. by Eq.(3): 4.3. Approximation of the Posterior with Particle Filter 2 M S i; j i; j i; j Once the deep neural network in Section 4.2 is L = 1 (3) P P 2 2 M + S trained, its weights are fixed during inference for i; j i; j i; j i; j computing the posterior p(x jz ) for new data. Ide- k 0:k where M denotes the ground truth binary catheter aly, the network detection output p(z jx ) should be a k k masks, S is the segmentation heatmap. The loss Gaussian PDF during inference, as it is trained with function for detection L is mean square error (MSE) labels of Gaussian PDFs. However, due to simi- given by Eq.(4): lar appearance of background structures or contrast residual, the detection output is unlikely to be a per- fect Gaussian (possibly non-Gaussian or having mul- L = jT D j (4) d i; j i; j w h tiple modes), which prevents the posterior p(x jz ) k 0:k iw; jh in Eq.(1) being solved with an analytical method. In where T denotes the ground truth PDF, w and h are practice, the posterior can be approximated using a the width and height of an image. The total training particle filter method ([3]). loss L is defined as Eq.(5): Particle filter methods approximate the posterior PDF by a set of N random samples with associated i i s weights fx ; w g ([3]). As N becomes very large, L = L + L (5) k k i=1 s d this discrete representation approaches the true pos- where  is a weight to balance L and L . terior. According to [3], the approximation of the s d 256 posterior p(x jz ) is given by Eq.(6): weights to be 1=N , so the number of e ective sam- k 0:k s ples which have actual contribution to approximate p(x jz ) is maximized ([3]). If the resampling is ap- k 0:k i i p(x jz )  w (x x ) (6) k 0:k k k k plied at every time step, the particle filter becomes a i=1 sampling importance resampling (SIR) filter, and the where () is the Dirac delta function. The weight weight update rule follows Eq.(9). w can be computed in a recursive manner as Eq.(7) i i once w is known ([3]): w / p(z jx ) (9) k1 k k k i i i The final decision on catheter tip location in frame k p(z jx )p(x jx ) k k k1 i i w / w (7) can then be computed as the expectation of x , x ˆ = k k1 i i k k q(x jx ; z ) k k1 x p(x jz )dx , which is in this case, the weighted k k 0:k k where q(x jx ; z ) is an importance density from sum of all samples: k k k1 which it should be possible to sample x easily. For simplicity, a good and convenient choice of i i x = w x : (10) i k k k the importance density is the prior p(x jx ) ([3]), k1 i=1 so that the weight update rule (7) becomes w / i i w p(z jx ). k 4.4. Summary k1 k A sample can be drawn from p(x jx ) in the fol- k1 The overall catheter tip tracking using a deep lowing way. First, a process noise sample v is k1 learning based Bayesian filtering method is summa- sampled from p (v ), the PDF of v ; then x is v k1 k1 rized in Algorithm 1. generated from x via the state transition model k1 i i i x = f (x ; v ). In this work, p (v ) is set to be a k v k1 k k1 k1 2 5. Experimental Setup GaussianN (0;  I). The choice of motion model for f is important for an accurate representation of the 5.1. Data true state transition prior p(x jx ). A random mo- k k1 Anonymized clinical imaging data were used for tion cannot characterize well the motion of catheter our experiments. The data were acquired with stan- tip in XA frames. In this work, we estimated the dard clinical protocol using Siemens AXIOM-Artis motion from adjacent frames using an optical flow system, and are from 55 patients who underwent a method, as this approach 1) takes into account of the PCI procedure at the Department of Cardiology at observation z , which results in a better guess of the Erasmus MC in Rotterdam, Netherlands. Out of catheter tip motion, and 2) enables estimation of a these data, we selected data from 37 patients which dense motion field where the motion of a sample x were acquired since the year 2014 to develop our can be eciently obtained. Therefore, f is defined method, and used the data from the other 18 patients as Eq.(8): acquired before the year 2013 for evaluation. The detailed information about the data is listed in Table x = x + u (x ) + v (8) k k1 k1 k1 k1 where u () is the motion from frame k1 to frame In order to evaluate the proposed roadmapping k1 k estimated with optical flow using the method by method, for which an o -line angiographic sequence [14], u (x ) is the motion from state x . is required for roadmap preparation and an on- k1 k1 k1 Once samples are drawn and their weights are line fluoroscopic sequence taken from the same C- updated, the so-called “resampling” of the samples arm position is needed for performing the actual should be performed to prevent the degenaracy prob- roadmapping (see Section 2), we selected the con- lem, where all but one sample will have negligi- trast frames from a real clinical sequence to simu- ble weight after a few iterations ([3]). The resam- late the o -line sequence, and chose the non-contrast pling step resamples the existing samples according frames from the same clinical sequence to simulate to their updated weights and then resets all sample the online sequence. The selected contrast sequence 9 Algorithm 1 Deep learning based Bayesian filtering for online tracking of catheter tip in X-ray fluoroscopy Require: fz ; : : : ; z g (sequentially observed frames), D (A trained network from Section 4.2), p(x ) (the 0 T 0 initial PDF),  (the variance of v ; k = 1; : : : ; T ), T (number of frames for tracking), N (number of k1 s samples) i i 1: Draw x  p(x ), set w = 1=N , 8 i = 1; : : : ; N 0 s s 0 0 2: for k = 1 to T do 3: Compute u from z to z using the optical flow method of [14] k1 k1 k 4: for i = 1 to N do i 2 5: Draw v  N (0;  I) k1 i i i 6: Compute the motion of x : u = u (x ) k1 k1 k1 k1 i i i i i i 7: Draw x  p(x jx ): x = x + u + v k k1 k k1 k1 k1 i i i 8: Update weight w = p(z jx ) = D (x ) k z k k k 9: end for i i s i 10: Normalize w w = w , 8 i = 1; : : : ; N k k i=1 k s i i 11: Prediciton in frame k: x ˆ = w x i=1 k k i i s i 12: Resample fx ; w g using the method of [3] (so all w are set to 1=N again) k k i=1 k 13: end for were ensured suciently long to cover at least one complete cardiac cycle. 5.2. Data Split for Catheter Tip Detection and Table 1: Basic information of the acquired X-ray image data Tracking for our experiments. The number in the parenthesis next to the To develop the catheter tip tracking method, 1086 pixel size indicates the possible image size. X-ray fluoroscopic images selected from 260 non- contrast sequences of 25 patients from the develop- Data Development Evaluation ment set were used for training the network from Fig- No. patients 37 18 ure 3; 404 images from 94 non-contrast sequences of No. sequences 354 34 another 12 patients from the development set were Frame rate (fps) 15 15 used as validation set for the network model and hy- Image size (px) 512  512 512  512 perparameter selection. In the training and validation 600  600 600  600 sets, 4-5 frames were randomly selected from each 776  776 776  776 sequence, which are not necessarily continuous. To 960  960 1024  1024 tune the parameters for tracking, 1583 images from 1024  1024 88 sequences out of the 94 from the same 12 patients Pixel size (mm) 0.108 (1024) 0.139 (1024) of the validation set were used (6 sequences were 0.139 (1024) 0.184 (600) not selected for this task due to very short sequence 0.184 (600) 0.184 (776) length not more than 5 frames). Finally, to evaluate 0.184 (776) 0.184 (1024) catheter tip tracking accuracy, 1355 images from 34 0.184 (960) 0.216 (512) non-contrast sequences of 18 patients from the eval- 0.184 (1024) 0.279 (512) uation set were used for testing. The frames selected 0.216 (512) for tracking from each sequence must be continuous; the number of selected frames for tracking might vary, depending on the number of the non-contrast frames in the sequences. Details of the datasets for training, validation and test are listed in Table 2. 10 Table 2: Dataset of training, validation and test for detection and tracking of catheter tip in X-ray fluoroscopic frames. Training Validation Validation Test (detection) (detection) (tracking) (tracking) No. patients 25 12 12 18 No. sequences 260 94 88 34 No. frames 1086 404 1583 1355 Continous frames? No No Yes Yes 5.3. Experimental Settings for Training the Deep e.g. coronary arteries in our case, is not directly vis- Network ible in the target image. One possible choice intro- duced by [50] is to use the guidewire as a surrogate This section describes the basic experimental set- of the target vessel centerline in non-contrast images, tings for training the deep neural network. Details of as guidewire is always inside vessels and commonly the training setup can be found in Appendix A. present in image sequences during interventions. In this work, we follow a similar strategy to evaluate the 5.3.1. Preprocessing accuracy of dynamic coronary roadmapping. As the image data have di erent size ranging from The first step is to select frames for roadmapping 512  512 to 1024  1024, all images were resam- evaluation. From each non-contrast sequence in the pled to a grid of 256 256 before being processed by test set for tracking in Section 5.2, we uniformly se- the neural network. In addition, the image intensities lect 8-20 frames to annotate guidewire. The num- were rescaled to the range from 0 to 1. ber of the selected frames from each sequence de- pends on the sequence length, the frame interval size 5.3.2. Training label and guidewire visibility. For some rare cases in our The standard deviation  of the Gaussian PDF for data where no guidewire is present in the image, we the training label of the detection network was set to discarded that non-contrast frame, and chose those 4 pixels in the resampled image space (256  256). frames with little vessel contrast from the same se- This choice corresponds to the estimation of the quence and annotated the vessel centerline. The se- maximal possible catheter tip radius. An example lection results in 409 frames from 34 sequences in to- of the Gaussian PDF is shown in Figure 2c. tal. Once the target non-contrast frames for evaluat- ing roadmapping are chosen, their corresponding an- 5.3.3. Evaluation Metric giographic frames were found using the ECG match- To select hyperparameters and model weights in ing method in Section 3. We then annotated the cen- training, an evaluation metric is required. As the terline of the vessel corresponding to the guidewire deep network is essentially a catheter tip detector, ac- in the non-contrast frames. curate detection of the tip location is desired. There- The next step is performing the transformation of fore, we chose the location with the highest value the labelled vessel centerline from the angiographic in the detection output, and computed the Euclidean frame to its corresponding target non-contrast frame distance between the chosen location and the ground via displacement of catheter tip in the two frames. truth tip coordinate as the evaluation metric to tune This step simulates the roadmapping transformation the deep network. in the last step in Figure 1. Finally, the distance between the guidewire anno- 5.4. Setup for Evaluating Dynamic Coronary tation in the target frame and the transformed vessel Roadmapping centerline is reported as the roadmapping accuracy. It is in general a challenge to evaluate the In order to compute the distance between two point roadmapping accuracy, as the structure of interest, sets of annotations (e.g. Figure 4a), point-point cor- 11 respondence between the two sets is required (Fig- 6. Experiments and Results ure 4b). The point sets were firstly resampled with the point interval being 1 mm. We then followed The following experiments are performed to as- the approach of [44] to find such correspondences sess the methods. First, In Section 6.1, the training which minimizes the sum of the Euclidean distance of the deep neural network is described. Then in Sec- of all valid point-point correspondence paths. This tion 6.2, the accuracy of catheter tip tracking using way guarantees no cross-over connection and each the optimized trained network and the tuned particle point in one set is connected to at least one point in filter is presented. Section 6.3 describes the accuracy the other set. As the annotated point sets may have evaluation of dynamic coronary roadmapping via the di erent size, the point correspondences to endpoints proposed catheter tip tracking method. Finally, in are excluded such that we only focused on the dis- Section 6.4, we measure the processing time of the tance between corresponding sections, not the entire proposed DCR approach. centerlines (Figure 4c). Once the point-point corre- spondence is available, the distance between the two 6.1. Training the Deep Neural Network points in a pair can be used for evaluating the accu- racy of DCR. The purpose of this experiment is to train the deep neural network to output reasonable likelihood prob- ability map. The network hyperparameters were tuned to optimize the detection performance. The training and validation data for detection men- tioned in Section 5.2 were used for training the deep neural network. The evaluation metric mentioned in Section 5.3, the mean Euclidean distance between the ground truth and the predicted tip location av- (a) (b) (c) eraged over all validation frames, was used as the validation criteria for selecting the optimal train- Figure 4: Correspondence between the labelled guidewire (green) and the transformed vessel centerline (red). The yellow ing epoch and the network hyperparameters. When lines connecting the two point sets illustrate the correspondence we evaluated hyperparameter settings, we firstly se- between red and green points. lected the training epoch with the lowest mean vali- dation error for each setting, then the settings were compared using the model weights (trainable net- work parameters) of their chosen epochs. 5.5. Implementation The network hyperparameters we investigated in the experiments include (1) the basic channel num- ber, i.e. the number of channels or feature maps in The proposed method was developed in Python. the first down block, (2) the network depth level, the The framework used for developing the deep learn- number of down or up blocks, and (3) the dropout ing approach for likelihood approximation is Py- probability. Torch. The major experiments of dynamic coronary roadmapping were performed on a computer with an The validation errors for di erent hyperparameter Intel Xeon E5-2620 v3 2.40 GHz CPU and 16 GB settings using the experimental settings in Section RAM running Ubuntu 16.04. The deep neural net- 5.3 are shown in Table 3. The table shows that the work and the optical flow method were running on hyperparameter setting with the lowest mean error, an NVIDIA GeForce GTX 1080 GPU. The approach which has 4 level in depth and 64 channels in the for evaluating dynamic coronary roadmapping was first down block, achieves a validation error of about developed and running in MeVisLab on a computer 2 mm. The table also shows other good choices of with an Intel Core i7-4800MQ 2.70 GHz CPU and network architecture that have a small validation er- 16 GB RAM running Windows 7. ror (shown in red in Table 3): 32 channels in the first 12 in Section 5.2 (see Appendix B for details). Then 1.2 training training 25 in Section 6.2.1, we evaluated the tracking accuracy 1.0 validation validation with the tuned optimal parameter setting (see Table 0.8 0.6 4) on the test dataset, and compared the proposed 0.4 5 tracking method with alternative approaches using 0.2 only the detection network in Section 4.2 or using 0 20 40 60 80 100 0 20 40 60 80 100 only optical flow. Finally, in Section 6.2.2, we inves- (a) Total loss (b) Detection error (mm) 0.5 tigated tracking accuracy with di erent ways of tip training 0.07 training 0.4 validation 0.06 initialization in the first frame. validation 0.05 0.3 0.04 0.2 0.03 6.2.1. Tracking Methods Evaluation 0.02 0.1 0.01 In this experiment, the proposed tracking method 0.0 0.00 0 20 40 60 80 100 0 20 40 60 80 100 in Algorithm 1 uses the ground truth tip probabil- (c) Segmentation loss (d) Detection loss ity map of the first frame as the initial PDF p(x ) to Figure 5: Learning curves for the chosen hyperparameter set- draw samples. This method is referred to as “Track- ting. ing”. In addition, we compared the proposed method with three alternatives. The first one tracks catheter down block with 4 or 5 levels in depth, or 64 chan- tip using only the detection network in Section 4.2 nels with 3 or 4 depth levels. The dropout regulariza- with the chosen network architecture and trained pa- tion improves the accuracy of the model, compared rameters in Section 6.1, therefore, no temporal infor- to the ones without dropout. mation is used. This method is referred to as “De- The learning curves of the training process with tection (Net)”. The other two methods in this ex- the chosen hyperparameter setting are illustrated in periment use only optical flow to track catheter tip Figure 5. The curves indicate that both segmentation starting from the ground truth tip position in the first and detection reach convergence after training 100 frame. The motion field towards the current frame, estimated by the two methods, was based on the de- epochs. We did not investigate a model with more than 64 formation from the previous frame or the first frame channels or 5 depth levels, because (1) it will fur- in the sequence, respectively. The same implemen- ther increase the processing time which makes on- tation setting as in Appendix B.1 was used for these line applications less feasible; (2) the results in Table two methods. They are called “Optical Flow (pre- 3 show that such a setting (64 channels, 5 depth lev- vious)” and “Optical Flow (first)”, or in short form, els) starts increasing the validation error compared to “OF (pre)” and “OF (1st)”. Additionally, we refer those less complex models. the interested readers to Appendix C.1 where the in- The subsequent experiments will be based on the fluence of catheter segmentation on the detection and tracking approaches is reported. network trained with the chosen hyperparameter set- ting indicated in Table 3 (64 channels, 4 depth levels, The tracking accuracies of all methods reported in dropout 0.2, also see Table 4). this section were obtained on the test set from Table 2. The mean, the median and the maximal tracking 6.2. Catheter Tip Tracking error between the predicted and the ground truth tip The purpose of this experiment is to assess the position of all test images are reported in Table 5. accuracy of catheter tip tracking with the proposed In addition, as the sequences in the test set have dif- method in Section 4. Guiding catheter tip is tracked ferent lengths, we also computed the mean and the in X-ray fluoroscopy using Algorithm 1 based on a median error per sequence, and report the the aver- trained network with the optimal hyperparameter set- age of the sequence mean and median errors, so that ting from Section 6.1. First, the parameters of the each sequence contributes equally in these metrics. optical flow method used in Algorithm 1 and particle Table 5 shows that the results from the detection net- filter were tuned on the validation data for tracking work have large average errors which are caused by 13 Table 3: Validation errors (mm) for di erent hyperparameter settings. Red cells show the settings with the 10 smallest validation errors. bold number indicates the setting with the lowest error. Basic Number Depth Dropout of Channels Level none 0.1 0.2 0.3 0.4 0.5 8 3 5.43 4.99 5.02 5.37 4.38 4.24 4 4.17 4.45 4.25 5.04 4.75 4.36 5 3 4.14 3.53 4.28 3.95 4.11 16 3 3.74 4.29 3.57 4.11 3.74 3.4 4 3.36 3.11 3.63 3.33 3.36 3.78 5 3.38 2.89 3.16 2.52 2.71 2.74 32 3 2.99 3.02 3.26 2.82 3.26 2.56 4 2.87 2.34 2.46 2.6 2.65 2.27 5 3.04 2.51 2.21 2.29 2.3 2.25 64 3 2.19 2.54 2.34 2.27 2.26 2.49 4 2.55 2.31 2.04 2.44 2.22 2.27 5 2.42 2.29 2.73 2.77 2.61 2.85 Table 4: The chosen (hyper-)parameters for di erent building blocks of the catheter tip tracking algorithm. The parameters of the optical flow method can be found in Appendix B.1. Building block (hyper-)parameters value Deep learning Basic channel number 64 Depth 4 Dropout 0.2 OF (pre) OF (1st) Detection Tracking OF (pre) OF (1st) Detection Tracking (a) Overall view of tracking er- (b) A zoom-in view of (a) Particle filter  (px) 5 rors N 1000 Figure 6: Tracking errors for the 4 methods on all test images. some completely failed cases. The proposed tracking racy when the target is on the track (row 4), they method has median errors of about 1 mm and mean errors of about 1.3 mm. It achieves the lowest errors present periodic error patterns in two sequences due to large cardiac motion. The detection method shows compared to the other 3 methods on all listed evalu- ation criteria. peaks of large errors, this is because temporal rela- tion between frames is not modeled by the approach, Figure 6 illustrates the boxplots of tracking er- thus the detection on di erent frames is independent rors made by the 4 methods on all test images. It of each other. The proposed tracking method over- shows that the proposed tracking approach outper- comes the problems that other methods have and forms the detection method by avoiding making ex- presents accurate detection on these 4 sequences. tremely large errors (Figure 6a); meanwhile, it main- The tracking results of the 4 methods on example tains as accurate as the detection method for cases frames from the 4 sequences are illustrated in Figure with small errors, and is more accurate than the methods based solely on optical flow (Figure 6b). Figure 7 shows longitudinal views of tracking er- Figure 9 illustrates how the proposed tracking rors of the 4 methods on 4 example sequences. Al- method works on the same 4 frames in Figure 8. though the optical flow methods show high accu- It shows that the prior hypotheses (samples) assists Ground truth / Prediction Distance (mm) Ground truth / Prediction Distance (mm) Table 5: Catheter tip tracking errors (mm) of the 4 methods on the test (tracking) dataset. y indicates that the di erence between that method and the “Tracking” method are statistically highly significant with the two-sided Wilcoxon signed-rank test (p < 0:001). Optical Flowy Optical Flowy Detection Nety Tracking Evaluation Metrics (previous) (first) (Section 4.2) Maximal error of all images 29.16 20.83 108.20 17.72 Median error of all images 1.78 1.22 0.96 0.96 Mean error of all images 3.74  4.93 3.05  4.05 5.62  15.91 1.29  1.76 Average of sequence median error 2.35  2.52 2.64  3.52 6.26  17.11 1.03  0.49 Average of sequence mean error 2.59  2.69 3.31  2.81 6.83  13.88 1.29  0.94 to focus on the correct target location and results views (Figure 7), the outliers for the tracking with au- in reliable posterior estimation, especially when the tomatic initialization are more consistent over time. detection produces ambiguity in cases of multiple Figure 11 shows the temporal change of tracking er- catheters or contrast residual presented in images. rors for the 6 sequences with outliers using the track- ing with automatic initialization. For the 3 sequences 6.2.2. Catheter Tip Initialization on the top row, the tracking with automatic initial- ization makes large errors at the beginning, but be- In this experiment, the initial PDF p(x ) from comes accurate very fast in a few frames; for the 3 which samples are drawn in the proposed tracking sequences on the bottom row, however, the tracking is investigated (Algorithm 1). In particular, we ex- errors remain large till the end of the sequences. plored and evaluated the tracking accuracy with an Figure 12 shows example frames to give an in- automatic initialization using the probability map ob- sight of the tracking with automatic initialization on tained from the trained detection network in Section the 6 sequences in Figure 11. For the 3 sequences 4.2 with the chosen setting in Section 6.1. on the top row (Figure 12a), although the initializa- Figure 10 shows the boxplot of tracking errors on tion on the first frame (frame 0) is overall not cor- all test images with automatic initialization (Auto) rect, the true tip positions are still covered by some and manual initialization (Manual) for which the samples; once the detection in subsequent frames is ground truth tip probability map of the first frame correct, the tracker can still converge to the right tar- was used. The tracking with automatic initialization get. For the 3 sequence on the bottom row (Figure presents an accuracy similar to the one with manual 12b), the initializations of samples are ambiguous in initialization for small tracking errors, but has more frame 0; the detection in subsequent frames focuses large tracking errors which influence the mean error on a wrong area also given by the initial samples due over all test images (Table 6). We, therefore, de- to residual of contrast agent or multiple catheters, the fined the tracking errors on the right side of the gap tracker then tends to find the wrong target. in the boxplot (> 40 mm) as outliers, and explored the statistics without those outliers. 6.3. Dynamic Coronary Roadmapping Table 6 indicates that, the mean and median er- ror of the tracking with automatic initialization ex- In this experiment, the accuracy of dynamic coro- cluding the outliers are only slightly higher than the nary roadmapping using the proposed method with tracking with manual initialization and the detection manual tip initialization was evaluated. For roadmap method. While the tracking with automatic initial- selection with ECG matching (Section 3), the num- ization has 100 outliers in total from 6 sequences, ber of online ECG signal points N was manu- ECG the detection method that has 10 sequences contain- ally determined so that the ECG signal stored in the ing 106 outliers. bu er corresponding to 12 X-ray frames (0.8 second Unlike the detection method for which the outliers in acquisition time). Following the setup in Section are mainly presented as the peaks in the longitudinal 5.4, we used the distance between the two points in 15 Optical Flow (previous) Optical Flow (first) Detection Tracking 80 80 80 60 60 60 40 40 40 20 20 0 0 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 80 80 80 80 60 60 60 60 40 40 40 20 20 20 20 0 0 0 0 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 80 80 80 60 60 60 60 40 40 40 20 20 20 20 0 0 0 0 0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25 80 80 60 60 60 40 40 40 40 20 20 20 20 0 0 0 0 0 20 40 60 80 100 120 0 20 40 60 80 100 120 0 20 40 60 80 100 120 0 20 40 60 80 100 120 Figure 7: Longitudinal view of tracking errors made by the 4 methods on 4 test sequences (one sequence per row). The x-axis denotes the time steps of a sequence, the y-axis is the tracking error (mm). each point pair as the evaluation metric for DCR (the and the frame mean point distances of all 409 evalua- length of a yellow line segment in Figure 4). As each tion frames are illustrated in Figure 13. The compar- frame may have di erent numbers of point pairs, de- ison of the three DCR approaches from Table 7 and pending on the length of the target guidewire, the av- Figure 13 indicates that the accuracy of the proposed erage point pair distance per frame was also com- DCR method has shown improvement over the DCR puted for evaluation. These distances were evaluated without tracking, and is only slightly less than the on 409 selected frames with manual annotation of DCR with manual tip tracking (although the di er- guidewires and vessel centerlines (Section 5.4). ence is statistically significant). Additionally, inter- ested readers are referred to Appendix C.2 where the In the experiment, we compared the DCR with the influence of catheter segmentation on the accuracy of proposed tracking method to those with manual tip DCR is investigated. tracking and without tracking. All three approaches were based on the same ECG matching method (Sec- Table 8 shows how the frame mean point distances tion 3) for selecting roadmaps. The accuracy of the of the 409 evaluation frames are distributed. The DCR without tracking in Table 7 shows that the mean DCR with the proposed method has similar error dis- distances are reduced to less than 3 mm by com- tribution as the one with manual tip tracking: they pensating only cardiac motion via roadmap selection both have about 1/3 of the distances less than 1 mm with ECG matching. Table 7 also shows that the and 1/3 of the distances between 1 and 2 mm. The DCR with the proposed method achieves median dis- proposed method has slightly more distances larger tances of about 1.4 mm and mean distances of about than 5 mm than manual tip tracking. Both methods 2 mm. The boxplots of the distances of all point pairs are more accurate than the DCR without tracking on 16 Sequence 1 Sequence 2 Sequence 3 Sequence 4 Optical Flow Optical Flow Detection Input Tracking Ground Truth (previous) (first) Net Figure 8: Tracking results on example frames from the same 4 sequences in Figure 7. The blue point indicates the predicted catheter tip location; the red point shows the ground truth location. (Best viewed in color) intervals of small errors (< 2 mm). network and the optical flow component of the track- ing method were running on the GPU. Figure 14 shows overlays of selected roadmaps on example frames of 4 sequences with the three DCR In the experiments, the runtimes for roadmap se- approaches. The DCR without tracking presents mis- lection (step 1) and roadmap transformation (step 3) match of catheters, guidewires or residual of contrast in Figure 1 were negligible (< 1 ms / frame). The agent in the images, whereas the other methods im- runtime of the proposed catheter tip tracking method prove the alignment and show good match between is shown in Table 9 and Figure 15. The average the structures in the original X-ray image and the time to compute the likelihood with the deep learning roadmaps. Compared to the DCR with manual tip setup (DL) is 31.5 ms / frame. The particle filtering tracking, the proposed method show similar visual (PF) step, which consists of the optical flow estima- alignment of the roadmaps to the original X-ray im- tion, sample propagation, sample weight update and ages. For a dynamic view of a roadmapping exam- normalization, prediction and resampling, takes on ple, we refer readers to the supplemental material. average 23 ms / frame. Therefore, the average track- ing time in total is 54.5 ms / frame. The total av- 6.4. Processing Time erage time of the proposed DCR including roadmap The processing time of all steps in the proposed selection, catheter tip tracking and roadmap transfor- DCR method was measured with the hardware and mation is still less than the acquisition time of our software setup in Section 5.5. The ECG matching data (66.7 ms / frame, 15 fps), indicating that the method for roadmap selection was running in Python proposed DCR method would run in real-time with on the CPU of the Linux machine; the deep neural our setup. 17 Sequence 1 Sequence 2 Sequence 3 Sequence 4 Manual Auto Segmentation Detection Particles Particles Prediction Input Ground Truth (heatmap) (likelihood) (prior) (posterior) (expectation) Figure 9: Workflow of the proposed tracking method on the same 4 frames in Figure 8. The high probability is shown with bright color in the detection map. Samples or particles are presented as green dots. The blue point indicates the predicted catheter tip location; the red point shows the ground truth location. (Best viewed in color) 50 60 40 60 0 0 0 0 20 40 60 80 100 0 5 10 15 20 25 0 10 20 30 40 50 60 70 60 80 60 60 40 30 40 20 20 0 20 40 60 80 100 0 0 0 Tracking error (mm) 0 3 6 9 12 15 18 0 10 20 30 40 50 60 70 0 1 2 3 4 5 6 7 Figure 10: Catheter tip tracking errors (mm) with manual and Figure 11: Longitudinal views of tracking errors (mm) for the automatic initialization. 6 sequences with outliers using automatic initialization. 7. Discussion using a proposed deep learning based Bayesian fil- We have presented a new approach to perform on- tering. The proposed tracking method represents and line dynamic coronary roadmapping on X-ray flu- tracks the posterior of catheter tip via a particle fil- oroscopic sequences for PCI procedures. The ap- ter, for which a likelihood probability map is com- proach compensates the cardiac-induced vessel mo- puted for updating the particle weights using a con- tion via selecting oine-stored roadmaps with ap- volutional neural network. In the experiments, the propriate cardiac phase using ECG matching, and proposed DCR approach has been trained and eval- corrects the respiratory motion of vessels by online uated on clinical X-ray sequences for both tracking tracking of guiding catheter tip in X-ray fluoroscopy and roadmapping tasks. 18 Table 6: Catheter tip tracking errors (mm) of detection and tracking with manual and automatic initialization Tracking Detection Manual init. Automatic init. Maximal error 108.20 17.23 98.58 Median error 0.96 0.96 0.96 Mean error 5.62  15.91 1.29  1.76 5.16  13.91 No. of outliers (> 40 mm) 106 0 100 No. of sequences with outliers 10 0 6 Maximal error of inliers 31.06 17.23 28.28 Median error of inliers 0.96 0.96 0.96 Mean error of inliers 1.17  1.78 1.29  1.76 1.34  2.15 Table 7: The statistics of DCR accuracy (mm) with three di erent tracking scenarios. With the two-sided Wilcoxon signed-rank test: y denotes that the di erence between the DCR without tracking and that with the proposed tracking method is statistically highly significant (p < 0:001); * indicates a statistically significantly di erence between the DCR using manual tip tracking and that with the proposed tracking approach (p < 0:05). Without Trackingy Proposed Tracking Method Manual Tip Tracking* All point pairs Maximal distance 27.19 20.24 13.12 Median distance 1.97 1.43 1.35 Mean distance 2.94  2.83 2.07  2.08 1.85  1.72 Frame mean distance Median distance 2.11 1.42 1.38 Average distance 2.76  2.08 1.91  1.52 1.75  1.30 One prerequisite of accurate tracking with the pro- formance does not increase much compared to those posed approach is to obtain a reasonably good like- with simpler settings, implying that the model starts lihood estimation, which requires to train the deep overfitting on our dataset. neural network to detect catheter tip well. In this In addition to the deep neural network, the other work, we have investigated the influence of three important component of the proposed tracking ap- network hyperparameters on the performance of the proach is the sampling in the particle filter that yields detection network (Section 6.1): the basic channel the samples for representing the prior and the poste- number and network depth level are model capac- rior of catheter tip position. First, a sucient number ity parameters, the dropout adds regularization to the of samples in the whole sample space are required model. The experiment showed that the detection ac- to well characterize the probability distributions (see curacy improves when the basic channel number and Appendix B.2). Second, the sample dynamics plays the network depth level increase (Table 3). This ob- an important role in tracking, in particular, as indi- servation matches the expectation that a more com- cated by Eq.(8), the process noise and the sample plex model has higher capacity to model the variation motion. The process noise influences the tracking in the data, hence results in better accuracy. When accuracy, according to Table 10 in Appendix B.2. the complexity reaches a certain level, e.g. 64 ba- Additionally, sample motion is another key aspect of sic channels and 5 level of depth, the network per- sample dynamics. Motion estimation has previouly 19 Table 8: Distribution of frame mean point distances of the 409 evaluation frames. Error Intervals (mm) Tracking Methods of DCR < 1 1-2 2-3 3-4 4-5  5 Without tracking 81 115 69 47 31 66 Proposed tracking method 131 145 61 32 17 23 Manual tip tracking 139 144 61 35 20 10 Table 9: Statistics of the runtime of catheter tip tracking (ms / frame) on the test (tracking) dataset. Deep Learning Particle Filtering Total Tracking Time Mean 31.5  10.3 23.0  8.7 54.5  12.3 Median 35.1 22.8 57.7 been incorporated in a motion-based particle filter, pecially when the CNN detector fails to find the cor- such as adaptive block matching ([6]). In our work, rect target area. The association of knowledge from optical flow was chosen for motion estimation, as two sources together improves the tracking accuracy its non-parametric nature allows to characterize the compared to each single source. complexity of motion in X-ray fluoroscopy well. In The initial state is a also key component of track- addition, the advantage of such approach from a the- ing approaches. In the context of Bayesian filtering, oretical point of view is that it takes into account of the initial state provides the prior knowledge of the the current observation, leading to a more optimal tracking target. Most tracking algorithms assume a importance density ([3]) compared to random mo- known initial state from which the target is tracked, tion. e.g. our proposed method with manual initialization The tracking results in Section 6.2.1 show that in Section 6.2.2. In this case, the prior knowledge is the proposed tracking approach is able to track the provided by human. In Section 6.2.2, we also investi- catheter tip in X-ray fluoroscopy accurately with an gated a scenario where the initial state is given by the average tracking error of about 1.3 mm. It also shows detection network, so that the complete tracking pro- advantages over methods based only on optical flow cess is fully automated. The results indicate that, the or the detection network. The OF (pre) method relies proposed tracking method with automatic initializa- heavily on tracking in the previous frame, hence the tion works reasonably well on most sequences even error could accumulate. The OF (first) method may when the initialization is sometimes incorrect (Fig- su er from large motion from the first frame to the ure 12a). This is because (1) the true position is cov- current frame. The detection method uses informa- ered by a few samples, and (2) the correct detection tion only from the current frame, no temporal rela- in later frames corrects the initial mistake in the first tion between frames is utilized; therefore, it results in frame. The automatic initialization fails when (1) spikes in the longitudinal view, as shown in Figure 7. a wrong position is covered by a few samples and The proposed tracking method has a CNN to provide (2) the wrong detection in subsequent frames con- an accurate observation on the current frame which firms the mistake in the initial frame (Figure 12b). improves the accuracy of optical flow tracking within This happens when there is contrast agent remaining the framework of Bayesian filtering. In the mean- in the image or there are multiple catheters, which time, the optical flow based particle filter maintains are the major sources causing ambiguity in detec- and propagates the prior knowledge from the initial tion. In practice, the automatic initialization would tip position to provide a constraint on searching for work well when contrast agent is washed out and the potentially correct positions, which is useful es- only one catheter is present in the field of view, oth- 20 Detection Initial Particles Detection Posterior Frame 0 Frame 10 Frame 0 Frame 10 Frame 0 Frame 0 Frame 5 Frame 5 Figure 13: Accuracy (mm) of DCR with three di erent tracking scenarios. Frame 20 Frame 0 Frame 20 Frame 0 (a) Sequence 1-3 on the top row in Figure 11 Detection Initial Particles Detection Posterior Original X-ray Without Catheter With The Proposed With Manual Tip Image Tip Tracking Method Tracking Frame 0 Frame 0 Frame 4 Frame 4 Frame 0 Frame 15 Frame 0 Frame 15 Frame 7 Frame 0 Frame 0 Frame 7 (b) Sequence 4-6 on the bottom row in Figure 11 Figure 12: Examples frames from the 6 sequences in Figure 11. The high probability in the detection heatmap is highlighted as bright color. Particles are presented as green dots. The red dots in the last column indicate the ground truth tip location. (Best viewed in color) Figure 14: Examples of superimposition of selected roadmaps (red) on X-ray fluoroscopic frames. (Best viewed in color) 21 Image acquisition data. Additionally, variation exists between di erent DL time 66.7 ms / frame cardiac cycles ([27]), therefore, choosing a roadmap from another cycle may cause inaccuracy for cardiac PF motion compensation. Finally, the way of DCR eval- Total uation in Section 5.4 might also introduce inaccura- 0 20 40 60 80 100 cies in the error measurement. Since guidewires of- Processing Time (ms) ten attach to the inner curves of a vessel to take the Figure 15: Runtime of catheter tip tracking (ms / frame) on all shortest path, the small di erence between the an- test frames. notation of guidewire and vessel centerlines was ig- nored in the evaluation. erwise manual initialization would be needed which In addition to accuracy, processing speed is also requires only one click to initiate tracking. critical for intraoperative applications. The results in Section 6.4 indicate that the total processing time of Dynamic coronary roadmapping is the direct ap- the proposed DCR approach is less than the image plication of the catheter tip tracking results. In our acquisition time, meaning that it runs in real-time on experiments, the DCR was performed with manual our setup. To build a real-time system for PCI in tip initialization to show the potential of the proposed practice, the overall latency of the complete system tracking method, and was compared with the DCR needs to be considered. It is also worth noticing that without tracking and with manual tracking. The the DL and PF steps of the proposed tracking method results indicate that using catheter tip tracking can are independent from each other. In practice, in case improve DCR accuracy, as the respiratory-induced more than one GPU are available, the proposed DCR vessel motion is corrected by the displacement of approach can be further accelerated by paralleling catheter tip in addtion to cardiac motion correction. the DL and PF steps, making them running on dif- The results also show that the proposed DCR reaches ferent GPUs. a good accuracy (mean error is about 2 mm) and per- Compared to the previous works on DCR, the pro- forms only slightly worse than its best case, the DCR with manual tip tracking which is not applicable for posed approach in this paper shows advancement in intraoperative use. Additionally, according to a pre- several aspects. First, our systems works on non- vious study by [10], the average lumen diameters of cardiac-gated sequences which does not require ad- human coronary arteries are between 1.9 mm (dis- ditional setups for cardiac motion gating that were tal left anterior descending artery) and 4.5 mm (left needed for some methods ([50, 26]). Second, our main artery). This means that the accuracy achieved approach compensates both respiratory- and cardiac- with the proposed approach is comparable with the induced vessel motion, which is more accurate than size of coronary arteries. systems that correct only cardiac motion ([12]). In addition, the proposed DCR approach follows a data- Apart from catheter tip tracking, several other pos- driven paradigm that learns target feature from se- sible factors in di erent steps of the experiments may quences acquired from di erent patients and various influence the final DCR accuracy. First, in the oine view angles, making it more robust than the method phase, the signal of contrast agent may become too that relies on traditional vesselness filtering ([20]) strong and completely cover the catheter tip, com- or methods that require specific tissue being present plicating the tip visibility in some cases. In this sit- ([50, 26]). These are the major advantages of the uation, the uncertainty in the manual tip annotation proposed DCR over the existing direct roadmapping may result in errors in roadmap transformation. Sec- systems. Compared to model-based motion compen- ond, in the roadmap selection step, the oine-stored roadmaps are only discrete samples of complete car- sation, our approach does not require the extraction of motion surrogate signals and train a motion model diac cycles which might not fully characterize every for each new patient, but can be directly run with a possible change in the cardiac motion. This problem trained model. could possibly be addressed in the future by inter- polating frames between the existing frames in the The proposed deep learning based Bayesian filter- 22 ing method has several advantages over the existing In the future, it may be worth investigating the instrument tracking methods. First, the deep learning following directions related to this work. As the component enables a more general framework to de- data used in this study was acquired from one hos- tect instruments in medical images than methods tai- pital using a machine from a single vendor, it would lored for specific tools ([25, 24]). Compared to the be interesting to evaluate the proposed approach on existing detection methods based on deep learning multi-center data acquired with machines from dif- ([5, 22, 11]), our approach takes into account of the ferent vendors. Next, since the ECG signals of our information between frames; the Bayesian filtering data appear to be regular, it may be necessary in a framework allows interaction between temporal in- future study to acquire data with irregular ECG that formation and the detection of a convolutional neural could be obtained in practice, and validate the pro- network, making the tracking more robust. Bayesian posed approach on those data. Besides, it would frameworks have been used in many previous tem- be also interesting to validate our approach during poral instrument tracking methods. Particularly, the PCI procedures in an environment simulating the real likelihood term in some works was designed based clinical settings. Additionally, from a methodolog- on registration or segmentation outcomes ([2, 38]) or ical point of view, although the proposed tracking traditional machine learning approaches with hand- method is invariant under di erent view angles, the crafted features ([45, 46, 29]). In our method, we whole DCR approach works only when the oine approximated the likelihood with a deep neural net- and online phase have the same view angle, i.e. it work learned from the clinical data which exempts is a 2D roadmapping system. Therefore, one future the need of feature engineering but yet possesses direction would be to develop a 3D DCR system that more discriminative power; the network directly out- would work with various view angles in the online puts the probability map, making it more straightfor- phase. ward to use. Finally, compared to the existing instru- ment tracking approaches based on Bayesian filter- 8. Conclusion ing ([2, 38, 39]), the state transition in our method was based on the motion estimated from two adja- We have developed and validated a novel approach cent frames, which is more reliable than totally ran- to perform dynamic coronary roadmapping for PCI dom motion or artificially-designed state transition image guidance. The approach compensates car- models. diac motion through ECG alignment and respiratory motion by guiding catheter tip tracking during fluo- From a practical point of view, the proposed DCR roscopy with a deep learning based Bayesian filter- approach could potentially fit into the clinical work- ing method. The proposed tracking and roadmap- flow of PCI. The oine phase of the method can ping approaches were trained and evaluated on clini- be done eciently by a person who assists the pro- cal X-ray image datasets and were proved to perform cedures: selecting and creating roadmaps from an accurately on both catheter tip tracking and dynamic angiography acquisition, annotating the catheter tip coronary roadmapping tasks. Our approach also runs (one point) in the images. This phase is typically in real-time on a setup with a modern GPU and thus done before a fluoroscopy acquisition during which has the potential to be integrated into routine PCI the guidewire advancement and stent placement are procedures, assisting the operator with real-time vi- performed. In the online phase, when a fluoroscopic sual image guidance. image is acquired, the proposed system selects the most suitable roadmap, tracks the catheter tip and transforms the roadmap to prospectively show a ves- Acknowledgment sel overlay on the fluorosocpic image. The online updated coronary roadmap can provide real-time vi- This work was supported by NWO-domain TTW sual guidance to cardiologists to manipulate inter- (The Division of Applied and Engineering Sciences ventional tools during the procedure without the need in The Netherlands Organisation for Scientific Re- of administering extra contrast agent. search), IMAGIC project under the iMIT program 23 (grant number 12703). Ries and Simon van Walsum [12] Elion, J.L., 1989. Dynamic coronary roadmapping. US Patent 4,878,115. are acknowledged for their contribution in the man- [13] Faranesh, A.Z., Kellman, P., Ratnayaka, K., Lederman, ual annotations. R.J., 2013. Integration of cardiac and respiratory motion into MRI roadmaps fused with X-ray. Medical physics References [14] Farneback, G., 2003. Two-frame motion estimation based on polynomial expansion, in: Scandinavian conference on [1] Ambrosini, P., Ruijters, D., Niessen, W.J., Moelker, A., Image analysis, Springer. pp. 363–370. van Walsum, T., 2017a. Fully automatic and real-time [15] Fischer, P., Faranesh, A., Pohl, T., Maier, A., Rogers, T., catheter segmentation in X-ray fluoroscopy, in: Inter- Ratnayaka, K., Lederman, R., Hornegger, J., 2018. An national Conference on Medical Image Computing and MR-based model for cardio-respiratory motion compen- Computer-Assisted Intervention, Springer. pp. 577–585. sation of overlays in X-ray fluoroscopy. IEEE transactions [2] Ambrosini, P., Smal, I., Ruijters, D., Niessen, W.J., on medical imaging 37, 47–60. Moelker, A., Van Walsum, T., 2017b. A hidden markov [16] Garc´ ıa-Peraza-Herrera, L.C., Li, W., Gruijthuijsen, C., model for 3D catheter tip tracking with 2D X-ray catheter- Devreker, A., Attilakos, G., Deprest, J., Vander Poorten, ization sequence and 3D rotational angiography. IEEE E., Stoyanov, D., Vercauteren, T., Ourselin, S., 2016. Trans. Med. Imaging 36, 757–768. Real-time segmentation of non-rigid surgical tools based [3] Arulampalam, M.S., Maskell, S., Gordon, N., Clapp, on deep learning and tracking, in: International Workshop T., 2002. A tutorial on particle filters for online on Computer-Assisted and Robotic Endoscopy, Springer. nonlinear/non-Gaussian Bayesian tracking. IEEE Trans- pp. 84–95. actions on signal processing 50, 174–188. [17] He, K., Zhang, X., Ren, S., Sun, J., 2016. Deep resid- [4] Baka, N., Lelieveldt, B., Schultz, C., Niessen, W., van ual learning for image recognition, in: Proceedings of the Walsum, T., 2015. Respiratory motion estimation in X- IEEE conference on computer vision and pattern recogni- ray angiography for improved guidance during coronary tion, pp. 770–778. interventions. Physics in Medicine & Biology 60, 3617. [18] Heibel, H., Glocker, B., Groher, M., Pfister, M., Navab, [5] Baur, C., Albarqouni, S., Demirci, S., Navab, N., N., 2013. Interventional tool tracking using discrete opti- Fallavollita, P., 2016. Cathnets: detection and single- mization. IEEE transactions on medical imaging 32, 544– view depth prediction of catheter electrodes, in: Interna- tional Conference on Medical Imaging and Virtual Real- [19] Honnorat, N., Vaillant, R., Paragios, N., 2011. Graph- ity, Springer. pp. 38–49. based geometric-iconic guide-wire tracking, in: Inter- [6] Bouaynaya, N., Schonfeld, D., 2009. On the optimality national Conference on Medical Image Computing and of motion-based particle filtering. IEEE transactions on Computer-Assisted Intervention, Springer. pp. 9–16. circuits and systems for video technology 19, 1068–1072. [20] Kim, D., Park, S., Jeong, M.H., Ryu, J., 2018. Regis- [7] Brost, A., Liao, R., Strobel, N., Hornegger, J., 2010. Res- tration of angiographic image on real-time fluoroscopic piratory motion compensation by model-based catheter image for image-guided percutaneous coronary interven- tracking during EP procedures. Medical Image Analysis tion. International journal of computer assisted radiology 14, 695–706. and surgery 13, 203–213. [8] Chang, P.L., Rolls, A., De Praetere, H., Vander Poorten, [21] King, A.P., Boubertakh, R., Rhode, K.S., Ma, Y., Chin- E., Riga, C.V., Bicknell, C.D., Stoyanov, D., 2016. Ro- chapatnam, P., Gao, G., Tangcharoen, T., Ginks, M., bust catheter and guidewire tracking using B-spline tube Cooklin, M., Gill, J.S., et al., 2009. A subject-specific model and pixel-wise posteriors. IEEE Robotics and Au- technique for respiratory motion correction in image- tomation Letters 1, 303–308. guided cardiac catheterisation procedures. Medical Image [9] Dannenberg, L., Polzin, A., Bullens, R., Kelm, M., Zeus, Analysis 13, 419–431. T., 2016. On the road: First-in-man bifurcation percu- [22] Laina, I., Rieke, N., Rupprecht, C., Vizca´ ıno, J.P., Eslami, taneous coronary intervention with the use of a dynamic A., Tombari, F., Navab, N., 2017. Concurrent segmenta- coronary road map and stentboost live imaging system. tion and localization for tracking of surgical instruments, International journal of cardiology 215, 7–8. in: International conference on medical image computing [10] Dodge, J.T., Brown, B.G., Bolson, E.L., Dodge, H.T., and computer-assisted intervention, Springer. pp. 664– 1992. Lumen diameter of normal human coronary arter- ies. influence of age, sex, anatomic variation, and left ven- [23] Ma, H., Dibildox, G., Banerjee, J., Niessen, W., Schultz, tricular hypertrophy or dilation. Circulation 86, 232–246. C., Regar, E., van Walsum, T., 2015. Layer separation for [11] Du, X., Kurmann, T., Chang, P.L., Allan, M., Ourselin, S., vessel enhancement in interventional X-ray angiograms Sznitman, R., Kelly, J.D., Stoyanov, D., 2018. Articulated using morphological filtering and robust PCA, in: Work- multi-instrument 2D pose estimation using fully convolu- shop on Augmented Environments for Computer-Assisted tional networks. IEEE transactions on medical imaging Interventions, Springer. pp. 104–113. 24 [24] Ma, Y., Gogin, N., Cathier, P., Housden, R.J., Gijsbers, image-guided cardiac interventions, in: 2010 IEEE Com- G., Cooklin, M., O’Neill, M., Gill, J., Rinaldi, C.A., puter Society Conference on Computer Vision and Pattern Razavi, R., et al., 2013. Real-time X-ray fluoroscopy- Recognition, IEEE. pp. 2948–2954. based catheter detection and tracking for cardiac electro- [36] Shechter, G., Ozturk, C., Resar, J.R., McVeigh, E.R., physiology interventions. Medical physics 40. 2004. Respiratory motion of the heart from free breath- [25] Ma, Y., King, A.P., Gogin, N., Gijsbers, G., Rinaldi, C.A., ing coronary angiograms. IEEE transactions on medical Gill, J., Razavi, R., Rhode, K.S., 2012. Clinical evalua- imaging 23, 1046. tion of respiratory motion compensation for anatomical [37] Shechter, G., Shechter, B., Resar, J.R., Beyar, R., 2005. roadmap guided cardiac electrophysiology procedures. Prospective motion correction of X-ray images for coro- IEEE Transactions on Biomedical Engineering 59, 122– nary interventions. IEEE Transactions on Medical Imag- 131. ing 24, 441–450. [26] Manhart, M., Zhu, Y., Vitanovski, D., 2011. Self- [38] Speidel, S., Delles, M., Gutt, C., Dillmann, R., 2006. assessing image-based respiratory motion compensation Tracking of instruments in minimally invasive surgery for fluoroscopic coronary roadmapping, in: Biomedical for surgical skill analysis, in: International Workshop on Imaging: From Nano to Macro, 2011 IEEE International Medical Imaging and Virtual Reality, Springer. pp. 148– Symposium on, IEEE. pp. 1065–1069. 155. [27] McClelland, J.R., Hawkes, D.J., Schae ter, T., King, [39] Speidel, S., Kuhn, E., Bodenstedt, S., Rohl, ¨ S., Ken- A.P., 2013. Respiratory motion models: a review. Medi- ngott, H., Muller ¨ -Stich, B., Dillmann, R., 2014. Visual cal image analysis 17, 19–42. tracking of da vinci instruments for laparoscopic surgery, [28] Milletari, F., Navab, N., Ahmadi, S.A., 2016. V-net: Fully in: Medical Imaging 2014: Image-Guided Procedures, convolutional neural networks for volumetric medical im- Robotic Interventions, and Modeling, International Soci- age segmentation, in: 3D Vision (3DV), 2016 Fourth In- ety for Optics and Photonics. p. 903608. ternational Conference on, IEEE. pp. 565–571. [40] Takimura, H., Muramatsu, T., Tsukahara, R., Nakano, M., [29] Pauly, O., Heibel, H., Navab, N., 2010. A machine learn- Nishio, S., Takimura, Y., Yabe, T., Kawano, M., Hada, T., ing approach for deformable guide-wire tracking in fluo- 2018. Usefulness of novel roadmap guided percutaneous roscopic sequences, in: International Conference on Med- coronary intervention for bifurcation lesions (the world’s ical Image Computing and Computer-Assisted Interven- first cases). Journal of the American College of Cardiol- tion, Springer. pp. 343–350. ogy 71, A1365. [30] Peressutti, D., Penney, G.P., Housden, R.J., Kolbitsch, C., [41] Tehrani, S., Laing, C., Yellon, D.M., Hausenloy, D.J., Gomez, A., Rijkhorst, E.J., Barratt, D.C., Rhode, K.S., 2013. Contrast-induced acute kidney injury following pci. King, A.P., 2013. A novel Bayesian respiratory mo- European journal of clinical investigation 43, 483–490. tion model to estimate and resolve uncertainty in image- [42] Timinger, H., Krueger, S., Dietmayer, K., Borgert, J., guided cardiac interventions. Medical image analysis 17, 2005. Motion compensated coronary interventional navi- 488–502. gation by means of diaphragm tracking and elastic motion [31] Petkovic, ´ T., Loncari ˇ c, ´ S., 2010. Guidewire tracking with models. Physics in Medicine & Biology 50, 491. projected thickness estimation, in: Biomedical Imaging: [43] Vandini, A., Glocker, B., Hamady, M., Yang, G.Z., 2017. From Nano to Macro, 2010 IEEE International Sympo- Robust guidewire tracking under large deformations com- sium on, IEEE. pp. 1253–1256. bining segment-like features (seglets). Medical image [32] Piayda, K., Kleinebrecht, L., Afzal, S., Bullens, R., ter analysis 38, 150–164. Horst, I., Polzin, A., Veulemans, V., Dannenberg, L., [44] van Walsum, T., Schaap, M., Metz, C.T., van der Giessen, Wimmer, A.C., Jung, C., et al., 2018. Dynamic coronary A.G., Niessen, W.J., 2008. Averaging centerlines: mean roadmapping during percutaneous coronary intervention: shift on paths, in: International Conference on Medical a feasibility study. European journal of medical research Image Computing and Computer-Assisted Intervention, 23, 36. Springer. pp. 900–907. [33] Rieke, N., Tan, D.J., di San Filippo, C.A., Tombari, F., [45] Wang, P., Chen, T., Zhu, Y., Zhang, W., Zhou, S.K., Co- Alsheakhali, M., Belagiannis, V., Eslami, A., Navab, N., maniciu, D., 2009. Robust guidewire tracking in fluo- 2016. Real-time localization of articulated surgical instru- roscopy, in: 2009 IEEE Conference on Computer Vision ments in retinal microsurgery. Medical image analysis 34, and Pattern Recognition, IEEE. pp. 691–698. 82–100. [46] Wu, W., Chen, T., Strobel, N., Comaniciu, D., 2012. Fast [34] Ronneberger, O., Fischer, P., Brox, T., 2015. U-net: tracking of catheters in 2D fluoroscopic images using an Convolutional networks for biomedical image segmenta- integrated CPU-GPU framework, in: Biomedical Imag- tion, in: International Conference on Medical image com- ing (ISBI), 2012 9th IEEE International Symposium on, puting and computer-assisted intervention, Springer. pp. IEEE. pp. 1184–1187. 234–241. [47] Wu, X., Housden, J., Ma, Y., Razavi, B., Rhode, K., [35] Schneider, M., Sundar, H., Liao, R., Hornegger, J., Xu, C., Rueckert, D., 2015. Fast catheter segmentation from 2010. Model-based respiratory motion compensation for echocardiographic sequences based on segmentation from 25 corresponding X-ray fluoroscopy for cardiac catheteriza- B. Parameter Tuning for Catheter Tip Tracking tion interventions. IEEE transactions on medical imaging This section describes the details of tuning the pa- 34, 861–876. rameters of optical flow and particle filter for catheter [48] Yabe, T., Muramatsu, T., Tsukahara, R., Nakano, M., tip tracking. Takimura, Y., Takimura, H., Kawano, M., Hada, T., 2018. The impact of percutaneous coronary intervention using the novel dynamic coronary roadmap system. Journal of B.1. Tuning Optical Flow Parameters the American College of Cardiology 71, A1103. The approach of [14] was used as the optical flow [49] Yatziv, L., Chartouni, M., Datta, S., Sapiro, G., 2012. To- implementation in Algorithm 1. A grid search to find ward multiple catheters detection in fluoroscopic image the optimal parameter setting was done on the fol- guided interventions. IEEE Transactions on Information Technology in Biomedicine 16, 770–781. lowing parameters of the method: (1) the image scale [50] Zhu, Y., Tsin, Y., Sundar, H., Sauer, F., 2010. Image- to build the pyramids, (2) the number pyramid levels, based respiratory motion compensation for fluoroscopic (3) the averaging window size, (4) the number of it- coronary roadmapping, in: International Conference on erations, (5) the size of the pixel neighborhood used Medical Image Computing and Computer-Assisted Inter- to find polynomial expansion in each pixel, and fi- vention, Springer. pp. 287–294. nally (6) the standard deviation of the Gaussian that is used to smooth derivatives used as a basis for the Appendix polynomial expansion. The above parameters were tuned independently A. Details of the Training Setup of the deep neural network, as optical flow di- A.1. Data Augmentation rectly estimates the catheter tip motion between two To increase the number of training samples and frames. To tune the parameters, we tracked the their diversity, data augmentation was used. The aug- catheter tip in X-ray fluoroscopy starting from the mentation includes geometric transformation such as ground truth tip position in the first frame using the flipping (left-right, up-down), rotation of multiple of motion field between two adjacent frames estimated 90 degrees, random ane transformation (translation with optical flow. The average and median distance -10 to 10 px, scaling 0.9 to 1.1, rotation -5 to 5 de- between the tracked tip position and the ground truth grees, shear -5 to 5 px), random elastic deformation were used as the evaluation criteria for the tuning. (deformation range -4 to 4 px, grid size of control The method of [14] was implemented by using the points 64 px). A training sample has 0.5 probability OpenCV function calcOpticalFlowFarneback. of being processed with one of the transformations. With consideration of the suggested parameter val- The probability for applying each transformation is: ues from the documentation, the parameter setting flipping left-right (1/24), flipping up-down (1/24), chosen for optical flow from the grid search is rotation of multiple of 90 degrees (1/12), ane trans- pyr scale = 0:5, levels = 3, winsize = 10, formation (1/6), elastic deformation (1/6), no trans- iterations = 30, poly n = 5, poly sigma = 1:1. formation (1/2). To make the trained model robust to Details of the parameters can be found on the func- noise, in addition to the geometric transformations, tion documentation page . we also augmented data by adding Gaussian noise to the pixel value with a zero mean and a standard de- B.2. Tuning Particle Filter Parameters viation between 0.01 and 0.03. The probability of The parameters to tune for the particle filter are adding the noise is 0.5. the number of samples N and the variance of pro- cess noise  . When tuning them, we fixed the pa- A.2. Training Settings v rameters of the trained network and the optical flow The  value in the training loss Eq. (5) was set method, and used their optimal parameter settings to 10 to make the scale of the two terms similar. during this experiment. Following Algorithm 1, we Adam optimizer was used to minimize the loss func- tion with a learning rate 0.0001. The number of train- ing samples in a batch is 4. The network was trained https://docs.opencv.org/2.4/modules/video/ with 100 epochs to ensure convergence. doc/motion_analysis_and_object_tracking.html? 26 tracked the catheter tip from the ground truth posi- Detection 6 Detection 120 Tracking Tracking tion (probability map) in the first frame, and used the mean and median distance between the tracked and the true position as the validation metric. The tracking results on the validation (tracking) 20 1 set are shown in Table 10. The table shows that 100 With Segmentation Without Segmentation With Segmentation Without Segmentation samples are suboptimal, while 1000 samples seem (a) Overall view of tracking er- (b) A zoom-in view of (a) sucient, as 10000 samples result in tracking accu- rors racies similar to 1000 samples. It also shows that the optimal choices of the standard deviation of the Figure 16: Tracking errors on all test images with and without process noise are 4 or 5 px for the downsampled im- catheter segmentation. ages. One possible reason for such choices may be that they are similar to the size of guiding catheters. C.1. Catheter tip detection and tracking In general, good choices for N are 1000 and 10000, The same metrics in Table 5 are used to report for  are 4 and 5. By considering the mean, the the accuracy of catheter tip detection and tracking standard deviation and the median of tracking errors, without catheter segmentation. Table 11 and Fig- the parameter setting  = 5, N = 1000 was chosen. v s ure 16 both manifest that the segmentation sub-task improves the accuracy of catheter tip detection and tracking. Although the improvement on the track- C. Detection, tracking and roadmapping without ing task is marginal and not statistically significant catheter segmentation (p = 0:06), the segmentation helps to reduce the magnitude and amount of outliers (large errors). Training of the network in Figure 3 requires catheter labels for detection and segmentation. As C.2. Dynamic coronary roadmapping the segmentation labels are often more expensive to In this experiment, the same setup in Section acquire than the detection label in practice, we also 6.3 was used to assess the accuracy of DCR using investigated the performance of catheter tip detec- catheter tip tracking without segmenting the catheter. tion, tracking and dynamic coronary roadmapping Table 12 indicate that tracking the catheter tip with without segmenting the catheter. To this end, we catheter segmentation improves the DCR accuracy used a similar encoder-decoder network architecture compared to that without catheter segmentation. Al- as Figure 3 which computes only the detection output though the improvement is not statistically signifi- directly after the last up block of the decoder with a cant (p = 0:43), the approach with segmentation is 1 1 convolution followed by a spatial softmax layer, more robust by making less large roadmapping er- instead of having two outputs. We then followed rors (Figure 17). the same way as the approach using the network with segmentation to search for (hyper-)parameters for the approach without segmentation. The follow- ing parameter setting was chosen for the experiments in this section: for deep learning, the basic channel number is 64, the depth is 5, the dropout rate is zero; for particle filtering,  = 3, N = 10000. With this v s setup, we compared the performance of the approach without catheter segmentation to the proposed ap- proach with segmentation on catheter tip detection and tracking (Appendix C.1) and dynamic coronary roadmapping (Appendix C.2) on the test set from Ta- ble 2. Ground truth / Prediction Distance (mm) Ground truth / Prediction Distance (mm) Table 10: Catheter tip tracking errors (mm) on the validation (tracking) dataset of di erent parameter settings for particle filter. The tracked tip point was rounded to the pixel center. The error of all images (mean  std / median) are presented. Red cells show the good choices of parameters; bold number indicates the chosen setting. v s (px) 100 1000 10000 3 1.52  2.19 / 0.79 1.49  2.18 / 0.79 1.48  2.18 / 0.79 4 1.50  2.17 / 0.79 1.46  2.17 / 0.79 1.47  2.18 / 0.79 5 1.52  2.21 / 0.79 1.47  2.17 / 0.74 1.47  2.19 / 0.74 6 1.53  2.39 / 0.79 1.49  2.33 / 0.79 1.48  2.29 / 0.74 7 1.56  2.42 / 0.79 1.50  2.29 / 0.74 1.50  2.39 / 0.74 8 1.58  2.41 / 0.79 1.51  2.40 / 0.74 1.51  2.42 / 0.74 9 1.56  2.22 / 0.79 1.53  2.43 / 0.79 1.52  2.45 / 0.61 10 2.25  6.18 / 0.79 1.54  2.46 / 0.79 1.53  2.47 / 0.61 Table 11: Catheter tip tracking errors (mm) with and without catheter segmentation on the test (tracking) dataset. y indicates that the di erence between the detection with and without segmentation is statistically highly significant with the two-sided Wilcoxon signed-rank test (p < 0:001). No statistically significantly di erence is observed between the tracking with and without segmenta- tion using the two-sided Wilcoxon signed-rank test (p = 0:06). With Segmentation Without Segmentation Evaluation Metrics Detectiony Tracking Detection Tracking Maximal error of all images 108.20 17.72 133.94 23.2 Median error of all images 0.96 0.96 0.96 0.96 Mean error of all images 5.62  15.91 1.29  1.76 9.32  19.73 1.75  3.17 Average of sequence median error 6.26  17.11 1.03  0.49 9.34  18.64 1.42  2.14 Average of sequence mean error 6.83  13.88 1.29  0.94 10.41  15.94 1.69  1.97 Figure 17: Accuracy (mm) of DCR via catheter tip tracking with and without catheter segmentation. 28 Table 12: The statistics of DCR accuracy (mm) via catheter tip tracking with and without catheter segmentation. With the two-sided Wilcoxon signed-rank test, no statistically significantly di erence is observed between the DCR with and without segmentation (p = 0:43). With Segmentation Without Segmentation All point pairs Maximal distance 20.24 25.20 Median distance 1.43 1.43 Mean distance 2.07  2.08 2.44  3.15 Frame mean distance Median distance 1.42 1.40 Average distance 1.91  1.52 2.23  2.59

Journal

Electrical Engineering and Systems SciencearXiv (Cornell University)

Published: Jan 11, 2020

References