Cutting Simulation in Unity 3D Using Position Based Dynamics with Various Reﬁnement Levels

: Augmented and Virtual Reality-based surgical simulations have become some of the fastest-developing areas, due to the recent technological advances and changes, in surgical education. Cutting simulation is a crucial part of the virtual surgery simulation in which an incision operation is performed. It is a complex process that includes three main tasks: soft body simulation, collision detection and handling, and topological deformation of the soft body. In this paper, considering the content developer’s convenience, the deformable object simulation, using position-based dynamics (PBD), was applied in the Unity 3D environment. The proposed algorithm for fast collision detection and handling between the cutting tool and the deformable object uses a sweep surface. In case of incision, the algorithm updates the mesh topology by deleting intersected triangles, re-triangulation, and reﬁnement. In the reﬁnement part, the boundary edges threshold was used to match the resolution of new triangles to the existing mesh triangles. Additionally, current research is focused on triangle surface meshes, which help to reduce the computational costs of the topology modiﬁcations. It was found that the algorithm can successfully handle arbitrary cuts, keeping the framerate within interactive and, in some cases, in the real-time.


Introduction
Virtual simulation is a reproduction of a real-world environment obtained with the help of a computer system, which allows user to interact with the simulated objects. Nowadays, due to the progress in computational science, virtual simulation became a rapidly growing area. Virtual simulation is used in many different fields, such as engineering, education, natural sciences, medicine, and so on. Many scientific papers were dedicated to study how the virtual simulation impacts education and teaching. Researchers in [1] provide information on a web-based simulation platform named EDISON, which can be used to simulate different physical phenomena. Ref. [2] proposes a Mixed Reality Chemistry Lab to simulate chemistry experiments and help students to obtain needed knowledge before starting to work with actual chemicals. The research in [3] summarizes the most important features of virtual simulations in education, such as visualization, inclusiveness, unlimited access to information, and increased engagement. Ref. [4] examines the application of virtual simulations in higher education. The study concludes that VR is a promising technology, and it is often used to obtain skills, which require declarative knowledge and procedural-practical knowledge. Nonetheless, there are proofs that AR and VR technologies are highly suitable in medicine, and there are many different applications of it in treatment and rehabilitation. Research in [5] effectively uses VR simulation in rehabilitation purposes to cure vertigo conditions. Ref. [6] proposes VR-based exposure therapy to treat mental illness. The research shows promising results and proves that blow-ups. The next physical model is FEM, which discretizes objects into smaller elements and applies governing differential equations to them [17]. Due to its precise simulation results, FEM is often used in physical phenomena analysis and engineering. However, the results of a FEM-based simulation are highly dependent on the chosen mathematical model, and to achieve higher accuracy the approach may require more refined mesh with bigger number of elements. Therefore, in many cases, because of high computational costs, FEM produces simulations in interactive framerate or offline mode [14]. The last numerical discretization model is PBD, which has the main feature of directly manipulating positions during the simulation [18]. The positions are controlled by applying restraining forces in the form of constraints. The most generally used constraints types are stretching and bending ones. The further discussion on PBD-based simulations will proceed in the Methodology part of the current article.
Active research on the cutting simulation is done based on the abovementioned geometry representation and numerical discretization models. The research in [19] proposed a tetrahedral volume mesh-based interactive cutting simulation, which exploits the look-up table in topology update and avoids cracks during tetrahedron subdivision. Additionally, the researchers used the axis-aligned bounding box technique to enhance collision detection between the cutting tool and the mesh. In [20], the researchers simulated an esophagus by using surface triangular mesh and the boundary element method. Here, a new approach was proposed in which binary tree is used to create a hierarchical database of triangles and trace each triangle subdivision. Ref. [21] designed an extended PBD(XPBD)-based hybrid cutting simulation. In the research, the refined surface triangular mesh is combined with coarser tetrahedral interior mesh to improve the performance of the simulation. Ref. [22] proposes an interesting solution to simulate the cutting and tearing of soft bodies in AR. The authors utilize a spatio-temporal registration technique to capture the deformation of tracked objects and apply it to the computational FEM model. In [23], the researchers proposed an adaptive octree mesh model with embedded structures to simulate the cutting of human organs with nerves and vessels inside them. Ref. [24] used PBD to execute a simulation of electrosurgery on a mesh-free model. The researchers proposed a heat-response threshold to make the cutting process more plausible. Additionally, the framework was able to simulate several human organs in the same scene in real time. In [25], a new hybrid method is proposed, which combines surface mesh and internal meshless point elements through the usage of virtual points. The algorithm simulates cutting in two steps: first, constructing the Bezier curve according to surface collision detection; second, calculating the deformation displacement for the point elements. Ref. [26] proposes a Mass Spring-based three-stage cutting simulation, which includes creation of a swept volume, composition of Bezier curve and retriangulation by using shortest distance node matching algorithm. Additionally, the proposed method calculates residual stress, which helps to simulate tissue shrinking during the cutting.
The current study is focused on the designing of a new method to simulate cutting during virtual surgery simulation basing on PBD. The numerical discretization model was chosen considering instability of MSM and high computational costs of FEM approaches. Additionally, compared to the simulation approaches in [21], the current model is applied only to the surface triangle mesh of the simulated soft body. We intentionally avoid using tetrahedral mesh or other volumetric geometry representation in order to reduce computational costs and achieve cutting simulation in real-time. Furthermore, current research proposes a method to patch and refine large holes, which may appear during cutting. Here, we apply the algorithm from [27], which achieves hole patching using the ear-clipping approach. The refinement level may be set to original surface mesh, with coarser or more refine triangulation depending on user input.
The remainder of this paper is organized as follows. Section 2 explains the overall methodology of the proposed cutting simulation. Section 3 shows the cutting experiments. Section 4 contains the performance evaluation data, and the last section summarizes the main findings of the current research, its limitations, and future works.

