**1. Introduction**

Cardiovascular disease (CVD) is one of the primary causes of death worldwide, with 22.2 million deaths expected by 2030. According to National Health and Nutrition Examination Survey (NHANES) data from 2013 to 2016, the prevalence of CVD was 48.0% in adults over the age of 20. The prevalence of CVD has a positive correlation with an increase with age [1]. The resulting social cost is estimated to have been 351.3 billion dollars in the U.S. alone, from 2014 to 2015. In particular, cardiovascular disease accounted for 14% of total medical spending in U.S., the highest rate among other major diagnostic groups—even higher than cancer. In the global population, the burden of expenditure is even more serious [1].

Providing sufficient information through image analysis acquired in the pre-operative diagnosis stage eliminates unnecessary examination and helps in developing patient-specific treatment plans.

As the heart is a continuously beating organ and there may be unexpected movements in patients (e.g., arrhythmia), if information on this movement can be obtained in advance, coronary artery and heart procedures may be more efficient.

As shown in Figure 1, the ECG signal is a change in potential that is correlated with the movement of the heart muscle. Motion of the heart produces changes in volume and pressure in the cardiac chambers; therefore, ECG provides important information about the movement of the heart. When CT is reconstructed by performing retrospective ECG synchronization, the movements of the heart and coronary arteries (according to the cardiac cycle) can be obtained geometrically, and movement information (e.g., coronary artery distortion), as well as characteristics of the coronary stenosis, can be obtained.

**Figure 1.** The ECG signal and volume of the ventricle during the different phases of a cardiac cycle.

Patient-specific 4D heart shape information facilitates the following applications: accurate coronary artery structure acquisition [2,3], analysis of 4D blood flow and stenosis [4–6], removal of motion artifacts in the vascular region [7–9], surgical simulation for each patient [10], postoperative evaluation and analysis [5,6], atrial motion analysis [11,12], vascular motion analysis [13,14], and real-time deformation prediction during surgery combined with 2D images [8].

In particular, cardiac CT angiography (CTA) has an isotropic spatial resolution of less than 0.5 mm and, so, can be used to observe the movement of the coronary artery and trabeculae of the ventricle. In addition to grading the degree of calcification of the coronary artery and the total amount of plaque from the CT image, it is also possible to measure the torsion of the coronary artery at a high resolution [15]. This high-resolution spatial information can help the operator perform a procedure appropriate for each patient before and after surgery. However, despite the high spatial resolution of the CTA, its temporal resolution is 50–200 ms, which is lower than that of 4D echocardiography (30–100 Hz) and cardiac MRI (30–50 ms). Therefore, proper shape interpolation for restoring high time resolution information from CTA imagery is essential for co-operation with the other applications while preserving high-precision anatomical details.

One of the essential requirements of cardiac modeling in many of these applications is that the topology of the mesh model constituting the cardiac model must be consistently preserved. In particular, in the case of a 4D heart model, the mesh should not cause new problems (e.g., self-intersection or mesh degeneration), even if the positions of the vertices constituting the shape change over time. To address these requirements, many researchers have employed the template-based registration scheme.

The registration process matches a template model to a target model through geometric shape deformation. The most popular registration algorithm is the iterative closest point (ICP) method, which consists of finding the correspondence between two models and finding the optimal transformation [16]. However, the ICP method is sensitive to the initial position and noise, while the shape registration is limited up to rigid transformation. Non-rigid registration deals with the deformation of the shape in addition to rigid transformations. Non-rigid registration is more challenging, however, as non-rigid transformations not only require more correspondences to be defined, but the solution space is much more extensive [17]. Research into the registration of 3D non-rigid shapes has been actively conducted using the methods of delineating shape deformation and shape correspondence.

