surgery, including tracked ultrasound [1], [2], cone-beam CT (CBCT) [3], and optical methods [4], [5], [6], [7], [8].

To integrate subsurface information from preoperative CT/MRI into the surgeon's intraoperative view, registration methods are employed to map critical information from the preoperative images into the intraoperative scene. Rigid registration is often used to align data collected in different coordinates. However, the initial configuration of the liver at the start of the intervention can significantly deviate from its preoperative state. Studies indicate that the deformation of the anterior surface of the liver can exceed 10 mm during laparoscopic liver surgery (LLS) and 7 mm during open liver surgery (OLS) [9]. In clinical contexts, the target registration error is expected to be below 5 mm [4]. Therefore, nonrigid registration may become essential to correct for such deformations.

Registration methods are broadly categorized into intensitybased [10] and geometric-based [11], [12], [13] approaches based on the type of information utilized to establish correspondences between the pre-and intraoperative data. Intensity-based methods employ pixel or voxel intensities for aligning images. However, challenges arise from intensity variations across different imaging modalities, such as CT to ultrasound, impacting the efficacy of these methods.

Geometric-based registration methods can be further classified based on the availability of 3D information extracted from intraoperative imaging. For scenarios where only 2D information is available, as with intraoperative monocular laparoscopes, 3D-2D registration methods are typically employed [11], [12]. In contrast, when intraoperative imaging modalities provide 3D information, usually presented as point clouds, 3D-3D registration methods [1], [5], [6] are preferred, as including 3D information introduces additional constraints to enhance registration accuracy [13].

This study focuses on the 3D-3D non-rigid registration problem in IGLS. Unlike the typical 3D-3D non-rigid registration problem in the computer vision [14], where the primary focus is achieving surface matching accuracy, the problem in IGLS emphasizes the necessity for the estimated organ deformation to be realistic, with a crucial consideration for the deformation beneath the tissue surface. Biomechanical models, specifically represented as finite element models (FEMs), have been extensively employed in modeling tissue deformation [15] and are key to solving the 3D-3D nonrigid registration problem [4], [5], [16]. The biomechanical model not only constraints displacements to enforce realistic deformations but also provides insights into volumetric deformation, capturing the complexities of organ deformations beyond surface-level changes.

We focus on registration techniques that utilize a FEM-based biomechanical model to align a preoperative model with an intraoperative point cloud acquired during the surgical procedure. This alignment subsequently aids in determining the locations of subsurface anatomical structures within the intraoperatively deformed liver.

In energy minimization approaches, it is common to integrate a FEM as a strain/deformation energy term to regularize a data term that ensures surface alignment between the preoperative model and the observed intraoperative point cloud data. The strain energy arises within the preoperative model due to the nonrigid displacement field. Rucker et al. [4] assumed the posterior side of the organ drove the deformation, then iteratively updated parameters that described the boundary condition and rigid transformation. Based on the similar parameters formation, the framework was further extended to impose constraints from the ligament [9] of the liver and vascular features [1]. The method proposed by Modrzejewski et al. [6] includes an iterative closest points (ICP)-based data term, a collision term preventing self-collision, and a strain energy term. Khallaghi et al. [17] proposed a framework referred to as the GMM-FEM method in the context of MRI and transrectal ultrasound (TRUS) prostate registration. This method incorporates the FEM into a coherent point drift (CPD) algorithm [18] as the regularization term and uses a Gaussian-mixture model to represent the preoperative surface. The strain energy term can also be minimized without the data term, as shown in the work presented by Peterlik et al. [3], where strain energy is minimized while satisfying the geometric constraints searched from ICP.

A limited number of methods directly incorporate the FEM into the data term. Mestdagh and Cotin [16], Mestdagh [19] integrated the FEM into the data term and solved it in an optimal control formulation. However, this approach necessitates manually identifying zero boundary conditions (ZBCs) and force locations.

Instead of incorporating the FEM as a data term or a strain/deformation term, Suwelack et al. [5] formulated the movement of the preoperative model to the intraoperative surface as an electrostatic-elastic problem, which the FEM solves. The preoperative model is assumed to be electrically charged and is attracted to the oppositely charged intraoperative surface by its electric potential field.

Despite being in the era of deep learning, learning-based methods for liver deformation registration are still constrained by several factors: the scarcity of large datasets for training and testing, issues with interpretability, and suboptimal accuracy. Additionally, effectively combining learning-based methods with biomechanical model-based approaches to address registration problems still remains a challenge. In learning-based techniques, the biomechanical model is often employed to simulate various deformations for training the neural network. Pfeiffer et al. [20] proposed V2Snet (volume-to-surface registration network), designed to estimate the deformation of a volume mesh to an intraoperative surface. Tagliabue et al. [21] introduced BA-Net (binary-attachment network) to predict the locations of the attachment points during tissue dissection. The predicted attachment points can be used as ZBCs to update a patient-specific biomechanical model intraoperatively [7].