Methodology
The main parts of the proposed cutting simulation method include composition of PBD model, intersection check, hole patching, and refinement. These steps are described in the following subsections.

Composition of PBD Model
In current research, PBD is used to physically simulate soft bodies during virtual cutting. Comparing it to many classical force-based approaches PBD influences positions directly, which helps to avoid overshooting and energy loss [18]. Hence, this method can achieve fast, stable, and rather accurate simulation both of rigid and soft bodies. A PBD model consists of N number of particles and M constraints. By tuning the stiffness value k of each constraint, we can adjust elasticity of simulated bodies and achieve different visual effects. The algorithm of PBD is described in Algorithm 1. Here, after receiving the initial positions x, velocities v, and inverted masses w of each vertex, the system solves Euler integration in lines 5-8 to obtain temporal positions p. After that, in line 9, the algorithm iteratively resolves all constraints in a Gauss Seidel-like manner to adjust the predicted positions p. Finally, in lines 10-14, velocity v and position x of all vertices are set based on corrected temporal positions p. To achieve the plausible simulation of soft bodies, two types of constraints were applied: stretching and bending. More information about these constraints can be found in [18]. for all vertices i do 6: v i ← v i + w i f e (x i )∆t 7: end for 9: for n times solveConstraints() 10: for all vertices i do 11: v i ← (p i − x i )/ ∆t 12: x i ← p i 13: end for 14: end while

