Surface Reconstruction for Planning and Navigation of Liver Resections

Computer-assisted systems for planning and navigation of liver resection procedures rely on the use of patient-specific 3D geometric models obtained from computed tomography. In this work, we propose the application of Poisson surface reconstruction (PSR) to obtain 3D models of the liver surface with applications to planning and navigation of liver surgery. In order to apply PSR, the introduction of an efficient transformation of the segmentation data, based on computation of gradient fields, is proposed. One of the advantages of PSR is that it requires only one control parameter, allowing the process to be fully automatic once the optimal value is estimated. Validation of our results is performed via comparison with 3D models obtained by state-of-art Marching Cubes incorporating Laplacian smoothing and decimation (MCSD). Our results show that PSR provides smooth liver models with better accuracy/complexity trade-off than those obtained by MCSD. After estimating the optimal parameter, automatic reconstruction of liver surfaces using PSR is achieved keeping similar processing time as MCSD. Models from this automatic approach show an average reduction of 79.59% of the polygons compared to the MCSD models presenting similar smoothness properties. Concerning visual quality, on one hand, and despite this reduction in polygons, clinicians perceive the quality of automatic PSR models to be the same as complex ∗Corresponding author Email address: rafael.palomar@hig.no (Rafael Palomar) Preprint submitted to Journal of LATEX Templates March 3, 2016 MCSD models. On the other hand, clinicians perceive a significant improvement on visual quality for automatic PSR models compared to optimal (obtained in terms of accuracy/complexity) MCSD models. The median reconstruction error using automatic PSR was as low as 1.03 ± 0.23 mm, which makes the method suitable for clinical applications. Automatic PSR is currently employed at Oslo University Hospital to obtain patient-specific liver models in selected patients undergoing laparoscopic liver resection.


Introduction
Liver cancer is the second most common cause of cancer death worldwide [1].
Hepatocellular carcinoma (HCC), which accounts for 70% to 80% of the cases [2], presents a 5-year survival rate below 12% [3]. Colorectal cancer metastatic to the liver, on the other hand, develops in 50% of the cases of colorectal cancer, and presents 5 5-year survival rates up to 58% for selected patients undergoing liver resection [4].
Liver resection is the treatment of choice for patients with localized HCC [5] and can potentially be a curative therapy. A successful surgical resection requires the complete removal of the tumoral cells including a safety margin, preserving as much healthy tissue as possible [6]. Resections with adequate tumor-free margins lead to a 10 better prognosis [7].
For more than a decade, computer-assisted surgical systems have been helping surgeons and other clinicians in the decision making process for planning and guiding surgical interventions. In the case of liver resection, these systems have recently found their way into the clinical practice, providing different patient-specific models: geomet- 15 ric, mechanical and functional as well as simulations. These models, ultimately rely on pre-operative computed tomography (CT) or magnetic resonance imaging (MRI).
The use of patient-specific 3D geometric models improve the capacity of surgeons to understand the liver and the underlying vascular structures. Using three-dimensional (3D) reconstructions have shown not only improvements in tumor localization and pre- 20 cision of surgery planning [8,9,10], but also an improved orientation and confidence of the surgeon while operating [11].
Marescaux et al. [12] were the first in implementing a system for planning and visualization applied to liver resection. In their approach, 3D models of the liver surface and vessels (tumors were artificially introduced) were built from MRI. The 3D models, 25 based on simplex meshes [13], could be manipulated for visualization and allowed deformations for simulation purposes. Later, Selle et al. [14] introduced the analysis of the hepatic vascular structures, which enabled the identification of the liver segments.
More recent approaches allowed the definition of a virtual resection as well as new visualization techniques. Reitinger et al. [15] presented a virtual reality environment in 30 which the planning was performed by direct manipulation of 3D structures. Lamata et al. [16] also introduced a progressive clipping visualization based on the advancement of resection during surgery. Hansen et al. [17] encoded distance to critical structures in the surface representing the virtual resection.
Most of the visualization systems for planning and navigation of liver resection 35 interventions, including all the aforementioned works, rely on the construction of 3D geometric models. These models have traditionally been obtained by isosurface extraction from segmented anatomical structures (parenchyma and venous vasculature) and tumors. Three-dimensional modeling of liver vasculature has been widely studied in the scientific literature [18,19,20,21,22,23], however, very little attention has been 40 given to the modeling of parenchyma and tumors.
Marching cubes (MC) is the most widespread technique to obtain 3D models of anatomical structures in the liver. The method was originally developed by Lorensen and Cline [24] as an isosurface extraction method from scalar volumetric data. A number of extensions to the original MC have been proposed along the scientific literature 45 [25]. MC provides with accurate isosurfaces which strictly adheres to the underlying data. In the context of isosurface extraction from medical images, MC leads to surfaces with a high number of polygons and vertices, and the appearance of staircase artifacts.
Different mechanisms like mesh smoothing [26,27,28] and decimation [29] have been employed to palliate these effects at the expense of adding more processing stages and 50 control parameters. and stereo cameras [33]. The noise resilience of PSR makes it a remarkable method for coping with the inherent noise acquired by these technologies. In the medical imaging context, PSR has found very little application perhaps due to the fact that, in this domain, the input data consist of 3D scalar fields rather than oriented cloud of points. In 60 order to apply PSR to medical images, Leonardi et al. [34] performed an estimation of the oriented clouds of points by fitting a plane to a segmented surface point and its knearest neighbors. Other approaches like [35] propose an extension to PSR to perform adaptive polygonization with application to vessel modeling.