Although utilizing the FEM as strain energy might restrict deformation, it fails to ensure coherent movement of volumetric vertices to maintain the original geometry. On the other hand, the direct solution of the FEM-based data term requires the identification of Zero-Boundary Conditions (ZBCs) and force locations, as demonstrated in the existing method by Mestdagh and Cotin [16].

Despite the numerous 3D-3D non-rigid registration methods that have been proposed, the liver registration community has limited access to their implementations, with only a few options available [16], [17], [20]. Furthermore, there is a scarcity of publicly available liver registration datasets [5], [6], [22].

In this paper, we introduce a 3D-3D non-rigid registration method based upon incorporating a modified FEM into the surface matching/data term. Most non-rigid registration methods incorporate the FEM model into the strain energy/regularization term. Only a few methods [16], [19] include the FEM in the surface matching data term, which, to our best knowledge, requires the manual identification and specification of boundary conditions.

Here, we alleviate the need to prescribe boundary conditions by using a modified stiffness matrix of the FEM, which is incorporated into the data term. We use a well-known numerical method (diagonal loading) to stabilize the stiffness matrix. As a result, the modified data term enables us to not need to manually identify and prescribe boundary conditions, and also to not need to use an additional strain energy term to be minimized.

To solve our formulation, we employ the Nesterov accelerated gradient algorithm and derive the necessary gradients. Additionally, we propose a novel optimal step size strategy that eliminates the need for manual step size tuning.

We validate our method using our both simulated and experimentally collected phantom datasets, as well as two publicly available datasets. Comparative assessments are conducted against two open-source methods -a learning-based approach and an optimization-based method utilizing FEM for regularization, along with state-of-the-art closed-source methods. This straightforward yet effective approach consistently demonstrates its effectiveness relative to the benchmarks.

Given the limited availability of public datasets and accessible non-rigid registration methods, we aim to contribute to the field by publicly disseminating our constructed liver dataset and the developed registration algorithm, thereby advancing liver registration benchmarks.

Before liver surgery, patient-specific biomechanical finite element models are created using preoperative CT and/or MRI scans. These models are based on the patient's liver geometry derived from imaging scans, along with the material properties assigned to the tissue.

The preoperative liver geometry is represented by a volumetric tetrahedral mesh model comprising n nodes x = {x 1 , . . . , x n |x i ∈ R 3 } ∈ R 3n . The preoperative liver surface is represented by a triangular mesh model ∂ composed of n s nodes x s ∈ R 3n s , which is a subset of the volumetric mesh nodes.

During surgery, the observed intraoperative surface is represented by a point cloud y = {y 1 , . . . , y m |y i ∈ R 3 } ∈ R 3m . The intraoperative surface can be obtained from various sources, such as CBCT, an RGBD depth camera, structure light, stereo laparoscope, or manual swabbing with an optically tracked probe. The visible intraoperative surface varies in surgery settings [9], [23]: for open surgery, the visible surface could reach over 50%, while for laparoscopic surgery, typically, the extent of visible surface is much more limited, varying from 20 to 30%.

The objective is to register a preoperative model with the deformed intraoperative surface data to map the sub-surface anatomical structures identified in the preoperative images into the intraoperative scene. A preliminary alignment of the preoperative model with the intraoperative surface is necessary before the non-rigid step, typically achieved through a rigid transformation. Existing IGLS systems commonly perform this alignment manually [24], [25]; nevertheless, the manual alignment can be further improved with a variant of the ICP [26], [27].

When surface forces f are applied to the liver boundary, the resulting displacements u = {u 1 , . . . , u n |u i ∈ R 3 } ∈ R 3n are a unique solution of the static equilibrium equation:

where P(u) = {P 1 (u), . . . , P n (u)} is a vector of nonlinear functions of u and are determined by the FEM model of the liver geometry and material properties. For a liver modeled by a linear elastic material, the forces are now a linear function of u where P(u) = Ku with K is a constant 3n x 3n stiffness matrix. For a nonlinear hyperelastic material, P(u) is a nonlinear function of the displacements u. The tangent stiffness matrix K T defined as the derivative of P(u) with respect u

is a matrix used in each iteration of the solution of a nonlinear problem, and it changes with each iteration. For a linear elastic material, the tangent stiffness matrix is simply the stiffness matrix. A more detailed description of the construction of the FEM models can be in the following book

In the forward problem, the displacement vector u throughout the domain can be determined if displacements and/or forces on the boundary ∂ are known. The inverse registration problem solves for the boundary forces that generate a displacement field u * that approximates the displacement field u between the partially observed intraoperative surface and preoperative model.

The liver geometry and material properties are encoded in the stiffness matrix of the finite element model. This matrix characterizes the object's resistance to nodal deformation when subjected to nodal forces. In the case of a linearly elastic biomechanical model [28], Young's modulus (E), and Poisson's ratio (¿) account for the material's elastic properties.

The stiffness matrix, as formulated, is singular and cannot be inverted. This singularity indicates a scenario where the structure lacks defined boundary conditions or constraints, allowing the structure to undergo rigid motion without internal deformations. In such a state, applied forces can result in indeterminate displacements.

The imposition of boundary conditions, specifying the solution's values on the boundary, constrain rigid-body motions. Consequently, the stiffness matrix becomes invertible, enabling the determination of a unique solution for displacements resulting from applied forces. In the context of assessing liver deformation during surgery, it is crucial to identify the locations featuring ZBC, which typically correspond to the attachment points of the liver. This identification is vital for a meaningful and accurate organ behavior simulation.

To construct a FEM in linear elastic version as an example, where f = Ku, the values of E and ¿ are used to construct the constant 3n x 3n stiffness matrix K [28], [29]. Given forces and ZBC, displacement u can be solved.

Our method employs a force-driven approach to deform the source model for optimal alignment with the target point cloud. Fig. 1 provides an overview of our registration method. The objective function comprises a surface matching term. The surface matching term utilizes a correspondence matrix generated by a closest point operator and a FEM. Minimization of the objective function is achieved through an accelerated gradient algorithm. During the force optimization phase, we initially determine the optimal step for calculating Nesterov's accelerated gradient [30]. Subsequently, the gradient is employed to update forces. In each iteration, we first update the correspondence matrix via the closest point operator, then update the force through force optimization; finally, we update the displacements using the FEM. The FEM ensures that the displacement or deformation preserves the geometry of the preoperative model, maintaining the integrity of the nodeto-node connections. After a set number of iterations, the predicted volumetric deformation enables the preoperative volume model to align with the intraoperative point cloud. In the following section, we will detail each of these steps.

We perform the registration by determining a surface force vector f that results in a displacement vector u f , which minimizes the objective function:

where J is a data term that measures the discrepancy between the predicted deformed preoperative surface ∂ u and the observed intraoperative point cloud data. The data term is given by

where C ∈ R 3m×3n is a soft correspondence matrix, x ∈ R 3n are spatial location of the preoperative mesh nodes, u ∈ R 3n are displacements associated with the nodes, and y ∈ R 3m is the intraoperative point cloud. The correspondence matrix C maps the deformed model to points on the intraop surface will be detailed in Sec. III-B. The liver biomechanical model can be incorporated into Eq. 4 by using the relationship between u and f given by Eq.2:

The resulting data term now becomes a function of the forces applied to the surface of the liver. Our data term is a reformulation of the data term from Mestdagh and Cotin [16], Mestdagh [19] where we have explicitly incorporated the correspondence matrix and the dependence of J (u) on u.

The correspondence matrix C contains the correspondences between the points in the intraoperative point cloud and preoperative model. Ideally, if corresponding fiducial points or anatomical landmarks could be consistently identified in the pre-and intraoperative spaces, the correspondence matrix would be a one-to-one binary matrix. However, obtaining an accurate binary correspondence matrix is non-trivial unless manual annotation is used, which is highly user-dependent and prone to error.

Instead of a binary correspondence matrix, a soft correspondence matrix is often constructed, indicating the probabilistic correspondence between a point in one space and another. The elements C i j of the correspondence matrix C denote the probability that the intraoperative point y i corresponds to the preoperative mesh point x j + u j . Each row of the correspondence matrix should sum to one.

Surface matching methods such as the closest point operator and Gaussian mixture models [17], [18] use the Euclidian distance between points in the pre-and intra-operative space as an alignment metric. Feature-based learning [31] or nonlearning [32] methods are used to measure the similarity between features in both spaces and can be used to construct a soft correspondence matrix.

Here, we use the closest point operator to construct a soft correspondence matrix. For a given intraoperative point y i , the closest point operator determines its closest point ỹi :

on the preoperative surface, which is defined by the three nodes x i, j,k , displacements vectors u i, j,k , and barycentric coordinates ¼ i, j,k associated with the surface triangle that contains the closest point. The ¼ i, j,k sums to one, representing confidence values C i,(i, j,k) . Thus, the 3m x 3n correspondence matrix is extremely sparse, where each row contains only three non-zero entries corresponding to the barycentric coordinates. The closest point operator and the resulting correspondence matrix are a function of the undeformed mesh x, intraoperative point cloud y, and the nodal displacements u.

