**4. Interpolation of Shape Control Points**

In this section, we introduce the shape interpolation and restoration for real-time usage of the 3D+t coronary artery model. According to Equation 3, the shape of the coronary artery relies on the locations of the control points. Therefore, as we derive the intermediate positions of control points among cardiac phases, the corresponding shape is restored. To interpolate the positions of control points, we consider a set of control point at the *k* th phase as a vector *P<sup>k</sup>* , such that *P<sup>k</sup>* = {*p*0, *p*1, . . . , *pn*−2, *pn*−1}, where *<sup>p</sup><sup>i</sup>* <sup>∈</sup> <sup>R</sup><sup>3</sup> and *<sup>n</sup>* is the number of cage control points. From the registration results, we interpolate the given sets of control points using periodic cubic spline interpolation [43], due to the (cyclic) nature of the heart's motion. The number of knots is the same as the number of reconstructions from 4D CTA.

Let a phase-varying vector *S*(*t*) = {*s*0(*t*), . . . ,*sn*−1(*t*)} be the set of interpolated control points, where *si*(*t*) is the *i*th control spline for the cardiac phase *t*. The spline vector *S*(*t*) has *C* 2 continuity with respect to the phase *t*. The spline function *S*(*t*) maps phase *t* to the set of cage control points, such that *<sup>S</sup>* : <sup>R</sup> <sup>→</sup> <sup>R</sup>3×*<sup>n</sup>* . Then, the vertices of shape are restored using the following equation:

$$v(t) = F(v; S(t)) = \sum\_{i \in I\_V} \varphi\_i(v) s\_i(t). \tag{26}$$

#### **5. Evaluations and Results**

We evaluated the proposed method both qualitatively and quantitatively on data from eight patients. The proposed method was tested on an Intel (R) Xeon (R) W-2133 workstation with CPU@3.60 GHz and 32 GB ram. We partially multi-threaded the computation of cost function measurements using OpenMp [44] and Thread building block [45] during the optimization process. The proposed method and comparison target methods were written in C++.

## *5.1. Quantitative Evaluations*

In the quantitative evaluation of non-rigid registration, we used metrics considering: (1) the closest point-mesh Euclidean distance (ED) from the target model to the matching result and (2) the dice coefficient (DC) obtained from the mesh boolean operation. As we set the number of iterations to 300/maxDepth for each depth, the total number of iterations for different max depths was set to be the same.

#### 5.1.1. Trade-off between Deformation Depth and Computation Time

First, we observed the trade-off between the degrees of freedom of deformation and computation time. As shown in Table 1, we compared the different depths of deformation incrementally, from 1 to 5. As shown in Figure 6, the ED and DC worsened, as the phases were far from the template phase. As the shape of the blood vessel was a thin tube shape, the ED and DC values noticeably deteriorated with slight movement. As the deformation depth increased, the ED values gradually decreased and the DC values increased more prominently; both metrics flattened for the other cardiac

phases. The metrics converged after deformation depth 4. The comparison results for the other patients are given in Appendix A.


**Table 1.** Trade-off between computation time and accuracy.

**Figure 6.** Effect of cage deformation depth for patient 1 in different cardiac phases: (**a**) dice coefficients; and (**b**) average distance from target to deformed model.

#### 5.1.2. Comparison with Other Methods

In the second experiment, the proposed registration method's performance was compared with that of other non-rigid matching algorithms. As the comparison target of non-rigid registration, we selected the variations of GMM methods, which are combinations of a deformation model and a cost function. The deformation models were thin-plate spline (TPS) and generalized radial basis function (GRBF), while the cost functions were kernel correlation (KC) or *L*<sup>2</sup> distance. The comparison targets used L-BFGS-B as an optimization method.

For comparison, each deformation model had the same number of deformation control points. Considering the convergence of accuracy from the previous analysis, we set the number of grid and control points as 16 × 16 × 16 and 4913, respectively. As shown in Figure 7, the proposed method had higher DCs than the comparison targets at the interval [15%, 35%] of cardiac phases, where the interval had a DC value of less than 0.2 before registration. Although the ED metrics showed a similar trend, compared with other methods, a significant improvement in DC was observed. The comparison results for the other patients are given in Appendix A.

**Figure 7.** Comparison with other algorithms for patient 1 at the different cardiac phases: (**a**) dice coefficients; and (**b**) average distance from target to deformed model.

#### 5.1.3. Interpolation Accuracy