65
The goal of this work is to address two common problems in 3D modeling techniques in the context of surgery planning and navigation. On one hand, integration of medical imaging, segmentation and 3D modeling methods into clinical workflows requires an elevated degree of automation. As an alternative to automation, expert users can perform manual adjustments of the reconstruction parameters. On the other 70 hand, models employed in planning and navigation are increasingly becoming more complex, both in terms of resolution and dimensionality of the underlying data. Processing these models is very demanding in terms of computing power, especially with the introduction of real-time constraints.
To overcome these problems and with the aim of facilitating the integration of 3D 75 modeling techniques in clinical workflows, we propose the application of PSR as a technique to model the liver surface (parenchyma). In order to adapt the segmented images to clouds of points (as required by PSR), a method based on gradient fields computation is described (Section 2.1).
Unlike traditional analysis of 3D modeling techniques, where different reconstruc-80 tion parameters are considered separately, we propose a multi-objective optimization framework consisting of a bi-dimensional accuracy/complexity optimization space (Section 3) . The aim is to obtain optimal parameters which can lead to automatic 3D model reconstructions.

85
The approach presented in this work is used to obtain a 3D geometric model of the liver surface from a segmented volumetric image. Fig. 1 shows the processing stages of both our PSR approach and state-of-the-art MCSD.
Isosurface extraction is considered to be an automatic process, however, adding processing stages either after or before, introduces additional parameters which make 90 the process more complicated and thus limiting its applicability and integration into clinical workflows. As opposed to MCSD, our strategy based on PSR requires only one parameter.
In the following, we describe our proposal. We start first by describing the computation of the oriented cloud of points from segmented images. Then we briefly describe 95 the application of PSR to obtain the 3D models of the liver parenchyma.

Oriented Cloud of Points from Liver Segmentation
As depicted in Fig. 1, medical images, obtained from CT, are represented by the scalar field defined as F : R 3 → R in which the point p i = (x i , y i , z i ) with i = 1, 2, ..., N is given a value representing an intensity level F(p i ) = v. Through a segmentation 100 process, the different points of the image, are assigned a class value (label) l i according to either anatomical or functional criteria. The segmentation is then defined as a new scalar field S : R 3 → {l 1 , ..., l k } with k the number of classes. In our case, since we are only interested in the liver parenchyma, our label map is constrained to the set {l p = 0, l b = 1} which only considers the classes parenchyma and background. (1) On the other hand, the points p i laying on the surface ∂M, can be discriminated in terms of the gradient of the segmented image: Using the properties established in Eq. 1 and Eq. 2, the definition of a oriented cloud of points V can be entirely expressed in terms of the gradient of the binary image: with i = 1, 2, ..., N.
In contrast to Leonardi et al. [34], in which the authors estimate the normals by local plane fitting of the k-nearest neighbors, we propose the obtention of oriented clouds of points through the approximation of the gradient of segmented images. Our implementation obtains the gradient of a 3D binary image through the application of the 3D Sobel operator. This operator approximates the partial derivatives individually through convolution of the segmented parenchyma with kernel functions: where K x , K y , K z ∈ R 3×3×3 are the third-order tensors representing the 3D Sobel kernel. The discrete nature of the Sobel operator makes that the condition stated in Eq. 2 becomes true for voxels near the surface of the liver (both inside and outside). Therefore, our implementation considers Eq. 3 only for those voxels belonging to the liver, this is, S (p i ) = l p .

Poisson Surface Reconstruction
with . the Euclidean norm.
Khazdan et al. [30] makes use of the divergence operator to transform this problem into a variational problem optimized by solving the Poisson equation: where X is an indicator function and V is a vector field. given by the functionX ∈ F O,F such that: Finally, in order to obtain the surface approximation ∂P, isosurface extraction is applied to the average value ofX at the sample positions:

Experimental Setup
In this work, we propose the application of PSR for planning and navigation of liver resection procedures compared to the state-of-the-art MCSD. For comparison purposes, we employ the implementation (C++) of MCSD included in 3D Slicer [37] in its version 4.4.0.
Our PSR implementation (C++), which incorporates the computation of oriented cloud of points described in 2.1, is based on Doria and Gelas [38], who adapted the original by Kazhdan et al. [30].

155
The experiments are performed on a data set consisting of six CT volumes acquired from the abdomen of patients undergoing liver resection under the Oslo-CoMet study (NCT01516710) [39]. These images were acquired with parameters normally employed for diagnostics and surgical purposes. The images, presenting different image spacing and different liver volumes (Table 1), were segmented (parenchyma, ves-160 sels, tumors) by biomedical engineers in a semi-automatic way using ITK-Snap [40] and reviewed by two laparoscopic surgeons. All the data sets present anisotropic image spacing, and therefore, are prone to generate staircase artifacts which may be distracting and confusing, since these do not have any anatomical foundation.

165
In this section, the evaluation criteria to compare the PSR approach described in the Section 2 to MCSD are established. There are different ways to evaluate and compare 3D geometric models, however, in this work, we focus on the desirable properties for planning an navigation of liver resection procedures, this is, accuracy, mesh complexity and smoothness. 170

Mesh Complexity and Accuracy
Precise planning and navigation of liver resection procedures requires high accuracy of the 3D models. Some works like Wu et al. [23] rely on the Hausdorff distance as a metric for accuracy, however, Hausdorff distance is known to be sensitive to noise, which is inherent to medical images. In the same line as Oeltze and Preim [18] and 175 Schumann et al. [19], we approximate the 3D reconstruction accuracy using surface distance based metrics, particularly the point-to-surface distanceδ(∂M) has been employed. This distance is computed from all points p i in a reference cloud of points V to all the points q in the model ∂M (either obtained by PSR or MCSD): where L is the number of points in ∂M.

180
In mesh modeling, mesh complexity and accuracy are variables related to each other. Decreasing the mesh complexity generally leads to a decreased accuracy (and vice versa). Therefore, obtaining an acceptable trade-off between complexity and accuracy can then be established in terms of a multi-objective optimization problem (i.e. obtaining low complexity and high accuracy).
units on different scales. To bring them into a unit-less objective space on the same scale, linear normalization is applied (Fig. 4), thus obtaining the normalized mean point-to-surface distanceδ ∈ [0, 1] and the normalized number of polygonsN ∈ [0, 1]: where N is the number of points in ∂M, andδ min ,δ max , N min and N max refer to the ex-190 treme values, considering all models (including MCSD and PSR) for a given parenchyma

M.
By using normalization, it is established that both objectives are equally important and so, the optimal reconstruction is obtained by minimizing the objective score σ defined as: Models presenting a score value closer to 0, exhibit better trade-off between accuracy and complexity. Then, the best model among all PSR models ∂P best , the best model among all MCSD models ∂C best and the absolute best model ∂M best can be computed as follows: where P and C are the sets of all PSR and MCSD respectively, for a given parenchyma.

