Three-Dimensional Atrial Wall Thickness Measurement Algorithm from a Segmented Atrial Wall using a Partial Differential Equation

Despite advancements in high-precision segmentation technology for computed tomographic angiography (CTA)-based cardiac wall segmentation, the accurate detection of the endocardial (Endo) and epicardial (Epi) boundaries remains a prerequisite for automated measurements of the cardiac wall thickness (WT). We proposed a novel algorithm for automated three-dimensional (3D) atrial WT (AWT) measurements, including an automatic Endo-Epi boundary detection. We detected the boundaries that were topologically indistinguishable due to an open geometry at the anatomical boundaries using the combined Convex hull and Poisson solver methods. The Laplace equation for the WT measurement was solved by a partial differential equation combining the two detected boundaries of the myocardial wall. We verified the robustness of our algorithm in mask images of the atrial wall that were separated from the CTA images of 20 patients and a phantom model. The accuracy of the automatically detected Endo-Epi boundaries was acceptable as compared to that manually extracted from the phantom model (Dice coefficient=0.979). The 3D AWTs calculated by the novel automated method from the CTA images obtained from 20 patients with atrial fibrillation had <10% error as compared to the conventional manual AWT measurement method for comparing the regional WT in all locations except the left antral area. The proposed algorithm’s AWT detection time was 27.15±6.99 s per patient, which was 1/40 that of the conventional method. Consequently, our results showed that the proposed automatic 3D AWT measurement algorithm had the potential to significantly improve the efficiency of calculating the AWT while maintaining the existing level of accuracy.


I. INTRODUCTION
Atrial fibrillation (AF) is a common, sustained, and progressive degenerative heart rhythm disorder with a prevalence of 2% [1,2]. Chronic atrial wall stress caused by repeated AF episodes or hemodynamic stress can lead to atrial dilation and heterogeneous changes in the cellular structure [3]. Thus, an altered atrial histology influences the electrophysiological activation patterns during arrhythmias by triggering reentry and vagally mediated AF along boundaries with a rapidly changing myofiber orientation and wall thickness (WT) [4,5]. Although radiofrequency (RF) catheter ablation has been introduced as an effective rhythm control method for AF treatment, there is a risk of procedurerelated complications and a hemopericardium [6]. Despite a successful ablation, the atrial WT (AWT) and atrial pressure increases associated with the progression of AF have been reported to increase the chances of a long-term recurrence [7]. Therefore, the AWT measurements are expected to reduce the risk of complications (such as cardiac perforation or pulmonary vein [PV] stenosis) with an RF energy titration and provide useful information for predicting the likelihood of an AF recurrence [7,8]. The improved resolution offered by a computed tomographic angiography (CTA) scan and the adoption of artificial intelligence (AI) for image segmentation are expected to generate a fully automatic, CTA-based, three-dimensional (3D)-cardiac anatomy reconstruction [9]. Despite advances in the CTA image analysis technology, the automation of the CTA-based AWT measurements that can be applied to the clinical workflow still faces two difficulties: the first is segmenting the atrial wall from the cardiac substructures, while the second is detecting the endocardial (Endo) and epicardial (Epi) boundaries of the atrial wall before the WT calculation [10]. The first issue (atrial wall segmentation) is expected to be partially resolved with the advent of commercial software (Segment CT, Medviso Inc., Lund, Sweden, and Simpleware AS Cardio, Synopsys Inc., CA, USA, etc.) with various AI-based tailored approaches to wall segmentation to identify the anatomical structures and obtain detailed geometric information [11][12][13]. However, the second issue (detecting the Endo and Epi boundaries of the atrial wall) still requires empirical manual work to identify the location information of an open hole, such as the mitral valve [14]. To the best of our knowledge, Wang et al. [8] was the first to present the automatic Endo and Epi surface extraction method with a Convex Hull-based method and reported the feasibility with one ex vivo post-mortem and two in vivo human atrial datasets. However, that approach can lead to an erroneous surface separation if the opening (especially the mitral valve) is concave relative to the atrium. Although the WT measurements have evolved to more robustly account for geometric changes throughout the atrial wall (from naïve methods such as surface-based methods [15,16] to Laplace-based solutions [14]), all of those methods operate under the premise that the Endo and Epi boundaries are already known. The purpose of this study was to propose a procedure to automate the process of detecting the AWT measurements from the segmented atrial walls. We developed an algorithm to automatically detect the Endo and Epi boundaries from the segmented walls based on a partial differential equation (PDE) approach that combined a Convex Hull and Poisson equation and process to automatically calculate the AWTs by applying the Laplacian equation to the detected Endo and Epi boundaries. The developed algorithm was validated using clinical data, and the measurement error and measurement time compared to a manual analysis were also evaluated.