To describe the deformation of an object, with respect to its dynamics and material properties, many researchers have assumed shape deformation to be a physical model, such as a linear elastic model [18], non-linear elastic model [19,20], viscous fluid [21], or diffusion model [22,23]. In particular, the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework provides robust deformation as a massive flow consisting of diffeomorphisms [24]. However, the physical models are computationally expensive and sensitive to mechanical properties. On the other hand, the statistical shape deformation model (SSM) uses a low-dimensional statistical model, in which shape deformation is inferred from a training data set [25,26]. Although SSM reduces the computational cost, the shape of variability is limited by the training data. Therefore, the deformation is hardly representative of the inter-variability of patients, such as the topological discontinuity of coronary arteries. Nora et al. described the motion modeling problem using the coronary arteries attached to the SSM of muscles [27]. Due to the representation of the coronary artery, the deformation poorly described the lumen diameter. Instead of modeling a priori physical and statistical information, there have been attempts to estimate the shape transformations with landmarks and coherent motions. Radial basis function methods express shape deformations as weighted sums of distance function for control point changes [28]. In particular, the thin-plate spline method minimizes the bending energy, which has a closed-form solution [29]. The most popular form of deformation is known as B-spline free-form deformation (FFD) [26,30,31]. Rueckert et al. [32] have proved the conditions for FFD to have diffeomorphic deformation. This method has disadvantages, however: the numerical cost increases with the number of control points and the degree of freedom of shape deformation is fixed. Therefore, the deformation has limited representation capacity.

In addition to deformation modeling, establishing correspondences between shapes is a critical problem in registration. The one-to-one correspondence of the ICP method is sensitive to the initial position and shape loss. To determine many-to-many correspondence points, Chui et al. [33] used the fuzzy correspondence between two shapes. The problem of selecting a robust point matching the correspondence has been interpreted as a combination of a Gaussian mixture model (GMM) and Expectation Maximization [34]. In the GMM model, one point is the centroid of the Gaussian distribution for the points constituting the shape, while the other point is regarded as the data to generate [35]. The variations of GMM have different deformation models, according to the obtained transformation parameters and the regularization term. For example, regularizing the second derivative of the transformation leads to a thin-plate spline transformation, while regularizing according to motion coherence theory leads to a coherent point drift transformation [36,37]. The variations of GMM have been generalized using the generalized Gaussian radial basis function [38]. To represent the local spatial representation, an *L*2*E* estimator has been proposed, which creates a robust sparse–dense correspondence [39,40].

However, the mixture model requires a high computational cost, as it generates a Gaussian distribution for each of the points. Furthermore, each Gaussian distribution shares its standard deviation among the points. The mixture model is thus sensitive to noise and shape loss. Especially for coronary arteries, narrow and tangled structures are very challenging to model in 3D+t. This is because the loss of blood vessel morphology can be observed in different heartbeats of the same patient, through motion artifacts and geometrical deformations.

In this paper, we propose a 3D+t coronary artery model that can be inferred in real-time, according to ECG signals. The overall structure of the proposed method is shown in Figure 2. The proposed patient-specific 3D+t coronary artery motion model is divided into two processing blocks, according to the timing of data processing; (1) a preoperative processing block, and (2) an intraoperative usage block.

**Figure 2.** General framework of the proposed method.

At the preprocessing step, we first perform the segmentation of 4D CTA volumes to generate artery models. After we have multiple coronary artery models, hierarchical cage-based registration is performed to construct a patient-specific 3D+t model with hyperelastic regularization. The proposed hierarchical cage deformation model more robustly/accurately registers coronary artery models in different cardiac phases. When updating the control points of the coronary artery model, we gradually increase the degrees of freedom of the deformation model. A modified hyper-elastic regularization term prevents mesh degeneration problems during the control point optimization step. After the optimal cage control point is obtained—which minimizes the shape dissimilarity of the source shape and the target shape—we interpolate the shape control point to build a continuous 3D+t model. The interpolated shape model provides fast shape-inference for intraoperative usage.

In the intraoperative usage step, the ECG phase can be assessed from the patient's real-time signal, which is then correlated to the geometric deformation of cardiac muscles, as shown in Figure 1. Therefore, the 3D shape of the coronary artery is provided by a continuous 3D+t model, according to change of the ECG signal and time.

The contributions of the proposed method consist of the following:

