Modeling algorithm for dynamic generalized body cavity

For the purpose of preoperative planning, this paper proposed a modeling algorithm for automatically creating the dynamic model of the generalized body cavity. The algorithm includes two phases. One is the static phase in which the initial cavity model is created from the cross sections; the other is the dynamic phase for modeling the dynamic cavity. In the static phase, the parallel cross sections are generated through slicing the organs enclosing the cavity; then construct the LBCs, and generate the uniform grids; and then create the hexahedron voxels according to the uniform grids; at last, post-processing is used to get the final cavity model represented by triangular meshes. In the dynamic phase, the cavity models at different moments are gained by applying the anchor connection method. We took the human articular cavity as an example to verify the proposed algorithm, and the experimental results showed that the proposed algorithm is effective. Subjects: Algorithms and Complexity; CAD CAE CAM-Computing and Information Technology; Computer Graphics and Visualization; Biomechanics


PUBLIC INTEREST STATEMENT
Minimally invasive surgery (MIS) is very difficult and has high risk because of the small field of view, the narrow cavity space and the way that surgeons perform an operation only through the endoscopic video. With the aid of a computer, visualization and digital analysis of the internal structures of surgical cavity will greatly improve the accuracy and reliability of MIS. Therefore, an algorithm of reconstructing the space model for operation according to patient's data was presented in this paper. Using this algorithm, the dynamic cavity models are reconstructed based on the results of numerical simulation. These models will be used for quantitative analysis to achieve the formulation of surgical program for individualized patient. This work is a part of the project that is designed to setup a systematic methodology for automatic preoperative planning.

Introduction
Nowadays, it is more efficient and more accurate to acquire patients' medical data as a result of the rapid development of medical imaging. In the past two decades, the 3D modeling of organs based on medical image processing has received great attention. The 3D models allow surgeons to gain insight into the anatomical structures and to discern the relative spatial position of the lesion in a more intuitive way. Therefore, they have been widely used in surgical planning, evaluation, virtual surgery system and prosthesis manufacturing.
In the field of surgical planning, however, 3D anatomical models were often limited to the visualization application. Due to the lack of in-depth quantitative analysis to the models, the current interactive software for surgical planning is low intelligence, which leads to a large workload of doctors, and the planning effect is over dependent on the surgeon's clinical experience.
More recently, automatic preoperative planning (De Momi et al., 2013, 2014Schoob et al., 2015) has received attention. However, previous studies focused mainly on specific application, and rare systematic research on automatic preoperative planning approach has been reported. The traditional path planning algorithms are often inefficient because they require intersection testing with the scattered organs. In addition, they are incapable of exploring the characteristics of the body cavity space. Different from the existing works, body cavity is the very environment to be analyzed quantitatively in our work. The body cavity in our work refers not only to the natural body cavity, but also to any anatomic structure which served as a candidate inside which artificial channel is to be created. So we named it generalized body cavity (GBC) (Mo, He, Li, & Wei, 2017).
We propose an algorithm for automatic creation of a GBC in this paper. An overview of the main features of the algorithm is as follows: (1) It can automatically create a static model of the body cavity which is enclosed by multiple organs known as enclosing body (EB).
(2) By simulating the motion of the EBs, it can automatically generate dynamic models of the body cavity at different moments.

Related work
This work focuses on the modeling algorithm for human body cavity, so we restricted the review of literatures in two aspects-the modeling of cavity and the 3D reconstruction algorithm from 2D crossing contours.
The previous studies on cavity modeling were common in computer assisted surgery and geological engineering. In the field of digital medicine, the 3D reconstruction methods for the human organs with cavities, such as blood vessels (Schumann, Neugebauer, Bade, Preim, & Peitgen, 2008), intestinal tract (Hong, Tavanapong, Wong, Oh, & De Groen, 2009), digestive tract (Feulner et al., 2009) and ear cavity (Rodt et al., 2002), have been widely studied. Similar to the solid organs, medical images (e.g., CT, MRI and ultrasound, etc.) were used to create the surface models or solid models of this kind of anatomical structures. However, they are different from the cavity model defined in this paper because they have wall thickness. To the best of our knowledge, the modeling algorithm for human body cavity, especially for the GBC, has not been reported. In the field of geological engineering, there were some researches on the modeling of the tunnel and the mined-out cavities. For example, Luo, Liu, Zhang, Lu, and Li (2008) created the 3D models of the mined-out cavities of copper mine by utilizing the laser scanning data. However, their works are different from our algorithm not only in the input data, but also the accuracy requirement and modeling method.
The 3D Boolean operation (Jiang et al., 2016) seems to be suitable for the cavity modeling, as it can create complex models efficiently from the existing models. Although the existing geometric models served as the input data in our algorithm, the cavity model cannot be achieved by Boolean operations. The reason lies in that the model obtained by Boolean operations may contain a part or the whole volume of the input models, while the body cavity in our application does not contain any solid objects.
Over the past three decades, much research (Holloway, Grimm, & Ju, 2016) has been done on the reconstruction of 3D objects from a finite set of planar cross sections. These studies have experienced a process from simple to complex, whether in the input data or in the complexity of objects to be rebuilt. The earlier studies were mainly focused on the reconstruction of 3D objects from parallel contours. Recently, the reconstruction of 3D objects from arbitrary cross sections has been paid more attention. Boissonnat and Memari (2007) extended the original Delaunay meshing method from parallel cross sections to address the issue of shape reconstruction from unorganized cross sections. Following the divide and conquer strategy of Boissonnat and Memari (2007), Liu, Bajaj, Deasy, Low, and Ju (2008) present an algorithm to handle curve networks of arbitrary shape and topology on cross section planes with arbitrary orientations. Bermano, Vaxman, and Gotsman (2011) proposed an algorithm to reconstruct the surface of smooth three dimensional multilabeled objects from sampled planar cross sections of arbitrary orientation. More recently, to achieve the correct topology of the reconstruction, Edwards and Bajaj (2011) present a method for reconstruction of spatially realistic and topologically correct models from planar cross sections of multiple objects. Zou, Holloway, Carr, and Ju (2015) proposed an algorithm that provides topological control during surface reconstruction from an input set of planar cross sections. Sharma and Agarwal (2017) utilized the Hermite mean value interpolation for triangular meshes based on split-andmerge method to surface reconstruction from unorganized planar cross sections.
A variety of methods (Bermano et al., 2011;McCrae, Singh, & Mitra, 2011) have been developed for surface reconstruction from cross sections, and they fall into two main categories, the contour stitching method (Zou et al., 2015) and the volumetric method (Curless & Levoy, 1996;Nießner, Zollhofer, Izadi, & Stamminger, 2013). The contour stitching method needs to solve the problem of contour correspondence, tiling, branching and fitting are needed to handle. Although this method has been well studied, there are still challenges to deal with the branching topology and the object with rich geometric features. The volumetric method usually needs to solve the isosurface which determined by the signed distances (Jones, Baerentzen, & Sramek, 2006;Sharma & Agarwal, 2017;Tang & Feng, 2018;Zollhöfer et al., 2015). The difference of the previous algorithms based on the volumetric method is mainly lies in the definition of the distance field (Gomes, Bellon, & Silva, 2014). Compared to the contour stitching method, the volumetric method is lower sensitive to the noise data, and no special consideration of the branching topology problem. However, the accuracy of the reconstruction result depends on the voxel resolution. Smaller voxels often lead to higher accuracy, but it will increase the computational cost and memory. In addition, the solution of the distances is time consuming.