Nesterov's accelerated gradient [30] is used to minimize the objective function Eq. 3. Nesterov's accelerated gradient is a momentum-based optimizer that adds a fraction of the previous update to the current one, creating a momentum effect that accelerates the algorithm toward the minimum. The following steps

are iterated until a stopping condition is met. k represents the number of current iterations.

Nesterov's accelerated gradient is given by Eq. 7 and Eq. 8. p k carries momentum from the previous iteration by combining the force estimates of the current and previous iterations. Eq. 8 updates the forces, where ∇ J (p k ) is the gradient of the data term (Eq. 5) with respect to p, and ³ k is the step size. ∇ J (p k ) is given by

where for a linear elastic material K T is a constant and for a nonlinear material K T is a function of u k and can be updated as needed. The correspondence matrix is recalculated during the iterative registration using the current estimate of the displacements. To simplify the calculation of ∇ J (p k ) we employ an iteration-lagging method which updates the correspondence matrix for each iteration using nodal displacements u k computed from the previous iteration. The correspondence matrix is now considered fixed and its variation with respect to u can be neglected (i.e. ∂C T ∂u = 0). Now ∇ J (p k ) is given by

1) Optimal Step Size Estimation: In Nesterov's accelerated gradient, the step is usually fixed. To further accelerate the optimization, we determine the optimal step size ³ k by substituting the gradient descent term (Eq. 8) into the data term (Eq. 4) and minimizing it with respect to ³:

Setting the derivative of Eq. 11 for ³ to zero, we can obtain the optimal step size:

As mentioned in Section II-B, without boundary conditions, the stiffness matrix K is ill-conditioned and cannot be inverted to yield a unique, stable solution. However, determining ZBCs and applied force locations is challenging during liver surgery.

To address these unknowns, we modify the stiffness matrix to stabilize the model. Inspired by the concept of soft springs used by the commercial simulation software SOLIDWORKS ® [33], we stabilize the model by adding a soft spring to each node. This corresponds to adding a small stiffness term to the diagonal terms of the stiffness matrix:

where k ss is the stiffness of the soft spring and I is the identity matrix. This is analogous to the work presented by Peterlik et al. [3], where a damping term was added to the stiffness matrix, which was then integrated into a strain/deformation energy term to regularize a data term. In our work, the modified stiffness matrix (Eq. 13) is integrated directly into the date term via the FEM. The added small stiffness increases the inertia of the liver model being deformed, preventing rigid motion without internal deformations. This approach enables us to allow the forces to be applied at any spatial location on the liver surface.

As a consequence of these modifications, we no longer need to quantify the actual location of ZBC's and physical forces applied to the liver's surface; rather, we empower the optimization algorithm to identify both the spatial distribution and magnitude of applied forces that optimally deform the preoperative model to match the intraoperative point cloud. This modified tangent stiffness matrix is incorporated into the gradient of the our data term described in Eq. 9.

This formulation is a key contributor to our goal of developing an algorithm that allows us to register the preoperative and intraoperative without prior knowledge. As we shall show, these modifications yield accurate registration results.

The experiments are conducted on the following datasets: 1) In silico Phantom: As illustrated in Fig. 2, we simulated the deformation of the in silico phantom using a linear elastic model with Young's modulus E = 1 and a Poisson's ratio ¿ = 0.49. Forces are applied perpendicular to the x-y plane on a liver mesh model obtained from OpenHELP [34]. The resulting deformed surface is cropped to generate the simulated intraoperative point cloud, consisting of 934 vertices. The volumetric preoperative mesh comprises 4291 vertices and 1919 tetrahedra.

2) In vitro Phantoms: We created two liver phantoms, A and B, each embedded with 53 and 176 fiducial markers, respectively, as illustrated in Fig. 4, and fabricated using synthetic gelatin (Humimic Medical Gelatin #0) and a 3D-printed mold based on a liver model [35]. To induce deformation, wedges featuring different gradients were strategically inserted underneath parts of the posterior side of the undeformed phantom, as depicted in Fig. 3 (b). The resulting phantom, with the wedge underneath, is referred to as the deformed phantom. We used the undeformed data as the preoperative data and the deformed data as the intraoperative data. The undeformed liver meshes, featuring around 9000 faces and 4000 vertices, were processed using TetGen [36] to generate the source volumetric meshes consisting of approximately 9000 volumetric vertices and 80000 tetrahedra.

Phantom A and Phantom B underwent two deformations, respectively, generating four deformed phantom configurations. We number the deformed phantoms (No. 1-2 from phantom A and No. 3-4 from phantom B). CT scans of both undeformed and deformed phantoms were acquired, and their surfaces and fiducial landmarks were manually segmented. Undeformed and deformed phantoms were manually aligned and refined with an ICP registration implemented using the Open3D Python library [37]. Initial registration errors are detailed in Table I.

Additionally, partial anterior surfaces were cropped to achieve surface visibility ratios ranging from approximately 18% to 27%, calculated as the ratios of cropped surface areas to the full surface areas, as summarized in Table I.

3) Open-CAS Dataset: The Open-CAS dataset contains three in silico phantoms and one in vitro phantom, as described in detailed in [5]. Deformations of the in silico phantoms are generated using a non-linear biomechanical model. The in vitro phantom featuring 6 fiducial markers was fabricated from silicone, and two intraoperative surfaces from a CT and stereo-endoscopic imaging are included.