Mesh Smoothness
Liver parenchyma is inherently a smooth organ. Smoothness in the 3D model is essential for visualization since it contributes to the natural appearance of the organ. Evaluation of the smoothness is performed through the discrete mean Gaussian curvature as a metric of the roughness of the model ∂M: where γ j denote the angles between pairs of edges converging at p i and A(p i ) is the sum of the areas of triangles having p i as vertex. Values of the mean Gaussian curvature closer to 0 are then interpreted as smoother appearance of the model.
Despite the fact that other curvature metrics (e.g. mean curvature) are available [41], gaussian curvature seems to be the most widespread in comparisons similar to the 210 one presented in this work [42,21,23].
For further comparisons in terms of accuracy and complexity, it will be useful to consider also the MCSD models ∂C sim which present the most similar smoothness to the best PSR models ∂P best . These models ∂C sim can be computed as: 3.3. Processing Time 3D models supporting planning and navigation of liver resection procedures are computed pre-operatively (i.e. before the operation) and therefore, there is no need to attend to real-time constraints. However, it is desirable that new methods involved in 215 surgery planning and navigation, do not alter the clinical workflow significantly in a negative way.
Medical image processing is often performed in terms of a region of interest ROI, which represents a reduced set of the original data. In general, this greatly improves the performance of algorithms in terms of processing time, since the number of ele-220 ments to process decreases dramatically. Therefore, the evaluation of the performance (processing time) presented in this work is performed in the basis of minimal ROIs containing the liver.

Visual Quality
Perceived visual quality of the 3D models is an important evaluation criterion since 225 it can affect not only the interpretation of the anatomy of the patient but also the confort of the clinicians working with the 3D models. In the scientific literature, visual quality is often discussed in terms of appearance of artifacts in the model and smoothness properties [19,21,22,23].
Some works explicitly introduce visual quality criteria in the evaluation, like in 230 Oeltze et al. [18] for 3D modeling or Feng et al. [43] for 3D visualization in laparoscopic surgery. In

Mesh Complexity and Accuracy
PSR was applied to the 6 liver parenchyma of the data set. A total of 24 surface models were obtained considering PSR at depths d ∈ {5, 6, 7, 8}. From the application of Eq. 12, the objective score value for all PSR and MCSD models was computed. The results, presented in Table 2, show that PSR with d = 7 always present the absolute best score values among all PSR and MCSD models, and 275 therefore better behavior in terms of complexity and accuracy. PSR reconstructions different from d = 7 present lower score values than those of the best MCSD. Another interesting comparison concerning accuracy and mesh complexity can be derived from the best MCSD (∂C best ), the best PSR models (∂P best ) and those MCSD models similar smoothness to ∂P best (Eq. 14). The results reveal that the main contribu-  Fig. 6 illustrates the compared models derived from P 1 .

285
A broader view of the evaluation results is presented in Fig. 7. In this figure, the performance of all PSR and MCSD are shown in terms of complexity and accuracy.
We computed the Pareto frontier, this is, the set of models which are not dominated by   I  II  III  I  II  III  I  II Fig. 7 show that all PSR models ∂P, together with some MCSD, including the best MCSD ∂C best models are part of the Pareto frontier.
None of the most similar MCSD ∂C sim are part of the Pareto frontier.

Mesh Smoothness
In order to evaluate the smoothness of the models, the mean Gaussian curvatureK 295 was computed for all PSR and MCSD models. A comparison of smoothness of PSR models with that of the best MCSD ∂C best and most similar MCSD models ∂C sim in terms of accuracy and complexity was performed. The results, presented in Table 4, show that, on one hand, the best MCSD ∂C best generally presents much higher curvature values than most of PSR and the most similar MCSD ∂C sim . On the other hand,  relatively less smoothness as their correspondent best PSR ∂P best (0.03 × 10 −3 average curvature).

Processing time
Our performance results (shown in Table 5) were obtained using a CPU Intel R Mean execution time values show higher performance for similar MCSD models ∂C sim (t ∂C sim = 12.05s) over best MCSD models ∂C best (t ∂C best = 14.08s) and best PSR 310 models ∂P best (t ∂P best = 14.50s), while variability of execution time is significantly lower for ∂P best models (s ∂C best = 1.24s) compared to similar MCSD ∂C sim (s ∂C sim =  4.06s) and best MCSD ∂C best (s ∂C best = 4.45s).

