Computer Methods in Biomechanics and Biomedical Engineering
ISSN: 1025-5842 (Print) 1476-8259 (Online) Journal homepage: http://www.tandfonline.com/loi/gcmb20
Guidewire path determination for intravascular applications Fernando M. Cardoso & Sergio S. Furuie To cite this article: Fernando M. Cardoso & Sergio S. Furuie (2015): Guidewire path determination for intravascular applications, Computer Methods in Biomechanics and Biomedical Engineering, DOI: 10.1080/10255842.2015.1055732 To link to this article: http://dx.doi.org/10.1080/10255842.2015.1055732
Published online: 15 Jul 2015.
Submit your article to this journal
Article views: 11
View related articles
View Crossmark data
Full Terms & Conditions of access and use can be found at http://www.tandfonline.com/action/journalInformation?journalCode=gcmb20 Download by: [University of Lethbridge]
Date: 06 November 2015, At: 17:57
Computer Methods in Biomechanics and Biomedical Engineering, 2015 http://dx.doi.org/10.1080/10255842.2015.1055732
Guidewire path determination for intravascular applications Fernando M. Cardoso1* and Sergio S. Furuie1,2 Department of Telecommunication and Control Engineering, Biomedical Engineering Laboratory, School of Engineering, University of Sao Paulo, Sao Paulo, Brazil
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
(Received 27 December 2013; accepted 25 May 2015) Vascular diseases are among the major causes of death in developed countries and the treatment of those pathologies may require endovascular interventions, in which the physician utilizes guidewires and catheters through the vascular system to reach the injured vessel region. Several computational studies related to endovascular procedures are in constant development. Thus, predicting the guidewire path may be of great value for both physicians and researchers. However, attaining good accuracy and precision is still an important issue. We propose a method to simulate and predict the guidewire and catheter path inside a blood vessel based on equilibrium of a new set of forces, which leads, iteratively, to the minimum energy configuration. This technique was validated with phantoms using a B0.33 mm stainless steel guidewire and compared to other relevant methods in the literature. This method presented RMS error 0.30 mm and 0.97 mm, which represents less than 2% and 20% of the lumen diameter of the phantom, in 2D and 3D cases, respectively. The proposed technique presented better results than other methods from the literature, which were included in this work for comparison. Moreover, the algorithm presented low variation (s ¼ 0:03 mm) due to the variation of the input parameters. Therefore, even for a wide range of different parameters configuration, similar results are presented for the proposed approach, which is an important feature and makes this technique easier to work with. Since this method is based on basic physics, it is simple, intuitive, easy to learn and easy to adapt. Keywords: guidewire; path simulation; minimum energy; intravascular ultrasound
1. Introduction Vascular diseases lie among the major causes of death in industrialized countries (Fuster and Kelly 2010). Endovascular interventions may be required for treating these pathologies, which require the use of guidewires and catheters through the vascular system to reach the injured vessel region, with the guidance of fluoroscopic imaging. The guidewire is a thin and flexible wire that can be inserted into a confined or tortuous/winding space and, as its name suggests, it is used to guide the following insertion of a stiffer or bulkier instrument. The guidewire works as a pathway through which catheters, stents, and balloons must travel. Since many simulation tools involving IVUS/OCT imaging, elastography, and vessel reconstruction are in constant development (Wahle et al. 1999; Moraes and Furuie 2011; Doyley 2012), having a reliable path for the guidewire and catheter may be valuable. Moreover, the guidewire and catheter path determination may be used for animation purposes in order to make disease and treatment easier to be understood by both patients and students. Several investigators with different approaches have proposed methods for modeling the guidewire position. Cotin et al. (2005) proposed model based on finite beam elements with an optimization scheme, which used substructure decomposition. Although this model was
*Corresponding author. Email: [email protected] q 2015 Taylor & Francis
linear, it utilized an incremental approach that allowed nonlinear simulation. Based on this method, Lenoir et al. (2006) developed a simulator that considers both catheter and guidewire, as well as the interaction between them. Alderliesten et al. (2007) also modeled the guidewire as a series of rigid segments connected by joint positions and, based on Hooke’s law, used an explicit analytical solution (in the first-order approximation), which minimizes the energy per joint. Konings et al. (2003) divided the guidewire into connected beam elements, and proposed an energy function model based on Hooke’s law and theoretical mechanics on bending energy, which was minimized iteratively. Schafer et al. (2009) used a graphtheoretical method based on the principle of minimal total potential energy, which was formulated as the summation of the elastic energy of the guidewire and the energy due to the deformation of the vessel wall. Ganji and JanabiSharifi (2009) presented a nonphysical model composed of rigid links and joints to estimate the position of the catheter tip; however, it assumes that the catheter is deformed with constant curvature. Luboz et al. (2011) proposed a hybrid mass-spring model where mass particles were connected by non-bendable springs. Differently from other studies from the literature, the evaluation of their approach was not based on root mean squared error (RMSE), but in dice similarity coefficient. Wen Tang et al.
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
2
F.M. Cardoso and S.S. Furuie
(2012) presented a physically based hybrid modeling approach, in which the guidewire shaft is simulated using nonlinear elastic Cosserat rods and the shorter flexible tip is modeled using a generalized bending model. Luo et al. (2014) proposed a method similar to Wen Tang et al. (2012); however, their study emphasized the performance of their customized haptic system. We propose an iterative method based on active contours (Terzopoulos and Fleischer 1988; Lobregt and Viergever 1995) to simulate the guidewire final position. The simulated guidewire consists in a sequence of connected vertices. The algorithm relies on basic physics and applies forces on the vertices to deform the curve. After a number of iterations, the equilibrium of forces occurs and the minimum energy position is reached. This simple, effective, and fast model is presented as follows: ‘Materials and methods’ will detail the construction of the curve, the forces applied on the vertices, and the algorithm inputs. ‘Evaluations’ explains how the evaluation of accuracy and variability is performed. ‘Results’ will detail the results. The manuscript concludes with the section ‘Discussion and conclusion’, with a discussion of the outcomes.
2.
Materials and methods
The numeric guidewire is a discretized model, which consists of consecutive segments united by vertices (Figure 1). The guidewire deformation is performed by forces that are applied iteratively at the vertices. The forces may be divided into three groups according to their nature: (1) longitudinal forces, which account for maintaining the guidewire length; (2) angular forces, which are responsible for unbending the guidewire, leading to the minimum energy position; and (3) external forces, which are performed by the blood vessel that constrains the guidewire to its interior. It is easy to notice that the longitudinal and angular forces are internal; therefore, reaction forces should also be applied to the guidewire. Detailed explanation about the forces will be presented in the following sections. Although the algorithm was designed to work in 3D, the following illustrations were made planar for clarity.
2.1
First, the initial length L 0 of each segment is registered. Then, for each iteration, the instantaneous length L is computed and the difference between L and L 0 ! L governs the amplitude of f i , which has the same direction d^ i as the segment. L0i ¼ V 0i 2 V 0i21 ; Li ¼ jV i 2 V i21 j; V i 2 V i21 ; d^ i ¼ jV i 2 V i21 j ! f Li ¼ d^ i : Li 2 L0i :
Where V 0i and V i are the vertex original position and the instantaneous position, respectively. The user introduces the vertices initial position and the initial length of each segment is computed using the pair of consecutive vertices adjacent to it. During the curve deformation, the resulting force may stretch or compress a segment. If the instantaneous length of si is greater than the original length, then f Li pulls V i toward V i21 and vice versa (Figure 2(b)). On the other hand, if the instantaneous length of si is smaller than the original length, then f Li pushes V i away from V i21 and vice versa (Figure 2(c)).
2.2 Curvature forces The minimum energy position of the guidewire consists of smoothest possible curvatures. A straight structure, when bended, has internal forces that tend to restore the straight position. The curvature forces f C simulate those internal forces (Figure 3).
Longitudinal forces
The longitudinal force f L works as an elastic force that tends to keep the original length of each segment.
Figure 1.
Representation of the numeric guidewire model.
Figure 2. (a) Illustration of si with its original length and its two consecutive vertices V i21 and V i ; (b) case in which the length of si is greater than the original length and f Li is applied in V i toward V i21 and vice versa; (c) case in which the length of si is smaller than the original length and f Li is applied in V i away from V i21 and vice versa. The vertical dashed lines indicate V i and V i21 original position.
Computer Methods in Biomechanics and Biomedical Engineering Let us define curvature intensity ci by V i 2 V i21 ; d^ i ¼ jV i 2 V i21 j ^ diþ1 2 d^ i ui sin ; ¼ 2 2
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
ci ¼
ui ; p
where d^ i and d^ iþ1 are the direction (unitary vectors) of the two consecutive segments si and siþ1 and 0 # ci # 1, where 0 represents d^ i and d^ iþ1 in the same direction (straight position) and 1 represents d^ i and d^ iþ1 in the opposite direction. The angle ui may also be obtained through the inner product between the unitary vectors d^ i and d^ iþ1 ; however, the presented method showed to be less time-consuming. The force direction f^iþ , applied to vertices V i and V iþ1 , is defined by f^iþ d^ i ¼ 0 f^i d^ i £ d^ i21 ¼ 0 ; þ ^ fiþ ¼ 1 where the first equation ensures that f^iþ is perpendicular to d^ i , the second equation bounds v~ to belong to the plane defined by d^ i and d^ i21 and the last equation defines the norm of f^iþ . The system of equations above provides two solutions of f^iþ ¼ ff^a ; f^b } in opposite directions, f^a ¼ 2f^b . Hence, the angle between d^ i21 and f^a , and the angle between d^ i21 and f^b are computed. Then, the solution force related to the angle that is smaller than 908 is applied to V iþ1 and the other is applied to V i . f Cþ ¼ ~fi iþ1 þ ~fi ¼ ci f^i ) þ þ f Cþ ¼ 2~f i iþ
Figure 3. Illustration of f C , which forces the alignment of V i21 , V i and V iþ1 .
3
Analogously, the force direction f^i2 , applied to vertices V i21 and V i , is defined by f^i2 d^ i ¼ 0 f^i d^ i £ d^ iþ1 ¼ 0 2 f^i ¼ 1 2 Again, two vectors compose the solution and the one that provides an angle smaller than 908 with d^ iþ1 is applied to V i21 and the other is applied to V i . f C2 ¼ ~f i2 i21 ~fi ¼ ci f^i ) 2 2 f C2 ¼ 2~f i i2 Since the curvature of each vertex provides four forces – two at itself, one at the previous vertex, and one at the next, the total curvature force at vertex i is f Ci ¼ ci21 f^i21þ 2 ci f^iþ 2 ci f^i2 þ ciþ1 f^iþ12 :
2.3 External forces The external force is simply the force of the vessel wall bounding the guidewire path to its interior. A lumen mask is created, where the lumen is the object and its exterior is the background. Then, the distance transform is performed (Fabbri et al. 2008). If a vertex instantaneous position is outside the lumen, f Ei is calculated by E d^ ¼ pnearest 2V i i j pnearest 2V i j ; ! E E ^ f ¼ D : d i DT i where V i is the vertex position of the guidewire path that is outside the lumen, pnearest is the position of the nearest E voxel inside the lumen, d^ i is the unit vector that defines the direction of f Ei and DDT is the distance between V i and pnearest . Since the algorithm checks only the vertices position, it is possible that two neighboring vertices are inside the lumen but the connecting edge crosses the vessel wall. The user may work with a small distance between neighboring vertices, generating a large number of vertices, making the resolution of the curve adequate for checking the guidewire inside the blood vessel. Another option is to perform an interpolation to check if an edge crosses the vessel wall and to apply the forces on the connected vertices.
4
F.M. Cardoso and S.S. Furuie
2.4 Global force The user enters the guidewire initial curve (IC). Then, the ! global forces f i on each vertex and the new position of the vertex are calculated as follows: !
!
!
!
f i ¼ wL f Li þ wC f Ci þ wE f Ei !
vi ¼ vi ð1 2 dc Þ þ f i
Dt m
f Li DLi ¼ E: 0 ) f Li ¼ kL :DLi ¼ kL : Li 2 L0i ; A Li
V i ¼ V i þ vi :Dt; !
!
!
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
where wL , wC , and wE are the weights of f Li , f Ci and f Ei , respectively vi ; d c ; Dt, and m are the vertex instantaneous speed, drag coefficient, time interval between iterations, and vertex mass, respectively. The drag coefficient, 0 # dc # 1, controls the velocity conservation – inertia – from one iteration to the next. If dc ¼ 0, the curve moves faster, but may become unstable or enter a resonance. If dc ¼ 1, the system becomes quasi-static, the curve deformation is more steady and slower. Table 1 organizes the algorithm inputs with their abbreviations.
2.5 Boundary conditions The boundary conditions, bc, consist of a set of vertices that should remain static during the curve deformation. If the user defines only the first vertex as bc, the insertion position will remain unchanged during the curve deformation. If the user declares the first two vertices as bc, then the insertion position and direction of the first segment will remain constant. A curve with no boundary conditions would represent a loose wire inside the vessel and the user may select any set of vertices – besides the first two vertices – according to his needs.
3. Theoretical motivation The dynamic modeling of the path simulator is based on beam theory and was adapted to work with a set of consecutive straight segments. Table 1.
3.1 Longitudinal forces As explained in Section 2.1, longitudinal forces are generated by the variation of length. The force is based on Hooke’s law and opposes the variation of the segment length in order to restore it to its original value. Let si be a segment with initial length L0i , the amplitude of the longitudinal force was formulated as follows:
Algorithm parameters and respective abbreviations.
Symbol
Name
IC M Dt wL wC wE dc
Initial curve Mass Iteration time Length force weight Curvature force weight External force weight Drag coefficient
where A is the area of the cross-section of the wire, E is the Young’s modulus and Li is the instantaneous length of the segment si . ! Therefore, the longitudinal force f Li is linearly dependent on DLi and has the same direction of the segment si . Figure 4 (left) depicts the lateral view of a guidewire that is locally compressed (green region) and the internal expanding forces (green vectors) that are responsible to restore the initial length. Analogously, Figure 4 (right) illustrates the lateral view of a guidewire that is locally extended (red region) and the internal compressing forces (red vectors).
3.2 Curvature forces ! The curvature forces, f Ci , are also based on Hooke’s law. However, differently from the longitudinal case, where the deformation is pure compression or elongation throughout the whole cross-section of the wire, in the bending case, both compression and extension occur in the same crosssection (Figure 5(b)). Since, in this method, the guidewire is represented by a set of straight segments, the curvature is concentrated on the vertices. Therefore, the internal unbending forces due to a local ! curvature in the vertex V i are performed by forces (f Ci ) perpendicular to the segment in the!neighboring vertices (Figure 5(e)). The formulation of f Ci is based on the moment of the segment. As seen in Figure 5(c), both compression and extension forces are applied in the crosssection of the wire and this generates a moment (Figure 5(d)). In this method, the moment is performed by two forces applied at the distal ends of the two segments that form the curvature and their reaction forces applied at the curvature central node (Figure 5(f)). From beam theory, the flexure formula is given by
s¼2
M:y ; I
where s is the bending stress, M is the bending moment, y is the distance from the neutral axis, and I is the area moment of inertia about the neutral axis
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
Computer Methods in Biomechanics and Biomedical Engineering
5
Figure 4. Left, sequence that illustrates the lateral view of a beam with a local compression (green region) and the internal expanding forces (green vectors). Right, sequence that illustrates the lateral view of a beam with a local stretch (red region) and the internal compressing forces (red vectors).
The moment from the flexure formula can be rewritten as: M flexure ¼ 2
s:I 1ðyÞ:E:I ¼2 : y y
!
Therefore, the curvature force f Ci is linearly dependent on the local angle of curvature u. As seen in Figure 6(b), u is also the angle between the directions of two adjacent segments.
3.3 Figure 6(a) depicts a curved beam section that is related to the region between the midpoints of the two segments adjacent to V i .AB ¼ r:u is the neutral axis, which is neither compressed nor extended, and A0 B0 ¼ ðr 2 yÞ:u is an off-axis segment. Hence, the bending strain is given by 1ðyÞ ¼
DL A0 B0 2 AB ðr 2 yÞ:u 2 r:u y ¼ ¼2 ; ¼ L r:u r AB
where r is the radius of curvature and u is the angle of curvature. Then, M flexure can be rewritten as: M flexure ¼
E:I u:E:I ¼ : r AB !
The moment M f generated from the force f Ci applied at the distal ends of the adjacent segments is !
M f ¼ f Ci :Li : If we equate M flexure ¼ M f , we have !
f Ci ¼
! u:E:I ) f Ci ¼ kc :u: AB:Li
External forces
The external force simply represents the normal force due to the contact between the guidewire and the blood vessel. If a vertex is momentarily located outside the lumen, then the vessel wall ‘pushes’ the vertex back to the lumen interior. As seen in Section 2.3, the direction of the external force is estimated used distance transform.
3.4 Global force As seen in Section 2.4, the user defines the weight of each force. This enables him to control the trade-off between speed and stability. It is possible to determine the catheter path even if the mechanical properties of the guidewire – such as Young’s modulus and area moment of inertia – are unknown. Figure 7 shows a catheter that suffers deflection due to the presence of an obstacle. This situation is similar to a cantilevered beam with end load (Figure 7). From beam theory, the deflection at the end point is given by wend ¼
P :2:L 3 : 6:E:I
In this work, the deflection is determined by the obstacle; hence, the deflections of the points of the
6
F.M. Cardoso and S.S. Furuie
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
Figure 6. (a) Section of a bended beam with compressed and stretched regions in green and red, respectively. (b) The quadrilateral ABCD shows that the angle between the two adjacent segments is equal to the curvature of the correspondent bended section.
Figure 5. Sequence that illustrates (a) a curved guidewire and (b) the lateral view of the corresponding beam with a local curvature, which generates compression (green region) and expansion (red region) in the same cross-section. (c) The internal expanding and compressing forces (green and red vectors, respectively) generate (d) a moment, which is remodeled as (e) forces perpendicular to the longitudinal axis of the beam in the curvature plane.
4. Evaluation of accuracy Trajectories of physical models were utilized as a reference for accuracy assessment of the proposed path determination. The models consist of a stainless steel guidewire ðf 0:33 mmÞ performing different curves due to different paths or obstacles. In all the models, the guidewire insertion length was 15 cm. A photography of each model was taken using a digital still camera Sony DSC-HX100V (Figure 8). The numeric models were built using the images, whereby the obstacles were manually segmented. In models A, B, C, and D, the boundary conditions consisted of the first (proximal) two vertices, defining initial position and direction. On the other hand, in model E, the boundary condition was just the first vertex, defining solely the initial position. A physical 3D phantom (model F) was also considered as a reference for path comparison. A PVC clear tube with internal diameter (lumen) 5 mm was used to mimic the blood vessel. Two orthogonal bi-planar images were acquired with the guidewire previously inserted. Then, the guidewire and the central line of the PVC clear tube were manually segmented in both images in order to render the 3D phantom. Figure 9 depicts the computational rendering of a 3D phantom. Finally, numeric estimation of trajectory was carried out and, after the deformation, the consecutive vertices were
guidewire may be calculated as follows:
w ð xÞ ¼
P wend 2 :x 2 :ð3:L 2 xÞ ¼ :x :ð3:L 2 xÞ: 6:E:I 2:L 3
In other words, since the deflection is determined by the position of the obstacle, the guidewire path does not depend on the Young’s modulus nor area moment of inertia.
Figure 7. Relation between guidewire path (top) and cantilevered beam (bottom).
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
Computer Methods in Biomechanics and Biomedical Engineering
Figure 8.
7
Physical models with guidewire and obstacles in different positions.
converted into a full curve (Figure 10) through cubic spline. Hence, an image of the determined path was obtained. The metrics utilized to evaluate the simulation performance were RMSE, Hausdorff distance (HD), and Hausdorff mean (HM). The following sections will explain each metric.
is the pixel position of the simulated curve, where psim i pref is the position of closest pixel to psim that belongs to i i the real curve, and N is the number of pixels of the simulated curve.
4.2 4.1 Root mean squared error In order to measure how close the simulated curve is to the real one, we calculated the RMSE as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 PN sim 2 pref i i¼1 pi RMSE ¼ ; N
Figure 9. Computational rendering of the 3D phantom. The blue line represents the guidewire. Axes scale in voxels, where the voxel side is 0.2 mm.
Hausdorff distance and Hausdorff mean
HD and HM also measure the accuracy of the simulated curve. ; HD ¼ max psim 2 pref i i
HM ¼
i ¼ 1; 2; . . . ; N;
PN sim p 2 pref i¼1
i
N
i
;
Figure 10. Top, vertices of the final curve. Bottom, full curve.
8
F.M. Cardoso and S.S. Furuie
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
Table 2. Parameters compared in this study, abbreviations, and parameters used. Symbol Name
Basic value
Tested values
IC m
Initial curve Mass
Figure 11(a) 20
Dt
Iteration time
1
wL
Length force weight
1
wC
Curvature force weight
2
wE
External force weight
1
dc
Drag coefficient
Figure 11 5, 10, 20, 40, 80 0.25, 0.5, 1, 2, 4 0.25, 0.5, 1, 2, 4 0.5, 1, 2, 4, 8 0.25, 0.5, 1, 2, 4 0.01, 0.25, 0.5, 0.75, 1
0.01
where psim is the pixel position of the simulated curve, i pref is the position of closest pixel to psim that belongs to i i the real curve, and N is the number of pixels of the simulated curve.
5.1 Mean standard deviation The mean standard deviation, d, measures the sensitivity to a parameter considering the corresponding values in Table 2. : V I þ V IIi þ V III i þ · · · þ Vi V i ¼ i : sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi I 2 2 : 2 2 V i 2 V i þ V IIi 2 V i þ V III i 2 Vi þ· · · þ V i 2 Vi di ¼ :
d1 þ d2 þ d3 þ · · · þ dN ; d ¼ N
where V ðjÞ i is the position of the i-th vertex in the j-th configuration, : is the number of configurations, and V i and di are the mean and standard deviation, respectively, of the position of V i considering all configurations. N is the number of vertices of the curve and d is the mean of the standard deviation considering all vertices.
5.2 5. Evaluation of parameters dependency We analyzed the dependence of the final curve on the variation of several parameters (Table 2). A model (Figure 11) was chosen and the parameters of Table 2 were adjusted to a basic set, which was determined empirically based on physical considerations. Then, each parameter was varied to different values according to Table 2, while the others remained set to basic values, yielding : ¼ 5 configurations for each parameter. The parameter sensitivity was measured by calculating the mean standard deviation of the position of the vertices and the greatest distance, which will be explained in the following sections.
Greatest variation
First, for each vertex, the greatest distance between any two simulations is obtained, generating N distances. Then, the greatest distance among the resulting N distances is selected to represent the greatest variation, GV.
6.
Results
The proposed technique was performed in a personal computer equipped with Intelw Coree i7 CPU (@2.67 GHz), 8GB RAM memory, Windows 7, and Matlabw R2011a. The number of vertices considered for the 2D models (A– E) (Figure 8) was 50 and the computation time for 10,000 iterations was, on average, 5.95 ^ 0.21 s. As for the 3D model F (Figure 8), the curve also consisted
Figure 11. Different ICs in red, tested to calculate the dependence of the final result on curve initialization (variability). In all the cases, the curve length was 15 cm.
Computer Methods in Biomechanics and Biomedical Engineering Table 3. Model # A B C D E F
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
a
9
Accuracy results. RMSEa
HMa
HDa
0.57 0.60 1.08 0.52 0.30 0.97
0.42 0.44 0.90 0.29 0.20 0.92
1.20 1.60 2.15 1.83 1.78 2.00
All the measures are in mm.
of 50 vertices, but 10,000 iterations consumed 50.3 ^ 1.61 s. This algorithm may be optimized through parallel processing in order to reduce computation time. The stopping criterion used in this work was the sum of the norm of the displacements of all the vertices from one iteration to the next one, less than 1 mm. In other words, the displacement of each vertex from one iteration to the next was calculated, next the norm of each displacement was computed, then the sum of all the norms was calculated and, if the sum was less than 1 mm, the algorithm was stopped. Performing a path determination with the parameters set to the basic values listed in Table 2, the algorithm was finished after 40, 105 iterations (, 24 s), which is a fast performance even with nonoptimal parameters. Obviously, the total number of iterations is dependent on the choice of the parameters and, consequently, on the user’s ability and experience. If the software is implemented with the possibility of changing some parameters during the iterations, the user may be able to work with the trade-off speed –stability and to reduce computation time.
Figure 12.
Figure 13. Rendered guidewire (blue) and simulated guidewire path (green) for visual assessment of a 3D model.
Table 3 provides the results for the proposed path determination for all the six models previously described. It can be seen that the RMSE of the path was smaller than 1 mm. Figures 12 and 13 illustrate the simulated guidewire for visual assessment. The error is noticed to be small and the solution is consistent with physical models. Table 4 provides the results that measure the variation of the final curve according to the variation of the parameters configuration. It shows that the proposed solution sensitivity to these parameters is small. Specifically, the solution is robust regarding the IC.
Determined guidewire path (red) for the models mentioned in Figure 8. The curve in black was obtained experimentally.
10
F.M. Cardoso and S.S. Furuie
Table 4. Variation of the final curve according to the variation of the parameters configuration. Parameter IC m Dt wL wC wE dc
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
a
da
GVa
0.01 0.03 0.03 0.02 0.02 0.04 0.01
0.11 0.32 0.35 0.14 0.16 0.33 0.14
All the measures are in mm.
It can be seen that the variation of each parameter separately provided a mean standard deviation, d; smaller than 50 mm; and the greatest variation observed was smaller than 0.5 mm. When we analyzed the dependence of the final curve on the variation of iteration time, the value Dt ¼ 4 generated an unstable curve. Similar instability was observed with wC ¼ 8. Neither result was considered (Table 4) for the calculation of the dependence on iteration time and curvature force weight, respectively.
7.
Discussion and conclusions
As seen in Table 5, the proposed method provides low RMSE even with large lumen diameter. Unfortunately, Alderliesten et al. (2007) and Wen Tang et al. (2012) did not inform the lumen diameter in their methods, which is important to assess how far the simulation could be from the actual path. A narrow path would limit RMSE. Unlike the previous works listed in Table 5, we included theoretical models with one obstacle in our validation (Figure 8 A – D). Although the vessel-mimicking phantoms are more realistic, they constrain the Table 5.
simulated guidewire to their interior, which may limit and underestimate errors. A large number of inputs usually make a technique difficult to work with. Nonetheless, although the proposed algorithm has seven inputs, it is easy and fast to learn and to adapt. The influence of each input is quite intuitive since the algorithm involves only basic physics. And, as can be seen in Table 4, even if an input is changed in a wide range, a similar final result is presented. In Figure 12(A – D), after the obstacle, even a small difference between the predicted and gold-standard directions may cause a considerable distance between the tips of both paths. Since after the obstacles, the paths consist in straight lines, as the points get further from the obstacle, the distance between the two paths becomes larger. However, this phenomenon occurs only in situation with loose tips, which is not realistic. Such simulations are important because they enable the evaluation of the performance of the proposed technique in simple situations, with only one obstacle with different angulation. In clinical cases, the tip is constrained by the blood vessel wall, so Figures 12(E) and 13 represent a more realistic situation. In those cases, the predicted position of the catheter tip is very close to the observed with the real guidewire. In Table 5, a considerable difference can be noticed between the results from 2D and 3D phantoms. Among the reasons that may have worsened the 3D phantom results, we can point out the difficulties in rendering the phantom. Besides, the friction between the guidewire and the tube internal wall may have prevented the guidewire from assuming the minimum energy position. This algorithm may also be applied to non-static simulations with proper adaptations, e.g., guidewire insertion. Friction, rotation, and torsion, for example, are not in the
Comparison of experimental results from previous works.
Year of publication
2D (planar) or 3D (volumetric) phantom
Bgw (mm)
Blum (mm)
Konings Alderliesten
2003 2007
2D 3D
0.89 0.8
16 N/I
Schafer
2009
3D
0.33
Tang
2012
3D
Cardoso
2013
2D 3D
0.97 0.97 0.97 0.51 0.33
3.6 4.5 N/I
Author and reference
16 5
Inserted length (mm) 208 102 138 174 143 N/I
150 340
RMSE (mm)
RMSE Blum %
1.60 1.10 0.80 1.00 0.68 1.30 1.00 1.49 1.25 1.28 0.30 0.86
10 – – – 18.9 28.9 – – – – 1.9 17.4
Note: N/I: Not informed. Bcath and Blum are the diameters of the guidewire and lumen, respectively.
Guidewire material N/I Terumo Stainless steel Plastic-coated steel core Bentson straight Cook straight Cook curved Terumo stiff Stainless steel
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
Computer Methods in Biomechanics and Biomedical Engineering scope of this work. The effects of friction were studied by Alderliesten et al. (2007) and Wen Tang et al. (2012). Although both studies contain a numeric demonstration of the impact of the friction on the guidewire path, no experimental– clinical comparison or validation was performed. We proposed a technique that is able to determine the path of a guidewire inside a blood vessel. Since this method is based on basic physics, it is simple and easy to understand. The influence of the input parameters is intuitive and similar results are presented even for a wide range of different configurations (Table 4), which makes this technique easy to work with. More importantly, it was shown that this method presented high accuracy (Table 3) and outperformed previous related works (Table 5). As the technique presented in this work, the method proposed by Luboz et al. (2011) is also based on forces applied to mass particles connected by non-bendable segments. However, the formulation of the forces and the simulation of the vasculature are different. In future works, we may include length increment during iterations in order to simulate guidewire insertion. And we may also consider a curve tip with a rigid shape to simulate, for example, a guidewire with a J-tip. Consequently, guidewire rotation around its longitudinal axis would also be considered. With the mentioned future works, this technique could serve different purposes: planning the guidewire/catheter path and make intravascular treatment more predictable and safer, and being an instrument of practice for apprentice surgeons.
Disclosure statement No potential conflict of interest was reported by the author(s).
Funding This work was supported by the FAPESP [grant number 2011/01314-3].
Notes 1. 2.
Av Prof. Luciano Gualberto, Travessa 3, 158 – sala D2-06, 05508-970. Email: [email protected].
References Alderliesten T, Konings MK, Niessen WJ. 2007. Modeling friction, intrinsic curvature, and rotation of guide wires for simulation of minimally invasive vascular interventions. IEEE Trans Biomed Eng. 54(1):29– 38. doi:10.1109/TBME. 2006.886659.
11
Cotin S, Duriez C, Lenoir J, Neumann P, Dawson S. 2005. New approaches to catheter navigation for interventional radiology simulation. In: Proceedings of medical image computing and computer-assisted intervention. 2005 Oct 26 – 30, Palm Springs (USA), MICCAI; p. 534– 542. Doyley MM. 2012. Model-based elastography: a survey of approaches to the inverse elasticity problem. Phys Med Biol. 57(3):R35– R73. doi:10.1088/0031-9155/57/3/R35. Fabbri R, Costa LF, Torelli JC, Bruno OM. 2008. 2d euclidean distance transform algorithms. ACM Comp Surv. 40(1):1–44. doi:10.1145/1322432.1322434. Fuster V, Kelly BB. 2010. Promoting cardiovascular health in the developing world: a critical challenge to achieve global health. Institute of Medicine (US) Committee on Preventing the Global Epidemic of Cardiovascular Disease: Meeting the Challenges in Developing Countries. Washington (DC): National Academies Press. Ganji Y, Janabi-Sharifi F. 2009. Catheter kinematics for intracardiac navigation. IEEE Trans Biomed Eng. 56(3):621– 632. doi:10.1109/TBME.2009.2013134. Konings MK, Van de Kraats EB, Alderliesten T, Niessen WJ. 2003. Analytical guide wire motion algorithm for simulation of endovascular interventions. Med Biol Imag Comput. 41(6):689– 700. doi:10.1007/BF02349977. Lenoir J, Cotin S, Duriez C, Neumann P. 2006. Interactive physically-based simulation of catheter and guidewire. Comput Graph. 30(3):416–422. doi:10.1016/j.cag.2006.02.013. Lobregt S, Viergever MA. 1995. A discrete dynamic contour model. IEEE Trans Med Imag. 14(1):12 – 24. doi:10.1109/42. 370398. Luboz V, Zhai J, Littler P, Gould D, Hot T, Bello F. 2011. Simulation of endovascular guidewire behaviour and experimental validation. Comput Methods Biomec Biomed Eng. 14(6):515– 520. doi:10.1080/10255842.2010. 484003. Luo M, Xie H, Xie L, Cai P, Gu L. 2014. A robust and real-time vascular intervention simulation based on Kirchhoff elastic rod. Comput Med Imag Graph. 38(8):735– 743. doi:10.1016/ j.compmedimag.2014.08.002. Moraes MC, Furuie SS. 2011. Automatic coronary wall segmentation in intravascular ultrasound images using binary morphological reconstruction. Ultrasound Med Biol. 37(9):1486 – 1499. doi:10.1016/j.ultrasmedbio.2011. 05.018. Schafer S, Singh V, Noe¨l PB, Walczak AM, Xu J, Hoffmann KR. 2009. Real-time endovascular guidewire position simulation using shortest path algorithms. Int J Comput Assist Radiol Surg. 4(6):597 – 608. doi:10.1007/s11548-0090385-z. Tang W, Wan TR, How T, John NW. 2012. A stable and realtime nonlinear elastic approach to simulating guidewire and catheter insertions based on Cosserat rod. IEEE Trans Biomed Eng. 59(8):2211 –2218. doi:10.1109/TBME.2012. 2199319. Terzopoulos D, Fleischer K. 1988. Deformable models. Visual Comput. 4(6):306 – 331. doi:10.1007/BF01908877. Wahle A, Prause G, Dejong S, Sonka M. 1999. Geometrically correct 3-D reconstruction of intravascular ultrasound images by fusion with biplane angiography-methods and validation. IEEE Trans Med Imag. 18(8):686– 699. doi:10. 1109/42.796282.
ISSN: 1025-5842 (Print) 1476-8259 (Online) Journal homepage: http://www.tandfonline.com/loi/gcmb20
Guidewire path determination for intravascular applications Fernando M. Cardoso & Sergio S. Furuie To cite this article: Fernando M. Cardoso & Sergio S. Furuie (2015): Guidewire path determination for intravascular applications, Computer Methods in Biomechanics and Biomedical Engineering, DOI: 10.1080/10255842.2015.1055732 To link to this article: http://dx.doi.org/10.1080/10255842.2015.1055732
Published online: 15 Jul 2015.
Submit your article to this journal
Article views: 11
View related articles
View Crossmark data
Full Terms & Conditions of access and use can be found at http://www.tandfonline.com/action/journalInformation?journalCode=gcmb20 Download by: [University of Lethbridge]
Date: 06 November 2015, At: 17:57
Computer Methods in Biomechanics and Biomedical Engineering, 2015 http://dx.doi.org/10.1080/10255842.2015.1055732
Guidewire path determination for intravascular applications Fernando M. Cardoso1* and Sergio S. Furuie1,2 Department of Telecommunication and Control Engineering, Biomedical Engineering Laboratory, School of Engineering, University of Sao Paulo, Sao Paulo, Brazil
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
(Received 27 December 2013; accepted 25 May 2015) Vascular diseases are among the major causes of death in developed countries and the treatment of those pathologies may require endovascular interventions, in which the physician utilizes guidewires and catheters through the vascular system to reach the injured vessel region. Several computational studies related to endovascular procedures are in constant development. Thus, predicting the guidewire path may be of great value for both physicians and researchers. However, attaining good accuracy and precision is still an important issue. We propose a method to simulate and predict the guidewire and catheter path inside a blood vessel based on equilibrium of a new set of forces, which leads, iteratively, to the minimum energy configuration. This technique was validated with phantoms using a B0.33 mm stainless steel guidewire and compared to other relevant methods in the literature. This method presented RMS error 0.30 mm and 0.97 mm, which represents less than 2% and 20% of the lumen diameter of the phantom, in 2D and 3D cases, respectively. The proposed technique presented better results than other methods from the literature, which were included in this work for comparison. Moreover, the algorithm presented low variation (s ¼ 0:03 mm) due to the variation of the input parameters. Therefore, even for a wide range of different parameters configuration, similar results are presented for the proposed approach, which is an important feature and makes this technique easier to work with. Since this method is based on basic physics, it is simple, intuitive, easy to learn and easy to adapt. Keywords: guidewire; path simulation; minimum energy; intravascular ultrasound
1. Introduction Vascular diseases lie among the major causes of death in industrialized countries (Fuster and Kelly 2010). Endovascular interventions may be required for treating these pathologies, which require the use of guidewires and catheters through the vascular system to reach the injured vessel region, with the guidance of fluoroscopic imaging. The guidewire is a thin and flexible wire that can be inserted into a confined or tortuous/winding space and, as its name suggests, it is used to guide the following insertion of a stiffer or bulkier instrument. The guidewire works as a pathway through which catheters, stents, and balloons must travel. Since many simulation tools involving IVUS/OCT imaging, elastography, and vessel reconstruction are in constant development (Wahle et al. 1999; Moraes and Furuie 2011; Doyley 2012), having a reliable path for the guidewire and catheter may be valuable. Moreover, the guidewire and catheter path determination may be used for animation purposes in order to make disease and treatment easier to be understood by both patients and students. Several investigators with different approaches have proposed methods for modeling the guidewire position. Cotin et al. (2005) proposed model based on finite beam elements with an optimization scheme, which used substructure decomposition. Although this model was
*Corresponding author. Email: [email protected] q 2015 Taylor & Francis
linear, it utilized an incremental approach that allowed nonlinear simulation. Based on this method, Lenoir et al. (2006) developed a simulator that considers both catheter and guidewire, as well as the interaction between them. Alderliesten et al. (2007) also modeled the guidewire as a series of rigid segments connected by joint positions and, based on Hooke’s law, used an explicit analytical solution (in the first-order approximation), which minimizes the energy per joint. Konings et al. (2003) divided the guidewire into connected beam elements, and proposed an energy function model based on Hooke’s law and theoretical mechanics on bending energy, which was minimized iteratively. Schafer et al. (2009) used a graphtheoretical method based on the principle of minimal total potential energy, which was formulated as the summation of the elastic energy of the guidewire and the energy due to the deformation of the vessel wall. Ganji and JanabiSharifi (2009) presented a nonphysical model composed of rigid links and joints to estimate the position of the catheter tip; however, it assumes that the catheter is deformed with constant curvature. Luboz et al. (2011) proposed a hybrid mass-spring model where mass particles were connected by non-bendable springs. Differently from other studies from the literature, the evaluation of their approach was not based on root mean squared error (RMSE), but in dice similarity coefficient. Wen Tang et al.
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
2
F.M. Cardoso and S.S. Furuie
(2012) presented a physically based hybrid modeling approach, in which the guidewire shaft is simulated using nonlinear elastic Cosserat rods and the shorter flexible tip is modeled using a generalized bending model. Luo et al. (2014) proposed a method similar to Wen Tang et al. (2012); however, their study emphasized the performance of their customized haptic system. We propose an iterative method based on active contours (Terzopoulos and Fleischer 1988; Lobregt and Viergever 1995) to simulate the guidewire final position. The simulated guidewire consists in a sequence of connected vertices. The algorithm relies on basic physics and applies forces on the vertices to deform the curve. After a number of iterations, the equilibrium of forces occurs and the minimum energy position is reached. This simple, effective, and fast model is presented as follows: ‘Materials and methods’ will detail the construction of the curve, the forces applied on the vertices, and the algorithm inputs. ‘Evaluations’ explains how the evaluation of accuracy and variability is performed. ‘Results’ will detail the results. The manuscript concludes with the section ‘Discussion and conclusion’, with a discussion of the outcomes.
2.
Materials and methods
The numeric guidewire is a discretized model, which consists of consecutive segments united by vertices (Figure 1). The guidewire deformation is performed by forces that are applied iteratively at the vertices. The forces may be divided into three groups according to their nature: (1) longitudinal forces, which account for maintaining the guidewire length; (2) angular forces, which are responsible for unbending the guidewire, leading to the minimum energy position; and (3) external forces, which are performed by the blood vessel that constrains the guidewire to its interior. It is easy to notice that the longitudinal and angular forces are internal; therefore, reaction forces should also be applied to the guidewire. Detailed explanation about the forces will be presented in the following sections. Although the algorithm was designed to work in 3D, the following illustrations were made planar for clarity.
2.1
First, the initial length L 0 of each segment is registered. Then, for each iteration, the instantaneous length L is computed and the difference between L and L 0 ! L governs the amplitude of f i , which has the same direction d^ i as the segment. L0i ¼ V 0i 2 V 0i21 ; Li ¼ jV i 2 V i21 j; V i 2 V i21 ; d^ i ¼ jV i 2 V i21 j ! f Li ¼ d^ i : Li 2 L0i :
Where V 0i and V i are the vertex original position and the instantaneous position, respectively. The user introduces the vertices initial position and the initial length of each segment is computed using the pair of consecutive vertices adjacent to it. During the curve deformation, the resulting force may stretch or compress a segment. If the instantaneous length of si is greater than the original length, then f Li pulls V i toward V i21 and vice versa (Figure 2(b)). On the other hand, if the instantaneous length of si is smaller than the original length, then f Li pushes V i away from V i21 and vice versa (Figure 2(c)).
2.2 Curvature forces The minimum energy position of the guidewire consists of smoothest possible curvatures. A straight structure, when bended, has internal forces that tend to restore the straight position. The curvature forces f C simulate those internal forces (Figure 3).
Longitudinal forces
The longitudinal force f L works as an elastic force that tends to keep the original length of each segment.
Figure 1.
Representation of the numeric guidewire model.
Figure 2. (a) Illustration of si with its original length and its two consecutive vertices V i21 and V i ; (b) case in which the length of si is greater than the original length and f Li is applied in V i toward V i21 and vice versa; (c) case in which the length of si is smaller than the original length and f Li is applied in V i away from V i21 and vice versa. The vertical dashed lines indicate V i and V i21 original position.
Computer Methods in Biomechanics and Biomedical Engineering Let us define curvature intensity ci by V i 2 V i21 ; d^ i ¼ jV i 2 V i21 j ^ diþ1 2 d^ i ui sin ; ¼ 2 2
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
ci ¼
ui ; p
where d^ i and d^ iþ1 are the direction (unitary vectors) of the two consecutive segments si and siþ1 and 0 # ci # 1, where 0 represents d^ i and d^ iþ1 in the same direction (straight position) and 1 represents d^ i and d^ iþ1 in the opposite direction. The angle ui may also be obtained through the inner product between the unitary vectors d^ i and d^ iþ1 ; however, the presented method showed to be less time-consuming. The force direction f^iþ , applied to vertices V i and V iþ1 , is defined by f^iþ d^ i ¼ 0 f^i d^ i £ d^ i21 ¼ 0 ; þ ^ fiþ ¼ 1 where the first equation ensures that f^iþ is perpendicular to d^ i , the second equation bounds v~ to belong to the plane defined by d^ i and d^ i21 and the last equation defines the norm of f^iþ . The system of equations above provides two solutions of f^iþ ¼ ff^a ; f^b } in opposite directions, f^a ¼ 2f^b . Hence, the angle between d^ i21 and f^a , and the angle between d^ i21 and f^b are computed. Then, the solution force related to the angle that is smaller than 908 is applied to V iþ1 and the other is applied to V i . f Cþ ¼ ~fi iþ1 þ ~fi ¼ ci f^i ) þ þ f Cþ ¼ 2~f i iþ
Figure 3. Illustration of f C , which forces the alignment of V i21 , V i and V iþ1 .
3
Analogously, the force direction f^i2 , applied to vertices V i21 and V i , is defined by f^i2 d^ i ¼ 0 f^i d^ i £ d^ iþ1 ¼ 0 2 f^i ¼ 1 2 Again, two vectors compose the solution and the one that provides an angle smaller than 908 with d^ iþ1 is applied to V i21 and the other is applied to V i . f C2 ¼ ~f i2 i21 ~fi ¼ ci f^i ) 2 2 f C2 ¼ 2~f i i2 Since the curvature of each vertex provides four forces – two at itself, one at the previous vertex, and one at the next, the total curvature force at vertex i is f Ci ¼ ci21 f^i21þ 2 ci f^iþ 2 ci f^i2 þ ciþ1 f^iþ12 :
2.3 External forces The external force is simply the force of the vessel wall bounding the guidewire path to its interior. A lumen mask is created, where the lumen is the object and its exterior is the background. Then, the distance transform is performed (Fabbri et al. 2008). If a vertex instantaneous position is outside the lumen, f Ei is calculated by E d^ ¼ pnearest 2V i i j pnearest 2V i j ; ! E E ^ f ¼ D : d i DT i where V i is the vertex position of the guidewire path that is outside the lumen, pnearest is the position of the nearest E voxel inside the lumen, d^ i is the unit vector that defines the direction of f Ei and DDT is the distance between V i and pnearest . Since the algorithm checks only the vertices position, it is possible that two neighboring vertices are inside the lumen but the connecting edge crosses the vessel wall. The user may work with a small distance between neighboring vertices, generating a large number of vertices, making the resolution of the curve adequate for checking the guidewire inside the blood vessel. Another option is to perform an interpolation to check if an edge crosses the vessel wall and to apply the forces on the connected vertices.
4
F.M. Cardoso and S.S. Furuie
2.4 Global force The user enters the guidewire initial curve (IC). Then, the ! global forces f i on each vertex and the new position of the vertex are calculated as follows: !
!
!
!
f i ¼ wL f Li þ wC f Ci þ wE f Ei !
vi ¼ vi ð1 2 dc Þ þ f i
Dt m
f Li DLi ¼ E: 0 ) f Li ¼ kL :DLi ¼ kL : Li 2 L0i ; A Li
V i ¼ V i þ vi :Dt; !
!
!
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
where wL , wC , and wE are the weights of f Li , f Ci and f Ei , respectively vi ; d c ; Dt, and m are the vertex instantaneous speed, drag coefficient, time interval between iterations, and vertex mass, respectively. The drag coefficient, 0 # dc # 1, controls the velocity conservation – inertia – from one iteration to the next. If dc ¼ 0, the curve moves faster, but may become unstable or enter a resonance. If dc ¼ 1, the system becomes quasi-static, the curve deformation is more steady and slower. Table 1 organizes the algorithm inputs with their abbreviations.
2.5 Boundary conditions The boundary conditions, bc, consist of a set of vertices that should remain static during the curve deformation. If the user defines only the first vertex as bc, the insertion position will remain unchanged during the curve deformation. If the user declares the first two vertices as bc, then the insertion position and direction of the first segment will remain constant. A curve with no boundary conditions would represent a loose wire inside the vessel and the user may select any set of vertices – besides the first two vertices – according to his needs.
3. Theoretical motivation The dynamic modeling of the path simulator is based on beam theory and was adapted to work with a set of consecutive straight segments. Table 1.
3.1 Longitudinal forces As explained in Section 2.1, longitudinal forces are generated by the variation of length. The force is based on Hooke’s law and opposes the variation of the segment length in order to restore it to its original value. Let si be a segment with initial length L0i , the amplitude of the longitudinal force was formulated as follows:
Algorithm parameters and respective abbreviations.
Symbol
Name
IC M Dt wL wC wE dc
Initial curve Mass Iteration time Length force weight Curvature force weight External force weight Drag coefficient
where A is the area of the cross-section of the wire, E is the Young’s modulus and Li is the instantaneous length of the segment si . ! Therefore, the longitudinal force f Li is linearly dependent on DLi and has the same direction of the segment si . Figure 4 (left) depicts the lateral view of a guidewire that is locally compressed (green region) and the internal expanding forces (green vectors) that are responsible to restore the initial length. Analogously, Figure 4 (right) illustrates the lateral view of a guidewire that is locally extended (red region) and the internal compressing forces (red vectors).
3.2 Curvature forces ! The curvature forces, f Ci , are also based on Hooke’s law. However, differently from the longitudinal case, where the deformation is pure compression or elongation throughout the whole cross-section of the wire, in the bending case, both compression and extension occur in the same crosssection (Figure 5(b)). Since, in this method, the guidewire is represented by a set of straight segments, the curvature is concentrated on the vertices. Therefore, the internal unbending forces due to a local ! curvature in the vertex V i are performed by forces (f Ci ) perpendicular to the segment in the!neighboring vertices (Figure 5(e)). The formulation of f Ci is based on the moment of the segment. As seen in Figure 5(c), both compression and extension forces are applied in the crosssection of the wire and this generates a moment (Figure 5(d)). In this method, the moment is performed by two forces applied at the distal ends of the two segments that form the curvature and their reaction forces applied at the curvature central node (Figure 5(f)). From beam theory, the flexure formula is given by
s¼2
M:y ; I
where s is the bending stress, M is the bending moment, y is the distance from the neutral axis, and I is the area moment of inertia about the neutral axis
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
Computer Methods in Biomechanics and Biomedical Engineering
5
Figure 4. Left, sequence that illustrates the lateral view of a beam with a local compression (green region) and the internal expanding forces (green vectors). Right, sequence that illustrates the lateral view of a beam with a local stretch (red region) and the internal compressing forces (red vectors).
The moment from the flexure formula can be rewritten as: M flexure ¼ 2
s:I 1ðyÞ:E:I ¼2 : y y
!
Therefore, the curvature force f Ci is linearly dependent on the local angle of curvature u. As seen in Figure 6(b), u is also the angle between the directions of two adjacent segments.
3.3 Figure 6(a) depicts a curved beam section that is related to the region between the midpoints of the two segments adjacent to V i .AB ¼ r:u is the neutral axis, which is neither compressed nor extended, and A0 B0 ¼ ðr 2 yÞ:u is an off-axis segment. Hence, the bending strain is given by 1ðyÞ ¼
DL A0 B0 2 AB ðr 2 yÞ:u 2 r:u y ¼ ¼2 ; ¼ L r:u r AB
where r is the radius of curvature and u is the angle of curvature. Then, M flexure can be rewritten as: M flexure ¼
E:I u:E:I ¼ : r AB !
The moment M f generated from the force f Ci applied at the distal ends of the adjacent segments is !
M f ¼ f Ci :Li : If we equate M flexure ¼ M f , we have !
f Ci ¼
! u:E:I ) f Ci ¼ kc :u: AB:Li
External forces
The external force simply represents the normal force due to the contact between the guidewire and the blood vessel. If a vertex is momentarily located outside the lumen, then the vessel wall ‘pushes’ the vertex back to the lumen interior. As seen in Section 2.3, the direction of the external force is estimated used distance transform.
3.4 Global force As seen in Section 2.4, the user defines the weight of each force. This enables him to control the trade-off between speed and stability. It is possible to determine the catheter path even if the mechanical properties of the guidewire – such as Young’s modulus and area moment of inertia – are unknown. Figure 7 shows a catheter that suffers deflection due to the presence of an obstacle. This situation is similar to a cantilevered beam with end load (Figure 7). From beam theory, the deflection at the end point is given by wend ¼
P :2:L 3 : 6:E:I
In this work, the deflection is determined by the obstacle; hence, the deflections of the points of the
6
F.M. Cardoso and S.S. Furuie
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
Figure 6. (a) Section of a bended beam with compressed and stretched regions in green and red, respectively. (b) The quadrilateral ABCD shows that the angle between the two adjacent segments is equal to the curvature of the correspondent bended section.
Figure 5. Sequence that illustrates (a) a curved guidewire and (b) the lateral view of the corresponding beam with a local curvature, which generates compression (green region) and expansion (red region) in the same cross-section. (c) The internal expanding and compressing forces (green and red vectors, respectively) generate (d) a moment, which is remodeled as (e) forces perpendicular to the longitudinal axis of the beam in the curvature plane.
4. Evaluation of accuracy Trajectories of physical models were utilized as a reference for accuracy assessment of the proposed path determination. The models consist of a stainless steel guidewire ðf 0:33 mmÞ performing different curves due to different paths or obstacles. In all the models, the guidewire insertion length was 15 cm. A photography of each model was taken using a digital still camera Sony DSC-HX100V (Figure 8). The numeric models were built using the images, whereby the obstacles were manually segmented. In models A, B, C, and D, the boundary conditions consisted of the first (proximal) two vertices, defining initial position and direction. On the other hand, in model E, the boundary condition was just the first vertex, defining solely the initial position. A physical 3D phantom (model F) was also considered as a reference for path comparison. A PVC clear tube with internal diameter (lumen) 5 mm was used to mimic the blood vessel. Two orthogonal bi-planar images were acquired with the guidewire previously inserted. Then, the guidewire and the central line of the PVC clear tube were manually segmented in both images in order to render the 3D phantom. Figure 9 depicts the computational rendering of a 3D phantom. Finally, numeric estimation of trajectory was carried out and, after the deformation, the consecutive vertices were
guidewire may be calculated as follows:
w ð xÞ ¼
P wend 2 :x 2 :ð3:L 2 xÞ ¼ :x :ð3:L 2 xÞ: 6:E:I 2:L 3
In other words, since the deflection is determined by the position of the obstacle, the guidewire path does not depend on the Young’s modulus nor area moment of inertia.
Figure 7. Relation between guidewire path (top) and cantilevered beam (bottom).
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
Computer Methods in Biomechanics and Biomedical Engineering
Figure 8.
7
Physical models with guidewire and obstacles in different positions.
converted into a full curve (Figure 10) through cubic spline. Hence, an image of the determined path was obtained. The metrics utilized to evaluate the simulation performance were RMSE, Hausdorff distance (HD), and Hausdorff mean (HM). The following sections will explain each metric.
is the pixel position of the simulated curve, where psim i pref is the position of closest pixel to psim that belongs to i i the real curve, and N is the number of pixels of the simulated curve.
4.2 4.1 Root mean squared error In order to measure how close the simulated curve is to the real one, we calculated the RMSE as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 PN sim 2 pref i i¼1 pi RMSE ¼ ; N
Figure 9. Computational rendering of the 3D phantom. The blue line represents the guidewire. Axes scale in voxels, where the voxel side is 0.2 mm.
Hausdorff distance and Hausdorff mean
HD and HM also measure the accuracy of the simulated curve. ; HD ¼ max psim 2 pref i i
HM ¼
i ¼ 1; 2; . . . ; N;
PN sim p 2 pref i¼1
i
N
i
;
Figure 10. Top, vertices of the final curve. Bottom, full curve.
8
F.M. Cardoso and S.S. Furuie
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
Table 2. Parameters compared in this study, abbreviations, and parameters used. Symbol Name
Basic value
Tested values
IC m
Initial curve Mass
Figure 11(a) 20
Dt
Iteration time
1
wL
Length force weight
1
wC
Curvature force weight
2
wE
External force weight
1
dc
Drag coefficient
Figure 11 5, 10, 20, 40, 80 0.25, 0.5, 1, 2, 4 0.25, 0.5, 1, 2, 4 0.5, 1, 2, 4, 8 0.25, 0.5, 1, 2, 4 0.01, 0.25, 0.5, 0.75, 1
0.01
where psim is the pixel position of the simulated curve, i pref is the position of closest pixel to psim that belongs to i i the real curve, and N is the number of pixels of the simulated curve.
5.1 Mean standard deviation The mean standard deviation, d, measures the sensitivity to a parameter considering the corresponding values in Table 2. : V I þ V IIi þ V III i þ · · · þ Vi V i ¼ i : sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi I 2 2 : 2 2 V i 2 V i þ V IIi 2 V i þ V III i 2 Vi þ· · · þ V i 2 Vi di ¼ :
d1 þ d2 þ d3 þ · · · þ dN ; d ¼ N
where V ðjÞ i is the position of the i-th vertex in the j-th configuration, : is the number of configurations, and V i and di are the mean and standard deviation, respectively, of the position of V i considering all configurations. N is the number of vertices of the curve and d is the mean of the standard deviation considering all vertices.
5.2 5. Evaluation of parameters dependency We analyzed the dependence of the final curve on the variation of several parameters (Table 2). A model (Figure 11) was chosen and the parameters of Table 2 were adjusted to a basic set, which was determined empirically based on physical considerations. Then, each parameter was varied to different values according to Table 2, while the others remained set to basic values, yielding : ¼ 5 configurations for each parameter. The parameter sensitivity was measured by calculating the mean standard deviation of the position of the vertices and the greatest distance, which will be explained in the following sections.
Greatest variation
First, for each vertex, the greatest distance between any two simulations is obtained, generating N distances. Then, the greatest distance among the resulting N distances is selected to represent the greatest variation, GV.
6.
Results
The proposed technique was performed in a personal computer equipped with Intelw Coree i7 CPU (@2.67 GHz), 8GB RAM memory, Windows 7, and Matlabw R2011a. The number of vertices considered for the 2D models (A– E) (Figure 8) was 50 and the computation time for 10,000 iterations was, on average, 5.95 ^ 0.21 s. As for the 3D model F (Figure 8), the curve also consisted
Figure 11. Different ICs in red, tested to calculate the dependence of the final result on curve initialization (variability). In all the cases, the curve length was 15 cm.
Computer Methods in Biomechanics and Biomedical Engineering Table 3. Model # A B C D E F
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
a
9
Accuracy results. RMSEa
HMa
HDa
0.57 0.60 1.08 0.52 0.30 0.97
0.42 0.44 0.90 0.29 0.20 0.92
1.20 1.60 2.15 1.83 1.78 2.00
All the measures are in mm.
of 50 vertices, but 10,000 iterations consumed 50.3 ^ 1.61 s. This algorithm may be optimized through parallel processing in order to reduce computation time. The stopping criterion used in this work was the sum of the norm of the displacements of all the vertices from one iteration to the next one, less than 1 mm. In other words, the displacement of each vertex from one iteration to the next was calculated, next the norm of each displacement was computed, then the sum of all the norms was calculated and, if the sum was less than 1 mm, the algorithm was stopped. Performing a path determination with the parameters set to the basic values listed in Table 2, the algorithm was finished after 40, 105 iterations (, 24 s), which is a fast performance even with nonoptimal parameters. Obviously, the total number of iterations is dependent on the choice of the parameters and, consequently, on the user’s ability and experience. If the software is implemented with the possibility of changing some parameters during the iterations, the user may be able to work with the trade-off speed –stability and to reduce computation time.
Figure 12.
Figure 13. Rendered guidewire (blue) and simulated guidewire path (green) for visual assessment of a 3D model.
Table 3 provides the results for the proposed path determination for all the six models previously described. It can be seen that the RMSE of the path was smaller than 1 mm. Figures 12 and 13 illustrate the simulated guidewire for visual assessment. The error is noticed to be small and the solution is consistent with physical models. Table 4 provides the results that measure the variation of the final curve according to the variation of the parameters configuration. It shows that the proposed solution sensitivity to these parameters is small. Specifically, the solution is robust regarding the IC.
Determined guidewire path (red) for the models mentioned in Figure 8. The curve in black was obtained experimentally.
10
F.M. Cardoso and S.S. Furuie
Table 4. Variation of the final curve according to the variation of the parameters configuration. Parameter IC m Dt wL wC wE dc
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
a
da
GVa
0.01 0.03 0.03 0.02 0.02 0.04 0.01
0.11 0.32 0.35 0.14 0.16 0.33 0.14
All the measures are in mm.
It can be seen that the variation of each parameter separately provided a mean standard deviation, d; smaller than 50 mm; and the greatest variation observed was smaller than 0.5 mm. When we analyzed the dependence of the final curve on the variation of iteration time, the value Dt ¼ 4 generated an unstable curve. Similar instability was observed with wC ¼ 8. Neither result was considered (Table 4) for the calculation of the dependence on iteration time and curvature force weight, respectively.
7.
Discussion and conclusions
As seen in Table 5, the proposed method provides low RMSE even with large lumen diameter. Unfortunately, Alderliesten et al. (2007) and Wen Tang et al. (2012) did not inform the lumen diameter in their methods, which is important to assess how far the simulation could be from the actual path. A narrow path would limit RMSE. Unlike the previous works listed in Table 5, we included theoretical models with one obstacle in our validation (Figure 8 A – D). Although the vessel-mimicking phantoms are more realistic, they constrain the Table 5.
simulated guidewire to their interior, which may limit and underestimate errors. A large number of inputs usually make a technique difficult to work with. Nonetheless, although the proposed algorithm has seven inputs, it is easy and fast to learn and to adapt. The influence of each input is quite intuitive since the algorithm involves only basic physics. And, as can be seen in Table 4, even if an input is changed in a wide range, a similar final result is presented. In Figure 12(A – D), after the obstacle, even a small difference between the predicted and gold-standard directions may cause a considerable distance between the tips of both paths. Since after the obstacles, the paths consist in straight lines, as the points get further from the obstacle, the distance between the two paths becomes larger. However, this phenomenon occurs only in situation with loose tips, which is not realistic. Such simulations are important because they enable the evaluation of the performance of the proposed technique in simple situations, with only one obstacle with different angulation. In clinical cases, the tip is constrained by the blood vessel wall, so Figures 12(E) and 13 represent a more realistic situation. In those cases, the predicted position of the catheter tip is very close to the observed with the real guidewire. In Table 5, a considerable difference can be noticed between the results from 2D and 3D phantoms. Among the reasons that may have worsened the 3D phantom results, we can point out the difficulties in rendering the phantom. Besides, the friction between the guidewire and the tube internal wall may have prevented the guidewire from assuming the minimum energy position. This algorithm may also be applied to non-static simulations with proper adaptations, e.g., guidewire insertion. Friction, rotation, and torsion, for example, are not in the
Comparison of experimental results from previous works.
Year of publication
2D (planar) or 3D (volumetric) phantom
Bgw (mm)
Blum (mm)
Konings Alderliesten
2003 2007
2D 3D
0.89 0.8
16 N/I
Schafer
2009
3D
0.33
Tang
2012
3D
Cardoso
2013
2D 3D
0.97 0.97 0.97 0.51 0.33
3.6 4.5 N/I
Author and reference
16 5
Inserted length (mm) 208 102 138 174 143 N/I
150 340
RMSE (mm)
RMSE Blum %
1.60 1.10 0.80 1.00 0.68 1.30 1.00 1.49 1.25 1.28 0.30 0.86
10 – – – 18.9 28.9 – – – – 1.9 17.4
Note: N/I: Not informed. Bcath and Blum are the diameters of the guidewire and lumen, respectively.
Guidewire material N/I Terumo Stainless steel Plastic-coated steel core Bentson straight Cook straight Cook curved Terumo stiff Stainless steel
Downloaded by [University of Lethbridge] at 17:57 06 November 2015
Computer Methods in Biomechanics and Biomedical Engineering scope of this work. The effects of friction were studied by Alderliesten et al. (2007) and Wen Tang et al. (2012). Although both studies contain a numeric demonstration of the impact of the friction on the guidewire path, no experimental– clinical comparison or validation was performed. We proposed a technique that is able to determine the path of a guidewire inside a blood vessel. Since this method is based on basic physics, it is simple and easy to understand. The influence of the input parameters is intuitive and similar results are presented even for a wide range of different configurations (Table 4), which makes this technique easy to work with. More importantly, it was shown that this method presented high accuracy (Table 3) and outperformed previous related works (Table 5). As the technique presented in this work, the method proposed by Luboz et al. (2011) is also based on forces applied to mass particles connected by non-bendable segments. However, the formulation of the forces and the simulation of the vasculature are different. In future works, we may include length increment during iterations in order to simulate guidewire insertion. And we may also consider a curve tip with a rigid shape to simulate, for example, a guidewire with a J-tip. Consequently, guidewire rotation around its longitudinal axis would also be considered. With the mentioned future works, this technique could serve different purposes: planning the guidewire/catheter path and make intravascular treatment more predictable and safer, and being an instrument of practice for apprentice surgeons.
Disclosure statement No potential conflict of interest was reported by the author(s).
Funding This work was supported by the FAPESP [grant number 2011/01314-3].
Notes 1. 2.
Av Prof. Luciano Gualberto, Travessa 3, 158 – sala D2-06, 05508-970. Email: [email protected].
References Alderliesten T, Konings MK, Niessen WJ. 2007. Modeling friction, intrinsic curvature, and rotation of guide wires for simulation of minimally invasive vascular interventions. IEEE Trans Biomed Eng. 54(1):29– 38. doi:10.1109/TBME. 2006.886659.
11
Cotin S, Duriez C, Lenoir J, Neumann P, Dawson S. 2005. New approaches to catheter navigation for interventional radiology simulation. In: Proceedings of medical image computing and computer-assisted intervention. 2005 Oct 26 – 30, Palm Springs (USA), MICCAI; p. 534– 542. Doyley MM. 2012. Model-based elastography: a survey of approaches to the inverse elasticity problem. Phys Med Biol. 57(3):R35– R73. doi:10.1088/0031-9155/57/3/R35. Fabbri R, Costa LF, Torelli JC, Bruno OM. 2008. 2d euclidean distance transform algorithms. ACM Comp Surv. 40(1):1–44. doi:10.1145/1322432.1322434. Fuster V, Kelly BB. 2010. Promoting cardiovascular health in the developing world: a critical challenge to achieve global health. Institute of Medicine (US) Committee on Preventing the Global Epidemic of Cardiovascular Disease: Meeting the Challenges in Developing Countries. Washington (DC): National Academies Press. Ganji Y, Janabi-Sharifi F. 2009. Catheter kinematics for intracardiac navigation. IEEE Trans Biomed Eng. 56(3):621– 632. doi:10.1109/TBME.2009.2013134. Konings MK, Van de Kraats EB, Alderliesten T, Niessen WJ. 2003. Analytical guide wire motion algorithm for simulation of endovascular interventions. Med Biol Imag Comput. 41(6):689– 700. doi:10.1007/BF02349977. Lenoir J, Cotin S, Duriez C, Neumann P. 2006. Interactive physically-based simulation of catheter and guidewire. Comput Graph. 30(3):416–422. doi:10.1016/j.cag.2006.02.013. Lobregt S, Viergever MA. 1995. A discrete dynamic contour model. IEEE Trans Med Imag. 14(1):12 – 24. doi:10.1109/42. 370398. Luboz V, Zhai J, Littler P, Gould D, Hot T, Bello F. 2011. Simulation of endovascular guidewire behaviour and experimental validation. Comput Methods Biomec Biomed Eng. 14(6):515– 520. doi:10.1080/10255842.2010. 484003. Luo M, Xie H, Xie L, Cai P, Gu L. 2014. A robust and real-time vascular intervention simulation based on Kirchhoff elastic rod. Comput Med Imag Graph. 38(8):735– 743. doi:10.1016/ j.compmedimag.2014.08.002. Moraes MC, Furuie SS. 2011. Automatic coronary wall segmentation in intravascular ultrasound images using binary morphological reconstruction. Ultrasound Med Biol. 37(9):1486 – 1499. doi:10.1016/j.ultrasmedbio.2011. 05.018. Schafer S, Singh V, Noe¨l PB, Walczak AM, Xu J, Hoffmann KR. 2009. Real-time endovascular guidewire position simulation using shortest path algorithms. Int J Comput Assist Radiol Surg. 4(6):597 – 608. doi:10.1007/s11548-0090385-z. Tang W, Wan TR, How T, John NW. 2012. A stable and realtime nonlinear elastic approach to simulating guidewire and catheter insertions based on Cosserat rod. IEEE Trans Biomed Eng. 59(8):2211 –2218. doi:10.1109/TBME.2012. 2199319. Terzopoulos D, Fleischer K. 1988. Deformable models. Visual Comput. 4(6):306 – 331. doi:10.1007/BF01908877. Wahle A, Prause G, Dejong S, Sonka M. 1999. Geometrically correct 3-D reconstruction of intravascular ultrasound images by fusion with biplane angiography-methods and validation. IEEE Trans Med Imag. 18(8):686– 699. doi:10. 1109/42.796282.