4) Sparse Dataset: We use the Sparse Non-rigid Registration Challenge [22], [38], [39], where 112 samples with ground truth are released [40]. Intraoperative surfaces are presented with sparse 3D point clouds collected using an optically tracked stylus. The dataset also provides initial rigid transformations to align the preoperative and intraoperative datasets. Lastly, the dataset also includes sub-dataset splits featuring different visibility ratios and noise magnitudes, as described in detail in [39].

The registration accuracy was estimated by calculating nodal displacement errors for the simulated cases and Target Registration Error (TRE) for the phantom cases. Their calculations share a common form:

where X i denotes a fiducial marker location or volumetric vertex position in the undeformed data, Y i represents a fiducial marker or vertex position in the deformed data, W i signifies the estimated transformation from a registration method for X i , and ∥∥ 2 denotes the Euclidean distance. Considering the four nearest points, the nearest neighbor interpolation is employed to propagate the estimated volumetric deformation to the fiducial markers.

The proposed method was implemented in C++ using the Eigen library [41], running on a 3.60 GHZ Intel i9-9900K CPU and 64 GB of RAM. The Cholesky decomposition implemented in the Eigen library, named SimplicialLLT, is used to invert the stiffness matrix from Eq. 13. An in-house implementation of a finite element model was used for the biomechanical calculations. We use a linear elastic biomechanical liver model for all experiments. Since a linear biomechanical model is employed, Young's modulus E solely influences the scale of the estimated forces and does not affect the displacement solution. For this study, we scale the force by setting E to 1. Poisson's coefficient is set as ¿ = 0.49, a typical value for modeling incompressible tissue [42]. The soft spring constant was set to k ss = 0.01, and the optimization algorithm was stopped at 200 iterations. The implemented procedure is detailed in Algorithm 1.

We use the official implementations of V2Snet [20] and GMM-FEM [17] as benchmarks for comparison. For V2Snet, we use the official weight included in the implementation. GMM-FEM [17] is an optimization-based method that incorporates the FEM as a strain energy term to regularize the CPD [18]. We implemented it using its suggested parameter values, which were tested in our in silico phantom described in Section V-A and gave the lowest errors. Furthermore, we include the Rigid Procrustes method that uses RANSAC ICP implemented in the Open3D library [37] with ground truth correspondences to estimate the lowest errors that may be achieved via the rigid registration.

We initiate our study using a noise-free and fully controlled in silico liver phantom to understand how our algorithm performs under various settings.

We first examine the algorithm's ability to predict ground truth forces when zero boundary conditions and force locations TABLE I COMPARISON OF REGISTRATION ERRORS ON THE PHANTOM DATASET, PRESENTED AS MEAN ± STANDARD DEVIATION (MAX) ERROR IN MILLIMETERS. THE LOWEST MEAN VALUES AND THE LOWEST MAXIMUM VALUES ARE HIGHLIGHTED IN BOLD Algorithm 1 Linear Elastic Surface Matching Registration Inputs: x the preoperative (source) nodes of liver mesh y the intraoperative (target) surface point cloud k ss the soft spring constant maxIters the maximum number of iterations Outputs:

u displacements for registered liver mesh 1: initialize displacements u 0 ← 0; 2: initialize forces f 0 ← 0; 3: calculate stiffness matrix K ← K + k ss ; 4: for k = 0 to maxIters do compute gradient of data term ∇ J (p k ) (Eq.10); 10: compute optimal step size ³ k (Eq.12); 11: update forces f k+1 (Eq.8); 12: end for are known. When provided with precise material parameters, exact force locations, and zero boundary conditions (ZBC), our method effectively recovers forces that closely approximate the ground truth forces, as illustrated in Fig. 5 (a). This accurate force recovery translates to minimal registration errors, lower than 2 mm, as depicted in Fig. 5 (b).

