Geometrical force constraint method for vessel and x-ray angiogram simulation

Abstract

This study proposes a novel geometrical force constraint method for 3-D vasculature modeling and angiographic image simulation. For this method, space filling force, gravitational force, and topological preserving force are proposed and combined for the optimization of the topology of the vascular structure. The surface covering force and surface adhesion force are constructed to drive the growth of the vasculature on any surface. According to the combination effects of the topological and surface adhering forces, a realistic vasculature can be effectively simulated on any surface. The image projection of the generated 3-D vascular structures is simulated according to the perspective projection and energy attenuation principles of X-rays. Finally, the simulated projection vasculature is fused with a predefined angiographic mask image to generate a realistic angiogram. The proposed method is evaluated on a CT image and three generally utilized surfaces. The results fully demonstrate the effectiveness and robustness of the proposed method.

Keywords

Vascular tree simulation force constraint angiogram

1 Introduction

In recent years, coronary artery disease (CAD) has become a major threat to human health because of its high mortality and abruptness [1, 2]. In 2012, the World Health Organization announced that CAD had become the major cause of death and hospital admissions worldwide [3]. CAD is mainly caused by the long-term accumulation of plaque in the inner walls of vascular structures. This accumulation leads to the narrowing of the arteries and to the reduction of the blood flow and oxygen supply of the heart. Given its high resolution and fast imaging speed, X-ray angiography remains the most widely adopted imaging modality in clinical practice for the interventional therapy and diagnosis of CAD [4]. However, this imaging method generates a perspective projection of the whole chest anatomy to the imaging plane. In this case, extensive information is lost in the acquired image. Moreover, catheters, bones, and myocardium chambers are superimposed in the grayscale image and thus interfere with vascular extraction and identification. The foreshortening of the vascular network as a result of the perspective projection makes vascular topology identification difficult [5, 6]. Hence, during the clinical diagnosis of CAD, several angiograms acquired at different imaging angles are necessary for the physicians to virtually reconstruct the 3-D volume image of the coronary arteries that run through their brain. An accurate vascular network extraction and identification in 2-D angiograms are obviously critical for the observation and treatment of CAD.

Numerous algorithms for the computer-aided diagnosis of CAD have been proposed, and much progress has been made in the fields of vasculature identification, optimal viewing angle determination, motion analysis, and 3-D reconstruction [7 –10]. Given that the centerline is the topological representation of vascular networks, the accuracy of centerline extraction seriously affects the effectiveness and robustness of the relative applications. However, the uneven perfusion of contrast agents, the distortion of imaging devices, the anatomical structure cover, the cardiac motion, and the complex vascular network make the centerline extraction of coronary arteries in angiograms the most challenging task in the computer-aided diagnosis of CAD [11 –13]. Several methods have been proposed to improve the accuracy of the enhancement, structure identification, and 3-D reconstruction of coronary arteries from angiographic images [14 –17]. Although great progress has been achieved, quantitative evaluation methods are still lacking [18 –20]. Existing algorithms can be evaluated by comparing their performance on a simulated phantom. With the accurate manipulation of the imaging parameters for vascular structures, the effectiveness of algorithms can be quantitatively evaluated.

3-D vascular network simulation, as well as dynamics calculation for corresponding vasculatures, has been extensively studied in recent years. Such simulation approaches are mainly based on stochastic algorithms and reconstruction theories. In 1993, Schreiner et al. [21] proposed the constrained constructive optimization (CCO) method to simulate vessels in a circular perfusion region. Based on the optimization principle of CCO, modeling methods for a spherical perfusion region and the surface of a hollow organ have also been proposed [22, 23]. Schwen et al. [24] extended the CCO model by constraining the generation of vascular structures with similar statistical information on different vascular segments. Yang et al. [25] optimized the angles, radii, and lengths of vascular branches to ensure the maximum space filling of the whole vasculature. Mitta et al. [26] reconstructed 3-D vessels according to the topological information of vascular trees. Although vascular tree simulation has been extensively studied, obtaining a high-fidelity image of the vascular network that satisfies the physiological features of the soft tissue of a specified human organ remains difficult. This difficulty arises from the fact that the anatomical parameters of vasculatures, especially for vessels in the capillary bed, are extremely difficult to measure. For many of the existing methods, the time consumed by the simulation increases as the 3-D vessels grow because constraints and terminal endings are continuously included in the calculation. To enhance calculation efficiency, Jurczuk et al. [27] simulated 3-D vessels by dividing the perfusion region into a number of subdomains. Hence, the generation of vessels in different subdomains can be achieved in different processors. The speed of the calculation thus experiences a tenfold increase. Although the structure and functional relationship of vascular trees have been investigated for many years, the anatomical parameters that statistically determine the space filling of vascular structures in a specified organ is still not well studied [28].

