Pulmonary Vessel Tree Matching for Quantifying Changes in Vascular Morphology

Abstract

Invasive right-sided heart catheterization (RHC) is currently the gold standard for assessing treatment effects in pulmonary vascular diseases, such as chronic thromboembolic pulmonary hypertension (CTEPH). Quantifying morphological changes by matching vascular trees (pre- and post-treatment) may provide a non-invasive alternative for assessing hemodynamic changes. In this work, we propose a method for quantifying morphological changes, consisting of three steps: constructing vascular trees from the detected pulmonary vessels, matching vascular trees with preserving local tree topology, and quantifying local morphological changes based on Poiseuille’s law (changes in $radius^{-4}$, $\triangle r^{-4}$). Subsequently, median and interquartile range (IQR) of all local $\triangle r^{-4}$ were calculated as global measurements for assessing morphological changes. The vascular tree matching method was validated with 10 synthetic trees and the relation between clinical RHC parameters and quantifications of morphological changes was investigated in 14 CTEPH patients, pre- and post-treatment. In the evaluation with synthetic trees, the proposed method achieved an average residual distance of $3.09 \pm 1.28$ mm, which is a substantial improvement over the coherent point drift method ($4.32 \pm 1.89$ mm) and a method with global-local topology preservation ($3.92 \pm 1.59$ mm). In the clinical evaluation, the morphological changes (IQR of $\triangle r^{-4}$) was significantly correlated with the changes in RHC examinations, $\triangle \text {sPAP}$ ($\mathrm{R}=-0.62$, p-value = 0.019) and $\triangle \text {mPAP}$ ($\mathrm{R}=-0.56$, p-value = 0.038). Quantifying morphological changes may provide a non-invasive assessment of treatment effects in CTEPH patients, consistent with hemodynamic changes from invasive RHC.

1 Introduction

Computed tomography (CT) pulmonary angiography (CTPA) is an important modality for assessing the severity and treatment effects of pulmonary vascular diseases, such as chronic thromboembolic pulmonary hypertension (CTEPH) [1]. Quantifying density changes in pulmonary vessels, by automatically comparing CTPA scans of pre- and post-treatment with image registration, can assess treatment effects of CTEPH [2]. CT measurements of pulmonary vascular morphology could reflect the severity of CTEPH disease [3]. However, invasive right-sided heart catheterization (RHC) serves as the gold standard for assessing disease severity and treatment effects of CTEPH [4], since it directly measures blood pressure at the main pulmonary artery. Quantifying morphological changes by matching pulmonary vessel trees of pre- and post-treatment CT scans may provide a non-invasive assessment of treatment effects.

Vascular tree matching can be treated as a point set registration task, in which the point sets represent the vessel trees. Myronenko et al. [5] proposed a coherent point drift (CPD) method for point sets registration based on a Gaussian mixture model (GMM), and with a regularization term for enforcing the motion coherence and preserving the global topology. The regularization is useful to constrain the global topology, however, its capacity to handle local deformation is low. Ge et al. [6] proposed a method with global-local topology preservation (GLTP), where a local topology term was used for regularization, based on a local linear embedding of the K nearest neighbors of each point. The main idea of local topology preservation is that local neighbors in the original point set should be preserved after transformation. The method works well with dense point sets in computer vision, such as data obtained from the Kinect depth sensor, however, the local topology constraint may induce errors in tree-like structures. This is because leaf points, that belong to different sub-trees, may still be located closely to each other. This method would then consider them as genuine neighbors, therefore over-regularizing the deformations of sub-trees.

In this paper, we propose, therefore, a method that preserves the local tree topology during vascular tree matching, and apply this method to quantify morphological changes of pulmonary vessel trees between pre- and post-treatment. The proposed method consists of three steps: (1) pre-processing for converting the detected vessels into tree structures; (2) vascular tree matching with geodesic paths for local tree topology preservation; and (3) quantification of vascular morphological changes on the basis of Poiseuille’s law [7]. The vascular tree matching method was validated with a synthetic data set, and a clinical data set consisting of 14 CTEPH patients, with CT scans and invasive RHC examinations before and after treatment.