Next, we examine a more realistic situation where prior knowledge of the location of zero boundary conditions and forces is unknown. We eliminate ZBCs and include the soft spring constant k ss . We explore two scenarios: the former assumes that the forces exist solely on the posterior surface, aligning with the assumption in [4]; the latter assumes that forces can be distributed across the entire surface. As illustrated in Fig. 5 (b), allowing the forces to be applied across the full liver surface yields superior registration results compared to confining the applied forces to the posterior surface. The results also indicate accurate registration results can be obtained without prior knowledge of force locations and ZBCs. As noted before, we are no longer computing the actual physical forces applied to the liver's surface, but we rather empower the optimization algorithm to determine the fictitious forces that minimize the registration error. Subsequent experiments are conducted, assuming that forces can exist on the full surface.

Furthermore, we show that incorporating Nesterov's acceleration and updating ³ enhances the optimization process, resulting in faster convergence and improved accuracy, as depicted in Fig. 5 (c).

The computation depends on the number of FEM nodes and intraoperative point cloud, as our method relies on the FEM model and the closest point operator. The overall computation time for this dataset is 22 seconds over the 200 iterations.

Table I summarizes the results for the phantom dataset, comparing the whole and partial surface registration errors of our proposed method relative to those achieved by V2S and GMM-FEM. Our method features a mean TRE below 5 mm across all testing samples. Notably, the differences in mean TREs between our predictions for partial and full surfaces are negligible, typically less than 1 mm. In contrast, GMM-FEM tends to exhibit lower errors on full surfaces, as illustrated by the mean and max TRE on the full surface of No. 4, which are 1.15 ± 1.62 (8.14) mm, compared to the partial surface, featuring an error of 17.54 ± 3.68 (24.10) mm. V2Snet encounters challenges in accurately predicting results for No. 2, resulting in a mean TRE exceeding 5 mm and a max TRE larger than 15 mm. Fig. 6 illustrates several qualitative results that provide the reader with a deeper and more tangible interpretation of the registration performance achieved by our proposed and comparative methods. Notably, blade-like structures are evident in the deformed meshes generated by V2Snet, as indicated by arrows in both No. 2 and No. 4. Specifically, for both the full and partial surface cases of the No. 2 phantom, the estimated deformed mesh of V2Snet exhibits discrepancies in fully aligning with the ridge of the target surface. This misalignment is also noticeable in the estimation by GMM-FEM, particularly in the case of the partial surface for the No. 2 phantom. In contrast, our method consistently produces reasonable deformations for both the full and partial cases of the No. 2 phantom.

The results are summarized in Table II. In the simulation dataset, featuring over 50% visibility, both GMM-FEM and our method exhibit comparable performance, with differences in mean TRE of less than 1 mm. Notably, for Simulation 2, our method achieves a maximum TRE of around 6 mm, significantly lower than the 10 mm achieved by GMM-FEM. V2Snet struggles to generalize on Simulation 1 and demonstrates a maximum TRE exceeding 20 mm for Simulations 2 and 3. The partial CT and stereo-reconstructed surfaces present approximately 20% partial visibility, posing a challenge. In this scenario, our method outperforms all others, yielding the lowest mean and max TRE.

The qualitative results in Fig. 7 support these findings. As observed previously, V2Snet generates deformed meshes featured with blade-like faces, evident in both cases in Fig. 7. GMM-FEM exhibits a similar phenomenon for the deformed mesh on the stereo-reconstructed surface. In contrast, our method consistently produces deformed meshes that maintain a realistic shape.

The registration errors are detailed in Tables III, IV, and V. Except for GMM-FEM and our proposed method, results from other approaches align with previous reports [39], [40], [45]. The other methods are either closed-source [1], [43], [44], require additional data processing [20], or involve manual input [16]. The most similar method to our proposed method is that of Mestdagh and Cotin [16], which also uses a FEM as the data term but requires manual boundary condition identification, which is known to pose significant challenges. Lastly, all biomechanical-based methods employ the linear elastic model.

The sparse dataset includes three initial transformations for the rigid alignment. Table III and Table V report results using the manually initialized ICP transformation, consistent with the performance shown by Heiselsman et al. [1]. Table IV presents results from our proposed method using all initial transformations, demonstrating its robustness across different starting conditions. The proposed method consistently achieves the lowest mean registration errors across various settings, including the entire sparse dataset (Table III), different visibility conditions (Table III), multiple initial transformations (Table IV), and both noise-free and noise-affected scenarios (Table V), demonstrating its effectiveness and robustness. As shown in Table III, our proposed method outperforms the method developed by Mestdagh and Cotin [16], highlighting its effectiveness in automatically identifying boundary conditions. Compared to other FEM-based methods [1], [17], our approach also achieves lower registration errors (Table III and Table V). Lastly, the learning-based methods [20], [43] show less satisfactory performance than the best biomechanical approaches.