Furthermore, we evaluated the accuracy of the interpolated 3D+t coronary artery models by comparing them with the segmented models. The proposed method created a smooth and non-degenerate 3D model by interpolating the cage control points over sampled cardiac phases, as shown in Figure 8. Our data sets were evenly reconstructed from 4D CT within the R-R peak with 5% sampling interval. Thus, we had 20 keyframes (*V*0%, *V*5%, . . . , *V*95%). To evaluate the effect of sampling the cardiac phases, we chose phase sets from the given 20 keyframes as follows: (1) Odd 10: [5, 15, 25, . . . , 85, 95]; (2) Even 10: [0, 10, 20, . . . , 80, 90]; (3) Odd 5: [5, 25, 45, 65, 85]; and (4) Even 5: [0, 20, 40, 60, 80]. Figure 9 shows the differences among the phase selections. Although the sampled phase became sparse, the DC of the interpolated result showed that the results had a lower bound. The ED still showed a flattened value, when compared to that before registration. The comparison results for the other patients are given in Appendix A.

**Figure 8.** Comparison of the registered model and interpolated model for patient 5: (**a**) Template (green) and 40% coronary artery (red); (**b**) registered model (blue); (**c**) interpolated model (yellow); and (**d**) comparison of registered model and interpolated model.

**Figure 9.** Comparison of interpolation sampling for patient 1 in the different cardiac phases: (**a**) dice coefficients; and (**b**) average distance from target to deformed model.

#### *5.2. Qualitative Evaluations*

In the qualitative evaluation, the role of the visualization effect of hyper-elastic regularization, geometrical comparison with other algorithms, comparison of the matching result and interpolation result model, and the limitations of the interpolation model were investigated.

#### 5.2.1. The Effect of Hyper-Elastic Regularization and Hierarchical Deformation

High-order deformation models often converge to a local minimum, which may look visually implausible. Figure 10a,b show examples of shape shrinkage when the target model contains loss of shape. The hyper-elastic regularization constraints lead to shape preservation, thus providing a plausible result, as shown in Figure 10c,d.

When compared with the other algorithms, Figure 11 shows an example where shape registration is defective at the excessively deformed and twisted parts. As the registration process converged to a local minimum, the deformation model represents the further details of local deformation. On the other hand, hierarchical cage deformation gradually acquired an optimal solution, passing from coarse to dense resolution, to avoid local minima, as shown in Figure 12. In this process, the low degree of freedom deformation serves as the initial value for the deformation in the next step. Therefore, local minima can be avoided more efficiently.

**Figure 10.** The effect of modified hyper-elastic regularization. We aligned the source model to the target model (red), which contains a loss of branch. The figures show effects: (**a**,**b**) without regularization (yellow) and (**c**,**d**) with regularization (blue).

**Figure 11.** Qualitative comparison of GMM and the proposed method: (**a**) Initial template model (green) and target coronary artery (red); (**b**) result of the proposed method (blue); and (**c**) result of Gaussian mixture modeling with TPS+*L*<sup>2</sup> (yellow).

**Figure 12.** Qualitative comparison of the proposed method while changing the deformation resolution: (**a**) Initial template model (blue) and target model (red); (**b**–**f**) the results of registration (blue) at different cage resolutions from [1, 1, 1] to [5, 5, 5], respectively.

## 5.2.2. The Representation Power of Interpolated Model

The proposed shape interpolation method may have limited representation ability for intermediate shapes. This limited shape representation is due to the recurring shape of the adjacent phases. We observed that the shape interpolation method restrictively delineates the intermediate shape. The shapes of the neighboring phases to the target phase resemble each other, but the target shape and the neighbor shapes are considerably different, as shown in Figure 13e.

**Figure 13.** Comparison of the registered model and interpolated model for patient 8: (a) Template (green) and 95% coronary artery (red edges); (**b**) registered model (blue); (**c**) interpolated model (yellow); (**d**) comparison of registration interpolation; and (**e**) comparison of neighboring coronary artery models (90%, dark blue; 0%, light green).

#### **6. Discussion and Conclusions**

In this paper, we proposed a method for generating a 3D + time vessel model from 4D CT images that can be used in real-time. Our purpose was to create a 4D vascular model without mesh degeneration, interpolate the model at high speed, and express a more precise shape.

To create a 4D vascular model, we matched the diastolic coronary artery model with the coronary artery model in other phases through hierarchical cage deformation. During the registration process, hyper-elastic regularization was used as a shape preservation constraint. The shape control points obtained as a result of registration were interpolated into a cyclic cubic spline to create a 3D+t model. The shape change depends only on the control points of the cage. The rapid deformation application and the preservation features of the local information are beneficial in the shape registration process.

To evaluate the precision of the proposed method, quantitative and qualitative evaluation was performed on 160 CTA volumes acquired from eight patients. In the quantitative evaluation, we assessed:


In the step of measuring the shape matching precision according to the hierarchical deformation, we observed that the matching precision converged in the fourth step of the hierarchical deformation, where the calculation time was 23.53 s on average. In the fifth deformation depth of hierarchical transformation, the matching accuracy slightly increased, but the required time increased by 28.70% (to 33.00 s). In the hierarchical transformation of the cage creation stage, the control points constituting the cage increased exponentially, as a regular grid was used. We obtained the mean distance with precision of a 0.543 mm and standard deviation 0.222 in step 4 of the hierarchical transformation, where the Dice coefficient obtained an average of 0.754 and a standard deviation of 0.064.