A. ATRIAL WALL SEGMENTATION
A flowchart of the proposed algorithm for fully automated 3D AWT measurements is shown in Figure 1. The input data were defined as a mask image (marked as 0 for the segmented atrial wall and 1 for the remaining regions) on the CTA image to run independently from the segmentation software ( Figure 2a). The format of the mask image was adopted with a common image format (such as a bitmap or jpeg), and the pixel resolution consisted of the same number of pixels as the CTA image. Here, the CTA images were acquired via an electrocardiogram (ECG)-gated acquisition considering the systolic and diastolic phases, and the mask images were manually acquired via AMBER software (Laonmed Inc., Seoul, Korea), which was clinically validated in previous studies [6,7]. Briefly, the Endo of the atrial wall was semi-automatically segmented on the CTA images using an edge detector and then was separated from the other tissue using the multi-Otsu threshold algorithm that is based on the histogram of the Hounsfield unit. Then, the atrial wall mask was extracted from the Endo as an overlapping area through morphology operations. This approach was used to match the segmentation information used in the manual method for the quantitative assessment of the WT rather than the difference due to segmentation. This protocol was designed for application of any segmentation software that generates binary mask images.

B. ENDOCARDIAL AND EPICARDIAL BOUNDARY DETECTION USING POISSON SOLVER
The thickness was defined as the distance connecting the starting and ending points. Since the Endo and Epi surfaces that corresponded to the starting and ending points were not explicitly differentiated in the wall mask images, we developed a Poisson-based algorithm to obtain the two surfaces (Figures 2b-f). In addition, for a high-speed calculation, all procedures were implemented with a graphics processing unit(GPU)-friendly architecture based on computed unified device architecture (CUDA, NVIDIA, Santa Clara, CA, USA).

1) CALCULATING THE BOUNDARY CONTOUR WITH CONVEX-HULL AND POISSON SOLVERS
Determination of the starting and ending points is required to define the boundary conditions of the Laplace equation for the AWT calculation. However, there were two problems that prevented defining the starting and ending points using only the atrial wall mask images. First, since the left atrial (LA) is not geometrically convex due to the PVs and LA appendage (LAA), the convex hull algorithm alone does not solve this problem (Figure 3a). In addition, the LA wall has a structure in which the inner chamber and outer surface are not discriminated topologically (homeomorphic to sphere) due to the open holes of the PVs and mitral valve. Since the wall mesh is a two-dimensional (2D) manifold, user intervention is required to apply conventional geometric algorithms, such as hole-filling algorithms [17]. Therefore, we represented the LA wall as an insulation wall and hypothesized a boundary condition in which the inner temperature no longer changes due to the outer temperature as a condition that can separate the inside from the outside.
In this research, we divided the Endo and Epi surfaces that preserved the shape of the chamber while satisfying the constraint that both surfaces were strictly inside the envelope surface ε in the open polygonal domain of the atrial wall. First, to obtain the envelope surface ε, we obtained the Convex Hull slice-by-slice in each 2D plane of the XYZaxes using the 2D quick hull algorithm [18]. After that, we merged and inverted the intersections of each Convex Hull computed separately along each orthogonal axis to generate an inverse envelope surface εInv ( Figure 2b). Next, we defined the Poisson equation by applying εInv as the Dirichlet boundary condition (ГD) and the atrial wall as the Neumann boundary condition (ГN) (left in Figure 2c). The Poisson equation with mixed boundary conditions was as follows: where 3 () fL  and ∆ denote the Laplacian operator. The Dirichlet boundary condition ГD (the envelope as the merged hull) was set to 10, and the Neumann boundary condition ГN (the atrial wall area) to 0; the calculation was stopped after 100 iterations. The Poisson equation was iteratively solved using the Jacobi method, and the Neumann boundary was multiplied by the binarized atrial wall mask at each iteration step (top right in Figure 2c). Finally, the boundary contour Г at which the inner temperature was no longer eroded was separated with a cutoff = 0.1. Then, binarization was performed to obtain the mask image (i.e., with a value of zero outside the surface and a value of one inside) from the final boundary contour Г (bottom right in Figure 2c).