TABLE II COMPARISON OF REGISTRATION ERRORS ON THE OPEN-CAS DATASET. MEAN ± STANDARD DEVIATION (MAX) ERROR IN MILLIMETERS IS PRESENTED. THE LOWEST MEAN VALUES AND THE LOWEST MAXIMUM VALUES ARE HIGHLIGHTED IN BOLD

TABLE III COMPARISON OF REGISTRATION ERRORS ON THE SPARSE DATASET. MEAN ± STANDARD DEVIATION (MEDIAN) OF MEAN ERROR IN MILLIMETERS IS PRESENTED. THE LOWEST MEAN VALUES AND THE LOWEST MEDIAN VALUES ARE HIGHLIGHTED IN BOLD TABLE IV SENSITIVITY ANALYSIS OF THE PROPOSED METHOD TO THE INITIAL ALIGNMENT ANALYSIS PERFORMED ON THE SPARSE DATASET. MEAN ± STANDARD DEVIATION (MEDIAN) OF MEAN ERROR IN MILLIMETERS IS PRESENTED

TABLE V PERFORMANCE OF THE PROPOSED METHOD TO NOISE-FREE AND NOISE-AFFECTED SPARSE DATASETS. MEAN ± STANDARD DEVIATION OF MEAN ERROR IN MILLIMETERS IS PRESENTED

Qualitative results for Set100 in Fig. 9 reveal challenges for GMM-FEM in aligning the right inferior ridges of the liver, where our method excels in accurately fitting the intraoperative point cloud.

The proposed method relies on three parameters: (1) k ss , a stabilizing parameter for the solution; (2) Poisson's ratio ¿, defining the compressibility of the tissue; and (3) iteration number. To understand the sensitivity of our method to these parameters, we vary each parameter individually while keeping the others constant. We examined their impact on registration error and conducted statistical Kolmogorov-Smirnov tests at a significance level of ³ = 0.05, using the entire sparse dataset. The results are shown in Fig. 8.

The registration exhibit a degree of sensitivity to the parameter k ss . The default setting of k ss = 0.01 outperforms other settings ( p < 0.001). The proposed method is not sensitive to Poisson's ratio ¿; hence, no statistical significance is found. The default setting of Poisson's ratio ¿ = 0.49 corresponds to nearly incompressible tissue behavior. While some prior studies [1], [4] adopt ¿ = 0.45, our experiments involving variations from 0.45 to 0.49 show stable performance. Lastly, when we reduce the iteration number to 50, statistical significance is found ( p < 0.001), suggesting the default iteration number 200 could be further tuned to reduce computation.

We first introduce a soft spring modification to the stiffness matrix to eliminate the need for priors. This modification has been demonstrated to substantially reduce registration errors, as evidenced in Fig 8 . Secondly, the method is designed to search for optimal forces autonomously. In the absence of manually identified ZBC and force locations, this strategy is more beneficial than the previous approaches [4] of constraining force locations to predefined areas, such as the posterior surface, shown in Fig. 5 (b). In addition, our proposed optimal step determination, as demonstrated in Fig. 5 (c), eliminates the need for manual tuning of the step size parameter and leads to quicker convergence and lower error.

We have evaluated our method across various datasets, encompassing our phantom datasets, the Open-CAS dataset, and the Sparse dataset. Moreover, in comparisons with two alternative open source methods (V2Snet and GMM-FEM) and closed source methods reported in the Sparse Data Non-rigid Registration Challenge [22], [39], our approach consistently outperforms or shows comparable performance, as elaborated in Sections V-B, V-C, and V-D.

In comparison with open source methods, while GMM-FEM incorporates a FEM into the strain energy, regularizing the deformation field estimated by the CPD registration method, it does not guarantee faithful preservation of the original geometry (i.e., vertex connectivity) of the deformed surface. This becomes apparent in situations with noisy intraoperative point cloud data, as illustrated in Fig. 7. In contrast, our approach, where FEM inherently guides the deformation, ensures the maintenance of volumetric vertex connectivity, preserving the geometry during deformation. Furthermore, as demonstrated in Section V-B, GMM-FEM exhibits sensitivity to the visibility of the intraoperative surface, as highlighted in its original paper [17]. V2Snet utilizes FEM to construct a training dataset [20]. While it shows competitive results in specific instances, its overall performance and generalization capabilities require improvement. Nevertheless, the use of learning-based methods, offering prior information, has the potential to assist optimization methods in achieving faster and more accurate results.

The sparse dataset enables comparison of our proposed method with other closed-source methods, particularly regarding surface visibility, initial alignment, and noise magnitude. Our method demonstrates robustness to these factors within reasonable variations and consistently achieves the lowest registration errors, as shown in Table III, IV, and V.