Compared with other algorithms, the GMM method requires the creation of a mixture model for each point and, so, even with the same degree of freedom of transformation, the calculation time was as high as 40 s for the GRBF model and 33 s for the TPS model. In addition, in the average distance index, TPS\_L2 was 0.530 mm, which had an error lower than that (0.543) of the proposed method; however, when comparing the Dice coefficient, TPS\_L2 was observed to be 0.690, 0.045 points lower than the index of the proposed method (0.735).

If the indicators of DC and AD values conflict with each other, it is necessary to determine which indicator is better to express the accuracy of the matching result; for example, (1) high DC value and high AD value or (2) low DC Value and low AD value. At this time, we decided case (1) was a better indicator. This was because, in the vascular model between 15% and 35%, which showed a great difference with the 75% phase, a difference in AD values between the algorithms was not significantly observed, but the difference in DC values was noticeable. This was not only consistently observed in the quantitative evaluation of patient data, but also explains the differences arising from inappropriate deformations occurring in excessively twisted blood vessels during the qualitative evaluation.

The precision of the shape interpolation model was measured with different phase-sampling sets, where the accuracy was worsened with a larger phase-sampling interval. However, this limitation may be resolved by increasing the temporal resolution and by co-operating with the other real-time imaging systems, such as X-ray angiography or 4D US.

A qualitative evaluation was performed to support the quantitative evaluation mentioned above, as well as to observe the effect of the proposed model on the visualization stage. In the qualitative evaluation, (1) the effect of the modified hyper-elastic regularization and hierarchical transformation, and (2) the limitations in shape interpolation were observed.

When we observed the effect of hyper-elastic regularization, the shape was transformed to be visually plausible when hyper-elastic regularization was used. This phenomenon was particularly well seen in the coronary artery model with loss of shape. This is a characteristic obtained by minimizing the excessive deformation by robustly obtaining face points and volume points against the rapid degeneration of the cage mesh due to incorrect correspondence pairs. Compared with other non-rigid algorithms, the proposed method was able to cope with the local minima that occur during the optimization process by hierarchically performing the transformation effectively. In particular, it was shown that, in the vascular model with an excessive twist, optimal deformation was gradually obtained from low-dimensional deformation.

In summary, the proposed 3D+t vascular modeling method utilized hierarchical deformation for robust shape registration, while interpolation of the registered vascular structure enabled the restoration of small and complex geometries due to the cardiac cycle.

The electric potential of the myocardium generates the ECG signal, and the electric potential of the myocardium is related to movement and contraction of the heart. As an observation tool of the myocardium, the proposed method can provide an alternative to real-time imaging by using the ECG signal and 4D CT. In intraoperative situations, invasive coronary angiography is a method for monitoring the movement of the coronary artery, which provides limited deformation information (up to the 2D plane). With the proposed method, we expect that 3D contraction and strain of the myocardium, according to the ECG signal, can be observed, as the coronary arteries are attached to the epicardium. The ability of motion monitoring is directly related to the evaluation of the physiological function of the myocardium.

In addition, the modified hyper-elastic regularization prevented implausible deformation and mesh degeneration, which must be avoided in the analysis of 4D blood flow. We demonstrated the accuracy of the proposed method by presenting qualitative and quantitative evaluations using data from eight patients.

Therefore, we expect that the proposed 3D+t vascular model can be utilized in real-time applications such as (1) pre-operative blood flow analysis, which requires rapid shape creation without mesh degeneration; and (2) 2D invasive coronary angiography-3D shape registration during the intervention, where it can be used as an image guidance tool that provides real-time shape information according to the ECG signal during percutaneous coronary intervention.

As a limitation of this study, since the deformed model of the proposed model is limited to a 3D mesh, blood vessel segmentation is required in phases other than the template model. In a future study, we will conduct a study on a volume–template mesh model matching method which can be applied to volumetric data, including shape loss, in order to eliminate unnecessary repetitive processes.

**Author Contributions:** Conceptualization, S.Y., C.Y., E.J.C. and D.L.; Data curation, S.Y., C.Y. and E.J.C.; Formal analysis, S.Y. and D.L.; Funding acquisition, D.L.; Investigation, S.Y.; Methodology, S.Y. and D.L.; Project administration, D.L.; Resources, S.Y., C.Y., E.J.C. and D.L.; Software, S.Y.; Supervision, D.L.; Validation, S.Y.; Visualization, S.Y.; Writing—original draft, S.Y.; Writing—review & editing, S.Y. and D.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by KIST intramural grants (2E30260) and the Ministry of Trade Industry & Energy (MOTIE, Korea), Ministry of Science & ICT (MSIT, Korea), and Ministry of Health & Welfare (MOHW, Korea) under Technology Development Program for AI-Bio-Robot-Medicine Convergence (20001655).

**Conflicts of Interest:** The authors declare no conflict of interest.