2) REMOVING OUTSIDE RESIDUAL VOXELS
Despite the erosion induced by the Poisson solver, residual voxels may still remain depending on the outer structure of the LA (Figure 3b). Therefore, two additional procedures were required to remove the outside residual voxels (Figures 2d-e). The first procedure subtracted the wall mask from the binarized final boundary contour mask to physically divide the outer residual voxels from inside the LA chamber ( Figure  2d). This approach involved a simple separation because the open holes that were likely to have connected with the

①). (f) The epicardium and endocardium are detected by identifying adjacent labels with a 3D connectivity filter at each voxel of the wall area. (g) Finally, the wall thickness is measured by solving the Laplace equation as the boundary condition of the separated wall area. (f) The wall thickness map is projected on the endocardium.
adjacent voxels were eroded by the Poisson solver ( Figure  2d). The next procedure removed the smaller fragments other than the one with the largest volume from the voxel fragmented by the first step ( Figure 2e). This step involved a GPU-based connected component labeling algorithm [19] to assign a label to each voxel and then remove the components other than the label of the largest volume to determine the inside of the LA chamber. That procedure could simply remove the residual voxels outside the atrial wall ( Figure 3c).

3) DETERMINING THE ENDOCARDIAL AND EPICARDIAL BOUNDARIES FROM THE ATRIAL WALL VOXELS
The following step was performed to determine the Endo and Epi boundaries at each voxel of the LA wall. For each voxel (the yellow dot in Figure 2f) with a 3D-connected filter (26connectivity), if the label of the adjacent voxel was inside the LA chamber, it was determined to be the Endo (the blue arrow in Figure 2f). When the LA wall label was detected in the adjacent voxel, it was classified as myocardial (the white arrow in Figure 2f), and the remaining labels were classified as Epi (the red arrow in Figure 2f), respectively. If the adjacent voxels contained more than one type, the boundary was determined as a type with a greater number of fetched voxels (Figure 2f).

C. ATRIAL WALL THICKNESS MEASUREMENTS USING LAPLACE EUQATION
The Laplace equation was computed to solve the gradient field used to measure the 3D AWT (Figure 2g). In the previous step, the boundaries detected as Endo, Epi, and myocardial were set to initial conditions of 100, 200, and 300, respectively. The Laplace equation was the second-order PDE for the scalar field  enclosed between two surfaces (Endo surface S and Epi surface S'), as shown in Equation (2): This Laplace equation was iteratively solved with the Jacobi method (3): where ( , , ) i x y z  was the value of the potential at each voxel during the i-th iteration. Convergence was calculated as the sum of the total energies of the scalar field as in Equation (4), and the solver stopped calculating when the threshold value or 400 iterations were reached. The sum of the total energies was calculated using the CUDA atomic function. The gradient of the scalar field was calculated with the following Equation (5): Finally, the integral length of the streamline was calculated using Euler's method (6) from each Endo voxel P to reach the Epi voxel P'.
(0) Euler's method was chosen because it was faster than the Runge-Kutta method and was also expected to have no large error due to the thin myocardial wall thickness; the length of each streamline calculated at all Endo voxels was defined as the myocardial thickness. The PDE stopped calculating when E=10 -5 or there were 400 iterations. Finally, the field line was calculated as the connection between the endocardium and epicardium using Euler's method (∆t=0.001) at all endocardium voxels. Therefore, to satisfy the mutuality between the points P on the endocardium and P' on the epicardium, the two points were connected by a curve.

D. VALIDATION PROCESSES
The study protocol adhered to the Declaration of Helsinki of 2013 and was approved by the Institutional Review Board of Yonsei University Health System (ClinicalTrials.gov; NCT 02138695). All patients provided written informed consent. To verify our algorithms, we performed two validation procedures: (1) Evaluate the accuracy of the Endo and Epi boundary detection methods with the phantom model, and (2) Compare the AWT using the previously proven manual measurement software, AMBER. In addition, the segmentation procedure was performed by two experts. The DSC ranged from zero (no common voxels between the closed and open geometries) to one (perfect match). The MAE and MAPE calculated the error point-by-point for the thickness obtained from the open and closed geometries, respectively, and the formulas were as follows: and are the WT measured using a manual or proposed algorithm, respectively. Second, the accuracy of the AWT was compared to previously proven manual software using the CTA images of 20 AF patients. For a quantitative regional analysis, two AWT maps generated by each software application were projected onto the same LA surface mesh, and the AWP map was divided into eight regions based on the embryological origin and clinical catheter ablation approaches [21,22]. The regional areas were segmented on the mesh surface using the surface segmentation tool in the customized software (CUVIA, Model: SH01, ver. 2.5; Laonmed Inc., Seoul, Korea). A statistical analysis was performed using the R package (3.1.0, R Foundation for Statistical Computing, Boston, MA, USA). The user interactions for previously proven manual measurement software were monitored using the Mousotron (12.1, Blacksun Software, Turnhout, Belgium). To evaluate the computational performance, the hardware system consisted of an Intel Core i9-9900KF CPU, 64GB RAM, and NVIDIA TITAN RTX. The operating system was Windows 10.0 (Microsoft Inc., Redmond, WA, USA).

A. EVALUATION OF THE ENDOCARDIAL AND EPICARDIAL BOUNDARY DETECTION ACCURACY
After a computer-aided design, the CTA spatial resolution of the 3D phantom model (Boston Scientific, MA, USA) generated by 3D printing was 0.246 mm on both the x-and yaxes and 0.5 mm on the z-axis.  Figure 4d).    (Table Ⅲ). Except for the left antral area, the regional MAPEs were <10%. It was a reasonable difference considering the morphology of the most complex LA appendage in the LA geometry. The correlation was 0.918 between the AWT maps generated using the two software methods. The regional AWT map measurements that were acquired by the software methods had a high correlation between R=0.89 to 0.99 (P<0.001, Figure 5).