Visual Quality
Subjective data obtained by questionnaires (Fig. 8a) show lower mean quality score 315 for ∂C best (s C best = 2.729) versus ∂C sim (s C sim = 3.291) and ∂P best (s P best = 3.333).
In the same line as Feng et al. [43], we test the statistical difference between groups of methods regarding the perceived quality. In our case, the comparison between methods is peformed by means of Welch two sample t-test using the R statistical environment (Fig. 8b). Statistical significance was obtained in the comparison between 320 ∂C best and ∂P best (p = 0.0006661), and in the comparison between ∂C sim and ∂C best (p = 0.001382). No statistical significance was found in the comparison between ∂P best and ∂C sim (p = 0.777).

325
The application of PSR to reconstruction of segmented images requires the transformation of the binary image (segmentation) to a cloud of oriented points. Perhaps due Our results suggest that PSR outperforms state-of-the-art MCSD for modeling liver parenchyma in different aspects. Overall, ∂P best present reconstructions with better ac-345 curacy/complexity trade-off than state-of-the-art MCSD. Comparing models with similar smoothness and visual quality (∂C sim and ∂P best ), ∂P best models present a dramatic reduction of complexity (73% less triangles on average) at the cost of decreasing the accuracy slightly. Accuracy values for PSR, however, remain within clinically acceptable limits. Despite ∂C best models can achieve similar complexity as ∂P best in some 350 cases, subjective evaluation experiments reveal that perceived visual quality of ∂C best models is lower with statistical significance than those of ∂P best and ∂C sim .
Based on the accuracy/complexity space ( Figure 5) two important observations can be made. First, PSR as method, defines the Pareto frontier, which suggests the goodness of the method over MCSD. Secondly, parameters to obtain ∂C best depend on the input 355 data while PSR obtains ∂P best using the same parameter (d = 7), therefore indicating that PSR, as method, is more stable than MCSD. One implication of the stability of PSR is the possibility to perform optimal parameter estimation which can be used to produce automatic reconstructions (see Section 5.4).
As it is shown earlier in Figure 1, state-of-art MCSD requires two parameters

PSR Compared to Other Reconstruction Methods
The wide variety of reconstruction techniques in the literature makes the choice of the method a relevant question. The answer to this question is often driven and limited by the scope of the application. There are a number of surface reconstruction methods other than MCSD or PSR. In particular, multi-level partition of unity implicits 375 (MPUI) [45] has been previously adopted for vessel modeling from clouds of points [19,42,20,46,23] and can be applied to liver modeling as well. This and other reconstruction methods are compared and discussed in [30]. In this comparison, MPUI shows generation of spureous surface sheets under noisy conditions, which are inherent in medical imaging, while PSR exhibits good noise resilience.

380
More recent approaches like dynamic particles [47,48] provide high-quality meshes suitable for visualization and simulation purposes. These methods show more flexibility (regularity of triangulation and accuracy of reconstructions can be controlled) at the cost of a more complex parameter space. Although these approaches can be better suited for simulation, widening the parameter space complicates the integration into 385 clinical workflows.

Accuracy, Complexity and Depth Parameter d
The choice of the depth parameter has an impact over the smoothness, accuracy and complexity of the reconstruction. Increasing the depth parameter value generally produces similar or superior smoothness, higher accuracy and higher complexity of the 390 model. Similar works as this study accuracy and complexity as independent dimensions of the problem [19,42,21,23]. and visualization as in [49].

