Haptic Rendering of Thin, Deformable Objects with Spatially Varying Stiffness
HHaptic Rendering of Thin, Deformable Objectswith Spatially Varying Stiffness
Priyadarshini Kumari and Subhasis Chaudhuri
Vision and Image Processing Laboratory, Department of Electrical Engineering.Indian Institute of Technology Bombay, Powai, Mumbai-400076 { priydarshini,sc } @ee.iitb.ac.in Abstract.
In the real world, we often come across soft objects havingspatially varying stiffness, such as human palm or a wart on the skin.In this paper, we propose a novel approach to render thin, deformableobjects having spatially varying stiffness (inhomogeneous material). Weuse the classical Kirchhoff thin plate theory to compute the deformation.In general, the physics-based rendering of an arbitrary 3D surface iscomplex and time-consuming. Therefore, we approximate the 3D surfacelocally by a 2D plane using an area-preserving mapping technique - Gall-Peters mapping. Once the deformation is computed by solving a fourth-order partial differential equation, we project the points back onto theoriginal object for proper haptic rendering. The method was validatedthrough user experiments and was found to be realistic.
Keywords:
Haptic rendering · deformable thin surface · spatially vary-ing stiffness · Kirchhoffs thin plate theory · Gall-Peters projection
Haptic perception becomes more realistic when we incorporate the physical andmaterial properties of an object. For the kinesthetic rendering of rigid objects, itis sufficient to calculate only the force as they do not undergo any deformationunder the applied force. However, the rendering of deformable objects is morecomplex since we need to calculate the force as well as the deformation. If thestiffness of the object varies over the surface, there is further complications asthe rendered force is a function of both local deformation and local variation instiffness. The continuity of deformation and a smooth movement of the proxy onthe deforming surface are some of the key concerns while rendering a deformableobject having a variable stiffness.In the last two decades, haptic rendering has found several applications oninteractions with the virtual environment, such as virtual museum[17], dentalsurgery simulation[10], and force modeling for needle insertion[14]. However,these interactions cannot handle skin or bowel simulator, which requires theobject to be both deformable and have variable stiffness. The proposed methodtries to achieve this special type of rendering. a r X i v : . [ c s . G R ] M a y Priyadarshini Kumari and Subhasis Chaudhuri
We propose a method to render a point cloud data representing thin de-formable objects having spatially varying stiffness such as skin and elastic sheets.First, the point cloud in the neighborhood of the proxy is locally fit into a hemi-spherical surface and using the parameter of the fitted sphere these points arethen projected on a plane using Gall-Peters projection[7]. This being an area-preserving mapping, preserves the stiffness variation on the projected 2D plane.Subsequently, we use the Kirchhoff’s thin plate theory[18] to calculate the de-formation due to an applied force at the proxy and for a given stiffness pattern.The solution is obtained by solving a fourth-order non-homogeneous Partial Dif-ferential Equation (PDE). Once the deformation map is estimated, it is backprojected on the original 3D surface for proper rendering.Most of the work done in the field of deformable object rendering can beclassified into two categories: geometry-based deformation and physics-baseddeformation. Geometry based models are fast but do not necessarily focus onthe physics involved in the deformation. On the other hand, physics-based ren-dering model is computationally expensive as the simulation of physical andmaterial properties of the object is complex. In geometry-based modeling mostof the methods use Gaussian[12] or a polynomial function to move the surround-ing vertices around Haptic Interaction Point (HIP). Another widely practicedmethod is spline-based in which control points are assigned on the object, andthese control points are manipulated to effect a smooth deformation [5], [11].A FEM based method, such as [3] provides the most accurate rendering. How-ever, this is computationally very demanding for haptic applications, requiringspecialized hardware. It has been shown in [9] that by sacrificing the accuracymarginally, the computations can be speeded up by a factor of 10 by using amass-spring model. However, the choice of dampers is very critical to ascertainthe desired dynamic behavior, and usually, they are chosen heuristically. Thepurpose of the proposed method is to avoid these heuristics by solving a properPDE over a 2D plate.In smoothed-particle-hydrodynamics (SPH) approach[4], authors have intro-duced a technique that renders the object in different possible states (fluid,elastic, and rigid) in the same scene. However, in this technique, the stiffness ofthe entire object is varied as opposed to spatially varying the stiffness. Moreoverin [4], a volume model of the object is used, while we use a 3D surface modelwhich is comparatively much faster than the volume model.The key contribution of our work is the development of an alternate, stable,physics-based haptic rendering method for an inhomogeneous, elastic deformableobject.
The input for our method is the point cloud representation of the 3D object.Since there is a one-to-one relationship between the point cloud and mesh data[6](albeit the rendering process is a bit different), the algorithm can be extendedto other forms of data representation. The point cloud is given by 4-tuple point aptic Rendering under Stiffness Variation 3
Fig. 1: Various steps involved in haptic rendering of a deformable object withvariable stiffness. P i = { ( x i , y i , z i ) , d i } , where ( x i , y i , z i ) is the location of the i th point and d i is the corresponding stiffness. The overall algorithm can be broken into severalstages as shown in Fig. 1. The various stages of algorithm are described in thefollowing subsections. With point cloud data, we do not have surface normal defined at a point. Inorder to provide the appropriate force feedback to the user, we need the normalinformation at the collision point. We follow a proxy-based technique proposedin [16] for computation of normals. In this method, the normal is computed fromthe radial overshoot of each point inside the spherical proxy. Once the collisionoccurs, HIP penetrates the object, and the proxy is constrained to lie on thesurface by using a dynamic function as defined in [16]. A collision is detectedwhen the condition ( v n . v h ) < v n is the computed surfacenormal at the point of collision and v h is the vector joining the proxy to HIP asshown in Fig. 2. Once the collision point is detected, we select the point cloud v n v h v h v n α < ◦ α > ◦ No collision Collision
Fig. 2: Illustration of collision detection. The proxy and the HIP are shown ingreen and red colors, respectively. The point cloud is shown in blue. Blue arrowsshow the computed outward surface normals.in the neighborhood of the point of collision for fitting a sphere locally.
We use the algebraic distance-based method proposed in [15] to fit the surfacepoints in the neighborhood of the proxy. Let f ( x, y, z ) be an equation of spherein the algebraic form Priyadarshini Kumari and Subhasis Chaudhuri f ( x, y, z ) = u ( x + y + z ) + u x + u y + u z + u = 0subject to u + u + u − u u = 1 (1)where u = [ u u u u u ] T represents the coefficients describing the sphere.The above mentioned constraint can be written in the matrix form as u T Cu = 1with C as given in (2). So any point ( x i , y i , z i ) lying on the surface of spherewill have an algebraic distance equal to zero. We minimize the sum of squaredalgebraic distances to estimate the best possible spherical fit u . Although severalpossible constraints have been suggested in [2],[1], we use the constraint u T Cu =1, as Pratt [15] has shown that such a constraint has only one point of singularityof being a zero radius sphere. The resulting cost function can be expressed byintroducing a Lagrangian multiplier λ as:minimize u ( (cid:107) Du (cid:107) − λ ( u T Cu − D = x + y + z x y z x n + y n + z n x n y n z n , C = −
20 1 0 0 00 0 1 0 00 0 0 1 0 − where D is an n × n points from the point cloud data,selected in the neighborhood of the collision point.The solution is given by solving the generalized eigenvalue problem D T Du = λ Cu , which can be solved using QZ algorithm developed by Moler and Stew-art [13]. The center and the radius of the sphere can be calculated from thecomponents of u as given in [8] c = − u [ u u u ] , r = (cid:114) c T c − u u . (3)Once we get the center c and the radius r of the sphere, all the points are locallyprojected on to a 2D plane as discussed in the next subsection.All points, however, need not lie on the surface of the best fit sphere, as shownin Fig. 3a. In order to approximate the object surface by an area-preservingplanar patch, we need to project all points on the surface of the fitted sphere.We use a projection technique in which we join all the points (which are notlying on the surface of the sphere) from the center. The joining line intersectsthe sphere at two points. The intersecting point nearest to the object point ischosen to be the projection of the corresponding point on the sphere. Figure 3bshows the best fit sphere to a set of object points near the proxy.As it was mentioned earlier that deformation due to an arbitrary variationin material properties is easier to compute for a thin sheet, the fitted sphereis mapped onto a 2D plane. For an equi-curvature surface (like a sphere), theamount of energy required for deformation is proportional to the surface area aptic Rendering under Stiffness Variation 5 P i P s A B D C O ProxyEastern hemisphereWestern hemisphere (a) (b)
Fig. 3: (a) Illustration of projection of points on the surface of the best fit sphereand employing the constraint of hemisphere fitting. The object points P i areshown in red, and the projected points P s are shown in blue. (b)Illustrationof the best fit sphere to the set of object points near the proxy. The proxy isshown with a red ball, and the green line points to the direction of normal at thecollision point. (Data Courtesy: http://graphics.stanford.edu/data/3Dscanrep/)when the stiffness is constant. Hence we require a projection operator g : IR → IR , which is area-preserving so that deformation computation is a function ofstiffness alone. Gall-Peters projection [7] used in cartography is one such operatorthat preserves the area. Hence we use Gall-Peters projection in a way so that theproxy position lies at the center of the unfolded globe.While fitting the sphere,if there are points on the western hemisphere, these points are not visible at theproxy, and hence they are rejected, and only a hemispherical facet is considered.This is illustrated in Fig. 3a, where D is the proxy on the surface, AB correspondsto the plane joining the north pole ( A ) and the south pole ( B ) and the point D corresponds to (0 ,
0) latitude and longitude of the sphere. Points lying on thewestern hemisphere (shown as the arc
ACB ) are discarded.
Since Gall-Peters projection describes a technique to project the globe onto aplane, we have used it for a sphere to plane projection. The longitude and latitudeof each point on the sphere is calculated from polar coordinates of the points.Considering the z-axis to be perpendicular to the computer screen, each pointon the sphere is projected onto the x-z plane using (4). x = rα √ z = r √ sin ( β ) , (4)where α and β are longitude and latitude of the point in radians respectivelyand r is radius of the sphere. While using the Gall-Peters projection some issuesemerge, when points lie either near the pole region or near the 180 ◦ longitude. Priyadarshini Kumari and Subhasis Chaudhuri
1. A point lying on the north/south pole gets stretched to a line on the planeafter projection.2. Points on 180 ◦ longitude get mapped on both the extreme of the 2D plane.We avert such problem by shifting the coordinate axis of sphere so that ourproxy (point D in Fig. 3a) falls in the region of zero latitude and zero longitude.After shifting the axes, we discard the points on the other half of the sphereopposite to the proxy (arc ACB in Fig. 3a). In this way, all neighborhood pointsaround the proxy fall near the center of the plane, which would be very usefulwhile employing an infinite horizon boundary condition in section 2.5. Oncethe points are projected on a plane, we compute the deformation of each pointusing Kirchhoff thin plate theory [18]. As mentioned earlier, computation ofdeformation on the planar sheet is easier than directly computing on the objectsurface. Therefore we compute the deformation of each point on 2D plane andsubsequently project the points back on the original 3D object.
In Kirchhoff theory of plate [18], the thickness of the plate is assumed to be verysmall as compared to the planar dimensions. In addition, some more assumptionsare made to reduce the three dimensional problem to two dimensions. Theseassumptions are summarized as follow1. The line normal to the neutral axis remain straight after deformation.2. The normal stress in the direction of thickness is neglected.3. Thickness of the plate does not change during bending.Let s , w and t be the displacements of a point along x, y and z directions,respectively, and | F | is a distributed load on the planar sheet. Since the planarsheet is on the x-z plane, we can write the variation of s and t across thicknessin terms of displacement w as s = − y ∂w∂x , t = − y ∂w∂z .Taking all the assumptions into account, the deformation of a point is gov-erned by the following PDE as explained in [18]. ∂ w∂x + 2 ∂ w∂x ∂z + ∂ w∂z = | F | D , or ∇ w = | F | D , (5)where ∇ is the biharmonic operator and D is flexural rigidity of the object.Interestingly, this equation was first obtained by Lagrange in 1811. The flexuralrigidity is a material property of the object and is defined as D = Eh − v ) , (6)where E , h , v are the modulus of elasticity, thickness of the plate and Poisson’sratio (usually 0 < v < . D of the object is constant. In this case, computation of deformation aptic Rendering under Stiffness Variation 7 also becomes easier as the governing differential equation (5) is a standard bihar-monic equation. However in real world, there are numerous objects which havespatially varying material properties when E and v are functions of x and z .Equation (5) cannot be applied to compute deformation for such objects. Hencewe move to the next stage where stiffness of the object is allowed to vary overthe surface.In order to compute the deformation in variable stiffness object we use theextended Kirchhoff thin plate theory[18]. The deformation of each point in suchcases is governed by the following equation D ∇ w + 2 ∂D∂x ∂∂x ( ∇ w ) + 2 ∂D∂z ∂∂z ( ∇ w ) + ∇ D ( ∇ w ) − (1 − v ) (cid:18) ∂ D∂x ∂ w∂z − ∂ D∂x∂z ∂ w∂x∂z + ∂ D∂z ∂ w∂x (cid:19) = | F | . (7)The detail derivation of (5) and (7) is explained in [18]. Unfortunately finding ananalytical solution of this form of PDE is extremely difficult. Hence we solve theequation in discrete domain by using Jacobi iterative method for a given bound-ary condition. As the Gall-Peters projection uses non-linear mapping function(4), points are not projected evenly on a plane. Therefore we first sample theprojected points on the uniform grid and subsequently we discretize all the vari-ables and parameters of (7) by using central difference. The value of D at agrid node is taken as the value of D of the closest projected point from the gridnode. Our initial condition on deformation is w ( i, j ) = 0. The final expressionfor iterative updating can be written in the form of a function L (): w n +1 ( i, j ) = L ( | F | , D ( i, j ) , w n N ( i, j )) (8)where the superscript n denotes the n th iteration and w n N ( i, j ) refers to thevarious lattice entries in the neighborhood of the central lattice ( i, j ), while D ( i, j ) is the given flexural rigidity map and its all derivatives (upto secondorder) at the lattice ( i, j ). In each iteration, we update the w matrix and compareit with w matrix of previous iteration. We stop the iteration when the Frobeniusnorm of difference between current and previous w matrices becomes very small.The final w matrix gives the deformation at the grid point. However, one needsthe deformation value at the location where the 3D point was originally projectedusing (4). The deformations at projected points are computed using bi-linearinterpolation. As we mentioned earlier, in order to obtain the solution of the PDE given in(7), we need suitable boundary conditions. For simplicity, we assume edges ofthe projected plane to be parallel to the coordinate axis X and Z .We need two boundary conditions at each edge. In our work, we assume allfour edges of projected plate to be fixed (i.e deformation at infinite horizon is Priyadarshini Kumari and Subhasis Chaudhuri D ( x )=2 . x (a) D ( x ) = x (b) Fig. 4: Computed deformations on the object surface with different flexural rigid-ity functions under the same amount of applied force: (a) Linear, (b) Hyperbolic( x in cm). Deformation at various locations are superimposed together for easeof visualization.zero). Hence the deflection and the first order derivative become zero at all fouredges. After getting the deformation locally on the projected plane, we projectback the points to compute the deformation on the original object for properrendering. Points on the object are displaced in the direction of the force F usingthe governing equation (9) (9)where ˆ P i is the final position of the i th point on the deformed surface and P i corresponds to the position before deformation with w i being the deformationat the point. In order to haptically render the object, we first detect the collision of HIPwith the object as discussed in section 2.1 and compute the force if a collisionis detected. Once the collision is detected, the force needs to be fed back bythe haptic device to provide the sensation of touch. Two important factors forforce computation are its magnitude and its direction. The magnitude of forceshould be proportional to the penetration depth of HIP from the surface and itsdirection should be in the direction of the normal at the point of collision. Hencethe reaction force is calculated using the expression, F = − EAh ( | X h − X p | ). Here X h is the HIP position, X p is the proxy position, A is the area of the plate and E is the modulus of elasticity at that location. The proposed method was implemented in Visual Studio 2010 in a Windows 7platform with a Core(TM2) Quad CPU Q8400 processor @ 2.66 GHz clock speedwith 8 GB RAM. We use OpenGL 2.0 for graphic rendering and HAPI libraryfor haptic rendering. We experimented with 3D point cloud model using Falconhaptic device from NOVINT. In order to make the rendering faster we segregate aptic Rendering under Stiffness Variation 9
Proxy
Fig. 5: Deformation shown on the bunny model.the tasks of deformation computation, proxy updation and force computationthrough three different threads while programming. We observed that the defor-mation thread, as expected, runs relatively slowly due to iterative solution forequation (7). The average time required to run the deformation thread once isaround 160 ms. We have set v = 0 . × ×
200 with inter-node space of 0.025cm. Figure 4 shows the deformations at several points onthe object surface with different flexural rigidity functions. In Fig. 4a, flexuralrigidity is varied gradually along the horizontal direction while in Fig. 4b it isvaried rapidly. We applied the same amount of force on the object surface andobserved that deformation decreases as flexural rigidity increases as shown in Fig.4a. It may be noted that the deformation is restricted to the neighborhood onlydue to the choice of boundary conditions. To explain this further, we experimenton a real object model and show the deformation at the ear lobe of a bunny inFig. 5. One may expect the ear lobe to be bent due to the applied force like acantilever. Instead one observes a dip within the ear lobe, exact shape of whichwould depend on the size of the chosen neighborhood.We also verify the performance of the proposed algorithm in terms of userexperience while interacting with various object models. We set up an experimentin which users were asked to interact with an object model (say Fig. 5) and ratethe experience on a scale of 1 −
5, with 1 being very poor to 5 being excellent.Most of the users (80%) rated the experience as very good (rating 4) remaining20% gave the rating of 3. When asked if they are able to perceive the change instiffness, 90% users responded positively.
We have proposed a numerical method to solve the haptic rendering problem ofelastic deformable objects with spatially varying stiffness. It is a physics-based rendering technique and gives a realistic feeling of touching the object. We havenot observed any divergence issue in the computation of the deformation mapor the proxy movement over the deforming surface. In order to mathematicallyhandle any object, an area-preserving map projection technique is used to trans-form an arbitrary 3D surface to a planar one following which Kirchhoff’s thinplate deformation propagation model involving a fourth-order PDE is used tocompute deformation in case of spatially varying stiffness. We also tested withsome subjects and observed that the hapto-visual rendering technique using theproposed algorithm greatly augmented the user’s experience. Our future workwill involve the choice of alternate boundary conditions and extension to dealingwith volumetric data.