Once the vasculature in space is generated, X-ray angiograms can be simulated by imitating the principle of X-ray passing through a space object with predefined absorption coefficients. Generally, the simulation of X-ray propagation can be conducted by two methods: Voxel-driven and pixel-driven methods [29, 30]. Voxel-driven approaches usually calculate the projection of every voxel in a predefined image plane. Then, the angiogram is obtained by interpolating the neighboring pixels. As a large number of voxels are projected from different X-ray sources, the simulated angiogram may contain a considerable number of artifacts. Pixel-driven methods first calculate the penetrating voxels in the X-ray propagating path, which connects the X-ray source and the calculating pixel in the angiogram. Then, the energy of the calculating pixel in the angiogram can be effectively determined by integrating the energy attenuation in the X-ray propagating path. As the pixel-driven method needs to calculate every interaction between the X-ray and the volume data in space, several coordinates must be determined. In this case, progress becomes relatively time consuming.

The present study proposes a novel geometric constraint method for vascular modeling and angiographic image simulation. The flowchart of the framework is presented in Fig. 1. First of all, the geometric constraint needs to be obtained by vascular topological constraint and surface constraint in the phase of vascular modeling. The former constraint is used to control the topological distribution of the vascular structure, where the constraint is established by the fusion of the space filling force, gravitational force and topological preserving force. The latter constraint is employed to determine the distribution of 3-D vessels on any predefined surface. And the surface covering force and surface adhesion force are introduced to construct the surface constraint. In addition, the radius and length of each vascular segments are calculated to construct the constraints. And then the space distribution of the vasculature can be globally improved based on structure refinement which integrates the pruning function and overlapping function. An angiogram of the generated 3-D vascular structure is simulated based on the X-ray projection principle. The projection and overlapping effects of the vascular structure are simulated according to the energy attenuation principle of X-rays. The angiogram imaged at any angle can be simulated since the location of the 3D model is able to be manipulated without restraint. We also evaluate the performance of our method in terms of 3D vessel simulation, and angiogram simulation in terms of clinical CT data and simulated vasculature. The outstanding results show that the proposed method is not only capable to generate 3-D vessels on any space shape effectively, but also has the ability to visualize the shortening, overlapping, and bifurcation of the vessel on both clinical CT data and simulated 3D vessel precisely. Moreover, the angiograms of the 3-D vessels in any view angles can be simulated.

2 3-D vessel simulation

2.1 Radius and length determination

A typical vascular structure is usually composed of numerous vascular branches that bifurcate continuously from the root segment to the capillary bed. In 1926, Murray [31] proposed that in a typical vascular system, the fluid volume of the parent branch should be equal to the sum of the fluid volume of the derived daughter branches. Murray constructed a cost function to describe the energy dissipation of the blood flow that passes through the bifurcation. Hence, the following minimal energy function can be obtained as follows: $r_{0}^{k} = r_{1}^{k} + r_{2}^{k}$ (1) where k is the bifurcation exponent. Murray’s law proposes that k should be equal to 3.0 for an ideal vasculature system. While Uylings et al. [32] demonstrated that the bifurcation exponent varies between 2.33 to 3.0, depending on the model utilized to calculate the blood flow.

To simulate a realistic vasculature that follows the minimal energy dissipation principle, we first construct a vessel tree that continuously bifurcates from a root segment. The diameters of the bifurcating daughter branches are expected to satisfy Equation (1). A realistic vasculature is usually presented as an asymmetrical structure. Hence, we add a random function to simulate the asymmetry of vascular structure, and the radii can be calculated as follows: $\begin{matrix} r_{1} = 2^{- \frac{1}{k}} * G_{1} (μ_{r}, σ_{r}) * r_{0} \\ r_{2} = {(1 - \frac{1}{2} * G_{1} (μ_{r}, σ_{r})^{k})}^{- \frac{1}{k}} * r_{0} \end{matrix}$ (2) where r ₀, r ₁ and r ₂ are the diameters of the parent and the two daughter branches. G ₁ (*) varies in the form of a Gaussian function whose mean and standard deviation are μ _r and σ _r separately.

The length of vascular branches can be simply expressed by the linear function of the radius [25, 33]. Without loss of generality, we introduce a random function to combine with the radius of vascular segment to control the length of vascular segment. Thus we arrive at the following equation: $l (r) = G_{2} (μ_{l}, σ_{l}) * r$ (3) where G ₂ is the random ratio of the vascular length, it is determined by a Gaussian function whose mean and standard deviation is μ _l and σ _l respectively. For the generated vascular structure, the variation of the branch lengths reflects the variation of the radii.