The proposed method is straightforward and effective, in contrast to registration methods driven by both surface-matching and strain energy terms. Our approach relies solely on the surface-matching term embedded within the FEM, eliminating the need to iteratively determine boundary conditions [1], [4] or manually specify them [16]. Additionally, our method depends on only a few parameters, making it easy to adapt to other datasets. Section V-E provides a sensitivity analysis of these parameters. The stabilizing parameter k ss should be set carefully, though we recommend against extensive tuning, as the default parameter value performs well across the samples included in this study. Potential future work includes automating the process of determining the k ss value.

We have shown that our algorithm's linear elastic model version is parameter robust and leads to accurate registrations for datasets with deformations typically observed in laparoscopic [9] and open liver surgery [23]. The linear elastic model is generally employed by most state-of-the-art methods [1], [16], [39]. In fact, all biomechanical-based methods reported in the Sparse Challenge dataset adopted the linear elastic model. The maximum TRE reported in the sparse data is 11.9 ± 3.7, consistent with clinically observed deformation magnitudes for open liver surgery reported in [23], [40].

Liver deformation in laparoscopic interventions can be larger due to the creation of pneumoperitoneum. The Open-CAS dataset features even more significant deformations than the Sparse Challenge dataset. Even following the initial rigid alignment correction, mean TRE still ranges from 5.88 mm to 13.95 mm, which reflects the typical deformations observed in laparoscopic liver surgeries reported in [9], where surface deformations ranged from 7.7 mm to 12.6 mm at 14 mmHg and from 6.2 mm to 11.9 mm at 7 mmHg. Our proposed registration method achieved mean TRE values below 5 mm across all these samples despite featuring large deformations, with the exception of one sample where the intraoperative surface was highly compromised due to a highly inaccurate stereo reconstruction. Hence, the widely adopted linear elastic biomechanical approach can yield reasonable registration accuracy, even within these typically large deformation ranges. This observation is also shared by previous works [4], [40], [46] that also suggested that a linear elastic model is perhaps adequate to achieve sufficiently accurate non-rigid organ registration.

The liver can also be represented using a nonlinear biomechanical model. While nonlinear models may be considered superior to linear elastic models at capturing large deformations, the material properties required for such models to be sufficiently realistic and accurate vary significantly among patients and are difficult to measure, especially during surgery. Additionally, using a nonlinear model increases computational complexity, as the tangent stiffness matrix must be recalculated with each algorithm iteration. To incorporate a nonlinear model into our method, line 5 of Algorithm 1 can be modified to utilize the FEM for calculating deformations and the tangent stiffness matrix. Various nonlinear models can be selected for this purpose. For further details, we refer the reader to the works described in [3], [15], [46], [47], leaving this aspect for future exploration.

The intraoperative point cloud, collected using an optically tracked stylus, is very sparse compared to other imaging modalities. However, despite this sparsity, the proposed method still achieved accurate results. To address extremely sparse intraoperative point clouds, a grid-fitting interpolation method similar to the one used by Collins et al. [38] could be employed to recover denser surface points. However, it is recommended to avoid such sparse conditions during data collection, and hence, an extremely sparse data analysis would become essentially unrealistic.

Before performing non-rigid registration, it is essential to (1) minimize intraoperative liver surface reconstruction errors, (2) ensure a high-quality surface with sufficient convergence, and (3) achieve accurate rigid alignment. These steps are critical prerequisites; otherwise, our method, like others, may yield inaccurate results. Integrating real-time imaging modalities, such as ultrasound or cone-beam CT, could potentially relax these requirements by providing additional constraints, such as vascular features (i.e., ligaments visible in ultrasound) or a complete intraoperative surface from cone-beam CT.

We have presented a unique non-rigid 3D-3D registration method that seamlessly integrates a biomechanical model into the surface-matching term without relying on priors. By embedding the biomechanical model in the surfacematching term, the estimated deformation is forced to preserve the geometry of the preoperative model without an additional strain energy term. Experiments on our in silico phantom, in vitro phantom datasets, and two publicly available datasets demonstrate the method's efficacy in achieving accurate registration results and preserving original geometry. Furthermore, our constructed phantom dataset and registration method represent significant contributions to liver registration benchmarking and, hence, valuable additions considering the limited availability of open-source datasets and methods.

As part of our future work, we intend to research the effects of nonlinear hyperelastic materials, incorporating learningbased methods, automating the process of determining the stabilizing parameter and the potential combination of this proposed method with other real-time intraoperative imaging modalities. We also note that the 3D-3D non-rigid registration problem tackled here shares many similarities with those used for other organs, such as the prostate, kidney, brain, and breast, and hence, we would be interested in investigating the performance of the proposed method on other organs.