C. MEASURING THE DICE SIMILARITY OF THE DETECTED SURFACES
We performed a similarity evaluation of the Endo and Epi surfaces reconstructed by two methods. For a fair comparison, we used the same atrial wall mask images for each patient. The reconstruction of the Endo and Epi surfaces was performed by our PDE-based method and the Convex hull-based method (by Wang et al.) for comparison. For Endo and Epi voxels detected by two methods, the similarity was evaluated by DSC, respectively. The similarity confirmed by the overlap of the two boundary voxels was a mean DSC of 0.673 for Endo and a mean DSC of 0.541 for Epi ( Figure 6). The PDE-based method effectively removes overflow voxels that the convex hull-based method failed to remove, showing improved results around open holes (Figure 7).

D. COMPUTATIONAL TIME AND NUMBER OF USER INTERACTIONS
The Endo and Epi boundary detection was within 10.15±2.92 s, and the 3D Laplace AWT calculation was within 17.00±6.85 s (Figure 1) in 20 patients. In our algorithm, the mask image with a segmented LA wall was an input, the number of user interactions was zero, and only monitoring of the result was needed. Conversely, the previously proven manual software took about 20 minutes to calculate the AWT   and required 687.2±46.79 user interactions (including mouse and keyboard events) for the LA segmentation that was conducted by two experts.