2.2 Space location optimization

This study proposes to generate vascular structures that adaptively adhere to any regular surface. When the parent branch splits into two daughter branches, the locations of the daughter branches are decided by the radii and lengths of the vascular branch, the intersection angles between the daughter branches and the parent branch, and the in-plane rotational angles. Let P _{i - 1,j}, P _i, P _i,1 and P _i,2 represent the four adjacent bifurcating points (Fig. 2). P _{i - 1,j} represents a daughter branch derived from a previous parent branch, j = 1, 2. θ ₁ is the intersectional angle of branches $\vec{P_{i - 1, j} P_{i}}$ and $\vec{P_{i} P_{i, 1}}$ in the plane generated by points P _{i - 1,j}, P _i and P _i,1. θ ₂ is the intersectional angle of branches $\vec{P_{i - 1, j} P_{i}}$ and $\vec{P_{i} P_{i, 2}}$ in the plane generated by points P _{i - 1,j}, P _i, and P _i,2. η ₁ is the rotational angle of branch $\vec{P_{i} P_{i, 1}}$ around branch $\vec{P_{i - 1, j} P_{i}}$ . η ₂ is the rotational angle of branch $\vec{P_{i} P_{i, 2}}$ around branch $\vec{P_{i - 1, j} P_{i}}$ . Hence, a local coordinate system is defined based on points P _{i - 1,j}, P _i, P _i,1 and P _i,2. In this system, axis x _i is perpendicular to vector $\vec{P_{i, 1} P_{i, 2}}$ while y _i is perpendicular to the plane comprising points P _{i - 1,j}, P _i,1 and P _i,2. Hence, z _I can be computed by the cross product of the unit vectors of axes x _i and y _i. Given that the angles to be optimized are made up of the intersection angles (θ ₁, θ ₂) between the parent branch and the daughter branches and the rotational angles (η ₁, η ₂) of the daughter branches around the parent branch, the lengths and radii of the vascular branches are kept constant. Hence, the corresponding 3-D location can be achieved by optimizing the four predefined angles. Figure 2 presents the optimization model of the vascular bifurcation.

To generate the vasculature that develops in any tissue shape, this study introduces two types of constraints for the structural optimization of the vascular tree: Vascular topological constraint and surface constraint. The first constraint involves three forces, namely, space filling force, gravitational force, and topological preserving force. The second constraint involves two kinds of forces, namely, surface covering force and surface adhesion force. The space filling force is designed to control the volume that the generated vasculature fills in space. The gravitational force is devised to control the heterogeneity direction of the perfusion model. The topological preserving force is utilized to constrain the topological consistency of the generated vascular structure. The surface covering force is designed to control the range of the extension of the generated vasculature on the surface of the predefined shape. The surface adhesion force is designed to control the degree of dispersion of the generated model to the surface of the predefined surface. Detailed descriptions of the forces are presented below.

(a) Space filling force. This force pulls each newly added daughter segment away from the existing segments. For the simulation process, the space filling force needs to be maximized so that the distance between the newly added branch and all the existing branches can be as large as possible. The energy function is given by: $E_{1} = E_{1}^{'} + E_{1}^{″} = (\sum_{m = 1}^{N} D (P_{i, 1}, p_{m}) + \sum_{m = 1}^{N} D (P_{i, 2}, p_{m})) + w * D (P_{i, 1}, P_{i, 2})$ (4) where D (* , *) is the Euclidean distance between two branching points in 3-D space. N represents the number of generated branches. P _i,1 and P _i,2 represent the newly added daughter branch pairs. w is the weighted factor to balance energy $E_{1}^{'}$ and E″. Here w is taken as 2 * N.

(b) Gravitational force. This force pushes the whole vasculature to the predefined direction, which is set as parallel to the direction of the root segment. By this force, the vasculature is gradually pushed away from the root segment and then distributed discretely along the predefined direction. We have: $E_{2} = D (P_{g}, ϖ (P_{i}, P_{i, 1}, P_{i, 2}))$ (5) where ϖ (*) is the plane determined by the three bifurcation points P _i, P _i,1 and P _i,2. P _g represents the bottom of the predefined shape.

(c) Topological preserving force. This force is defined to preserve the topology continuity of the newly generated segments. At each bifurcation, the daughter segment with a large radius has smaller angle differences compared with branch with smaller radius. Hence, we have $E_{3} = arccos (\frac{\vec{P_{i - 1, j} P_{i}} * \vec{P_{i} P_{i, l}}}{| \vec{P_{i - 1, j} P_{i}} | * | \vec{P_{i} P_{i, l}} |})$ (6) where $\vec{P_{i} P_{i, l}}$ is the branch which has a larger intersectional angle with the branch $\vec{P_{i - 1, j} P_{i}}$ .

