An All-In-One Geometric Algorithm for Cutting, Tearing, and Drilling Deformable Models
AAn all-in-one geometric algorithm for cut-ting, tearing, drilling deformable models
Dr. Manos Kamarianakis and Prof. George Papagiannakis
Abstract.
Conformal Geometric Algebra (CGA) is a framework thatallows the representation of objects, such as points, planes and spheres,and deformations, such as translations, rotations and dilations as uniformvectors, called multivectors . In this work, we demonstrate the merits ofmultivector usage with a novel, integrated rigged character simulationframework based on CGA. In such a framework, and for the first time,one may perform real-time cuts and tears as well as drill holes on arigged 3D model. These operations can be performed before and/or aftermodel animation, while maintaining deformation topology. Moreover,our framework permits generation of intermediate keyframes on-the-flybased on user input, apart from the frames provided in the model data.We are motivated to use CGA as it is the lowest-dimension extensionof dual-quaternion algebra that amends the shortcomings of the major-ity of existing animation & deformation techniques. Specifically, we nolonger need to maintain objects of multiple algebras and constantly trans-mute between them, such as matrices, quaternions and dual-quaternions,and we can effortlessly apply dilations. Using such an all-in-one geo-metric framework allows for better maintenance and optimization andenables easier interpolation and application of all native deformations.Furthermore, we present these three novel algorithms in a single CGArepresentation which enables cutting, tearing and drilling of the inputrigged model, where the output model can be further re-deformed ininteractive frame rates. These close to real-time cut,tear and drill algo-rithms can enable a new suite of applications, especially under the scopeof a medical VR simulation.
Mathematics Subject Classification (2010).
Primary 68U05.
Keywords.
Conformal Geometric Algebra (CGA), Skinning, Inter-polation, Cutting Algorithm, Tearing Algorithm, Drilling Algorithm,Keyframe Generation.
The authors are affiliated with the University of Crete, Greece and the ORamaVR company( ).This is an extended version of work originally presented in the CGI 2020 conference, on theENGAGE workshop. [1]. a r X i v : . [ c s . G R ] F e b
1. Introduction
This is an extended version of the work originally presented in the CGI2020 ENGAGE workshop[1]. In this extended work, we introduce a novelalgorithm to perform drill in the rigged model and provide further backgroundknowledge for representing and applying translations, rotations and dilationsin multivector form. Furthermore, we give better insight regarding multivectorinterpolation and provide the updated performance results of our optimizedcutting algorithm.Rigged models and their animation and deformation techniques havebeen among the most studied topics in computer graphics since their inception,and especially in the past few years due to the rapid growth of the industryof Virtual/Augmented Reality and computer games.Although the linear-blend skinning algorithm for rigged models ([2]) hasnot radically changed over the years, the demand for more robust and efficientreal-time implementations of the animation, led researchers into developingmore complex mathematical frameworks to enhance the overall performanceand decrease running times. Originally [3], the animation techniques werebased on matrix representation of the three basic deformations: translation,rotation and dilation. The core idea was to be able to apply these defor-mations to 3D point by simply multiplying the respective matrices, in thedesired order, with the homogeneous coordinates of the point. Since matrixmultiplications are extremely fast to perform due to GPUs’ ability of parallelprocessing, matrices became and still remain the favorite representation classof deformations for the majority of current state-of-the-art skeletal animationframeworks.The major drawback of using matrices was discovered when the need ofcreating interpolated keyframes highlighted the fact that the interpolationresult of two rotation matrices does not correspond to a rotation matrix. Theidea of using the original Euler angles instead of the derived rotation matrixdid not solve the problem as it yielded an even greater one; the famous gimballock . Modern implementations tackle the issue, using quaternions ; an algebraof 4 dimensions, originally introduced by Hamilton in 1843. Quaternions,often denoted by H , are an extension of the complex numbers, using two morenegative dimensions, i.e., they include, besides i , two more distinct imaginarybasis elements j, k such that j = k = − . It was proved that a certainsubset of quaternions, called unit quaternions , could encapsulate the essenceof a rotation and also support interpolation.The idea of using unit quaternions to store rotations provided a solutionto the matrix interpolation problem and remains until today the world stan-dard in computer graphic’s modern engines. However, it also introduced theneed to constantly transmute rotations from quaternion to matrix form andvice versa in every intermediate step, adding an extra performance burden.Matrices are still needed in such implementations to store translation anddilation data, while vertices are kept in homogeneous coordinates. geometric algorithm for cutting, tearing, drilling deformable models 3As in improvement to this situation, an algebraic extension of quaternionscalled dual quaternions was used [4]. A specific subset of these 8-dimensionalobjects, called the unit dual quaternions , was proved to be able to uphold bothrotation and translation data and still allow for effortless and inexpensive linearblending. Nevertheless, this advance did not solve the uniformity problem,however it reduced artifacts appearing during animation [5], while furtherpost-processing can be used to further minimize them [6].Our approach utilizes the Conformal Geometric Algebra (CGA) frame-work to perform both model animation and more complex techniques such ascutting, tearing and drilling. CGA is an algebra containing of dual-quaternions,where all entities such as vertices, spheres, planes as well as rotations, trans-lations and dilation are uniformly expressed as multivectors [7, 8, 9]. Theusage of multivectors allows model animation without the need to constantlytransmute between matrices and (dual) quaternions, enabling dilation to beproperly applied along with translation [10, 11]. Furthermore, the interpolationof two multivectors of the same type correctly produces the expected inter-mediate result [12], which makes creation of keyframes trivial to implement.Finally, usage of the proposed framework demands a single representationtype for all data and results, which is the current trend in computer graphics[13]. The use of Conformal Geometric Algebra and multivector representationallows the creation of simpler algorithms to perform complicated tasks, asfundamental geometric predicates are baked in the framework. For example,the intersection of two planes can be determined by simply evaluating theygeometric product.Therefore, complex operations such as cutting, tearing and drilling amodel are now easier to be accomplished, with near real-time results. Suchoperations have become a major research topic as they appear in increasingfrequency in real-time simulation applications, for both academic as well asindustrial purposes. Current algorithms [14, 15] handle such deformationsusing tetrahedral mesh representations of the model, which demands a heavypre-processing to be performed. Since originally introduced, cutting andtearing methods have been upgraded and extended to allow almost real-timeresults, using mostly finite element methods, intuitive optimization and heavypreprocessing [16, 17, 18]. To make the final results even more realistic, physicsengines utilizing position-based dynamics are used to simulate soft-tissue cutsat the expense of performance [19, 20, 21]. Our contribution:
The novelty of our work involves the complete im-plementation of rigged model animation in terms of CGA, extending the workof Papaefthymiou et al. [11] in a python-based implementation that enableskeyframe generation on-the-fly . The original animation equation involvingmatrices is translated to its equivalent multivector form (see Section 4.1) andall information required to apply the linear blend skinning algorithm (vertices,animation data) is obtained from the model and translated as multivector.This enables us to have future animation models in CGA representation only, M. Kamarianakis and G. Papagiannakiswhich, in combination with an optimized GPU multivector implementation,produces faster results under a single framework. Another major noveltyof our work is the cutting, tearing and drilling algorithms that are beingapplied on top of the previous framework; given the input animated model,we perform cuts, tears and drills on the model surface with the ability tofurther re-deform the newly processed model. The subpredicates used inthese algorithms utilize the multivector form of their input, so they can beimplemented in a CGA-only framework. Their design was made in such away that little to no pre-processing of the input model is required whileallowing a future integration with a physics engine. Furthermore, using ourmethod, we can generate our own keyframes instead of just interpolatingbetween pre-defined ones. Our all-in-one CPU python implementation is ableto process an existing animation model (provided in .dae or .fbx format)and translate the existing animation in the desired CGA form while furthertweaks or linear-blend deformations are available in a simple way to perform.Such an implementation is optimal as far as rapid prototyping, teaching andfuture connection to deep learning is concerned. It also constitutes the basefor interactive cutting, tearing and drilling presented in Section 4.2. Thesimplicity and robustness of our algorithms design promise real-time results ifrun in a compiled programming language such as C++ or C
2. Introduction to Conformal Geometric Algebra
Three dimensional Conformal Geometric Algebra (CGA) can be seen asanother algebra containing dual-quaternions which allows round elementssuch as spheres to be represented as objects of this algebra, i.e., as multivectors .To be more precise, CGA is the lowest possible extension where this is possible.Being able to represent round elements in conjunction with the ability toreflect on objects using the so-called sandwich operation presented in thefollowing sections, CGA is also able to represent dilators (uniform scaling) asmultivectors. Therefore, CGA is a geometric algebra where round elements(points, spheres, circles), flat elements (lines, planes, point pairs) and all basicdeformations (translations, rotations, dilations) can be expressed explicitly inmultivector form.In order to create the model of 3D Conformal Geometric Algebra (CGA),we extend the basis { e , e , e } of the original Euclidean space R by twoelements e + and e − . These elements have positive and negative signaturerespectively, i.e., it holds that e = − e − = 1 . The resulting non-Euclideanspace is usually denoted as R , or G (4 , has 1 negative ( e − ) and 4 positive( e , e , e , e + ) basis elements.It is convenient to define a null basis given by the original basis vectors e , e , e of R and e o = 12 ( e − − e + ) , e ∞ = e − + e + . (1) geometric algorithm for cutting, tearing, drilling deformable models 5The elements e o and e ∞ are called null vectors because e o = e ∞ = 0 , wherethe operation implied is the geometric product described in the followingsections. R , A generic vector Y of R , is a linear combination of the basis elements { e , e , e , e ∞ , e o } , i.e., Y = y e + y e + y e + y ∞ e ∞ + y o e o , y i ∈ R . (2)Note that CGA is a projection space where the elements Y and Z areequivalent if and only if there is a λ ∈ R such that Y = λZ . Due to thisequivalence, we usually assume, without loss of generality, that the coordinateof e o is either 0 or 1. In this algebra, points, spheres and planes are easilyrepresented as vector objects of the space, as described below. Points.
A point x = ( x , x , x ) = x e + x e + x e of R is up-projected into the conformal vector X = x + 12 x e ∞ + e o (3) = x e + x e + x e + 12 ( x + x + x ) e ∞ + e o (4) Spheres.
A sphere s of the R , centered at x = ( x , x , x ) with radius r is up-projected into the conformal vector S = X − r e ∞ (5) = x e + x e + x e + 12 ( x + x + x − r ) e ∞ + e o , (6)where X is the image of x in R , . Planes.
A plane π of the original space, perpendicular to the normal vector (cid:126)n = ( n , n , n ) , such that d ( π, O ) = d (where d ( · , · ) is the Euclideandistance) is up-projected into the conformal vector Π = (cid:126)n + de ∞ (7) = n e + n e + n e + de ∞ . (8) R , There are 3 major products in R , : the inner, the outer and the geometric.Each of these products is initially defined among the vectors of either baseof the space, i.e., { e , e , e , e − , e + , e o , e ∞ } . The respective definition is thenextended to included any vector (or multi-vector ) of the space. Below wepresent some of the basic properties of these products; further informationcan be found in [8, 7]. Inner.
The inner product (denoted by · ) of the basis elements is definedas follows: • e i · e j := δ ij for i, j ∈ { , , , + }• e − · e − := − • e − · e j := 0 for j ∈ { , , , + } M. Kamarianakis and G. Papagiannakis • e o · e o := e ∞ · e ∞ = 0 • e o · e ∞ := − • e i · e j := 0 for i ∈ { , , , + } and j ∈ { o, ∞} Outer.
The outer product of the basis elements e i and e j is denoted as e i ∧ e j . The outer product is an associative operation that can be appliedto more than two elements, e.g., e i ∧ e j ∧ e k and e i ∧ e j ∧ e k ∧ e ∞ areproperly defined. The outer product of k basis vectors is called a k - blade and k is usually referred to as the grade of this blade. A sum of k -bladesis called a k -vector and the addition of k -vectors of different grades is a multivector .The importance of the outer product derives from the fact that itallows us, in certain cases, to obtain the intersection of two objects bysimply evaluating their outer product. Specifically, a circle (resp. line)can be seen as the intersection - outer product of two spheres (resp.planes). The outer product of a circle with an intersecting sphere orequivalently, the outer product of three intersecting spheres represent aset of two points, usually referred to as a point pair . Geometric.
The most important product in R , is the so-called geometric product. For basis vectors e i and e j , their inner product e i e j is definedas the addition of the outer and inner product of the elements, i.e., e i e j := e i ∧ e j + e i · e j . (9)Note that, by definition, e i e j = e i ∧ e j for every i, j ∈ { , , , ∞ , o } suchthat { i, j } (cid:54) = {∞ , o } . The inner product of random multivectors canbe obtained by evaluating the sum of their respective inner and outerproducts. First, let us denote the pseudoscalar I of R , , I := e ∧ e ∧ e ∧ e + ∧ e − = e ∧ e ∧ e ∧ e ∞ ∧ e o . (10)Using I , we may define the dual object O (cid:63) of a multivector O is to be O (cid:63) := − OI, (11)where the operation between O and I is the geometric product. Notice that itholds that ( O (cid:63) ) (cid:63) = − O and therefore we can easily obtain the normal form O of an object from it’s dual form O (cid:63) and vice versa.The dual form of certain objects holds strong geometric meaning, asdescribed below. • The outer product of 4 non-coplanar points yields the dual form of thesphere defined by these points. • The outer product of 3 non-collinear points and e ∞ yields the dual formof the plane defined by these points. • The outer product of 3 non-coplanar points yields the dual form of thecircle defined by these points. geometric algorithm for cutting, tearing, drilling deformable models 7 • The outer product of 2 points and e ∞ yields the dual form of the linedefined by these points. So far we have shown that objects (or their duals) such as points, planes, circles,spheres, lines and point pairs are represented as multivectors. However, thebeauty and versatility of this algebra comes from its ability to also representrotations, translations and dilations as multivectors as described below.
Rotation.
A rotation in CGA is encapsulated in a multivector R := exp( − b φ − I u φ φ − uI sin( φ , (12)where φ is the angle of the rotation, b is the normalized plane of therotation, u is the normalized axis of the rotation and I := e e e . Alloperation are geometric products and exp( · ) denotes the exponentialfunction. Translation .
The multivector T := exp( − te ∞ ) = 1 − te ∞ , (13)where t = t e + t e + t e is a euclidean vector, represents a translationby t in CGA. Dilation .
The multivector D = 1 + 1 − d d e ∞ ∧ e o , (14)corresponds to a dilation of scale factor d with respect to the origin.The conformal space model allows us to apply any or multiple of theoperations above not only to a point but also to any object O that waspreviously defined. To apply the transformations M , M , . . . M n (in thisorder), to an object O , we first define the multivector M := M n M n − · · · M ,where all in-between products are geometric. Let M − denote the inverse of M , i.e. the resulting multivector after all elements of the form e i ∧ e j ∧ · · · , iswritten in the reverse order · · · ∧ e j ∧ e j . Then, the object O (cid:48) := M OM − (15)represents the final form of O after all instead transformations are applied. Interpolation of data is an essential part for Computer Graphics as it isneeded in every animation procedure of a rigged model. The poses of themodel with respect to time are not stored in a continuous manner but ratherat discrete time-steps. If additional intermediate frames are demanded, wehave to interpolate the animation data between two provided keyframes.As in the case of matrix ([3]) or (dual) quaternion quaternion interpola-tion [5], a blending of two multivectors can be accomplished in various ways,yielding different results [12, 22]. Choosing a proper interpolation techniqueis not a simple task as it may depend on the model or other factors. However, M. Kamarianakis and G. Papagiannakis (a) (b) (c)
Figure 1.
Linear Versus Logarithmic Interpolation. A singleface is interpolated between the bottom blue and top greykeyframe. Both translation ( t ), rotation ( r ) and scaling( s )data that have to be interpolated where in multivector form.In (A), we have used a classic linear interpolation, whereas in(B) we have used a logarithmic blending. In (C), the outcomesof these two methods are superimposed.two methods remain dominant in analogue with the quaternion case: thelinear and the logarithmic blending.Linear blending of the multivectors m and m , which, in a modelanimation context, may represent translations, rotations or dilations , is doneby evaluating (1 − α ) · m + α · m , for α ∈ [0 , . Another blending method is theso-called logarithmic interpolation where we evaluate m exp( α · log( m − · m )) ,for α ∈ [0 , , where the exponential and logarithmic function of a multivector m are approximated in our case using the respective Taylor series expansion.Using different blendings, we obtain different results, as shown in Figure 1.More information regarding the evaluation and properties of multivectorlogarithms/exponentials can be found in [23, 10, 24, 25].In our framework, linear blending is preferred when generating frameson-the-fly by the user, whereas logarithmic blending is used when reading themodels existing animation data.
3. State of the art
The current state of the art regarding skeletal model animation and defor-mation is based on the linear-blend skinning algorithm ([2]) and the repre-sentation of bones animation via transformation matrices and quaternions geometric algorithm for cutting, tearing, drilling deformable models 9or dual-quaternions. Such an implementation allows for efficient and robustinterpolation methods between keyframes. A major drawback of such animplementation is that a dilation can not be applied as a scaling matrixalways refers to the origin and not the parent bone [11].To be more precise regarding the mechanics of the deformation process,in the case of a simple rigged model, every bone b i amounts to an offset matrix O i and an original transformation matrix t i . The skin of the model is importedas a list of vertices v and a list of faces f . A bone hierarchy is also providedwhere { t i } are stored along with information regarding the animation of eachjoint. This information, usually referred to as TRS data , is provided in theform of a quaternion, a translation vector and a scaling vector that representrespectively the rotation, displacement and scaling of the joint with respectto the parent joint for each keyframe (see Section 3.1).In order to determine the position of the skin vertices at any given time k and therefore render the scene by triangulating them using the faces list,we follow the steps described below. Initially, a matrix G is evaluated asthe inverse of the transformation matrix that corresponds to the root node.Afterwards, we evaluate the global transformation matrix for every bone b i attime k and denote it as T i,k . To evaluate all T i,k , we recursively evaluate thematrix product T j,k t i,k where b j is the parent bone of b i , given that T r,k is theidentity matrix (of size 4), where b r denotes the root bone. The matrix t i,k isa transformation matrix equal to t i if there is no animation in the model; inthis case, our implementation allows to generate the keyframes ourselves inreal-time. Otherwise, t i,k is evaluated as t i,k = T R i,k
M R i,k S i,k (16)where T R i,k , M R i,k , S i,k are the interpolated matrices that correspond to thetranslation, rotation and scaling of the bone b i at a given time k .After evaluating the matrices { T i,k } for all bones { b i } , we can evaluatethe global position of all vertices at time k , using the rigged deformationequation : V k [ m ] = (cid:88) n ∈ I m w m,n GT n,k O n v [ m ] (17)where • V k [ m ] denotes the skin vertex of index m (in homogeneous coordinates)at the animation time k , • I m contains up to four indices of bones that affect the vertex v [ m ] , • w m,n denotes the “weight”, i.e., the amount of influence of the bone b n on the vertex v [ m ] , • O n denotes the offset matrix corresponding to bone b n , with respect tothe root bone, • G denotes the inverse of the transformation matrix that corresponds tothe root bone (usually equals the identity matrix) and • T n,k denotes the deformation of the bone b n at animation time k , withrespect to the root bone.0 M. Kamarianakis and G. Papagiannakis The modern way to represent the TRS data of a keyframe is to use matricesfor the translation and dilation data as well as quaternions for the rotationdata. Let { T R i , R i , S i } , denote such data at keyframe i ∈ { , } , where: • T R i = x i y i z i and S i = sx i sy i sz i
00 0 0 1 represent thetranslation by ( x i , y i , z i ) and the scale by ( sx i , sy i , sz i ) respectively and • R i is a quaternion representing the rotation.Before quaternions, Euler angles and the derived rotation matrices wereused to represent rotation data. However the usage of such matrices induced agreat problem: a weighted average of such matrices does not correspond to arotation matrix and therefore interpolating between two states would requireinterpolating the Euler angles and re-generate the corresponding matrix. Thisin turn would sometimes lead to a gimbal lock or to ‘candy-wrapper’ artifactssuch as the ones presented in [5].The usage of quaternions allowed for easier interpolation techniques whileeradicating such problems. Nevertheless, a transformation of the interpolatedquaternion to corresponding rotation matrix was introduced since the GPUcurrently handles only matrix multiplications in a sufficient way for skinningreasons. Therefore, the interpolation between the two keyframes mentionedabove follows the following pattern:1. the matrices T R a = (1 − a ) T R + aT R and S a = (1 − a ) S + aS areevaluated for a given a ∈ [0 , ,2. the quaternion R a = (1 − a ) R + aR is determined and finally,3. the rotation matrix M R a that corresponds to R a is calculated.The interpolated data T R a , M R a and S a are then imported to the GPUin order to determine the intermediate frame, based on the equation(17).Using the method proposed in this paper, all data are represented inmultivector form. A major implication of this change is that the interpolationbetween two states is done in a more clear and uniform way as presented inSection 4. This also makes the need to constantly transform a quaternionto a rotation matrix redundant, although we are now obliged to performmultivector additions and multiplications as well as down project pointsfrom R , to R to parse them to the GPU. However, since all our dataand intermediate results are in the same multivector form, we could (ideally)program the GPU to implement such operations and therefore greatly improveperformance. geometric algorithm for cutting, tearing, drilling deformable models 11 (a) (b) (c) Figure 2.
Skinning via multivectors versus skinning via dualquaternions. The original model is deformed using multivec-tors and depicted in magenta wireframe, superimposed withthe color-graded result (based on the z coordinate of eachvertex) of the quaternion method for the same deformation.It is qualitatively verified that linear blending of multivec-tors produces similar results with the current state-of-the-artmethod. Evaluating the vector differences of all vertices forthe two methods, we have evaluated the approximation errorassuming the quaternion method to be the correct, usingthe infinity ( (cid:96) ∞ ) norm. (A) Applying rotation on a bone,approximation error . . (B) Applying dilation on a bone,approximation error . . (C) Applying translation, ap-proximation error . The model used contains 1261 verticesand 1118 faces.
4. Our Algorithms and Results
The deformation equation (17), core of the animation algorithm, yields fastresults (especially when combined with a GPU implementation) but deniesus a robust way to dilate with respect to a bone. Our motivation is to extendand apply the animation equation for multivector input as proposed in [11].To be more specific regarding our method, we propose the replacementof all matrices appearing in (17) with multivectors for animation purposes.The transformation matrix of t i of each bone b i as well as all informationregarding translation and rotation for each keyframe can be easily convertedto multivectors [7, 8]. Consequently, we can evaluate the multivector M i,k which is equivalent to the matrix T i,k by following the same procedure of2 M. Kamarianakis and G. Papagiannakisdetermining the latter (described in Section 3) while substituting all involvedmatrices with the corresponding multivectors.Note that various techniques can be used to interpolate between twokeyframes to obtain M i,k ; for existing keyframes logarithmic blending ispreferred [12, 5], whereas for keyframe generation we use linear blending. Inboth scenarios, the intermediate results are multivectors of the correct type.Furthermore, each offset matrix O n and each skin vertex v [ m ] is trans-lated to their CGA form B n and c [ m ] respectively. Finally, G matrix isnormalized to identity and is omitted in the final equation.Our final task is to translate in CGA terms the matrix product T n,k O n v [ m ] ,where apparently each multiplication sequentially applies a deformation tovertex v [ m ] . To apply the respective deformations, encapsulated by M n,k and B n , to CGA vertex c [ m ] , we have to evaluate the sandwich geometric product ( M n,k B n ) c [ m ]( M n,k B n ) (cid:63) where V (cid:63) denotes the inverse multivector of V (see[4, 7] for details).Summarizing, if the multivector form of the vertex V k [ m ] , which corre-sponds to the final position of the m -th vertex at animation time k , is denotedby C k [ m ] , then the multivector deformation equation becomes C k [ m ] = (cid:88) n ∈ I m w m,n ( M n,k B n ) c [ m ]( M n,k B n ) (cid:63) (18)After the evaluation of C k [ m ] for all m , we can down-project all these conformalpoints to the respective euclidean ones in order to represent/visualize themand obtain the final result of the keyframe at time k .The replacement of matrices with multivectors enables the introductionof dilations in a simple way. The multivector M i,k that represents a rotationand translation with respect to the parent bone of b i can be replaced with M i,k D i,k where D i,k is the corresponding dilator and the operation betweenthem is the geometric product. The dilator corresponds to a scale factor withrespect to the parent bone, information that could not be easily interpretedvia matrices. However, since the application of a motor and/or a dilator to avertex is a sandwich operation, such a dilation becomes possible when usingmultivectors.A comparison between the results of our proposed method and thecurrent state-of-the-art is shown in Figure 2, where we successfully applydilation to different bones and obtain similar results. Rotations, dilationsand translations are obtained in our method using multivectors only, under asingle framework with simpler notation/implementation; linear blending isused to interpolate between keyframes. A novelty we present in this paper is the cutting, tearing and drilling algorithmson skinned triangulated models. As the name suggests, the first module enablesthe user to make a planar cut of the model whereas the second is used toperform smaller intersections on the skin. The last module can be utilized todrill holes in the skinned model. In the following sections, we provide a detailed geometric algorithm for cutting, tearing, drilling deformable models 13presentation of the algorithms involved as well as certain implementationdetails.
Cutting a skinned model is implemented incurrent bibliography in many forms [14, 26, 27, 28, 29]. The most commontechnique is via the usage of tetrahedral meshes ([16]) which require a heavypre-processing on the model and currently do not enable further animation ofthe model or scale to VR enviroments. Our work includes an algorithm forplanar model cut, where the final mesh is deformable, as we implemented afunction to calculate weights for all additional vertices that did not originallyexist (see Figure 3). Most of the subpredicates used in the cutting algorithmare implemented in terms of conformal geometry and therefore can be usedeven if the model is provided in multivector form.Our proposed planar cut implementation is summarized as Algorithm 1.A description of how we tackle the weight evaluation in step 2 is found inSection 4.3. Our algorithm does not require tetrahedral meshed models andrequires minimum to none pre-processing. It is GA-ready and the low numberof operations it demands make it suitable for VR implementations.
Algorithm 1
Cutting Algorithm
Input:
Triangulated Mesh M = ( v, f ) ( f is the face list), and a plane Π . Output:
Two meshes M = ( v , f ) and M = ( v , f ) , result of M gettingcut by Π Evaluate (using GA) and order the intersection points of Π with each faceof M . Evaluate the weights and bone indices that influence these points. Re-triangulate the faces that are cut using the intersection points. Separate faces in f and f , depending on which side of the plane they lie. From f and f , construct M and M . The purpose of this module is to enable partialcuts on the skinned model, in contrast with the cutting module where thecut is, in a sense, complete. The importance of this module derives from thefact that most of the surgical incisions are partial cuts and therefore they areworth replicating in the context of a virtual surgery. Towards that direction,our work involves an algorithm that both tears a skinned model and alsoenables animation of the final mesh (see Figures 4 and 5).To understand the philosophy behind the design of the tearing algorithmthat is described below, one must comprehend the differences between cuttingand tearing. In tearing, the movement of a scalpel defines the tear ratherthan a single plane. To capture such a tear in geometric terms, we have totake into consideration the location of the scalpel in either a continuous way(e.g. record the trail of both endpoints of the scalpel in terms of time) or adiscrete way (e.g. know the position of the scalpel at certain times t i ). For VRpurposes, the latter way is preferred as it yields results with better fps, since4 M. Kamarianakis and G. Papagiannakis (a) (b)(c) (d) Figure 3.
Cutting module intermediate steps. (A) The orig-inal animated model. (B) The model where the (red) intersec-tion points of the cutting plane and the mesh are calculatedand re-triangulated. (C) The model after the cut. (D) Themodel is deformed by a rotation (axis= (0 , , , . rad), atranslation (vector= (13 , , ) and a dilation (factor = 0.5)at joint 1 (elbow), as well as another rotation (axis= (0 , , , . rad) at joint 2 (wrist). Note that minimal artifacts occurin the final result. The vertices in (D) are colored dependingon the influence of joint 1 which is mostly deformed. Thevertices in (A)-(C) are colored based on their z coordinate.input is hard to be monitored and logged continuously in a naive way. For geometric algorithm for cutting, tearing, drilling deformable models 15these reasons, our implementation requires the scalpel position to be knownfor certain t i .The proposed tearing algorithm is summarized in Algorithm 2. A descrip-tion of how we tackle the weight evaluation in step 4 is found in Section 4.3. Algorithm 2
Tearing Algorithm
Input:
Triangulated Mesh M = ( v, f ) , and scalpel position at time steps t i and t i +1 Require:
Scalpel properly intersects M at these time steps Output:
The mesh M t = ( v t , f t ) resulting from M getting torn by the scalpel Determine the intersection points S i and S i +1 of M with the scalpel attime step t i and t i +1 respectively. Determine the plane Π , containing S i and the endpoints of scalpel at time t i +1 . Small time steps guarantee that Π is well-defined. Evaluate the intersection points Q j of Π and M , s.t. the points S i , Q , Q , . . . , Q m , S i +1 appear in this order on Π when traversing theskin from S i to S i +1 . Assign weights to points S i , S i +1 and all Q j . Re-triangulate the torn mesh, duplicating Q j vertices. Move the two copies of Q j away from each other to create a visible tear(optional).Our major assumption is that all intermediate intersection points lie onthis plane, which is equivalent to the assume that the tearing curve is smooth,given that t i and t i +1 are close enough. In our implementation, during step 6,the intermediate torn points are moved parallel to the direction of the normalof the plane Π and away from it, to replicate the opening of a cut humantissue. The usage of Virtual Reality by surgeons andtheir need to drill holes in a simulated 3D model motivated the creation ofthe drilling module. Given a triangulated mesh and finite cylindrical drill, wewould like to evaluate the mesh that corresponds to the drilled model.Designing the drilling predicate was more intriguing, compared to therespective cutting and tearing algorithms, as multiple ideas turned out tobe inadequate. The initial idea of substituting the cylinder with a prism of n -surfaces, for some suitable n , looked promising enough, as it would enableusing drill as a special case of tear. However, one would have to provide aneasy way to determine an n that would be sufficiently large to produce asmooth hole-like effect in the outcome mesh. On the other hand, choosing anarbitrary large n would result in many surfaces and therefore many costlytear operations had to be performed, hindering our chances of a real-timeimplementation. The prismatic approach also yielded the question of how tochoose the position the edges of the prism such that the intersection points ofthe prism and the mesh would be re-triangulated in a clever and robust way.6 M. Kamarianakis and G. Papagiannakis (a) (b) (c) Figure 4.
Tearing module intermediate steps. (A) The orig-inal animated model and the scalpel’s position at two con-secutive time steps. (B) The plane defined by the scalpels(depicted as a red triangle) intersects the skin in the yel-low points. (C) The intermediate points are used in there-triangulation, and are «pushed» away from the cuttingplane to form an open tear.Of course, if the edges of the prism were selected such that they intersectedthe mesh’s faces only on their boundaries, the re-triangulation would be moreefficient and not produce a lot of slither faces. However, if we had to decide theoptimal prism, that would be equivalent to specify the intersection points ofall mesh edges with our initial cylinder, which is the idea behind our proposedalgorithm.In the core of our drill module lies a point-versus-cylinder predicatethat allows us to determine the intersection point of every edge of the givenmesh with the cylindrical drill. Since the drill is described by its radius r and two endpoints A (the “tip” of the drill) and B that define its axis,we can easily determine the plane Π that is perpendicular to its axis andgoes through B . Given an edge e defined by the vertices v i and v j of themesh, we first determine if any of these two vertices lie inside the semi-finitecylinder (we ignore the existence of A for now and consider that the cylinderis only bounded by Π and goes indefinitely towards the direction of A ). Toaccomplish such task for the vertex v ∈ { v i , v j } , we project it to the plane Π and compare the distance of the projected point P ( v ) and B with r ; if it geometric algorithm for cutting, tearing, drilling deformable models 17 (a) (b) (c) (d) Figure 5.
Deformation of a torn model. (A) The originalmodel after applying the tear. (B) Two rotations are appliedto the torn model, one at elbow joint around y -axis by -1 rad, and another at wrist joint around y -axis by 1 rad.(C) A dilation of scale 1.5 is applied to the torn model, atelbow joint. (D) A translation is applied to the torn modelat elbow joint with translation vector (18 , , . In all cases,minor artifacts only arise, despite the great magnitude ofthe applied deformations. In (B),(C) and (D), vertices arecolored depending on the influence of elbow joint which ismostly deformed. In (A), vertices are colored based on their z coordinate.is smaller (respectively larger) then v lies strictly inside (resp. outside) thecylinder. In the case of equality, the vertex v lies on the cylinder.If the vertices v i and v j lie on different sides with respect to the cylinder,then we first evaluate the intersection point of the edge defined by µ := P ( v i ) and ν := P ( v j ) with the sphere centered at B with radius r . The coordinatesof µ and ν can be explicitly determined as they are the projections of thepoints v i and v j respectively on the plane going through B with normal (cid:126)n = (cid:126)AB/ || (cid:126)AB || . Therefore, the projected image of (cid:126)v ∈ { v i , v j } on the planeis (cid:126)v − (cid:104) (cid:126)n, (cid:126)v − (cid:126)OB (cid:105) (cid:126)n , where (cid:104)· , ·(cid:105) denotes the classic inner product.Since every point on the projected edge is of the form α · µ + (1 − α ) · ν forsome α ∈ [0 , , and the edge is intersected by the sphere - as it is intersectedby the cylinder- there exists an α that corresponds to the intersection point.For this α , the point ξ := α · µ + (1 − α ) · ν must have exact distance from B equal to r . If d ( · , · ) denotes the Euclidean distance, solving the equation d ( ξ, B ) = r in terms of α yields that α is a root of the quadratic equation8 M. Kamarianakis and G. Papagiannakis (a) (b) (c) Figure 6.
Drilling module intermediate steps. (A) The drillintersects the model skin in the yellow points. (B) The inter-section points are used in the re-triangulation. (C) The elbowjoint of the drilled model is translated by (1 , , , rotatedby . rad around all 3 axis and then dilated by a factor of2. The weight function ensures that minimal to no artifactsarise in the drilled area despite the deformation. Kα + Lα + N = 0 , where K = ( x µ − x ν ) + ( y µ − y ν ) + ( z µ − z ν ) (cid:54) = 0 (19) L = 2( x µ − x ν )( x ν − x B ) + 2( y µ − y ν )( y ν − y B )+ 2( z µ − z ν )( z ν − z B ) (20) N = ( x ν − x B ) + ( y ν − y B ) + ( z ν − z B ) − r (21)and ξ = ( x ξ , y ξ , z ξ ) , for ξ ∈ { µ, ν, B } .Therefore, we conclude that α is the only root of the quadratic thatbelongs in [0 , . Since for this α , the point α · P ( v i ) + (1 − α ) · P ( v j ) is theintersection of the cylinder with the projected edge, the point α · v i +(1 − α ) · v j is a good approximation of the intersection point of cylinder with the originaledge.Except of the basic edge-cylinder intersection where the endpoints v i and v j of the edge lie on different sides with respect to the cylinder, another twocases have to be taken into consideration. It is possible that both endpointslie outside of the cylinder but the edge intersects the cylinder in two pointsor is tangent to the cylinder in one point. These cases are equivalent to both P ( v i ) and P ( v j ) lying outside the spheres centered at B with radius r andthe quadratic Kα + Lα + N = 0 has one or two root(s) α ∈ [0 , . As before, geometric algorithm for cutting, tearing, drilling deformable models 19the intersection point(s) is(are) approximated by α · v i + (1 − α ) · v j for these α . After evaluating the intersection points of the drill with the model andsince all of them lie on some edge of the original mesh, a robust and efficienttriangulation can be easily applied. If the number of intersection points isbelow some threshold, e.g. , we can perform a “split” operation on all affectedfaces and drill again. To split a triangular face one may connect the middlepoints of all edges and therefore create four smaller sub-triangles similar tothe original. This operation will create a more “dense” triangulation in thespecific part of the model, resulting in more intersection points with the drilland hopefully in a more realistic result. Although generating more intersectionpoints when needed is not a difficult task, we have to take into considerationthat it has to be done in a clever way so as not to hinder the re-triangulationprocess in terms of performance or implementation complexity.The results of our drilling module when applied to our arm model aredemonstrated in Figure 6. As in the precious modules, we can assign weightsto the newly introduced intersection points, allowing re-deformation of thedrilled model. The outline of the proposed drilling algorithm is described inAlgorithm 3. A summary of how we address the weight evaluation in step 4 isfound in Section 4.3. Algorithm 3
Drilling Algorithm
Input:
Triangulated Mesh M = ( v, f ) and drill position via its endpoints A (“tip” of the drill), B and radius r . Require:
Drill properly intersects M Output:
The mesh M (cid:48) = ( v (cid:48) , f (cid:48) ) resulting from M being drilled Let Π denote the plane perpendicular to the drill axis going through theendpoint B . Determine the faces of M that are pierced by the drill and run a BFSalgorithm that checks, for this and all neighboring faces, if at least one ofits three vertices, when projected to Π , has distance from B less than r ,i.e., if it lies withing the drill. Mark such a face as “affected”. For all affected faces, determine in the plane Π the intersection pointsof the drill and the projected face and then down-project them to theoriginal mesh. Assign weights to the intersection points of the original mesh. Re-triangulate the drilled mesh by replacing the affected faces by appro-priate ones.
The main framework used for skinning and animation with the use of multi-vectors is Python’s PyAssimp and Clifford package for the evaluation of the PyAssimp Homepage: https://pypi.org/project/pyassimp/ Clifford Homepage: https://clifford.readthedocs.io/ R , for each vertex. There are two possibleways of achieving this task. The first way is to evaluate the sum and thendown project the final result to obtain each vertex in Euclidean form. Thesecond way is to down project each term and then add them to get the finalresult. Although not obvious, the second method yields faster results since theaddition of 4 multivectors (32-dimensional arrays) and one down-projection isslower than down-projecting (up to) 4 multivectors and adding 4 euclideanvectors of dimension 3.A final implementation detail regards the weight evaluation for newlyadded vertices in the cutting and tearing modules. In the former module, suchvertices necessarily lie on an edge of the original mesh, whose endpoints bothlie on different sides of the cutting plane. Another method is the one usedin the tearing module where the intersection point can also lie inside a face.Assuming the point X lie somewhere on the face ABC , we can explicitly write OX = pOA + qOB + rOC for some a, b, c ∈ [0 , such that p + q + r = 1 .The tuple ( p, q, r ) is called the barycentric coordinate of X with respect tothe triangle ABC . Each of the vertices
A, B, C are (usually) influenced byup to 4 bones, so let us consider that they are all influenced by a set of N ( ≤ vertices, where the bones beside the original 4 have weight 0. Let w A , w B , w C , w X denote the vectors containing the N weights that correspondto vertices A, B, C and X respectively, for the same ordering of the N involvedbones. To determine w X , we first evaluate w = pw A + qw B + rw C and considertwo cases. If w contains up to 4 non-zero weights, then w X = w . Otherwise,since each vertex can be influenced by up to 4 bones, we keep the 4 greatervalues of w , set the others to zero, and normalize the vector so that the sum ofthe 4 values add to 1; the final result is returned as w X . We denote this weightas weight of X via barycentric coordinates . Variations of this technique can beapplied in both modules to prioritize or neglect influences on vertices lying ona specific side of the cutting plane. Different variations of the weight functionallows for less artifacts ([30]), depending on the model and the deformationsubsequent to the cutting/tearing. Performance:
Running the cutting algorithm in the arm model (5037faces, 3069 vertices) took for a simple cylinder model a total of 898ms: 42msfor vertex separation, 757ms for re-triangulation of the 92 intersection points,87ms to split faces in two meshes and 12ms to update the weights. To cut thearm model, it took 4666ms as shown in Table 4.3, where most time (2205ms) geometric algorithm for cutting, tearing, drilling deformable models 21Time Spent Using Time Spent UsingSubroutine Euclidean Tools GA ToolsSubroutine 1 0,036439578 + sec 0,072434584 secSubroutine 2 2,050984303 sec 2,205727577 secSubroutine 3 0,061544961 secSubroutine 4 2,326937914 secCutting time 4,475906756 sec 4,666645036 sec
Table 1.
Running times of the four main subroutines ofthe cutting algorithm. In the 2nd column, point-versus-planerelative positions for subroutine 1 and segment-plane intersec-tions for subroutine 2 where determined using only Euclideansubpredicates. In the 3rd column, the same operations werecarried using Geometric Algebra equivalent subpredicates.The subroutines 3 and 4 are independent of the model datarepresentation. Subroutines: (1) Check vertices locations withrespect to the cutting plane, (2) Detect which faces are inter-sected by the cutting plane, evaluate the intersection pointsand triangulate them, (3) Evaluate weights for the intersec-tion points, (4) Split original model into submodels.was spent on the evaluation and triangulation of the intersection points of thecutting plane and the model. The offline pre-processing time of the model,i.e., the time required to translate the model skin or animation data fromEuclidean coordinates or matrices respectively to multivector form is nottaken into account in the above measurements.Applying the tearing algorithm to the arm model took 2437ms for thefinal output, for 34 intersection points. Most of this time (2411ms) were neededjust to determine which two faces were intersected by the scalpel. Tearing asimple cylinder model (758 faces, 634 vertices) took 362ms for 17 intersectionpoints. Again, most time (331ms) was spend for the scalpel intersection.The drilling algorithm for the arm model takes on average 274ms fora hole consisting of 17 intersection points on our arm model. For a hole ofthe same diameter consisting of 20 intersection points, the algorithm requires319 ms to return the final outcome whereas, the running time grows to 595ms when the diameter is increased from 2 to 3 and the intersection pointsbecome 33. As a rule of thumb, there is an average running time of 16-18 msper intersection point.These running times can be greatly improved as our current unoptimizedCPU-based Python implementation has to thoroughly search all faces forcuts/tears. A GPU implementation optimized for multivector operations wouldallow the comparison of our proposed method with the current state-of-the-artmethods, which however do not allow further deformation of the model. Therunning times of our algorithms indicate that there is only a small percentage2 M. Kamarianakis and G. Papagiannakisof performance load added when using Geometric Algebra representationforms instead of Euclidean ones to perform cuts/tears and drill holes.
5. Conclusions and Future Work
This work describes a novel way to perform model animation and deformationas well as cutting, tearing and drilling under a single geometric frameworkin Conformal Geometric Algebra. We focus towards a pure geometric-basedimplementation that can be applied to rigged models even in low-spec VRheadsets and ultimately enable real-time operations such as the ones presentedhere. Our current results were obtained using python but, since our goal isto have a full implementation in real-time virtual reality simulation, we willinevitably have to use more suitable programming languages and platformssuch as C
References [1] Manos N Kamarianakis and George Papagiannakis. Deform, Cut and Tear aSkinned Model Using Conformal Geometric Algebra.
CGI , 12221(3):434–446,2020.[2] Nadia Magnenat-Thalmann, Richard Laperrire, and Daniel Thalmann. Joint-dependent local deformations for hand animation and object grasping. In
InProceedings on Graphics interface’88 . Citeseer, 1988.[3] Marc Alexa. Linear combination of transformations.
ACM Trans. Graph. ,21(3):380–387, 2002.[4] Ben Kenwright. A beginners guide to dual-quaternions: What they are, howthey work, and how to use them for 3D character hierarchies. In
WSCG 2012- Conference Proceedings , pages 1–10. Newcastle University, United Kingdom,December 2012.[5] Ladislav Kavan, Steven Collins, Jiří Žára, and Carol O’Sullivan. Geometricskinning with approximate dual quaternion blending. dl.acm.org , 27(4), October2008.[6] Young Beom Kim and Jung Hyun Han. Bulging-free dual quaternion skinning.In
Computer Animation and Virtual Worlds , pages 321–329. Korea University,Seoul, South Korea, John Wiley & Sons, Ltd, January 2014.[7] D Hildenbrand.
Foundations of geometric algebra computing, 2013 . Springer.[8] Leo Dorst, Daniel Fontijne, and Stephen Mann.
Geometric algebra for computerscience - an object-oriented approach to geometry.
The Morgan Kaufmann seriesin computer graphics, 2007. geometric algorithm for cutting, tearing, drilling deformable models 23 [9] Rich Wareham, Jonathan Cameron, and Joan Lasenby. Applications of Con-formal Geometric Algebra in Computer Vision and Graphics.
IWMM/GIAE ,3519(1):329–349, 2004.[10] George Papagiannakis. Geometric algebra rotors for skinned character animationblending. In
SIGGRAPH Asia 2013 Technical Briefs, SA 2013 , December 2013.[11] Margarita Papaefthymiou, Dietmar Hildenbrand, and George Papagiannakis.An inclusive Conformal Geometric Algebra GPU animation interpolation anddeformation algorithm.
The Visual Computer , 32(6-8):751–759, June 2016.[12] H Hadfield and J Lasenby. Direct Linear Interpolation of Geometric Objects inConformal Geometric Algebra.
Advances in Applied Clifford Algebras , 2019.[13] Matthias Müller, Nuttapong Chentanez, and Miles Macklin. Simulating visualgeometry. In
Proceedings - Motion in Games 2016: 9th International Conferenceon Motion in Games, MIG 2016 , pages 31–38, 2016.[14] C D Bruyns, S Senger, A Menon, K Montgomery, S Wildermuth, and R Boyle.A survey of interactive mesh-cutting techniques and a new method for imple-menting generalized interactive mesh cutting using virtual tools‡.
The Journalof Visualization and Computer Animation , 13(1):21–42, February 2002.[15] Jun Wu, Rüdiger Westermann, and Christian Dick. A Survey of PhysicallyBased Simulation of Cuts in Deformable Bodies.
Computer Graphics Forum ,34(6):161–187, September 2015.[16] D Bielser, P Glardon, M Teschner, and M Gross. A state machine for real-timecutting of tetrahedral meshes. In , pages 377–386. IEEE Comput. Soc, 2004.[17] Andrew B Mor and Takeo Kanade. Modifying Soft Tissue Models: Progres-sive Cutting with Minimal New Element Creation. In
Advances in ComputerGraphics , pages 598–607. Springer Berlin Heidelberg, Berlin, Heidelberg, 2000.[18] Cynthia D Bruyns and Steven Senger. Interactive cutting of 3D surface meshes.
Computers & Graphics , 25(4):635–642, August 2001.[19] Daniel Bielser, Volker A Maiwald, and Markus H Gross. Interactive Cutsthrough 3-Dimensional Soft Tissue.
Computer Graphics Forum , 18(3):31–38,1999.[20] Jan Bender, Matthias Müller, Miguel A Otaduy, Matthias Teschner, and MilesMacklin. A survey on position-based simulation methods in computer graphics.
Computer Graphics Forum , 33(6):228–251, September 2014.[21] Iago U Berndt, Rafael P Torchelsen, and Anderson Maciel. Efficient Surgi-cal Cutting with Position-Based Dynamics.
IEEE Computer Graphics andApplications , 37(3):24–31, 2017.[22] Rich Wareham, Jonathan Cameron, and Joan Lasenby. Applications of confor-mal geometric algebra in computer vision and graphics. In
Computer algebraand geometric algebra with applications , pages 329–349. Springer, 2004.[23] Leo Dorst and Robert Valkenburg. Square root and logarithm of rotors in 3dconformal geometric algebra using polar decomposition, 2011.[24] Pablo Colapinto.
Articulating Space: Geometric Algebra for Parametric Design–Symmetry, Kinematics, and Curvature . PhD thesis, UC Santa Barbara, 2015.[25] Rich Wareham.
Computer Graphics Using Conformal Geometric Algebra . PhDthesis, University of Cambridge, 2007. [26] Xiufen Ye, Jianguo Zhang, and Yanyang Gu. An improved collision detection andcutting algorithm of the soft tissue in virtual surgical simulation.
InternationalJournal of Mechatronics and Automation , 4(4):236–247, 2014.[27] Daniel Wang, Yuru Zhang, Yuhui Wang, Y-S Lee, Peijun Lu, and Yong Wang.Cutting on triangle mesh: local model-based haptic display for dental prepa-ration surgery simulation.
IEEE Transactions on Visualization and ComputerGraphics , 11(6):671–683, 2005.[28] Zhongping Ji, Ligang Liu, Zhonggui Chen, and Guojin Wang. Easy mesh cutting.In
Computer Graphics Forum , volume 25, pages 283–291. Wiley Online Library,2006.[29] Xiufen Ye, Xi Ji’er, Ling Zhu, and Rui Yan. Research on soft tissue deformationand cutting in the virtual surgery. In
The 2011 IEEE/ICME InternationalConference on Complex Medical Engineering , pages 340–345. IEEE, 2011.[30] Rich Wareham and Joan Lasenby. Bone glow: An improved method for theassignment of weights for mesh deformation. In
International Conference onArticulated Motion and Deformable Objects , pages 63–71. Springer, 2008.[31] Hugo Hadfield, Dietmar Hildenbrand, and Alex Arsenovic. Gajit: SymbolicOptimisation and JIT Compilation of Geometric Algebra in Python withGAALOP and Numba. In
Advances in Computer Graphics , pages 499–510.Springer, 2019.Dr. Manos KamarianakisDepartment of Mathematics & Applied Mathematics,University of Crete,Voutes Campus, 70013 Heraklion, GreeceOrchid ID: 0000-0001-6577-0354e-mail: [email protected]