2 Methods

We aim to align trees $T^x$ and $T^y$, which can be treated as a point set registration, with reference point set $\text {X}=[x_1,...,x_N]^T$ corresponding to nodes in $T^x$ and template point set $\text {Y}=[y_1,...,y_M]^T$ corresponding to those in $T^y$, $x_n, y_m \in \mathbb {R}^3$.

2.1 Vascular Tree Construction

For each CT scan, pulmonary vessels were segmented with a graph-cuts method [8]. The skeletons of the pulmonary vessels were extracted with a skeletonization method based on a distance transform [9] (‘DtfSkeletonization’ of MeVisLab), and the radius was recorded at the corresponding voxels on the skeleton. The skeletons were converted into a directed graph g. In the directed graph, an edge e from a start-node a to end-node b, is called an out-edge of node a and an in-edge of node b. The graph g was processed by stripping cyclic edges, so that each node (except for root node) has only one in-edge, and was converted to a tree T. A node, then, represents a bifurcation point or a leaf point, and an edge represents a branch. For each node, the average radius of the in-edge was calculated by iterating along the voxels on that in-edge and was assigned to the corresponding node. The pseudo-code of the algorithm for constructing vascular trees is given in Algorithm 1.

2.2 Vascular Tree Matching

In GMM-based methods, point sets $\text {X}$ and $\text {Y}$ can be registered by maximizing the likelihood function and an additional regularization term $R(\varTheta )$ where $\varTheta$ represents the deformation parameters. This framework minimizes the energy function:

$$\begin{aligned} E = -\text {log}\left( p(\text {X} ) \right) + R(\varTheta ). \end{aligned}$$

(1)

$\text {X}$ is considered to be distributed from a GMM with centroids $\text {Y}'$ and all Gaussians are equally-weighted with the same isotropic variance $\sigma ^2$, where $\text {Y}'$ is deformed from $\text {Y}$, $\text {Y}'=\text {Y}+\text {GW}$, $\text {G}$ is a Gaussian kernel matrix with elements $g_{ij}=exp(-\frac{\parallel y_i-y_j \parallel ^2}{2\beta ^2})$ and $\text {W}$ is $M \times D$ weight matrix of the Gaussian kernel. $\text {W}$ can be calculated by minimizing Eq. (1) when fitting the GMM to $\text {X}$.

The CPD [5] and GLTP [6] use this GMM framework, with different regularization terms. CPD uses a global regularization term $R_{cpd}(\text {W})=\frac{\lambda }{2}\text {Tr}(\text {W}^T\text {G}\text {W})$ and GLTP [6] adapted the regularization term by adding a term for preserving local topology:

$$\begin{aligned} R(\text {W}) = \frac{\lambda }{2}\text {Tr}(\text {W}^T\text {G}\text {W}) + \frac{\alpha }{2}\text {Tr}\{(\text {Y}+\text {G}\text {W})^T\text {M}(\text {Y}+\text {G}\text {W})\}, \end{aligned}$$

(2)

where $\text {M}$ is an $M \times M$ kernel matrix for preserving local deformation obtained by minimizing local linear cost function embedded with K nearest neighbors [10]. Instead of using K nearest neighbors, we compute a geodesic path with K connected nodes $N^g$ for local topology preservation, which is more suitable for the vascular trees’ deformation. $\text {M}=(\text {I}-\text {H})^T(\text {I}-\text {H})$, where $\text {H}_{ij}$ is calculated by minimizing: $\varPhi (\text {Y}) = \sum _i|y_i-\sum _{j\in N^g_i} \text {H}_{ij}y_j|^2$ [10]. For each node, a geodesic path is generated by iteratively searching the parent node and child node with a depth-first strategy, the pseudo-code is described in Algorithm 2. E is optimized with the EM algorithm, by minimizing its upper bound is:

$$\begin{aligned} Q = \sum _{n=1}^{N}\sum _{m=1}^{M}p^{prev}(y'_m|x_n)\dfrac{\parallel x_n - y_m-\text {G}(m,.)\text {W} \parallel ^2}{2 \sigma ^2} +R(\text {W}), \end{aligned}$$

(3)

where $p^{prev}(y'_m|x_n)=p(y'_m)p(x_n|y'_m)/p(x_n)$ is the posterior probability computed with the parameters from the previous step. For optimizing Eq. (3), the derivative of Q with respect to $\text {W}$ is:

$$\begin{aligned} \dfrac{\partial Q}{\partial \text {W}}=\dfrac{1}{\sigma ^2}\text {G}(diag(\text {P}\cdot 1)(\text {Y}+\text {GW})-\text {P}\text {X})+\lambda \text {GW} + \alpha \text {GM(Y + GW)}, \end{aligned}$$

(4)

in which $\text {P}$ is an $M\times N$ matrix with elements $p^{prev}(y'_m|x_n)$. By setting the function 4 to zero and right multiplying it by $\sigma ^2\text {G}^{-1}$, we have:

$$\begin{aligned} \lbrace diag(\text {P}\cdot 1)+\sigma ^2\lambda \text {I} + \sigma ^2 \alpha \text {MG}\rbrace \text {W} = \text {PX}-diag(\text {P}\cdot 1)\text {Y} +\sigma ^2 \alpha \text {MY}. \end{aligned}$$

(5)

In the M-step of the EM algorithm, $\text {W}$ is calculated by solving Eq. (5). In the E-step, $\text {P}$ is updated with the weight $\text {W}$. After optimizing the energy function, the matching pair C between $\text {X}$ and $\text {Y}$ can be built by searching point $x_n$ that maximizes the posterior probability, $C(m)={\mathop {\text {argmax}}\limits _{n}} \{ p(y'_m|x_n) \}$.

2.3 Quantitative Analysis

The nodes of the vascular trees in pre-treatment CT scans can be compared with those in post-treatment CT scans based on C. As the average radius of a branch is assigned to its end-node, morphological changes in each branch can be quantified based on C between vascular trees. Poiseuille’s law [7] describes the relation between the resistance (ratio between pressure difference and flow rate, $\triangle P/F$) and the radius r in a tube:

$$\begin{aligned} \triangle P/F=\frac{8\eta L}{\pi r^4}, \end{aligned}$$

(6)

where L is the length of the tube and $\eta$ is the fluid viscosity. Assuming that L and $\eta$ of a local branch do not change after treatment, its resistance changes can be estimated by the changes in $r^{-4}$ ($\triangle r^{-4}$). Thus, the morphological changes of vascular trees are quantified based on $\triangle r^{-4}$ of matched branches. The median and interquartile range (IQR) of the $\triangle r^{-4}$ are calculated over all branches and are used as global assessments of morphological changes.

3 Experiment

The method for constructing pulmonary vascular trees was implemented as a module in MeVisLab 2.7.1, the methods for matching vascular trees and quantifying morphological changes were implemented in Matlab, which is benefiting from the open source tools of CPD [5]. The experiments were performed on a local PC, with a 2.67 GHz CPU, 24 GB memory and a 64-bit Windows 7 system.

To evaluate the performance of vascular tree matching, synthetic vascular trees were obtained with a tree editing method [11]. In short, an initial tree $T^0$ with 3176 nodes was obtained from the left lung of a clinical CT scan and 10 synthetic trees $T^i, i=1,...,10$ were generated by randomly removing $30*i$ leaf nodes and deformed with Elastix using different non-rigid transformation parameters [12]. To simulate both deletions and additions, the synthetic tree $T^5$ and $T^i$ were matched with the proposed method (settings: MaxIteration = 100, $\beta =1, \lambda =3, \text {outlier}=0.05, \alpha =100, K=5$), furthermore, CPD [5] (MaxIteration = 100, $\beta =1, \lambda =3, \text {outlier}=0.05$) and GLTP [6] (MaxIteration = 100, $\beta =1, \lambda =3, \text {outlier}=0.05, \alpha =100, K=5$) were adopted for comparison. The Euclidean distance between nodes in $T^5$ and $T^i$ were calculated, based on the corresponding point pairs. The average and standard deviation (STD) of the residual distances were used for evaluation.

The quantification of morphological changes was validated with 14 CTEPH patients [2], who were treated with balloon pulmonary angioplasty (BPA), referred to the Tohoku University Hospital. All patients underwent both CTPA scans and RHC examinations, pre- and post-BPA treatment. The invasive RHC examinations, including pulmonary artery pressure (PAP, systolic, diastolic and mean; sPAP, dPAP and mPAP) and pulmonary vascular resistance (PVR), are examined at the main pulmonary artery. The RHC parameters changes ($\triangle$PAP and $\triangle$PVR) were used as reference measurements for assessing treatment effects. The morphological changes in vascular trees were quantified with the proposed method. The relation between the quantifications of morphological changes and hemodynamic changes ($\triangle \text {sPAP}$, $\triangle \text {dPAP}$, $\triangle \text {mPAP}$, $\triangle \text {PVR}$) were validated with Pearson’s correlation.

4 Results

The proposed method obtained an average residual distance of $3.09 \pm 1.28$ mm, while CPD and GLTP obtained an average distance of $4.32 \pm 1.89$ mm and $3.92 \pm 1.59$ mm, respectively. In comparison with CPD and GLTP, the proposed method achieved a substantial improvement, as shown in Fig. 1. The 3D visualization of vascular tree matching can be found in the supplement.

The relation between morphological changes in pulmonary vascular trees and changes in RHC measurements were investigated with 14 CTEPH patients. The IQR of $\triangle r^{-4}$ significantly correlated with $\triangle \text {sPAP}$ ($\mathrm{R}=-0.62$, p-value = 0.019) and $\triangle \text {mPAP} (\mathrm{R}=-0.56$, p-value = 0.038), but the median of $\triangle r^{-4}$ did not have a significant correlation with hemodynamic changes. Quantitative analysis of vascular morphological changes in two selected patients are shown in Fig. 2. Pearson’s correlation results are given in Table 1, and scatter plots are shown in Fig. 3.

Fig. 1.

Evaluation for vascular tree matching, average and STD of distance.

Table 1. Pearson’s correlation R (p-value) between morphological changes and hemodynamic changes.

Fig. 2.

Morphological changes of pulmonary vessels for two patients, patient A in the first row and B in the second row. Left column, initial position of vascular trees; middle column, matched vascular trees; right column, color-coded vascular trees, based on morphological changes (red: a large increase in $r^{-4}$; blue: a large decrease; green small changes).

Fig. 3.

Scatter plot for IQR of $\triangle r^{-4}$ against $\triangle \text {sPAP}$ and $\triangle \text {mPAP}$ (A and B are corresponding to patient A and B in Fig. 2).

5 Discussion and Conclusion

We present a method for quantifying morphological changes in pulmonary vascular trees, pre- and post-treatment, using vascular tree matching. The vascular tree matching method with geodesic paths for local topology preservation showed a better performance, in comparison with methods of CPD and GLTP. The IQR of $\triangle r^{-4}$, calculated based on Poiseuille’s law, had a significant negative correlation with the $\triangle \text {sPAP}$ and $\triangle \text {mPAP}$, which implies that a higher variation in $\triangle r^{-4}$ corresponds to a bigger treatment effect of decreasing pulmonary arterial pressure. This finding is consistent with a previous observation that a higher variation in density changes was related to bigger drop in pressure. In future work, we will focus on a more detailed validation of the vascular tree matching with manually annotated corresponding point pairs. By applying methods of artery-vein separation, the quantification of morphological changes may become more specific for CTEPH, since that is an arterial disease.

In conclusion, morphological changes can reflect hemodynamic changes, and quantifying morphological changes by matching vascular trees can provide a non-invasive assessment of treatment effects in CTEPH patients.