(d) Surface covering force. This force pulls the newly generated segment to cover the surface of the predefined shape as large as possible. The force can be defined as follows: $E_{4} = \sum_{m = 1}^{(M - 1) / 2} s (P_{m - 1, j}, P_{m, 1}, P_{m, 2})$ (7) where s (*) is the triangle area that comprises three bifurcation points P _m-1,j, P _m,1 and P _m,2. M represents the number of the triangle areas composed by the bifurcation points of P _m-1,j, P _m,1 and P _m,2.

(e) Surface adhesion force. This force pulls the whole vasculature so that it adheres closely to the predefined surface. We have the following equation: $E_{5} = D (P_{i, 1}, ϖ (\sum_{k} p_{k})) + D (P_{i, 2}, ϖ (\sum_{k} p_{k}))$ (8) where p _k represents the point on the structure surface while ϖ (∑_k p _k) is the surface composed by a series of points ∑_k p _k on the predefined surface.

The proposed five energies can be utilized to control the topology and space distribution of the vessel development on the predefined shape of the tissue. The constraint energy function hence can be defined as follows: $\underset{α_{1}, β_{1}, α_{2}, β_{2}}{arg min} (F_{j} (α_{1}, β_{1}, α_{2}, β_{2})) = w_{1} * 1.0 / E_{1} + w_{2} * E_{2} + w_{3} * E_{3} + w_{4} * A_{s} / E_{4} + w_{5} * E_{5}$ (9) where A _s is the area of the vessel on the structure surface. (w ₁, w ₂, w ₃, w ₄, w ₅) are the weighting coefficient of the five forces. Given that each force in the cost function may have a different order of magnitude, the determination of the weighting values must ensure that all the energies have roughly similar effects. Hence, the ratio of the weighted factor w₁ , w₂ , w₃ , w₄ and w₅ are expressed in the form of 10^t. The value of t is decided by the first three vascular segments. By incorporating the weighting coefficients, all the energies can be transformed to roughly the same scale. Then the optimization procedure is performed by the well-known Powell algorithm [34].

2.3 Structure adjustment

To improve the fidelity of the generated vascular structure, we use a cardinal spline to smoothen the generated vessels. In this study, the intersection angle between the parent branch and the daughter branch is defined as the criterion for classifying the connecting branch. Generally, the daughter branch with a large intersection angle belongs to the main branch that follows the parent branch. In simulating the detectability of the imaging device, a threshold function is defined to remove the branches with diameters of less than the predefined threshold. We have: ${\vec{P_{i} P_{i, l}} | length (\vec{P_{i} P_{i, l}}) \leq min {{len}_{t} | t = 2^{T - 1}, 2^{T - 1} + 1, \dots, 2^{T} - 1}}$ (10) where len _t is the set of clustered vascular branches, $length (\vec{P_{i} P_{i, l}})$ is the length of the connected segments, and i = 2^T-2, 2^T-2 + 1, ⋯ , 2^T-1 - 1.

The number of newly added branches increases gradually with the increase of the branching level in the proposed model. Hence, some small segments may inevitably overlap. We introduce the following collision function to adjust the vascular structure further: $D (\vec{P_{m 1} P_{n 1}}, \vec{P_{m 2} P_{n 2}}) \leq min {r_{t} | t = 1, 2 \dots, 2^{T} - 1}$ (11) where $\vec{P_{m 1} P_{n 1}}$ is the segment defined by points P _m1 and P _n1 while $\vec{P_{m 2} P_{n 2}}$ is the branch that connects points P _m2 and P _n2. r _t is the radius of the generated branch. By calculating the distance between the intersection point and the terminal point, we determine the segment with a short length and cut it off. The procedure is repeated until all the segments are adjusted.

3 Angiogram simulation

A typical angiographic image is generated by capturing the passing of X-rays through an object to reach an energy receiver after a specific path. In the angiographic imaging system, the X-ray source and receiver are separately placed in the two sides of the C-arm. Generally, the source-to-object distance (SOD) and the source-to-image distance (SID) determine the amplification ratio of the imaging object. The relative imaging orientation of the C-arm is decided by the left/right anterior oblique angles (LAO/RAO) and the caudal/cranial angles (CAU/CRA). Fig. 3 shows the angiographic imaging system.