Cutting Simulation
The proposed cutting simulation method is shown in Algorithm 2. The cutting starts with the intersection check between mesh triangles and a cutting tool in line 5. Lines 6-12 describe the procedures taken in case of intersection. Firstly, we create a set of intersected triangles and a set of intersected edges. Then, the intersection points are calculated, and intersected triangles are removed from the initial surface mesh. After that, in line 14, to maintain the physical part of the simulation accurately, the algorithm removes constraints corresponding to the edges of deleted triangles. More details on constraint deletion can be found in Appendix A. Then, in lines 15-16, we define the cutting plane by calculating plane's normal n cp . The cutting plane is required to divide all intersected edges into two groups: ones that lie on the positive side of the plane, and ones on the negative side (lines [17][18][19]. The deletion of intersected triangles produces a hole in the mesh with uneven boundaries. Therefore, the next step in the cutting algorithm is to reconstruct the boundary edges of the hole. For the beginning, we delete intersection points which lie too close to each other in lines 20-25. This helps to avoid degenerate triangles in further steps. After adding intersection points to the vertices list, the algorithm creates new boundary for positive side edges in line 27. Then basing on the updated boundary edges, the hole is patched with initial triangles. The triangulation technique is described in Appendix A. However, to match the original triangles' resolution, we propose to further refine triangles in line 29. Additionally, we provide options to change the refinement level from coarse to fine based on user input. Finally, the algorithm adds new triangles from line 29 to the surface mesh. Then, lines 28-30 are repeated to triangulate negative side edges too. if(intersection occurs) 7: add t i to the intersectedTrianglesList 8: add edges e 1 , e 2 , e 3 of triangle t i to the intersectedEdgesList 9: find intersection points pint 10: remove In current research, we use sweep plane to track the cutting occurrences. The sweep plane is constructed as pictured in Figure 1. The dissecting part of the cutting tool is represented as cutting edge, which is shown in red in Figure 1b. By tracing the position of cutting edge in previous and current frames, we can find the sweep area and define the final sweep plane. The sweep plane is represented by two triangles, as in Figure 1c, which are used to detect triangle-edge intersections with element edges of surface mesh.
As it was mentioned above, after deleting the intersected triangles, there is a hole left in the surface mesh, which requires further patching. However, the boundary edges form an uneven shape, which can cause crenelations during triangulation and unplausible cutting effects. Therefore, in this study we propose a method described in Algorithm 3 to reconstruct boundaries to achieve smooth and fair cutting. The method consists of two steps: first, connect current boundary edges with the closest intersection points ( Figure 2c); second, triangulate edges, which form an acute angle (Figure 2d). In the first stage, the algorithm receives one of the edge sets defined in Algorithm 2 lines 17-19, and for each edge e i , it calculates its center point c, and distances from each edge e i vertex to the cutting plane (lines 6-8). Next, in line 9, we test whether distance from edge vertices to the plane is more than threshold β. On the occasion of opposite current edge is skipped, assuming that it is close enough to the cutting plane and will not provide distortion during triangulation. Whenever all vertices of edge e i lie far from the plane, the algorithm iterates through intersection points in lines 10-13 to find the closest point pint n . After that, a new triangle is created by triangulating vertices of edge e i and closest point pint n . The new edges are created and added into boundary edges list and current edge e i is deleted (lines [15][16][17]. reconstruct boundaries to achieve smooth and fair cutting. The method consists of two steps: first, connect current boundary edges with the closest intersection points ( Figure  2c); second, triangulate edges, which form an acute angle (Figure 2d). In the first stage, the algorithm receives one of the edge sets defined in Algorithm 2 lines 17-19, and for each edge ei, it calculates its center point c, and distances from each edge ei vertex to the cutting plane (lines 6-8). Next, in line 9, we test whether distance from edge vertices to the plane is more than threshold β. On the occasion of opposite current edge is skipped, assuming that it is close enough to the cutting plane and will not provide distortion during triangulation. Whenever all vertices of edge ei lie far from the plane, the algorithm iterates through intersection points in lines 10-13 to find the closest point pintn. After that, a new triangle is created by triangulating vertices of edge ei and closest point pintn. The new edges are created and added into boundary edges list and current edge ei is deleted (lines 15-17).  reconstruct boundaries to achieve smooth and fair cutting. The method consists of two steps: first, connect current boundary edges with the closest intersection points ( Figure  2c); second, triangulate edges, which form an acute angle ( Figure 2d). In the first stage, the algorithm receives one of the edge sets defined in Algorithm 2 lines 17-19, and for each edge ei, it calculates its center point c, and distances from each edge ei vertex to the cutting plane (lines 6-8). Next, in line 9, we test whether distance from edge vertices to the plane is more than threshold β. On the occasion of opposite current edge is skipped, assuming that it is close enough to the cutting plane and will not provide distortion during triangulation. Whenever all vertices of edge ei lie far from the plane, the algorithm iterates through intersection points in lines 10-13 to find the closest point pintn. After that, a new triangle is created by triangulating vertices of edge ei and closest point pintn. The new edges are created and added into boundary edges list and current edge ei is deleted (lines 15-17).  Yet, after executing lines 5-19 there might be duplicated edges in current boundary edge list. Hence, in lines 20-24 we remove all duplicates to proceed to the second part of the Algorithm 3. The edges should also be sorted for further hole patching as shown in lines 25-31. Finally, the algorithm iterates through a list of sorted edges and calculates the angle between current edge es i and next edge es i+1 in line 34. In case of positive acute angle, we create a new triangle from current edge es i vertices and end vertex of the next edge es i+1 as in line 36. A new edge e new is formed from start vertex of edge es i , and end vertex of edge es i+1 . e new is further added to the boundary. In case of insufficient angle, no new triangle is created, and the current edge is added to the boundary edges list in line 38. In the end of Algorithm 3, a complete even boundary is created, as it can be seen in Figure 2d.
In current research, we propose to refine initial patch triangulation as it is described in Algorithm 4. Here, we use dictionary structure edgeDictionary, where edges are keys and triangles' indices are values. The dictionary contains only interior edges of patch triangulation and adjacent triangles' information, which is used during the swapping procedure. Figure 3 shows the main steps of the refinement approach: firstly, it receives initial patch triangulation as in Figure 3b and, then, the method subdivides triangles, as in Figure 3c, and swaps edges to obtain final triangulation, as presented in Figure 3d. However, before starting subdivision, the pre-refinement edge swapping is performed in lines 6-7. The experimental findings showed that this technique helps to achieve more accurate and balanced triangulation. Figure 3e depicts the triangulation without pre-refinement edge swapping. It is clear to see that, without additional edge swapping, the triangulation is not fair and even, and there is a big difference between the smallest and the largest triangles. calculate distance dist1 from startPointe i to the plane cp 8: calculate distance dist2 from endPointe i to the plane cp 9: if(dist1 > β and dist2 > β) 10: for all intersection points pint j do 11: After pre-refinement edge swapping, the algorithm proceeds with triangles subdivision in lines 9-25. It calculates average edge length curLength for triangle tp i in patch triangulation. In case of the curLength is bigger than threshold α, the triangle tp i will be subdivided. Here 1-to-3 triangle subdivision is executed through division of triangle tp i into 3 new triangles. For that the center point vc is calculated in line 14. In case of subdivision triangle information in edgeDictionary must be updated the corresponding edges.
The edges of new triangles should also be entered into edgeDictionary as in lines 21-23. Finally, the algorithm iterates over all interior edges and swaps them one by one to achieve Delaunay-like fair triangulation. The detailed explanation on swapping algorithm can be found in Appendix A. As it can be seen in line 7, the subdivision and swapping are repeated for n times. Experiments showed that different models and refinement resolutions may require more or less iterations to achieve better triangulation results. In Algorithm 4, we introduced threshold α to match triangle resolution to the original one of the surface mesh. To find threshold α, firstly, the average boundary edge length t should be calculated as in Formula (1). Here, be i is a boundary edge. Next, α is obtained through the multiplication of average boundary edge length t with a value z, as in Formula (2). The value z is a user-defined input, which helps to manipulate the refinement level. In current research, we set five such levels: original, 2 times coarser, 1.5 times coarse, 2 times more refine, 1.5 time more refine. The examples of each refinement resolution are shown in Figure 4. for (i = 0; i < patch triangles number − 2; i++) 9: let edge 1 , edge 2 , edge 3 be the vertices of triangle tp i 10: find average length of edges of triangle tp i : curLength ← (edge 1 + edge 2 + edge 3 )/3 11: if(curLength >α) 12: let the algorithm iterates over all interior edges and swaps them one by one to achieve Delaunay-like fair triangulation. The detailed explanation on swapping algorithm can be found in Appendix A Algorithm A3. As it can be seen in line 7, the subdivision and swapping are repeated for n times. Experiments showed that different models and refinement resolutions may require more or less iterations to achieve better triangulation results. In Algorithm 4, we introduced threshold α to match triangle resolution to the original one of the surface mesh. To find threshold α, firstly, the average boundary edge length t should be calculated as in Formula (1). Here, bei is a boundary edge. Next, α is obtained through the multiplication of average boundary edge length t with a value z, as in Formula (2). The value z is a user-defined input, which helps to manipulate the refinement level. In current research, we set five such levels: original, 2 times coarser, 1.5 times coarse, 2 times more refine, 1.

Virtual Cut Simulation Using PBD Model
In the current section the information on cutting simulation experiments is provided. The experiment environment consisted of Intel i7-7700 3.60 GHz CPU with 8 cores, Ge-Force RTX 2080 SUPER GPU, 32 GB RAM and 11 GB V-RAM. The development and tests

Virtual Cut Simulation Using PBD Model
In the current section the information on cutting simulation experiments is provided. The experiment environment consisted of Intel i7-7700 3.60 GHz CPU with 8 cores, GeForce RTX 2080 SUPER GPU, 32 GB RAM and 11 GB V-RAM. The development and tests were executed in Unity 3D game engine version 2019.2.12f1 and Visual Studio 2019 IDE Community edition. In the current study, the PBD part of simulation was executed using GPGPU parallel processing; however, the cutting simulation was implemented purely on CPU.
For the experiments we use 7 models: Sphere20, Sphere80, Sphere 768, Bunny, Octopus, Banana, and Armadillo. The smallest model in the test consists only of 20 triangles, and the largest one has over 30,000 triangles.

Slicing Test
PBD was applied to all experimental models to simulate correct physical behavior. In the current experiment the models were cut with a large sweep plane, which divided them into two parts. According to Algorithm 2, we executed hole patching after deleting the intersected triangles. The initial triangulation was further refined using the Algorithm 4 in order to match the patch triangulation with the original surface mesh. Figures 5 and 6 show the results of slicing test. In Figure 5, we provide the cutting simulation of simple models: Sphere20, Sphere80, and Sphere768. Figure 6 presents the cutting results of more complex models, such as Bunny and Banana. In all cases, the algorithm provided fair and even triangulation matched to the original resolution. In the case of Sphere20 and Sphere80, the triangle subdivision and swapping were finished in 2-3 iterations. However, for bigger models it required more than 3 iterations of refinement to achieve the best results.

Refinement Level Test
In the current section we present the refinement test results. As it was mentioned in Section 2.2, the Algorithm 4 proposes five options of refinement level, which can be controlled through the user input. Figures 7 and 8 show different refinement levels of simple models Sphere20 and Sphere 768. Figures 9 and 10 provide cutting simulation results for more complicated Bunny model. In case of the smallest model, Sphere20, the 2 times coarser, 1.5 times coarser, and original refinement resolutions produce the same results in topology update. The reason for such behavior is that the model is too small to accommodate triangles with larger edges. In other cases, the refinement algorithm provides fair and even triangulation, which corresponds to the original resolution of surface mesh triangles. In Figures 7-10, the second parts of the models were intentionally deleted to show, in detail, the refinement levels.
the intersected triangles. The initial triangulation was further refined using the Algorithm 4 in order to match the patch triangulation with the original surface mesh. Figures 5 and 6 show the results of slicing test. In Figure 5, we provide the cutting simulation of simple models: Sphere20, Sphere80, and Sphere768. Figure 6 presents the cutting results of more complex models, such as Bunny and Banana. In all cases, the algorithm provided fair and even triangulation matched to the original resolution. In the case of Sphere20 and Sphere80, the triangle subdivision and swapping were finished in 2-3 iterations. However, for bigger models it required more than 3 iterations of refinement to achieve the best results.

Refinement Level Test
In the current section we present the refinement test results. As it was mentioned in Section 2.2, the Algorithm 4 proposes five options of refinement level, which can be controlled through the user input. Figures 7 and 8 show different refinement levels of simple models Sphere20 and Sphere 768. Figures 9 and 10 provide cutting simulation results for more complicated Bunny model. In case of the smallest model, Sphere20, the 2 times coarser, 1.5 times coarser, and original refinement resolutions produce the same results in topology update. The reason for such behavior is that the model is too small to accommodate triangles with larger edges. In other cases, the refinement algorithm provides fair and even triangulation, which corresponds to the original resolution of surface mesh triangles. In Figures 7-10, the second parts of the models were intentionally deleted to show, in detail, the refinement levels.    Overall, the proposed solution can handle different kinds of holes and produces satisfactory refinement results. In Figure 10, it can be observed that the hole patching and refinement algorithm manages even thin difficult shapes, such as the ear part on the bunny model. The results of cutting simulation for other models can be found in Appendix A, Figures A1-A4.  Overall, the proposed solution can handle different kinds of holes and produces satisfactory refinement results. In Figure 10, it can be observed that the hole patching and refinement algorithm manages even thin difficult shapes, such as the ear part on the bunny model. The results of cutting simulation for other models can be found in Appendix A, Figures A1-A4.  Overall, the proposed solution can handle different kinds of holes and produces satisfactory refinement results. In Figure 10, it can be observed that the hole patching and refinement algorithm manages even thin difficult shapes, such as the ear part on the bunny model. The results of cutting simulation for other models can be found in Appendix A, Figures A1-A4.

Results
To evaluate the proposed algorithm, firstly, we compared frames per second (FPS) rate of only PBD simulation of soft body and the case of PBD with cutting simulation (CS). The FPS measurements of the seven test models are presented in Figure 11. Here, the y axis stands for the FPS rate, and the x axis represents the time of simulation in seconds. In both cases first 15 s of simulation were evaluated.   Table 1 we show the FPS difference between the two cases. The results from Figure 11 and Table 1 reveal that, in case of only PBD, as well as in PBD with the cutting simulation all models, except for Armadillo, real-time FPS rate is achieved. The FPS difference and its percentage representation from Table 1 shows that the bigger models demonstrate greater FPS drops during cutting simulation. The reason for the FPS decrease is that the collision detection algorithm consumes a considerable amount of computational costs, which grow as the number of elements in the mesh increases.   Table 1 we show the FPS difference between the two cases. The results from Figure 11 and Table 1 reveal that, in case of only PBD, as well as in PBD with the cutting simulation all models, except for Armadillo, real-time FPS rate is achieved. The FPS difference and its percentage representation from Table 1 shows that the bigger models demonstrate greater FPS drops during cutting simulation. The reason for the FPS decrease is that the collision detection algorithm consumes a considerable amount of computational costs, which grow as the number of elements in the mesh increases. Beside the FPS comparison, current research also provides data on average time per cut in milliseconds, as shown in Table 1. The proposed method exploits surface triangle meshes as the geometry representation to reduce computational costs and avoid tetrahedron subdivision, such as in [19,21]. Therefore, the per-cutting time includes time spent on collision detection between surface mesh and the sweep plane, topology update and refinement of the triangulated area. In the case of three sphere models and Bunny model, the per-cut time demonstrates real-time rates, and Octopus model achieves interactive time. For the two largest models, Banana and Armadillo, it takes a great deal of time for cutting due to large number of triangles in the surface mesh. However, comparing to previous studies, we execute rather large cuts and provide different levels of triangle refinement, as were shown in Sections 3.1 and 3.2, accordingly. In case of large models, such as Banana and Armadillo, triangulation and refinement lead to a decrease in simulation performance, as they require more computational costs if a larger number of triangles must be subdivided.
Overall, it is shown that the proposed algorithm produces fast and plausible simulations of large cuts with various refinement levels. The implemented triangulation and refinement algorithms produce even and satisfactory triangulation of holes during cutting simulation, such as those shown in Figures 9 and 10. Such effects were achieved by introducing the pre-subdivision swapping method and through controlling the iteration number of the refinement algorithm. Another advantage of the proposed method is reducing computational costs through avoiding tetrahedral volumetric meshes. However, further research on volume preservation should be done to confirm the reliability of the simulation. One more drawback of the proposed model is that it is restricted to slicing simulations, which also induce high computational costs when the simulated mesh grows in number of elements. The future research should be focused on developing cost-effective progressive cutting method.

Conclusions
In current research, we proposed an algorithm to execute large cutting simulations and compared cutting FPS and average per cut time for seven different models. The results showed that, even executed purely on CPU, the cutting simulation can be done in real-time for small models and interactive time for the larger ones. More than that, we introduced the method to refine or coarsen the area of cutting. However, the current method ran into a few limitations. The first one is that the proposed approach can handle a cutting simulation, in real time and interactively, only for models with the number of elements under 5000. This is due to high computational costs of collision detection and refinement algorithms in large models. Another drawback of the proposed method is that it lacks volumetric representation, even though the PBD's stretching and bending constraints provide some volume preservation during simulation. The last limitation of the algorithm is that it can't be applied in a progressive cutting simulation and cracks simulation. Therefore, our future goal is to implement the proposed algorithm, using GPGPU parallel programming to improve cutting performance and expand the capacity, to process large models. Nonetheless, based on our method, we aim to update it and achieve an effective progressive cutting simulation. Besides, the future research should be focused on volume representation of the simulated soft body without involving computationally expensive tetrahedral meshes. Overall, current research is expected to be useful in surgery simulation, especially the simulation of comparatively large cuts. As the proposed method was implemented in Unity game engine, we anticipate it to be of help to other users when simulating cuts in games and other applications.  if(index equal to 0) 8: add new triangle t new with vertices startPointe last , endPointe index and startPointe index to the patchTrianglesList 9: remove e last and e index 10: add new edge e new with vertices startPointe last and startPointe index 11: else 12: add new triangle t new with vertices startPointe index-1 , endPointe index and startPointe index to the patchTrianglesList 13: remove e index and e index-1 14: add new edge e new with vertices startPointe index-1 and startPointe index 15: add new key e new to edgeDictionary 16: add current triangle's index to key e new in edgeDictionary 17: end for 18: add new triangle t new with vertices startPointe last , endPointe 0 and startPointe last to the patchTrianglesList 19: return patchTrianglesList Algorithm A3. Swap 1: Require: set of patch triangles patchTrianglesList 2: Require: set of vertices v 3: Require: triangles indices t 1 and t 2 of key ed x fromedgeDictionary 4: find the index v 1 of triangle t 1 , which is not equal to startPointed x and endPointed x 5: find the index v 2 of triangle t 2 , which is not equal to startPointed x and endPointed x 6: if(edgeDictionary contains edge with vertices v 1 and v 2 ) return false 7: let ac ← v 1 -startPointed x 8: let ab ← endPointed x -startPointed x 9: crossProduct ← ac x ab 10. toCenter ← ((crossProd x ab) * ||ac|| + (ac x crossProd) *||ab||)/2*||crossProd|| 11. radius ← ||toCenter|| 12: center ← startPointed x + toCenter 13: if(||v 2 -center|| < radius) 14: swap ed x with new edge e new with virtices v 1 and v 2 15: return true 16: end if 17: let ac ← v 2 -startPointed x 18: let ab ← endPointed x -startPointed x 19: crossProduct ← ac x ab 20: toCenter ← ((crossProd x ab) * ||ac|| + (ac x crossProd) *||ab||)/2*||crossProd|| 21: radius ← ||toCenter|| 22: center ← startPointed x + toCenter 23: if(||v 1 -center|| < radius) 24: swap ed x with new edge e new with virtices v 1 and v 2 25: return true 26: return false