Problem description and the algorithm flow
The problem of creating the model of GBC from a set of enclosing bodies can be formally described as follows.
Suppose that a cavity is surrounded by n organs (enclosing bodies) whose shapes may be changed dynamically. The geometric models of these EBs at t moment are denoted as M t k È É , M t k : ¼ ðV t k ; T t k Þ, k ¼ 1; 2; :::; n, where V t k is the vertex set, T t k is the triangle set of the k-th EB which represented as triangular mesh. At t moment, the surface model of GBC is defined as M t c : is composed of two parts. The first part inherits from the EB surfaces, and it is denoted as [ n k¼1 Γ t k , where Γ t k is the partial or complete boundary surface of the k-th EB; and the second part is denoted as [ 1 i;j , where 1 i;j is the patches suturing two adjacent surfaces Γ t i and Γ t j .
Construction of the surface patches Γ t k and 1 i;j is the key for the modeling of the GBC; however, it is a challenging issue to construct them directly in 3D space. To simplify the problem, we proposed a two-phase algorithm. Although the algorithm converts the 3D models to a series of parallel cross sections first, an alternative method is employed to enhance the robustness, instead of applying the contour stitching method.
The flow of the two-phase algorithm is detailed in Figure 1. As the input data of the algorithm, all the initial models M t 0 k ; k ¼ 1,n n o of EBs are reconstructed from medical images. In the static phase, the initial model M t 0 c of the GBC is extracted though five steps: (1) The parallel cross sections are generated through slicing the EBs; (2) Construct the layered bounding contours (LBCs) which can determine the space extent of the cavity model at each layer; (3) By performing the uniform rectangular grids generation, the space of each layer will be divided into boundary-grids, outside-grids and cavitygrids; (4) Create voxel model of the body cavity according to the grid type of two adjacent layers; (5) Post-processing is to be performed to get the triangular mesh model of body cavity.
In the dynamic phase, the cavity models at different moment are obtained by applying the anchor connection method (ACM).

Generation of cross sections
Cross sections are generated by cutting every enclosing bodies with a series of parallel planes which perpendicular to a given vector d. The crossing sections of the k-th EB form a set ; i ¼ 1; 2; :::; n, where n is the number of layers, L i k is the crossing sections by slicing the EB M t 0 k on the i-th layer. L i k consists of a certain number of closed contours, that is, where RðL i k Þ ! 0 is the cardinal number which represents the number of closed contours.

Construct the LBCs
The approach of constructing the LBCs is contours merging. At first, two LBCs are merged, and then add one LBC at a time until only one LBC left in a certain layer. The general strategy of the algorithm is to merge the contours of the same EB first, and then merge the contours between different enclosing bodies. When merging two section contours, connection is required because there is usually a gap between the enclosing bodies, and their section contours do not intersect at all. However, it is a challenging issue how to connect two contours since the actual connection points are not determined. In general, it depends on the application context to set the connection conditions, and to determine the connection scheme. We introduce two merging schemes, as shown in Figure 2(a,b), respectively. The first scheme connects two contours with B-spline curve passing through three points. As shown in Figure 2(a), the points e 1 and e 2 are the centers of the cross section Ω 1 and Ω 2 , respectively. Starting from e 1 , draw two lines l 11 and l 12 tangent to Ω 2 , and from e 2 draw the other two lines l 21 and l 22 tangent to Ω 1 , then construct two B-spline curves passing through one intersection point (I 1 orI 2 ) and two tangent points at same side, such as T 1 and T 2 . As the result of merging, only one contour Ω 3 is left at one layer. The second scheme also uses B-spline curves to connect cross sections. However, the points controlling B-spline curve are different from the first scheme. As shown in Figure 2(b), the lines l 1 , l 2 and l 3 are three control lines which include the three control points P 1 , P 2 and P 3 ; the points I 1 and I 2 are intersection points between the control lines. Construct two B-spline curves passing through the point set P 1 ; I 1 ; P 3 f gand P 2 ; I 2 ; P 3 f g , respectively. The control points and control lines depend on the application context, for example, in the case of section 4, the distribution ranges of the articular cartilage were used as the constraint condition.

Rectangular grids generation and voxels creation
The generation of the layered rectangular grids is the precondition for constructing the cavity model. To align the grids on the adjacent two layers, set the same size of the bounding rectangles (as shown in Figure 3) at each layer and align their centers. Each layer is divided into uniform rectangular grids with r rows and c columns. These grids will be marked with different types-outside-grid (OG), cavity-grid (CG) or boundary-grid (BG). As shown in Figure 3, the OG defines the space between the bounding rectangle and the layered bounding contour; the CG defines the space of the GBC; the BG is passed through by the layered bounding contour or by the cross sections of the EBs. The BGs are divided into three subcategories which are denoted as G 1 b , G 2 b and G 3 b (as shown in Table 1). The category G 1 b is passed through only by the cross section, and G 2 b only by the layered bounding contour, while the G 3 b passed through by the cross section and layered bounding contour as well. All the grid categories and their symbols are listed in Table 1.
Multisource region growing algorithm is applied to label the grid type. The labeling order is boundary-grids, outside-grids and cavity-grids. The boundary-grids are served as the sources for labeling cavity-grids, and the outside-grids provide cues to indicate the direction for region growing.
First, label the boundary-grids by intersection test between the cross section Å and the rectangular grid g i;j . If Å intersects with g i;j , then g i;j is classified as boundary-grid. Then label the outside-grids by region growing which restricted by the labeled boundary-grids. Select the grids locating at the first row, the last row, the first column and the last column as seeds, besides those grids labeled as boundary-grid. Starting from an unlabeled seed, recursively access and mark their unlabeled neighbors. Repeat the above process until there is no unlabeled seed.
At last, region growing is still used to label the cavity-grids. Select the neighbor grids of the boundary-grid G 2 b as candidate seeds, excluding those grids belong to the type G o , G 1 b or G 3 b .
After labeling all the cavity-grids by repeating the above steps for each layer, start from the first layer to build the voxel model of the body cavity. If both the grid m i j;k (i, j and k is the layer number, row number and column number, respectively) and m iþ1 j;k belong to cavity-grid, create a hexahedron voxel for this two grids. After this step, the cavity model presented by voxels is achieved.

Post-processing
The post-processing is aimed to remove small isolated blocks, generate the triangular mesh model and perform mesh optimization.
The algorithm for removing small isolated blocks is detailed as following: Step 1. Select a voxel that has not been grouped as a seed, mark it with the current group number, and then put it into the queue Q; The boundary-grid of layered bounding contour The hybrid boundary-grid G 3 b Step 2. If Q is empty, turn to Step 3; otherwise, delete the first element Q 0 from Q, pick out its neighbors which have not been grouped, and mark them with the current group number, and then put it into Q; Step 3. If all voxels have been grouped, stop grouping and turn to Step 4; else, update the current group number, and then turn to Step1; Step 4. Among all the groups, remove those voxels in the groups in which the number of voxels is less than a given threshold. Finally, generate the triangular mesh model and perform mesh optimization by applying the methods presented in (Yoo, 2011).

The modeling method of dynamic cavity
The shape variation of the body cavity is caused by the dynamic changes in shape or position of the EBs. Three factors may cause these dynamic changes of the EBs-the relative motion between rigid bodies, the deformation of flexible bodies and tissue resection. The relative motion of rigid bodies may transform the positions of the vertexes on the object surface, but it does not affect their topologies. As for the flexible bodies, e.g., the soft tissues, they may deform under the external forces, such as stretching, extrusion and shearing.
Updating the models of EBs is the basis for the construction of dynamic cavity models. According to different material characteristics of the organs and the force acting on them, the geometric model set M t k ; k ¼ 1; :::; n È É at t moment can be updated. For the rigid bodies, the displacements of vertexes can be obtained by solving kinetic equations, and carrying out translation or rotation transformation. For the flexible bodies, the finite element method is utilized to calculate their deformations.
After updating the model set, the cavity model at t moment is to be created automatically by applying the algorithm given in section 3. However, repeated contours merging, rectangular grids generation, voxel creation and post-processing at every time step will result in low efficiency. We present an alternative approach-the ACM-to automatically create the cavity model at each moment. Because there exist common vertices between the body cavity and its enclosing bodies, the corresponding points on the cavity model have the same displacements with their counterparts. The ACM updates the cavity model M t c with the corresponding relationship between the vertices. We defined the anchor pair at t moment as pðv t where v t j is the closest point to v t i . Then, we can update the vertices on the surface of cavity model at t þ 1 moment by It is not necessary to update all the anchor pairs at each time step. If large displacements occurred on the local surface of the enclosing bodies, it needs to compute the anchor pairs at the next moment to avoid large errors. Updating completely the anchor pairs is time consuming, hence local adaptive subdivision is employed to improve computational efficiency. Subdivide those triangular meshes which contain vertices with large displacement, and then only calculate the anchor pairs of the new added vertices, which can improve accuracy without significantly increasing the computational time.
In some extreme cases, for example, expanding the cavity space by tissue resection, directly solving the anchor pairs may cause error matching due to the large gap between the cavity surface and the enclosing bodies. In this case, union operation is to be performed to add the expanding space to the current cavity model, before recalculating the anchor pairs.

Experimental results and discussion
We took the human articular cavity as an example to verify the proposed algorithm.
As shown in Figure 4, the articular cavity of the knee joint is enclosed by femur, tibia plateau, ligaments, meniscus and muscles. It is very difficult for surgical instruments to reach the lesion due to the knee joint cavity is very narrow and complex. During knee arthroscopic surgery, surgeons need to constantly adjust the patient's lower limb to obtain necessary operating space. Therefore, only the experienced doctor can be qualified for this job.
The shape of the knee joint cavity can be changed by adjusting the degree of flexion of the tibia. In order to make surgeon understand thoroughly the variation of the cavity shape under different flexion angles, we created the dynamic models for the knee cavity by applying the proposed algorithm. The cavity model will be further applied to the automatic surgical planning and guidance for designing surgical instruments.

Create the initial cavity model
We reconstructed the models of the enclosing bodies (e.g., femur, tibia, cartilage, ligaments and menisci, etc.) basing on the medical images. Served as the input data of the proposed algorithm, these models were used to create the initial cavity model M t 0 c in the static phase of our algorithm. The results of the main steps in the static phase are displayed in Figures 5-8. As shown in Figure 6(b), when constructing the LBCs, the distribution ranges of the articular cartilage were used to position the control points and control lines. If the most common anterior portal is intend to be adopted in the knee arthroscopy surgery, surgical instruments  are sometimes allowed to access the cavity through the infrapatellar fat pad (as shown in Figure 6(b)), so we integrated it into the cavity model.

Create the dynamic models of the knee joint cavity
To obtain the accurate dynamic models of the knee joint cavity, it is necessary to acquire the patient-individual data under different flexion angles of the tibia before creating the dynamic cavity models. However, by applying this method, there are several disadvantages including high cost, exposing patients to excessive radiation and requiring designing a special auxiliary device to accurately adjust the posture of patient's lower limb. Therefore, we proposed an alternative approach that requires the acquisition of patient data just in a fixed pose, while the anatomical models under different poses will be updated through dynamics simulation.
The locomotion of the knee joint is very complicated, since it contains both the motion of rigid bodies and the soft tissue deformation. In general, it can be abstracted as the sixdegrees-of-freedom motions of the tibio-femoral joint (TFJ) and patello-femoral joint (PFJ). As shown in Figure 9(a), there are three rotations-the flexion-extension, abduction-adduction, internal-external rotation, and three transitions-anterior-posterior, medial-lateral and proximal-distal transition. In order to establish the kinetic equations, the local coordinate systems (as shown in Figure 9) of the femur and tibia were setup. The motion of the TFJ and PFJ are mainly affected by the inertia force, ligament force, quadriceps tendon force and contact force. Friction forces were neglected because of the extremely low coefficient of friction of the articular surfaces. In this paper, the kinetic model in Caruntu & Hefzy (2004) was used to solve the internal-external rotation angles (α 2 ), and the abduction-adduction rotation angles (α 3 ) under different tibial flexion angles (α 1 ). We listed a subset of data in Table 2.
Once the angles α 2 and α 3 were obtained for a given flexion angle α 1 , the position vector of any point on the tibia can be calculated by the following transformation (2) where P o is the position vector of the origin of the tibia coordinate system relative to the femoral coordinate system; r denotes the local coordinate of the tibia; and the rotational matrix R describing the orientation of the tibia relative to the femoral coordinate system is given as R ¼ s 2 c 3 s 2 s 3 c 2 Àc 1 s 3 À s 1 c 2 c 3 c 1 c 3 À s 1 c 2 s 3 s 1 s 2 s 1 s 3 À c 1 c 2 c 3 Às 1 c 3 À c 1 c 2 s 3 c 1 s 2 2 4 3 5

Flexion-Extension
(3) Among the organs enclosing the knee cavity, the tibia and the femur are regarded as rigid bodies. The femur remains static, while the tibia is allowed to the six-degrees-of-freedom movement. The position vectors of the vertices on the surface of tibia were evaluated by the Equation (2). As to the soft tissues, for example, the anterior cruciate ligament, posterior cruciate ligament and meniscus, the forces (such as the quadriceps tendon force, contact forces and ligamentous forces) acting on them were computed by applying the method in Caruntu & Hefzy (2004). Their deformations were computed by the finite element method. In the implementation phase, we calculated the deformation by using the finite element module of SOFA [10][11] which is an open source framework. Thus, we can update all the models of the enclosing bodies at different poses. If the pose p is regarded as a function of time, that is p ¼ f ðtÞ, then we get the dynamic models at different moments. Figure 10 displays the variation of the body cavity at different poses. When the flexion angle of the tibia is smaller, the upper surface of the tibial condyle and the lower surface of the femoral condyle are almost parallel, and the volume of the anterior part of the cavity between the tibia and the femur is almost the same as the volume of the posterior part. When the flexion angle increases, the contact point of the TFJ moves backward, the volume of the lower part of the anterior cavity increases, while the volume of the posterior part decreases, until the space of the anterior cavity and the posterior part is almost completely blocked (Figure 10(f)). Comparing with the experimental studies (Tranberg, Saari, Zügner, & Kärrholm, 2011), the shape variation of the cavity model is consistent with the measured results by qualitative analysis.

Discussion
We quantitatively analyzed the accuracy only for the static cavity model due to the limitation of the experimental conditions. As for the ideal static cavity model, its vertices on the surface should be located on the surface of a certain EB. However, the voxel resolution and post-processing may affect the accuracy of the initial cavity model. To quantitatively calculate the accuracy, we defined the mean distance deviation (MDD) of the anchor pairs, and take it as the evaluation index of the precision. The MDD is computed by the equation " δ ¼ ∑ n k¼1 d k =n, where n is the total number of the anchor pairs; d k is the distance between vertex v t 0 i and v t 0 j , and they constitute an anchor pair pðv t 0 i ; v t 0 j Þ. According to this definition, the values of MDD are 0.201, 0.137 and 0.085 mm with respect to the grid resolution 200 × 300 × 300, 500 × 300 × 300 and 500 × 500 × 500, respectively. This suggested that the accuracy would be improved as the resolution increases.
Deviations exist in the models of the enclosing bodies updated through simulation, so they should not be used as the basis to evaluate the precision of dynamic cavity. To accurately evaluate the accuracy of the dynamic cavity, the key lies in collecting the patient's medical image at different poses. The special device which can be able to accurately adjust the flexion angle of the tibia is needed to be designed in the future.
The time complexities of the main tasks of our algorithm are listed in Table 3. Because the updating of the models of enclosing bodies is achieved by dynamic simulation, the time complexity of this task depends on the kinetic model used and the numerical method. Therefore, the time complexity of our algorithm without consideration of the task of updating the models of enclosing bodies is Oðn t þ n v þ n vox Þ. Among them, n t is the total triangles of the EBs, n v is the number of vertexes in GBC, and n vox is the number of voxels in GBC.
Compared with the existing algorithms, our algorithm has the following characteristics because of the specific application context.
(1) The input data of our algorithm is quite different. As shown in Figure 3 and Figure 6(b), the LBCs are obtained by slicing multiple organs. These contours are not single and continuous closed curves, and they include not only the target space (the space of GBC), but also the space outside of the GBC. Therefore, it is hard to form the 2D cross sections utilized in the contour stitching algorithms or the volumetric algorithms.
(2) The purpose of modeling is different. For the purpose of preoperative planning, the surface model and voxel model of GBC are needed to reconstruct simultaneously. These models are served as the basis of topological analysis (Mo et al., 2017) and the quantitative analysis of spatial accessibility. Although post-processing was implemented in our reconstruction process, the smooth surface model is not necessary.
The above two characteristics determine that our application cannot be solved by simply employing the existing algorithms. And the different purpose of modelling makes our focus different.

Conclusion
An algorithm for automatically creating the dynamic body cavity with multi enclosing bodies is presented in this paper. In view of the difficulty of directly solving the boundary surface of the body cavity in 3D space, we converted the 3D models to a series of parallel cross sections. To improve the efficiency of creating dynamic model, the cavity model at different moment is updated with the relationship between the two vertices of the anchor pair. Finally, we took the knee joint cavity as an example to verify the proposed algorithm. The results of qualitative analysis showed that the dynamic cavity model is consistent with the law of the motion with respect to the knee joint. In addition, for the initial cavity model, the quantitative analysis indicated that the high precision model can be created by our algorithm, and the precision is related to the voxel resolution.
Several aspects of the proposed approach can be improved in future research: (1) Due to the limitations of the current experimental conditions, there is no quantitative analysis for the accuracy of the dynamic model. It is needed to design a special experimental device to verify the dynamic cavity models of the human knee joint.
(2) The vascular and neural tissues have not been taken into account in our experiments; the effects of these models need to be considered in future studies.
(3) To verify the algorithm in more fields, such as creating the abdominal cavity for the laparoscopic MIS.