IV. DISCUSSION
We developed a novel algorithm to measure the 3D AWT from the cardiac CTA images using automatic detection of the Endo and Epi boundaries with Convex Hull and Poisson solvers and distance measurements of the streamlines calculated with Laplace equations and then validated that it could reduce the time consumption by 1/40 without a significant difference when compared with the conventionally measured WT. A major contribution of our work was to independently separate the segmentation technique and WT measurement technique while maintaining the accuracy of the WT, which is summarized as follows: • Availability: Our algorithm (code repository: https://github.com/ohseokkwon/AutoAWT) was designed to work independently of the segmentation procedure for easy integration with commercialized atrial wall segmentation software. • Efficiency: The algorithm can reduce the time for measuring the WT. Indeed, the proposed algorithm can calculate the WT more than 40 times faster than the conventional manual method. • Generality: Our algorithm could be flexibly applied to various cardiac geometries because our algorithm was designed to automatically measure the thickness from mask images of chambers that are composed of a single label.

A. AUTOMATIC ENDOCARDIAL AND EPICARDIAL BOUNDARY DETECTION
The procedure for the AWT measurement consisted of extracting the myocardial region from the image volume and calculating the distance between the Endo and Epi boundaries. Karim et al. proposed a method for segmenting the Endo boundary and then identifying the Epi boundary as a regional growth method and measuring the WT [10]. However, this approach has the limitation that the WT measurement accuracy is highly dependent on the boundary segmentation quality. Bishop et al. generated spheres defining the mitral valve from the cropped segmentation images and then used a flood fill algorithm to identify the Endo, Epi, and remaining boundary surfaces by the size of the connected nodes [14]. Although their proposed method is independent of the segmentation procedure, it still requires an additional landmark detection algorithm to identify the location of the mitral valve. Wang et al. [8] proposed the first-ever approach for automatically extracting the inner and outer membrane surfaces based on the Convex Hull, which operates independently of the segmentation method. However, the naive Convex Hull-based method is error-prone when considering various anatomical shapes, such as the concave opening of the atrium. In contrast, our proposed method combined the Convex Hull and Poisson solvers and therefore has the advantage of being robust to local protrusions because the atrial wall is used as an adiabatic wall and the external and internal boundaries are calculated through PDEs. This method suggests that it is possible to distinguish the inner and outer surfaces of a porous structure that cannot be distinguished topologically and that it could be applied to various shapes, such as the ventricle.

B. ATRIAL WALL THICKNESS MEASUREMENTS
Previously measured AWT values reported from individual patient cohorts have a significant variability between investigators due to the different methodologies used, chronic progressive changes that take place in the AWT, complexity involving fibrous and adipose tissue, and ambiguity of the inner and outer boundary definitions [23][24][25]. An in-depth analysis of the inter-population variability of the WT requires high-resolution CTA imaging and 3D global measurements along with the acquisition of large numbers of samples. Despite this knowledge, in several previous studies, the AWT measurements were performed only at selected 2D or 3D positions using digital calipers. [10]. Although there are simple and intuitive methods (such as the nearest neighbor and surface normal approaches) for 3D global measurements, these methods have the disadvantage of being error-prone for significantly complex structures [15,16]. In contrast, the Laplace solution approach is a robust approach for thickness measurements because it considers both the Endo and Epi surfaces [8]. The accuracy of the Laplace solution approach has been demonstrated in atrial wall, ventricular wall, and cortical thickness studies [14,26,27]. However, a complex segmentation procedure is involved because the Endo and Epi boundaries must be obtained to compute the gradient field and streamline trajectories of the Laplace solution. Our method simplified this procedure by applying an algorithm to automatically detect the Endo and Epi boundaries. We reported the 3D AWT of 20 patients with AF with an error of <10% when compared with the manual WT measurements. However, the relatively high error observed between the left PV and LAA junction was a reasonable difference considering the morphology of the most complex LAA in the LA geometry [28] and the variability of the AWT in the anterior and posterior regions of the left PV [29,30]. Another potential cause for the high error percentage could be the presence of a poor segmentation between the left PV and LAA due to the practical limitations of cross-sectional imaging [31] and the subjectivity of manual segmentation. Alternatively, it is estimated as the effect of the accumulated error of the proposed procedure as the boundary conditions change due to an inaccurate boundary determination as the left PV and LAA junction.

C. LIMITATIONS AND FUTURE WORK
Our study results showed the potential to dramatically reduce the AWT detection time but had some limitations. First, although the acquisition timing of the CTA images was measured after the ECG gating, the AWT in an arrhythmia state was not standardized. Second, depending on the size and number of open holes, the Endo boundary could have been underestimated in the Endo and Epi detection procedure. Finally, depending on the spatial resolution of the CTA image, the accuracy of the WT could have been degraded near the area where the WT changed rapidly or around a salience structure. We suggest that this limitation could be overcome with improved segmentation techniques and better spatial resolution of the CTA. By overcoming those limitations, the currently proposed algorithm for automated measurements of the AWT is expected to be used in various fields of application, including human arrhythmia research or clinical practice. The AWT measurements for real-time monitoring of the RF energy titration during AF catheter ablation requires reliable, fast, and convenient software. We expect to be able to use this software to construct large-scale datasets with fully automated AWT measurement software. In the future, we plan to conduct large-scale clinical studies using a high-precision segmentation technology based on AI to improve the WT accuracy dependent on the segmentation.

V. CONCLUSION
We proposed and verified a novel automated algorithm for AWT measurements with the atrial wall binary masks segmented in the cardiac CTA images. Since our approach was designed to work independently of the segmentation procedure, it could be easily integrated with clinically proven commercial atrial segmentation software. It could also efficiently improve the labor-intensive AWT measurement procedures and drastically reduce the time required for manual AWT measurements. In addition, our algorithm can be flexibly applied to various geometries with chambers because it is designed to automatically distinguish between the interior and exterior without user intervention in the open geometry of a 2D manifold with chambers.