To simulate the angiographic X-ray image, the 3-D simulated vessels are placed in the angiographic system. Let o - x _W y _W z _W, c - x _C y _C z _C and o ^- uv represent the world coordinate, camera coordinate and image coordinate, respectively. The imaging of the 3-D volume data in the X-ray transmission path can be formulated as follows: $p_{I} = M_{I} M_{C} p_{W}$ (12) where M _C is the transformation matrix from system o - x _W y _W z _W to system c - x _C y _C z _C while M _I is the transformation matrix from system c - x _C y _C z _C to system o ^- uv. p _I is the homogeneous matrix in system o ^- uv, and p _W is the homogeneous matrix in o - x _W y _W z _W. Based on the pinhole camera model, the projection matrix can be reformulated into the following homogeneous form: $p_{I} = {KQp}_{W} / z_{C}$ (13) where K and Q are the internal and external matrix, respectively. z _C is the z coordinate in the system c - x _C y _C z _C. We can finally determine the homogenous pixel coordinates in the imaging plane. The internal matrix is formulated as follows: $[K = [\begin{matrix} SID / α_{u} & 0 & u_{0} & 0 \\ 0 & SID / α_{v} & v_{0} & 0 \\ 0 & 0 & 1 & 0 \end{matrix}]$ (14)

The external matrix can be formulated by the rotational angles of the C-arm: $Q = [\begin{matrix} - cos β & sin α sin β & - cos α sin β & SOD sin β \\ 0 & cos α & sin α & - SOD sin α \\ sin β & sin α cos β & - cos α cos β & SOD \sqrt{1 - {sin}^{2} α - {sin}^{2} β} \\ 0 & 0 & 0 & 1 \end{matrix}]$ (15) where (α _u,α _v) represent the resolution of the angiograms. (u ₀,v ₀) are the center of the image plane. α and β are the LAO/RAO and CAU/CRA rotational angles of the C-arm.

When X-ray passes through 3-D objects, the energy attenuates gradually in its propagating path according to the absorption ratio of the object. X-ray attenuation has been proven to come in the form of a negative exponent, which is decided by the attenuation coefficient and voxel size. As the angiographic image represents the energy integration of X-ray in its propagation path, the gray value of each pixel should be computed by the weighted integration of its absorptions. Thus, we have: $\begin{matrix} I_{i} = I_{0} e^{- \sum_{m = 1}^{G} κ_{m} λ_{m} d_{m}} \\ λ_{m} = [\frac{v_{m}}{1000} + 1] λ_{water} \\ κ_{i} = \frac{s_{m}}{s_{i}} \end{matrix}$ (16) where I ₀ is the initial energy of X-ray. λ _m and d _m are the attenuation coefficient and size of the mth voxel of the vascular structure, respectively. κ _i is the weight, G is the number of voxels in the propagation path. v _m is the attenuation coefficient of the mth voxel of the object that is determined by the attenuation coefficient of water λ _water. s _i is the area of the i th pixel. s _m is the overlapping area of the projection of m th voxel and the pixel. The calculation model of angiograms is shown in Fig. 4.

4 Experiment results

The proposed method is developed by the python programming language, and the experiments are performed on a relative low cost PC with CPU i7-4770, 16G computer memory and GeForce GTX 760 graphics card.

4.1 Evaluation of 3-D vessel

We simulate 15 vascular structures on ellipsoid, sphere, and heart surfaces. Each surface contains five different topological structures. The radius of the root segment is set to 0.384. Let M - EV _m (m = 1, …, 5), M - SV _m (m = 1, …, 5) and M - HV _m (m = 1, …, 5) represent the five different vasculatures on the sphere, ellipsoid, and heart surfaces, respectively. The sizes of the ellipsoid and sphere surfaces are set to [9, 4, 4] and [5, 5, 5], respectively. The heart data are obtained from the Biomed Town Heart Model Dataset 1 . The key parameters of the vasculatures are presented in Table 1.

Figure 5 shows the 15 vascular structures as presented in Table 1. Figure 5(A1) to 5(A5) represent the vessels on the ellipsoid surface, 5(B1) to 5(B5) represent the vessels on the sphere surface, and 5(C1) to 5(C5) represent the vessels on the heart surface. All the generated vascular structures regularly bifurcate from the root to the terminal branches, and they strictly adhere to the predefined surface. When the bifurcation exponent, asymmetric ratio, and length–diameter ratio vary, the shape of the vascular structure varies accordingly. Even on the same surface, the structures of the vessels differ to a large extent.

The calculation times of the vessel centerline generation on the three surfaces are 4 . 4 ± 0 . 0s, 6 . 1 ± 0 . 0s, and 35 . 0 ± 0 . 0s, respectively. The time efficiency of centerline generation is mainly depended on the shape and the size of the predefined surface.

4.2 Evaluation of clinical data

To simulate the X-ray angiograms of the vascular structure, the imaging parameters of the rotational angiographic system are defined according to the Siemens AXIOM-Artis imaging device. The detailed imaging parameters of the projection processes are provided in Table 2.

The algorithm is first validated on the coronary artery data sets obtained from the publicly available data repository released by the MICCAI Grand Challenge 2008 2 . The size of the CT data is 512 × 512 × 272 while the resolution is [0.36, 0.36, 0.40]. The data set includes the centerline and diameter of the coronary artery in the CT image. This information can thus be utilized as the golden standard for the evaluation process.

Figure 6 shows an example of the generated angiogram at the imaging angle of (α = 0°β = 0°). Figure 6(A) shows the 3-D rendering surface of the extracted vasculature, and Fig. 6(B) shows the simulated result of the whole anatomical structure of the CT data set after the coronary arteries are enhanced. Figure 6(B1) to 6(B4) demonstrate the magnified views of the marked region in Fig. 6(B). A realistic angiogram is clearly generated from the simulated 3-D vasculature. In the angiogram, the effects of overlapping and foreshortening can be conveniently visualized.

Figure 7 shows the simulated angiograms in different view angles. The first and third columns are the direct projected vasculatures while the second and fourth columns are the projected images of the whole anatomical structure. The first to the second rows show the simulated images at different imaging angles. The overlapping effects of the perspective projection are realistically simulated in the projection images. When integrated with the whole anatomical structure, practical angiograms can be simulated with high fidelity. From the projection images at different imaging angles, the perspective projection effects are vividly simulated. The average calculation time of X-ray angiogram simulation in different viewing angles is 55 . 0 ± 0 . 1s.

Figure 8 provides the accuracy evaluation of the proposed angiogram simulation method. The quantitative assessment of the angiogram simulation is realized by calculating the error of the gray histogram between the clinical angiogram and the simulation angiogram. The gray histogram of the clinical angiograms is obtained by computing the average histogram of sixty angiograms. Angles in fig. 8 are CAU 10° CAU 20° CAU 30° CRA 10° CRA 20° CRA 30° respectively. The gray histogram of the imaging angles is achieved by calculating the average histogram of six angiograms in the viewing angle which varies from LAO 15° to LAO 90° with interval of 15°. And the gray histogram of the remaining three angles is acquired by computing the average histogram of thirteen angiograms in the viewing angle from RAO 15° to RAO 195° with interval of 15°. From the error statistics in the range of the vascular gray distribution, the histograms of the clinical angiograms and the simulated angiograms present high similarity. Hence, the proposed method is very robust for the angiogram simulation.

4.3 Evaluation of the simulated vasculature

Figure 9 demonstrates the simulated angiograms through the projection of the 3-D vasculature at different imaging angles. The first column shows the volume rendering of the simulated vasculature in the 3-D space, and the second column shows the projected vasculature that overlaps with a predefined mask image. The third column provides the magnified regions of interest that correspond to the location of the image in the second column. The first, second, and third rows correspond to the structural models M-EV2, M-SV5, and M-HV1, respectively. The 3-D vascular models effectively preserve the shape of the attaching surfaces. Meanwhile, the projection effects of the vasculature can be effectively visualized in the simulated angiograms. In the magnified views of the regions of interest, the overlapping effects of the bifurcations and branches are realistically simulated.

Figures 10 to 12 demonstrate the simulated angiograms in the different imaging angles of the models M-EV2, M-SV5, and M-HV1. By projecting the 3-D model in the image plane and integrating it with the mask image, realistic angiograms in different imaging angles are simulated. The coronary artery is presented at a high contrast because the vasculatures are designed to have higher X-ray absorption ability than the surrounding anatomical structures. As branches in terminal endings have small radii, which make their X-ray attenuation ability lower than that of root branches, their appearance in the angiograms is of lower contrast compared with the appearance of the other branches. In addition, the overlapping or foreshortening effects of the vasculature in the angiograms are realistically simulated. The attaching model of the surface can be effectively visualized by comparing the vasculatures in different imaging angles. As the topological relationship of the vasculature can be exactly controlled, the projection effects of the angiogram can be vividly simulated, as shown in the images with little overlapping or shortening.

5 Conclusion and discussion

This study proposes a novel geometrical force constraint method for vascular structure modeling and angiogram simulation. Spacing filling force, gravitational force, and topological preserving force are constructed to form the topological constraint of the vasculature. Meanwhile, surface covering force and surface adhesion force are constructed to constrain the generated vasculature on any surface. According to the combination effects of the topological and surface adhering forces, a realistic vasculature can be effectively simulated on any surface. To construct a vasculature following the commonly accepted energy dissipation principle, Murray’s law is introduced in the construction of the diameter relationship between the parent and the two daughter branches at all bifurcations of the generated vascular structure. An angiographic image of the generated 3-D vascular structures is simulated according to the perspective projection principle of X-rays. The projection effect of the vasculature is simulated according to the X-ray attenuation principle, in which X-rays pass through an object with predefined densities. To obtain a realistic simulation effect for the angiogram, the simulated projection image of the vascular structure is fused with a predefined angiographic mask image. As the 3-D vasculature and imaging angle can be exactly controlled, the angiogram of the vasculature at any imaging angle can be vividly simulated.

3-D vessels on three different geometrical surfaces, namely, ellipsoid, sphere, and heart surfaces, are generated. Five different structures of vessels with different symmetries and length–diameter ratios are generated on each surface. Based on the simulated 3-D vasculature, angiograms at different imaging angles of the X-ray angiographic device can be effectively simulated. The proposed angiographic simulation method is evaluated on the CT image obtained from the MICCAI Grand Challenge 2008. As the coronary structure and the volume data of the chest are both known in the CT image, realistic angiograms are effectively simulated at any imaging angle. In addition, 3-D structure and projection angles are exactly known for the proposed method. The method thus provides a useful way to evaluate the performance of angiographic image analysis methods, such as vascular segmentation, 3-D reconstruction of coronary arteries, optimal view angle determination, and hemodynamics analysis.

Competing interests

The authors declare that they have no competing interests.

Footnotes

Acknowledgments

This research was supported by the National Basic Research Program of China (2013CB328806), the Key Projects in the National Science & Technology Pillar Program (2013BAI01B01), the NationalHi-Tech Research and Development Program (2015AA043203), the National Science Foundation Program of China (81430039, 61501030), the China Postdoctoral Science Foundation (2014M560050).

References

Lloyd-Jones

, Adams

R.J.

, Brown

T.M.

, Carnethon

, Dai

, De Simone

, et al., Heart disease and stroke statistics— update A report from the American Heart Association, Circulation 121 (2010), e46–e215.

Mathers

C.D.

, Fat

D.M.

and Boerma

, The Global Burden of Disease: 2004 update: World Health Organization, 2008.

Finegold

J.A.

, Asaria

and Francis

D.P.

, Mortality from ischaemic heart disease by country, region and age: Statistics from World Health Organisation and United Nations, International Journal of Cardiology 168 (2013), 934–945.

Cademartiri

, La Grutta

, Palumbo

, Malagutti

, Pugliese

, Meijboom

W.B.

, et al., Non-invasive visualization of coronary atherosclerosis: State-of-art, Journal of Cardiovascular Medicine 8 (2007), 129–137.

Tipper

, Detection and evaluation of intracranial aneurysms with 16-row multislice CT angiography, Clinical Radiology 60 (2005), 565–572.

Hyde

D.E.

, Fox

A.J.

, Gulka

, Kalapos

, Lee

D.H.

, Pelz

D.M.

, et al., Internal carotid artery stenosis measurement comparison of 3D computed rotational angiography and conventional digital subtraction angiography, Stroke 35 (2004), 2776–2781.

Shechter

, Devernay

, Coste-Manière

, Quyyumi

and McVeigh

E.R.

, Three-dimensional motion tracking of coronary arteries in biplane cineangiograms, Medical Imaging, IEEE Transactions on 22 (2003), 493–503.

Chen

, Hoffmann

and Carroll

, Computer assisted coronary intervention: 3D reconstruction and determination of optimal views, in Computers in Cardiology (1996), 117–120.

Chen

S.J.

and Carroll

J.D.

, 3-D reconstruction of coronary arterial tree to optimize angiographic visualization, Medical Imaging, IEEE Transactions on 19 (2000), 318–336.

10.

Waechter

, Bredno

, Hermans

, Weese

, Barratt

D.C.

and Hawkes

D.J.

, Model-based blood flow quantification from rotational angiography, Medical Image Analysis 12 (2008), 586–602.

11.

Schaap

, Metz

C.T.

, van Walsum

, van der Giessen

A.G.

, Weustink

A.C.

, Mollet

N.R.

, et al., Standardized evaluation methodology and reference database for evaluating coronary artery centerline extraction algorithms, Medical Image Analysis 13 (2009), 701–714.

12.

Kirbas

and Quek

, A review of vessel extraction techniques and algorithms, ACM Computing Surveys (CSUR) 36 (2004), 81–121.

13.

Duncan

J.S.

and Ayache

, Medical image analysis: Progress over two decades and the challenges ahead, Pattern Analysis and Machine Intelligence, IEEE Transactions on 22 (2000), 85–106.

14.

Frangi

A.F.

, Niessen

W.J.

, Vincken

K.L.

and Viergever

M.A.

, Multiscale vessel enhancement filtering, in Medical Image Computing and Computer-Assisted Interventation–MICCAI’98, ed: Springer, 1998, pp. 130–137.

15.

Aylward

S.R.

and Bullitt

, Initialization, noise, singularities and scale in height ridge traversal for tubular object centerline extraction, Medical Imaging, IEEE Transactions on 21 (2002), 61–75.

16.

Liao

, Luc

, Sun

and Kirchberg

, 3-D reconstruction of the coronary artery tree from multiple views of a rotational X-ray angiography, The International Journal of Cardiovascular imaging 26 (2010), 733–749.

17.

Yang

, Wang

, Liu

, Tang

and Chen

, Novel approach for 3-D reconstruction of coronary arteries from two uncalibrated angiographic images, Image Processing, IEEE Transactions on 18 (2009), 1563–1572.

18.

Yin

, Adel

and Bourennane

, Retinal vessel segmentation using a probabilistic tracking method, Pattern Recognition 45 (2012), 1235–1244.

19.

Bouix

, Siddiqi

and Tannenbaum

, Flux driven automatic centerline extraction, Medical Image Analysis 9 (2005), 209–221.

20.

Greenspan

, Laifenfeld

, Einav

and Barnea

, Evaluation of center-line extraction algorithms in quantitative coronary angiography, Medical Imaging, IEEE Transactions on 20 (2001), 928–941.

21.

Schreiner

, Computer generation of complex arterial tree models, Journal of Biomedical Engineering 15 (1993), 148–150.

22.

Karch

, Neumann

and Schreiner

, A three-dimensional model for arterial tree representation, generated by constrained constructive optimization, Computers in Biology and Medicine 29 (1999), 19–38.

23.

Schreiner

, Karch

, Neumann

, Szawlowski

and Roedler

, Optimized arterial trees supplying hollow organs, Medical Engineering & Physics 28 (2006), 416–429.

24.

Schwen

L.O.

and Preusser

, Analysis and algorithmic generation of hepatic vascular systems, International Journal of Hepatology 2012 (2012), Article ID 357687, 1–17.

25.

Yang

and Wang

, Design of vascular networks: A mathematical model approach, International Journal For Numerical Methods in Biomedical Engineering 29 (2013), 515–529.

26.

Mittal

, Zhou

, Ung

, Linares

, Molloi

and Kassab

, A computer reconstruction of the entire coronary arterial tree based on detailed morphometric data, Annals of Biomedical Engineering 33 (2005), 1015–1026.

27.

Jurczuk

, Kretowski

and Bezy-Wendling

, Vascular system modeling in parallel environment-distributed and shared memory approaches, Information Technology in Biomedicine, IEEE Transactions on 15 (2011), 668–672.

28.

Waters

S.L.

, Alastruey

, Beard

D.A.

, Bovendeerd

P.H.

, Davies

P.F.

, Jayaraman

, et al., Theoretical models for coronary vascular biomechanics: Progress & challenges, Progress in Biophysics and Molecular Biology 104 (2011), 49–76.

29.

Galigekere

R.R.

, Wiesent

and Holdsworth

D.W.

, Cone-beam reprojection using projection-matrices, Medical Imaging, IEEE Transactions on 22 (2003), 1202–1214.

30.

Riddell

and Trousset

, Rectification for cone-beam projection and backprojection, Medical Imaging, IEEE Transactions on 25 (2006), 950–962.

31.

Murray

C.D.

, The physiological principle of minimum work: I. The vascular system and the cost of blood volume, Proceedings of The National Academy of Sciences of The United States of America 12 (1926), 207.

32.

Uylings

, Optimization of diameters and bifurcation angles in lung and vascular tree structures, Bulletin of Mathematical Biology 39 (1977), 509–520.

33.

Suwa

and Takahashi

, Morphological and Morphometrical Analysis of Circulation in Hypertension and Ischemic Kidney, Urban & Schwarzenberg München–Wien, 1971.

34.

Powell

M.J.

, UOBYQA: unconstrained optimization by quadratic approximation, Mathematical Programming 92 (2002), 555–582.

35.

Staal

, Abramoff

M.D.

, Niemeijer

, Viergever

M.A.

and van Ginneken

, Ridge-based vessel segmentation in color images of the retina, Medical Imaging, IEEE Transactions on 23 (2004), 501–509.