- Email: [email protected]
An improved graph matching algorithm for the spatio-temporal matching of a coronary artery 3D tree sequence
Disponible en ligne sur
ScienceDirect www.sciencedirect.com IRBM 36 (2015) 329–334
RITS 2015
An improved graph matching algorithm for the spatio-temporal matching of a coronary artery 3D tree sequence H. Feuillâtre a,b , J.-C. Nunes a,b,∗ , C. Toumoulin a,b a University of Rennes 1, Laboratoire Traitement du Signal et de l’Image (LTSI), Rennes, France b INSERM U1099, Rennes, France
Received 18 January 2015; received in revised form 1 July 2015; accepted 4 September 2015 Available online 9 October 2015
Abstract The paper describes an inexact tree-matching algorithm to register non-isomorphic 3D coronary artery trees over time. This work is carried out in the frame of the determination of the optimal viewing angles on the C-arm acquisition system for coronary percutaneous procedure. The matching method is based on association graph and maximum clique. Different similarity measures are compared, which use tree characteristics and geometric features of vascular branches. In order to take into account the topology variation between 3D vascular trees and thus improve the performance of the algorithm, we propose to insert artificial nodes in the association graph. Results show that unmatched node rate significantly decreases with the insertion of artificial nodes. © 2015 AGBM. Published by Elsevier Masson SAS. All rights reserved.
1. Introduction Both diagnostic and percutaneous vascular procedures rely on the choice of angiographic views and operator’s 3-D mental perception of the coronary arterial tree. The objective is to view as much of the coronary tree as possible, while at the same time minimizing vessel overlap and foreshortening of vessel segments in areas of interest (i.e. stenotic lesions). Indeed, improper views may result in suboptimal visualization of the coronary vessel shape, thus errors in the assessment of the severity and characteristics of this lesion and in a certain number of cases, in the choice of the material (catheter size, type of stent, etc.). Several works have been published that exploit the 3D coronary tree to compute the corresponding viewing angles on the C-arm acquisition system [1,2]. This 3D tree is either reconstructed from a sequence of 2D angiographic images acquired from a C-arm Cone Beam CT [3] or gotten from a preliminary CTA exam. The angioplasty procedure being performed under 2D observations, the determination of these op-
timal viewing angles involves taking into account the motion and the deformation of the coronary artery segments over the set of the cardiac cycle. It necessitates to compute for each vascular segment, the optimal viewing directions or for multiple vessel segments (when considering the case of bifurcations), the optimal viewing regions. As a consequence, the trajectory of each vascular branch has to be computed over the sequence. We propose to do it by using a tree matching technique. It is now well established that the global matching results on nonrigid transform using only geometric features are worse than the results of combining both geometric and topological features (such graph structure of coronary tree). This is because the motion (or deformation) is highly local such as for the heart. Combining this useful information with the skeleton graph can potentially improve coronary matching. So to take into account the local or elastic deformation of the artery segments, artificial nodes have been inserted. This paper focuses on this matching process; Section 2 describes the algorithm while Section 3 reports the results. 2. Materials and methods
* Corresponding author.
E-mail addresses: [email protected] (H. Feuillâtre), [email protected] (J.-C. Nunes). http://dx.doi.org/10.1016/j.irbm.2015.09.002 1959-0318/© 2015 AGBM. Published by Elsevier Masson SAS. All rights reserved.
We propose an inexact tree-matching algorithm to find each vascular branch in each volume of the sequence [4–7]. The pro-
330
H. Feuillâtre et al. / IRBM 36 (2015) 329–334
istics. Metzen et al. [4] suggested for instance considering the normalized path between the root and a given node. This path is described as a normalized polyline. The similarity between the polylines is evaluated on the basis of criteria such as the area difference between the curves or the deviation of the tangents at each point of the polylines. The normalization aims to make the measure independent of the orientation, the length and the absolute position of the polylines. The association graph is created by introducing all the vertices that correspond to every possible pair of nodes in the two trees Ti and Ti+1 . To reduce the computation time at this early stage, node assignment is constrained within a given perimeter around the considered node in the reference graph. The second step inserts then edges in the association α ,β β graph between two graph nodes vi,i+1 = (viα , vi+1 ) ∈ VGi,i+1 Fig. 1. Main steps of the tree matching algorithm.
α ,β
β
posed method is an adaptation of Metzen et al. [4]. It is based on an association graph approach. Let consider, we have at our disposal a dynamic 3D coronary tree sequence (including both 3D centerlines and meshes) [8]. Each volume i is associated with a phase of the cardiac cycle, i.e. a position of the heart at a ti time. Based on centerline information, a 3D tree structure is deduced for each coronary tree. Then, each remarkable point such as bifurcation point, root point and terminal point is considered as a tree node or vertex. A centerline segment located between two nodes is an edge. Thus, let be Ti = (Vi , Ei , ri ) and Ti+1 = (Vi+1 , Ei+1 , ri+1 ) two rooted coronary trees at successive phases ti and ti+1 . n The node sets are represented by Vi = {vi1 , . . . , viα , . . . , vi i } 1 , . . . , v β , . . . , v ni+1 } with n and n and Vi+1 = {vi+1 i i+1 rei+1 i+1 spectively the number of node in Vi and Vi+1 , α = {1, . . . , ni } and β = {1, . . . , ni+1 }.Ei and Ei+1 represent the sets of edges and ri and ri+1 the root nodes of Ti and Ti+1 respectively. Gi,i+1 = (VGi,i+1 , EGi,i+1 ) is the association graph to be built. It represents an intermediate data structure created from the two graphs Ti and Ti+1 . As shown in Fig. 1 after the construction of the graph Gi,i+1 , the maximum clique is extracted and finally artificial node are introduced to refine the matching results. 2.1. Association graph The association graph Gi,i+1 of the trees Ti and Ti+1 is defined in two steps (Fig. 1a) and b)). α,β β The first step inserts nodes vi,i+1 = (viα , vi+1 ) with viα ∈ Vi β
and vi+1 ∈ Vi+1 in the association graph iff the pairing of the two nodes appears consistent. This one is established on the basis of a set of unary constraints CU, each defining a similarity measure between two nodes. Each γ unary constraint α,β cuγ (vi,i+1 ) is normalized in the range [0, 1] through a selectivity parameter s(cuγ ). The value 1 indicating a strong correspondence. These constraints may be based on local properties of the nodes (hierarchical level or spatial coordinates of a node, length between the current node and the root node, size of its subtree [4,5]) or can take advantage of the centerline character-
β
and vi,i+1 = (viα , vi+1 ) ∈ VGi,i+1 with viα , viα ∈ Vi and β
vi+1 , vi+1 ∈ Vi+1 and 1 ≤ α ≤ ni , 1 ≤ α ≤ ni and 1 ≤ β ≤ ni+1 , 1 ≤ β ≤ ni+1 . An edge is added in the set EGi,i+1 if
β
β
the two paths pi (viα , viα ) ∈ Ti and pi+1 (vi+1 , vi+1 ) ∈ Ti+1 are considered as similar. A distance measure between the two paths can be computed from a set of binary constraints CB. Similarly to unary constraints, each δ binary constraint α ,β α ,β cbδ (vi,i+1 , vi,i+1 ) is normalized in the range [0, 1] through the selectivity parameter s(cbδ ). Basic binary constraints CB include the Euclidean or topological distance between the two extracted paths [4,5], the normalized polylines, the path length and curvature [7]. We added from our side, a binary constraint based on the Dynamic Time Warping (DTW) to determine the degree of similarity between two paths [9]. This dynamic algorithm computes a distance matrix based on the Euclidean distance between each point of the two paths pi and pi+1 and a cost function representing the minimum cumulated distance of their values of the neighbor points. All points of the two paths are considered and the result of the algorithm is normalized by the minimum path length. 2.2. Maximum clique Once the association graph Gi,i+1 established, we look for the complete sub-graph with maximum cardinality. This problem corresponds to search for the maximum clique C into the graph e.g. the maximum number of nodes that are connected by an edge and corresponds to maximal sub-tree isomorphism. A node corresponds to a pair of matching nodes in the two trees. Several methods are available to solve the maximum clique problem, which used either heuristic (approximate) [4] or exact algorithms [10]. In our case, we applied the exact algorithm developed by Östergård [10]. This method guarantees an optimal solution and is quite fast for small graph (with less than thousand nodes). Indeed, the size of the association graph is reduced by the selection introduced with the unary constraints CU. The maximum clique is the one which displays the higher sum over the set of binary constraints CB for all the edges.
H. Feuillâtre et al. / IRBM 36 (2015) 329–334
331
φ
φ
Fig. 2. Example of artificial node ai+1 insertion for an unmatched bifurcation node vi .
2.3. Artificial nodes insertion Our contribution lies at this level. The trees to be matched are not topologically identical: the number of bifurcation points and the length of the vascular segments vary from one tree to another. It is mainly due to the upstream processing (segmentation, reconstruction, etc.). As a result, the number of paired nodes is reduced because it corresponds to the nodes common to all the trees of the sequence. It only covers a small subset of the available nodes and a certain number of bifurcation and leaf nodes are not matched. Our objective is to refine the matching process by adding artificial nodes to take into account the topology variation between the vascular trees. These artificial nodes correspond to bifurcation points and leaf nodes that only exist in some of the vascular trees. We proceed in two times. Firstly, we draw up the list of unmatched nodes over the sequence. Then, we look for matching unpaired bifurcation nodes φ vi (corresponding to the start of spurious branches in the refφ erence tree Ti ) with an artificial node ai+1 that corresponds to a point of one centerline in the tree Ti+1 (see Fig. 2). For all φ unmatched nodes vi , we look for the nearest matched parent φ φ node mp(vi ) and the nearest matched children node mc(vi ). Their matching nodes in the successive tree Ti+1 are respecφ φ tively mp(ai+1 ) and mc(ai+1 ). If they both exist, the branch φ having the point which corresponds to the artificial node ai+1 is φ found. This branch corresponds to the path pi+1 with mp(ai+1 ) φ as source node and mc(ai+1 ) as target node. In some case, the φ children node mc(ai+1 ) cannot exist; we replace it by one leaf φ of the subtree induced by the node mp(ai+1 ). If more than one
leaf exists, we have to compare all paths to determine the most φ similar path pi+1 to the path pi with its source node mp(vi ) and φ vi as its target node. There are as many leave as paths. To determine the path pi+1 , the DTW algorithm is applied between pi and each possible paths. The difference of path lengths are tak-
ing into account: the algorithm is stopped when the last point in φ pi , corresponding to the bifurcation node vi , is reached. With this step, we know the most similar path pi+1 and also the point φ in this path corresponding to the unmatched node vi . This point φ is inserted as an artificial node ai+1 in the tree Ti+1 , and this φ
tree is updated (an edge is split). So, the artificial node ai+1 is φ
matched with the bifurcation node vi . This new pair of nodes is inserted in the best matching obtained from the maximum clique. ψ Secondly, unmatched leaf nodes vi are matched with an arψ tificial node ai+1 . We consider only the leave in Ti which has a corresponding point in a branch of the tree Ti+1 . The method used for the unmatched bifurcation node is applied. For all leaf ψ nodes, the nearest matched parent mp(vi ) is found and its corψ responding node mp(ai+1 ) in Ti+1 . We search the most similar ψ
path pi+1 beginning with the node mp(ai+1 ) to the path pi , ψ
ψ
from the source node mp(vi ) to the target node vi . There are as many possible paths as leave in the subtree induced from ψ mp(vi ). If the most similar path pi+1 is found with the DTW algorithm, the point in this path corresponding to the leaf node ψ ψ vi is also known. The artificial node ai+1 is inserted in Ti+1 , ψ
ψ
and the couple (vi , ai+1 ) added in the best matching result. 3. Results We evaluated the algorithm with a 3D sequence of ten coronary trees of the same patient. These ten volumes (3D meshes and centerlines) were acquired on a 64-slice GE Light-Speed CT coronary angiography system. The tree structure data Ti was obtained from the centerlines and nine matching between two successive trees (i.e. two successive phases) are computed. We obtained 196 nodes that we classified according to the results of the matching: node correctly matched, correct node
332
H. Feuillâtre et al. / IRBM 36 (2015) 329–334
Table 1 Results for each unary constraint cuγ with the best value of selectivity parameter s(cuγ ). MN is the number of well-matched nodes and PN is the number of all extracted potential nodes in the graph Gi,i+1 . cuγ
s(cuγ )
MN (%)
PN (%)
Runtime (ms)
cu1 cu2 cu3 cu4
10 1.5 0.3 1.3
53 (100.0) 52 (98.1) 52 (98.1) 51 (96.2)
68 (7.2) 432 (46.3) 807 (86.4) 727 (77.9)
2 3 72 255 32
without assignment (e.g. leaf node of a spurious branch), leaf node correctly matched with artificial node, bifurcation node correctly matched with artificial node, unmatched node and node incorrectly matched. We present our result in three parts: 1) comparison between the unary constraints (with the best value of selectivity parameters), 2) comparison between the binary constraints and 3) the relevance of the proposed method with the insertion of artificial nodes. 3.1. Unary constraint comparison We evaluated four unary constraints cuγ : – cu1 : spatial coordinates of the two potential nodes viα ∈ Vi β and vi+1 ∈ Vi+1 , β
– cu2 : length of the paths pi (ri , viα ) and pi+1 (ri+1 , vi+1 ), – cu3 : area between two normalized polylines defined by the paths pi and pi+1 , – cu4 : global curvature between two paths pi and pi+1 . To determine the best configuration among the set of unary constraints cuγ and the selectivity parameter values s(cuγ ), we verified if all the nodes, which had to be matched, were selected and also the number of all potential nodes. Table 1 reports the result for the four constraints as the computation time. The best one was achieved with the first constraint cu1 based on the Euclidean distance between nodes: all paired nodes were selected and a maximal score was achieved in a fast execution time. 3.2. Binary constraint comparison We compared then each binary constraint by taking into account all possible pairs of nodes in the association graph: the selection of the nodes that was performed at the first step of the algorithm (i.e. based on unary constraints) is not applied. Also four binary constraints were tested:
β
β
– cb1 : path lengths pi (viα , viα ) ∈ Ti and pi+1 (vi+1 , vi+1 ) ∈ Ti+1 , – cb2 : curvature between the two paths pi and pi+1 , – cb3 : area between the two normalized polylines derived from the paths pi and pi+1 , – cb4 : distance between the paths pi and pi+1 derived from the DTW algorithm.
Fig. 3. Results for each binary constraint cbδ with the best value of selectivity parameter s(cbδ ).
In the same way as for unary constraints, each binary constraint with different values of selectivity parameters (empirically chosen) was tested. Results are displayed in Fig. 3 for each constraint. The best results are obtained with the fourth constraint cb4 (i.e. based on our DTW criterion) with a computing time equal to about 90 sec for the nine matches. The faster time was obtained for cb1 or cb2 with about 1 sec. 3.3. Artificial node insertion results As shown in the two previous sections, the best matching results were respectively obtained with the unary constraint cu1 (Euclidean distance between tree nodes) and the binary constraint cb4 (DTW algorithm). We then combined them to evaluate the relevance of the artificial node insertion stage. Fig. 4 provides a comparison of the matching algorithm without adding artificial nodes (A1) and with artificial nodes (A2). The number of correctly matched nodes is not modified between the two options, and only improperly paired nodes as unmatched nodes were affected by this third stage. Indeed, a better score is reached in 70% of the cases. To confirm the position of artificial nodes and pairs of matched nodes, we test our method with the same parameter values on each coronary tree Ti of the sequence with respectively their modified trees Ti . The synthetic coronary tree Ti represents the tree Ti with a reduction of 25% of terminal branch length. In some case, we have the removal of segments which causes the loss of bifurcations (for example trees T1 , T2 and T5 in the Fig. 5). The previous matching test is taken as reference to evaluate the branch location of artificial nodes. For these nine tests (188 nodes), we have as results only one artificial node incorrectly located and three unmatched nodes (Table 2). The major errors of the matching algorithm are located in terminal branches. In this region, the movements of coronary arteries are highly localized and therefore same centerlines are more deformed (different length, position, curvature, etc.) in the cardiac phases. When the centerline’s curvature differs significantly, the constraint based on Euclidean distance in the DTW’s algorithm shown its limit.
H. Feuillâtre et al. / IRBM 36 (2015) 329–334
333
Fig. 4. Classification of tree nodes after the matching between two successive cardiac phases without artificial node (A1) and with artificial nodes (A2).
Fig. 5. The sequence of ten 3D coronary trees which represent successive phases ti in the cardiac cycle. Both black and orange centerlines represent the coronary tree Ti extracted from a CT. Black centerlines represent the synthetic tree Ti which corresponds to a reduction of 25% of the terminal branch length of the tree Ti and the suppression of short branch. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 2 Matching results between Ti and Ti . Trees
Correctly matched nodes (%)
Artificial nodes: Bifurcation (%)
Artificial nodes: Leaf (%)
Wrong artificial nodes (%)
Unmatched nodes (%)
Incorrectly nodes (%)
T0 − T1 T1 − T2 T2 − T3 T3 − T4 T4 − T5
14 (70) 14 (78) 16 (80) 18 (81.8) 16 (80)
2 (10) 1 (5.5) 1 (5) 1 (4.6) 2 (10)
4 (20) 1 (5.5) 1 (5) 3 (13.6) 2 (10)
0 (0) 1 (5.5) 0 (0) 0 (0) 0 (0)
0 (0) 1 (5.5) 2 (10) 0 (0) 0 (0)
0 (0) 0 (0) 0 (0) 0 (0) 0 (0)
17 (85) 16 (72.8) 21 (87.5) 16 (72.8)
1 (5) 3 (13.6) 2 (8.3) 1 (4.5)
2 (10) 3 (13.6) 1 (4.2) 5 (22.7)
0 (0) 0 (0) 0 (0) 0 (0)
0 (0) 0 (0) 0 (0) 0 (0)
0 (0) 0 (0) 0 (0) 0 (0)
T5 − T6 T6 − T7 T7 − T8 T8 − T9
4. Conclusion We proposed a modified version of Metzen’s algorithm [4] to register 3D coronary artery trees represented by their centerlines and 3D meshes. This hierarchical algorithm exploits different constraints to evaluate the similarity between nodes and edges, which are based on both geometrical features and node
characteristics. We exploited a new binary constraint based on the DTW algorithm (Euclidean distance) for matching the edges between promising assignments, and introduced a new step in the association graph matching procedure to deal with topology variation between coronary trees and demonstrated that the artificial node insertion allowed to increase the correctly matched node number of 13.8% in real CT-scan data. The
334
H. Feuillâtre et al. / IRBM 36 (2015) 329–334
major matching errors are located in terminal segments. The deformation and movement of these terminal branches are more important due to cardiac cycle. Some of these centerlines correspond to a local deformation and the similarity measure based only on the Euclidean distance in the DTW’s algorithm shows its limits in these particular cases. Acknowledgements The authors thank Marc Bedossa for his contribution and his advices all along this research project. References [1] Chen S-YJ, Carroll JD. Computer assisted coronary intervention by use of on-line 3D reconstruction and optimal view strategy. In: MICCAI. Springer; 1998. p. 377–85. [2] Kitslaar PH, Marquering HA, Jukema WJ, Koning G, Nieber M, Vossepoel AM, et al. Automated determination of optimal angiographic viewing angles for coronary artery bifurcations from CTA data. In: Medical imaging 2008: visualization, image-guided procedures, and modeling. Proceedings of SPIE, vol. 6918. 2008. p. 69181J.
[3] Hu Y, Jung M, Oukili A, Yang G, Nunes J-C, Fehrenbach J, et al. Sparse reconstruction from a limited projection number of the coronary artery tree in X-ray rotational imaging. In: ISBI. 2012. p. 804–7. [4] Metzen JH, Kröger T, Schenk A, Zidowitz S, Peitgen H-O, Jiang X. Matching of anatomical tree structures for registration of medical images. Image Vis Comput 2009;27(7):923–33. [5] Tschirren J, McLennan G, Palagyi K, Hoffman EA, Sonka M. Matching and anatomical labeling of human airway tree. IEEE Trans Med Imaging 2005;24(12):1540–7. [6] Lohe T, Kröger T, Zidowitz S, Peitgen H-O, Jiang X. Hierarchical matching of anatomical trees for medical image registration. ICMB, vol. 4901. Springer; 2007. p. 224–31. [7] Kaftan JN, Kiraly AP, Naidich DP, Novak CL. A novel multipurpose tree and path matching algorithm with application to airway trees. In: Medical imaging 2006: physiology, function, and structure from medical images, vol. 6143. 2006. p. 61430N. [8] Velut J, Lentz P-A, Boulmier D, Coatrieux J-L, Toumoulin C. Assessment of qualitative and quantitative features in coronary artery MRA. IRBM 2011;32(4):229–42. [9] Keogh EJ, Pazzani MJ. Derivative dynamic time warping. In: SIAM international conference on data mining. 2001. [10] Östergård PRJ. A fast algorithm for the maximum clique problem. Discrete Appl Math 2002;120(1–3):197–207.
ScienceDirect www.sciencedirect.com IRBM 36 (2015) 329–334
RITS 2015
An improved graph matching algorithm for the spatio-temporal matching of a coronary artery 3D tree sequence H. Feuillâtre a,b , J.-C. Nunes a,b,∗ , C. Toumoulin a,b a University of Rennes 1, Laboratoire Traitement du Signal et de l’Image (LTSI), Rennes, France b INSERM U1099, Rennes, France
Received 18 January 2015; received in revised form 1 July 2015; accepted 4 September 2015 Available online 9 October 2015
Abstract The paper describes an inexact tree-matching algorithm to register non-isomorphic 3D coronary artery trees over time. This work is carried out in the frame of the determination of the optimal viewing angles on the C-arm acquisition system for coronary percutaneous procedure. The matching method is based on association graph and maximum clique. Different similarity measures are compared, which use tree characteristics and geometric features of vascular branches. In order to take into account the topology variation between 3D vascular trees and thus improve the performance of the algorithm, we propose to insert artificial nodes in the association graph. Results show that unmatched node rate significantly decreases with the insertion of artificial nodes. © 2015 AGBM. Published by Elsevier Masson SAS. All rights reserved.
1. Introduction Both diagnostic and percutaneous vascular procedures rely on the choice of angiographic views and operator’s 3-D mental perception of the coronary arterial tree. The objective is to view as much of the coronary tree as possible, while at the same time minimizing vessel overlap and foreshortening of vessel segments in areas of interest (i.e. stenotic lesions). Indeed, improper views may result in suboptimal visualization of the coronary vessel shape, thus errors in the assessment of the severity and characteristics of this lesion and in a certain number of cases, in the choice of the material (catheter size, type of stent, etc.). Several works have been published that exploit the 3D coronary tree to compute the corresponding viewing angles on the C-arm acquisition system [1,2]. This 3D tree is either reconstructed from a sequence of 2D angiographic images acquired from a C-arm Cone Beam CT [3] or gotten from a preliminary CTA exam. The angioplasty procedure being performed under 2D observations, the determination of these op-
timal viewing angles involves taking into account the motion and the deformation of the coronary artery segments over the set of the cardiac cycle. It necessitates to compute for each vascular segment, the optimal viewing directions or for multiple vessel segments (when considering the case of bifurcations), the optimal viewing regions. As a consequence, the trajectory of each vascular branch has to be computed over the sequence. We propose to do it by using a tree matching technique. It is now well established that the global matching results on nonrigid transform using only geometric features are worse than the results of combining both geometric and topological features (such graph structure of coronary tree). This is because the motion (or deformation) is highly local such as for the heart. Combining this useful information with the skeleton graph can potentially improve coronary matching. So to take into account the local or elastic deformation of the artery segments, artificial nodes have been inserted. This paper focuses on this matching process; Section 2 describes the algorithm while Section 3 reports the results. 2. Materials and methods
* Corresponding author.
E-mail addresses: [email protected] (H. Feuillâtre), [email protected] (J.-C. Nunes). http://dx.doi.org/10.1016/j.irbm.2015.09.002 1959-0318/© 2015 AGBM. Published by Elsevier Masson SAS. All rights reserved.
We propose an inexact tree-matching algorithm to find each vascular branch in each volume of the sequence [4–7]. The pro-
330
H. Feuillâtre et al. / IRBM 36 (2015) 329–334
istics. Metzen et al. [4] suggested for instance considering the normalized path between the root and a given node. This path is described as a normalized polyline. The similarity between the polylines is evaluated on the basis of criteria such as the area difference between the curves or the deviation of the tangents at each point of the polylines. The normalization aims to make the measure independent of the orientation, the length and the absolute position of the polylines. The association graph is created by introducing all the vertices that correspond to every possible pair of nodes in the two trees Ti and Ti+1 . To reduce the computation time at this early stage, node assignment is constrained within a given perimeter around the considered node in the reference graph. The second step inserts then edges in the association α ,β β graph between two graph nodes vi,i+1 = (viα , vi+1 ) ∈ VGi,i+1 Fig. 1. Main steps of the tree matching algorithm.
α ,β
β
posed method is an adaptation of Metzen et al. [4]. It is based on an association graph approach. Let consider, we have at our disposal a dynamic 3D coronary tree sequence (including both 3D centerlines and meshes) [8]. Each volume i is associated with a phase of the cardiac cycle, i.e. a position of the heart at a ti time. Based on centerline information, a 3D tree structure is deduced for each coronary tree. Then, each remarkable point such as bifurcation point, root point and terminal point is considered as a tree node or vertex. A centerline segment located between two nodes is an edge. Thus, let be Ti = (Vi , Ei , ri ) and Ti+1 = (Vi+1 , Ei+1 , ri+1 ) two rooted coronary trees at successive phases ti and ti+1 . n The node sets are represented by Vi = {vi1 , . . . , viα , . . . , vi i } 1 , . . . , v β , . . . , v ni+1 } with n and n and Vi+1 = {vi+1 i i+1 rei+1 i+1 spectively the number of node in Vi and Vi+1 , α = {1, . . . , ni } and β = {1, . . . , ni+1 }.Ei and Ei+1 represent the sets of edges and ri and ri+1 the root nodes of Ti and Ti+1 respectively. Gi,i+1 = (VGi,i+1 , EGi,i+1 ) is the association graph to be built. It represents an intermediate data structure created from the two graphs Ti and Ti+1 . As shown in Fig. 1 after the construction of the graph Gi,i+1 , the maximum clique is extracted and finally artificial node are introduced to refine the matching results. 2.1. Association graph The association graph Gi,i+1 of the trees Ti and Ti+1 is defined in two steps (Fig. 1a) and b)). α,β β The first step inserts nodes vi,i+1 = (viα , vi+1 ) with viα ∈ Vi β
and vi+1 ∈ Vi+1 in the association graph iff the pairing of the two nodes appears consistent. This one is established on the basis of a set of unary constraints CU, each defining a similarity measure between two nodes. Each γ unary constraint α,β cuγ (vi,i+1 ) is normalized in the range [0, 1] through a selectivity parameter s(cuγ ). The value 1 indicating a strong correspondence. These constraints may be based on local properties of the nodes (hierarchical level or spatial coordinates of a node, length between the current node and the root node, size of its subtree [4,5]) or can take advantage of the centerline character-
β
and vi,i+1 = (viα , vi+1 ) ∈ VGi,i+1 with viα , viα ∈ Vi and β
vi+1 , vi+1 ∈ Vi+1 and 1 ≤ α ≤ ni , 1 ≤ α ≤ ni and 1 ≤ β ≤ ni+1 , 1 ≤ β ≤ ni+1 . An edge is added in the set EGi,i+1 if
β
β
the two paths pi (viα , viα ) ∈ Ti and pi+1 (vi+1 , vi+1 ) ∈ Ti+1 are considered as similar. A distance measure between the two paths can be computed from a set of binary constraints CB. Similarly to unary constraints, each δ binary constraint α ,β α ,β cbδ (vi,i+1 , vi,i+1 ) is normalized in the range [0, 1] through the selectivity parameter s(cbδ ). Basic binary constraints CB include the Euclidean or topological distance between the two extracted paths [4,5], the normalized polylines, the path length and curvature [7]. We added from our side, a binary constraint based on the Dynamic Time Warping (DTW) to determine the degree of similarity between two paths [9]. This dynamic algorithm computes a distance matrix based on the Euclidean distance between each point of the two paths pi and pi+1 and a cost function representing the minimum cumulated distance of their values of the neighbor points. All points of the two paths are considered and the result of the algorithm is normalized by the minimum path length. 2.2. Maximum clique Once the association graph Gi,i+1 established, we look for the complete sub-graph with maximum cardinality. This problem corresponds to search for the maximum clique C into the graph e.g. the maximum number of nodes that are connected by an edge and corresponds to maximal sub-tree isomorphism. A node corresponds to a pair of matching nodes in the two trees. Several methods are available to solve the maximum clique problem, which used either heuristic (approximate) [4] or exact algorithms [10]. In our case, we applied the exact algorithm developed by Östergård [10]. This method guarantees an optimal solution and is quite fast for small graph (with less than thousand nodes). Indeed, the size of the association graph is reduced by the selection introduced with the unary constraints CU. The maximum clique is the one which displays the higher sum over the set of binary constraints CB for all the edges.
H. Feuillâtre et al. / IRBM 36 (2015) 329–334
331
φ
φ
Fig. 2. Example of artificial node ai+1 insertion for an unmatched bifurcation node vi .
2.3. Artificial nodes insertion Our contribution lies at this level. The trees to be matched are not topologically identical: the number of bifurcation points and the length of the vascular segments vary from one tree to another. It is mainly due to the upstream processing (segmentation, reconstruction, etc.). As a result, the number of paired nodes is reduced because it corresponds to the nodes common to all the trees of the sequence. It only covers a small subset of the available nodes and a certain number of bifurcation and leaf nodes are not matched. Our objective is to refine the matching process by adding artificial nodes to take into account the topology variation between the vascular trees. These artificial nodes correspond to bifurcation points and leaf nodes that only exist in some of the vascular trees. We proceed in two times. Firstly, we draw up the list of unmatched nodes over the sequence. Then, we look for matching unpaired bifurcation nodes φ vi (corresponding to the start of spurious branches in the refφ erence tree Ti ) with an artificial node ai+1 that corresponds to a point of one centerline in the tree Ti+1 (see Fig. 2). For all φ unmatched nodes vi , we look for the nearest matched parent φ φ node mp(vi ) and the nearest matched children node mc(vi ). Their matching nodes in the successive tree Ti+1 are respecφ φ tively mp(ai+1 ) and mc(ai+1 ). If they both exist, the branch φ having the point which corresponds to the artificial node ai+1 is φ found. This branch corresponds to the path pi+1 with mp(ai+1 ) φ as source node and mc(ai+1 ) as target node. In some case, the φ children node mc(ai+1 ) cannot exist; we replace it by one leaf φ of the subtree induced by the node mp(ai+1 ). If more than one
leaf exists, we have to compare all paths to determine the most φ similar path pi+1 to the path pi with its source node mp(vi ) and φ vi as its target node. There are as many leave as paths. To determine the path pi+1 , the DTW algorithm is applied between pi and each possible paths. The difference of path lengths are tak-
ing into account: the algorithm is stopped when the last point in φ pi , corresponding to the bifurcation node vi , is reached. With this step, we know the most similar path pi+1 and also the point φ in this path corresponding to the unmatched node vi . This point φ is inserted as an artificial node ai+1 in the tree Ti+1 , and this φ
tree is updated (an edge is split). So, the artificial node ai+1 is φ
matched with the bifurcation node vi . This new pair of nodes is inserted in the best matching obtained from the maximum clique. ψ Secondly, unmatched leaf nodes vi are matched with an arψ tificial node ai+1 . We consider only the leave in Ti which has a corresponding point in a branch of the tree Ti+1 . The method used for the unmatched bifurcation node is applied. For all leaf ψ nodes, the nearest matched parent mp(vi ) is found and its corψ responding node mp(ai+1 ) in Ti+1 . We search the most similar ψ
path pi+1 beginning with the node mp(ai+1 ) to the path pi , ψ
ψ
from the source node mp(vi ) to the target node vi . There are as many possible paths as leave in the subtree induced from ψ mp(vi ). If the most similar path pi+1 is found with the DTW algorithm, the point in this path corresponding to the leaf node ψ ψ vi is also known. The artificial node ai+1 is inserted in Ti+1 , ψ
ψ
and the couple (vi , ai+1 ) added in the best matching result. 3. Results We evaluated the algorithm with a 3D sequence of ten coronary trees of the same patient. These ten volumes (3D meshes and centerlines) were acquired on a 64-slice GE Light-Speed CT coronary angiography system. The tree structure data Ti was obtained from the centerlines and nine matching between two successive trees (i.e. two successive phases) are computed. We obtained 196 nodes that we classified according to the results of the matching: node correctly matched, correct node
332
H. Feuillâtre et al. / IRBM 36 (2015) 329–334
Table 1 Results for each unary constraint cuγ with the best value of selectivity parameter s(cuγ ). MN is the number of well-matched nodes and PN is the number of all extracted potential nodes in the graph Gi,i+1 . cuγ
s(cuγ )
MN (%)
PN (%)
Runtime (ms)
cu1 cu2 cu3 cu4
10 1.5 0.3 1.3
53 (100.0) 52 (98.1) 52 (98.1) 51 (96.2)
68 (7.2) 432 (46.3) 807 (86.4) 727 (77.9)
2 3 72 255 32
without assignment (e.g. leaf node of a spurious branch), leaf node correctly matched with artificial node, bifurcation node correctly matched with artificial node, unmatched node and node incorrectly matched. We present our result in three parts: 1) comparison between the unary constraints (with the best value of selectivity parameters), 2) comparison between the binary constraints and 3) the relevance of the proposed method with the insertion of artificial nodes. 3.1. Unary constraint comparison We evaluated four unary constraints cuγ : – cu1 : spatial coordinates of the two potential nodes viα ∈ Vi β and vi+1 ∈ Vi+1 , β
– cu2 : length of the paths pi (ri , viα ) and pi+1 (ri+1 , vi+1 ), – cu3 : area between two normalized polylines defined by the paths pi and pi+1 , – cu4 : global curvature between two paths pi and pi+1 . To determine the best configuration among the set of unary constraints cuγ and the selectivity parameter values s(cuγ ), we verified if all the nodes, which had to be matched, were selected and also the number of all potential nodes. Table 1 reports the result for the four constraints as the computation time. The best one was achieved with the first constraint cu1 based on the Euclidean distance between nodes: all paired nodes were selected and a maximal score was achieved in a fast execution time. 3.2. Binary constraint comparison We compared then each binary constraint by taking into account all possible pairs of nodes in the association graph: the selection of the nodes that was performed at the first step of the algorithm (i.e. based on unary constraints) is not applied. Also four binary constraints were tested:
β
β
– cb1 : path lengths pi (viα , viα ) ∈ Ti and pi+1 (vi+1 , vi+1 ) ∈ Ti+1 , – cb2 : curvature between the two paths pi and pi+1 , – cb3 : area between the two normalized polylines derived from the paths pi and pi+1 , – cb4 : distance between the paths pi and pi+1 derived from the DTW algorithm.
Fig. 3. Results for each binary constraint cbδ with the best value of selectivity parameter s(cbδ ).
In the same way as for unary constraints, each binary constraint with different values of selectivity parameters (empirically chosen) was tested. Results are displayed in Fig. 3 for each constraint. The best results are obtained with the fourth constraint cb4 (i.e. based on our DTW criterion) with a computing time equal to about 90 sec for the nine matches. The faster time was obtained for cb1 or cb2 with about 1 sec. 3.3. Artificial node insertion results As shown in the two previous sections, the best matching results were respectively obtained with the unary constraint cu1 (Euclidean distance between tree nodes) and the binary constraint cb4 (DTW algorithm). We then combined them to evaluate the relevance of the artificial node insertion stage. Fig. 4 provides a comparison of the matching algorithm without adding artificial nodes (A1) and with artificial nodes (A2). The number of correctly matched nodes is not modified between the two options, and only improperly paired nodes as unmatched nodes were affected by this third stage. Indeed, a better score is reached in 70% of the cases. To confirm the position of artificial nodes and pairs of matched nodes, we test our method with the same parameter values on each coronary tree Ti of the sequence with respectively their modified trees Ti . The synthetic coronary tree Ti represents the tree Ti with a reduction of 25% of terminal branch length. In some case, we have the removal of segments which causes the loss of bifurcations (for example trees T1 , T2 and T5 in the Fig. 5). The previous matching test is taken as reference to evaluate the branch location of artificial nodes. For these nine tests (188 nodes), we have as results only one artificial node incorrectly located and three unmatched nodes (Table 2). The major errors of the matching algorithm are located in terminal branches. In this region, the movements of coronary arteries are highly localized and therefore same centerlines are more deformed (different length, position, curvature, etc.) in the cardiac phases. When the centerline’s curvature differs significantly, the constraint based on Euclidean distance in the DTW’s algorithm shown its limit.
H. Feuillâtre et al. / IRBM 36 (2015) 329–334
333
Fig. 4. Classification of tree nodes after the matching between two successive cardiac phases without artificial node (A1) and with artificial nodes (A2).
Fig. 5. The sequence of ten 3D coronary trees which represent successive phases ti in the cardiac cycle. Both black and orange centerlines represent the coronary tree Ti extracted from a CT. Black centerlines represent the synthetic tree Ti which corresponds to a reduction of 25% of the terminal branch length of the tree Ti and the suppression of short branch. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 2 Matching results between Ti and Ti . Trees
Correctly matched nodes (%)
Artificial nodes: Bifurcation (%)
Artificial nodes: Leaf (%)
Wrong artificial nodes (%)
Unmatched nodes (%)
Incorrectly nodes (%)
T0 − T1 T1 − T2 T2 − T3 T3 − T4 T4 − T5
14 (70) 14 (78) 16 (80) 18 (81.8) 16 (80)
2 (10) 1 (5.5) 1 (5) 1 (4.6) 2 (10)
4 (20) 1 (5.5) 1 (5) 3 (13.6) 2 (10)
0 (0) 1 (5.5) 0 (0) 0 (0) 0 (0)
0 (0) 1 (5.5) 2 (10) 0 (0) 0 (0)
0 (0) 0 (0) 0 (0) 0 (0) 0 (0)
17 (85) 16 (72.8) 21 (87.5) 16 (72.8)
1 (5) 3 (13.6) 2 (8.3) 1 (4.5)
2 (10) 3 (13.6) 1 (4.2) 5 (22.7)
0 (0) 0 (0) 0 (0) 0 (0)
0 (0) 0 (0) 0 (0) 0 (0)
0 (0) 0 (0) 0 (0) 0 (0)
T5 − T6 T6 − T7 T7 − T8 T8 − T9
4. Conclusion We proposed a modified version of Metzen’s algorithm [4] to register 3D coronary artery trees represented by their centerlines and 3D meshes. This hierarchical algorithm exploits different constraints to evaluate the similarity between nodes and edges, which are based on both geometrical features and node
characteristics. We exploited a new binary constraint based on the DTW algorithm (Euclidean distance) for matching the edges between promising assignments, and introduced a new step in the association graph matching procedure to deal with topology variation between coronary trees and demonstrated that the artificial node insertion allowed to increase the correctly matched node number of 13.8% in real CT-scan data. The
334
H. Feuillâtre et al. / IRBM 36 (2015) 329–334
major matching errors are located in terminal segments. The deformation and movement of these terminal branches are more important due to cardiac cycle. Some of these centerlines correspond to a local deformation and the similarity measure based only on the Euclidean distance in the DTW’s algorithm shows its limits in these particular cases. Acknowledgements The authors thank Marc Bedossa for his contribution and his advices all along this research project. References [1] Chen S-YJ, Carroll JD. Computer assisted coronary intervention by use of on-line 3D reconstruction and optimal view strategy. In: MICCAI. Springer; 1998. p. 377–85. [2] Kitslaar PH, Marquering HA, Jukema WJ, Koning G, Nieber M, Vossepoel AM, et al. Automated determination of optimal angiographic viewing angles for coronary artery bifurcations from CTA data. In: Medical imaging 2008: visualization, image-guided procedures, and modeling. Proceedings of SPIE, vol. 6918. 2008. p. 69181J.
[3] Hu Y, Jung M, Oukili A, Yang G, Nunes J-C, Fehrenbach J, et al. Sparse reconstruction from a limited projection number of the coronary artery tree in X-ray rotational imaging. In: ISBI. 2012. p. 804–7. [4] Metzen JH, Kröger T, Schenk A, Zidowitz S, Peitgen H-O, Jiang X. Matching of anatomical tree structures for registration of medical images. Image Vis Comput 2009;27(7):923–33. [5] Tschirren J, McLennan G, Palagyi K, Hoffman EA, Sonka M. Matching and anatomical labeling of human airway tree. IEEE Trans Med Imaging 2005;24(12):1540–7. [6] Lohe T, Kröger T, Zidowitz S, Peitgen H-O, Jiang X. Hierarchical matching of anatomical trees for medical image registration. ICMB, vol. 4901. Springer; 2007. p. 224–31. [7] Kaftan JN, Kiraly AP, Naidich DP, Novak CL. A novel multipurpose tree and path matching algorithm with application to airway trees. In: Medical imaging 2006: physiology, function, and structure from medical images, vol. 6143. 2006. p. 61430N. [8] Velut J, Lentz P-A, Boulmier D, Coatrieux J-L, Toumoulin C. Assessment of qualitative and quantitative features in coronary artery MRA. IRBM 2011;32(4):229–42. [9] Keogh EJ, Pazzani MJ. Derivative dynamic time warping. In: SIAM international conference on data mining. 2001. [10] Östergård PRJ. A fast algorithm for the maximum clique problem. Discrete Appl Math 2002;120(1–3):197–207.