Integration of PSR in Clinical Workflows
The application domain of our study is planning and navigation of liver resection procedures. In this domain, visual realism and accuracy are of paramount importance.
Visual realism is achieved through smoothness, which removes staircase artifacts and The difference in computing time between PSR and MCSD is, for pre-operative purposes, neglectable (Table 5) Works like [31] and [50] show the feasibility of using GPUs for real-time reconstruction of complex scenes using PSR. As intra-operative systems (e.g. surgical navigation) move towards the use of deformable models, aspects related to algorithm performance 435 and parallelization capabilities will get more relevance in medial systems design.
Process automation during imaging, segmentation and 3D modeling is the key for improving the adoption of 3D patient-specific models in clinical workflows. to obtain segmented images (Figure 1). Open source software like 3D Slicer [37] or ITK-Snap [40] can also make use of this work not only for the reconstruction of liver surfaces, but also to investigate further applications.
Liver surgery is moving towards minimally invasive surgery (laparoscopy). The clinical advantages of this approach must take into consideration the increased com-450 plexity of the surgical procedure (reduced maneouvrability and visual field). In this context, the use of 3D patient-specific models, together with intra-operative imaging (ultrasound) are becoming increasingly relevant. Recent studies highlight the advantages of integrating these models as part of augmented reality guidance systems in laparoscopy [55] and open surgery [53]. The combination of 3D patient-specific model, 455 together with the latest trends in surgery, like robotic surgery, have the potential to make surgical interventions easier, faster and probably safer [56]. Some of these works highlight the importance of accuracy, and reduction of human interaction (automation of processes), topics which are widely studied and discussed in this work.

460
Our application of PSR to liver modeling (parenchyma) from scalar volumes inevitably raises the question of the adequacy of PSR to obtain 3D models of other anatomical structures. In the context of planning and navigation of liver resection pro- another set of techniques (known as model-based methods) [20], widely used in the context of vessel modeling for surgery planning.

Future Work
Leonardi et al.
[34] suggest the use of PSR to construct geometric models of other organs than kidney. In the same line, and despite the little attention PSR has been 485 given in the medical domain, we believe that its use can be extended to other organs outperforming state-of-the-art methods.
Smooth organs absent of sharp features are, in principle, good candidates to undergo PSR. Evaluations similar to Wu et al. [23] could also consider the intrinsic relationship between number of polygons and error according to our multi-objective 490 optimization framework. To the best of our knowledge, modeling of tumors has not been subject to an exhaustive evaluation like the one presented in this work or in Wu et al. [23], which can be of great interest.
New reconstruction methods that may arise, can be evaluated using this work as guideline. Despite the more complex parameter space of methods based on dynamic 495 particles [47,48], these can support an interesting comparison with PSR for modeling of anatomical structures for different purposes.

Conclusion
In this work, we propose the application of PSR to obtain patient-specific models of liver parenchyma for planning and navigation of liver resection procedures. For 500 the application of PSR to medical images, we propose an efficient transformation of the segmented images to oriented cloud of points based on computing gradient fields.
In order to make an automatic PSR, we found the PSR parameter obtaining the best accuracy/complexity trade-off (d = 7).
Comparing PSR with the state-of-the-art (MCSD) in terms of accuracy, complexity 505 and smoothness, PSR shows better reconstruction performance and stability of results.
This study also reveals that PSR liver models using the optimal parameter d = 7 not only are smooth, but also present better accuracy/complexity trade-off than MCSD models. Reconstructions obtained through automatic PSR (d = 7), presents median errors within 1.03 ± 0.23 mm, which makes them suitable for clinical applications. On 510 average, these models have 79.59% less polygons compared to MCSD models with similar smoothness, while clinicians do not perceive a significant quality difference.
Optimal PSR models d = 7, exhibit a significant improvement of visual quality compared to optimal MCSD in terms of accuracy/complexity trade-off. Automatic PSR can be seamlessly integrated in clinical workflows already using MCSD, since the process-515 ing time is similar to that of MCSD.
The contribution of this work, is therefore, a step towards the automation and quality needed for a wide adoption of 3D patient-specific models in the medical community.
Currently, at Oslo University Hospital (The Intervention Centre), PSR is employed in a fully automatic way (after segmentation, which takes place in a semi-automated way) 520 to obtain patient-specific models of liver parenchyma in selected patients undergoing laparoscopic liver resection. Fig. 9a shows a complete patient-specific liver model in which the parenchyma was obtained through PSR (d = 7) while MCSD was applied for vessels and tumors. During operation, (Fig 9b) the patient-specific model, which includes the resection path, helps the surgeons to perform the resection according to 525 the pre-operative plan.