Direct Computation of Skeleton Points for 3D CAD Models

 

D. Brujic, M. Ristic, J. G. Lambourne
Imperial College London, Exhibition Road, London SW7 2AZ Z, England
d.brujic@imperial.ac.uk

Z. Djuric
Silvaco , Silvaco Technology Centre  St.Ives, Cambridge PE27 5JL, England

Contents


Abstracts:

We present an algorithm for computing the skeleton points of a trimmed Non-Uniform Rational B-Spline (NURBS) model, by finding the centre and the radius of maximal spheres contained within the modelled object.  This information can then be visualised as a colour coded thickness map on the object's surface, allowing intuitive interpretation by designers and engineers. By making use of surface normal information the task of finding maximal spheres can be reduced to a simple minimisation problem, avoiding the calculation of Delaunay tetrahedra or Voronoi diagrams. Also our method allows the skeleton points to be evaluated one point at a time, so computation can be focused on a localised region of the model.  The simplicity of this method makes implementation straightforward allowing it to be rapidly integrated into existing CAD/CAM packages.

 

Key words: Medial axis, Skeleton, NURBS

 

1. Introduction

 

Thickness information is intuitively useful when considering the mechanical, physical, electrical and thermal properties of solid objects. As such the rapid calculation and visualisation of thickness information would be of great use to designers and engineers.  Thickness is of particular interest when designing parts intended to be manufactured by plastic injection molding.  During the injection molding process thin regions can restrict the flow of molten plastic [13], causing the mold to be only partially filled. Although this process can be accurately simulated [24], experienced designers often have a set of heuristic rules, based on cavity thickness, which can be used at the design stage to cut out the first iteration on the design/analysis cycle.

For an arbitrary 3D solid, the thickness related to a point on the surface is defined as the radius of the largest sphere which can be placed inside the object, touching that point [14]. Thickness is related to the well known Medial Axis Transform (MAT) which is the locus of the centres of all maximal spheres, along with their corresponding radii [5].  Over the years a number of approaches have been developed for computing medial axis transforms.  These include voxel based techniques [23], algorithms using distance and vector distance functions [11, 8] and tracing methods [21].  When working with CAD models almost all research has focused on the use of Voronoi diagrams and Delaunay tetrahedral [6, 1, 4, 19, 20, 18].   In this paper we present an algorithm which, given a point on a parametric surface, returns the radius of the maximal sphere which touches that point tangentially.  A second touching point for the same maximal sphere, useful for a variety of automatic meshing applications [4, 18], is also recovered.  The algorithm can be used to find the thickness associated with points individually, which is not possible using methods based on Delaunay tetrahedralisation. Further, when evaluated for a mesh of points covering the object, the thickness can be represented graphically, either as a colour coded thickness map on the objects surface or in the form of the object's medial axis transform.  Ultimately our goal is to construct a parametric representation of this information such as that described in [25], using the parameterisation of the original boundary surfaces as a starting point.

The precision of our method comes from its use of iterations on the objects surface.  Ang et. al [3] have also suggested a method which uses an iterative procedure on parametric surfaces to optimise the touching points of maximal spheres found using the conventional Delaunay tetrahedralisation approach. Their paper contains details of how differential geometry can be employed for the adaptive meshing of the medial axis.  However, our method is somewhat simpler than theirs, requiring the solution of a smaller system of equations and providing its own initial guesses.

 

2. Tangential and Maximal Spheres

 

Let us consider a CAD model which defines a solid object . The model comprises of a set of trimmed NURBS surfaces which define the object's boundary, ∂Λ, and have their normal directions oriented such that they point into the inertia of Λ.  As defined above, the thickness of Λ at a point p on the objects surface is the radius of the largest sphere contained within Λ and touching p.  This will sometimes be called the maximal sphere associated with p.

Sherbrooke [22] has shown that any sphere contained within Λ must touch the boundary surface tangentially. For the algorithm described here, spheres which touch the boundary surface tangentially, even when not fully contained in Λ, are of great importance.  If we consider a sphere, touching the surface tangentially at the point p and touching one other point q  ∂Λ, this sphere will in general intersect with ∂Λ and so cannot be contained in Λ. The maximal sphere associated with p is simply the smallest sphere tangentially touching p and at least one other point on the boundary ∂Λ.

By definition the centre of a sphere, x, touching the surface tangentially at a point p, must lie somewhere on the line

where  is the unit surface normal at the point p. 

 

Fig. 1. A sphere touching the point p tangentially and also touching one other surface point q

 

When this sphere also touches a point q, as shown in figure 1, the value of R can be calculated by equating expressions for cos(α)

giving

The reader should notice that the term  can be negative, corresponding to spheres lying outside the object.  If we are concerned only with spheres contained within the object it is convenient to define

As q approaches p this equation become susceptible to small perturbations in the relative positions of p and q. Since CAD models are only watertight to a modelling tolerance ε, small patch overlaps can cause problems close to patch boundaries.  Consequentially R(q) should be set to infinity when the distance between p and q are smaller than the modelling tolerance, with the exception of sharp convex corners where R(q) should be set to zero.

The task of finding maximal spheres can now be restated as a simple minimisation problem.  For a set of parametric surfaces ∂Λ and the point p, we wish to find the point q* for which equations 3 and 4 give the smallest value of R(q) of any q  ∂Λ.  The thickness is then the radius of this smallest sphere.

The solution framework for this task is identical to that used for point projection [15, 12], where given an arbitrary `datum point' in space, the closest point on the surface must be found.  Algorithms for point projection are well established in the literature and implemented as an integral part of most software which deals with CAD models.  Only two adaptations are needed to convert the algorithm described in [2] to find the point q*.  First, the metric used when searching for an initial guess [7] must be changed as described in section 3.  Second, the `cost function' to be minimised [17], must be changed to the criteria stated in equations 3 and 4.  This methodology is now described in detail.

 

3. A simple algorithm for the computation of maximal spheres

 

As a first step the surface ∂Λ must be sampled with a set of points q0, q1, ... qn.  This sampling will be referred to as the `search set' for the continuation of this paper.  The density of these sample point should be such that the greatest distance between a sample qi and any of its neighbours along the surface is smaller than some density parameter ω [6].  Methods for generating such a sample set are described in [15, 10].  Brandt [ 6] has shown that for a sampling with density parameter ω any region of the surface for which the associated maximal spheres have radii greater than ω is adequately sampled. In principle ω should be smaller than the minimal thickness present at any point on the object, but sampling using this criteria is often impossible as the thickness falls to zero at sharp edges.  In practice we find that sharp edges rarely cause problems and so this criteria can be relaxed.  A reasonable heuristic is to set ω small enough such that more than 90\% of the surface has associated maximal spheres with radii greater than ω.

To find the radius of the maximal sphere associated with a point p  ∂Λ the following procedure should be used.  For each point qi in the search set, calculate R(qi), using equations 3 and 4 and select the point qj for which the corresponding value of R(qj) is smallest. A 3D tree system [7] can be used to accelerate this process, avoiding extensive searching of the set and finding the initial guess in expected O(log(n)) time.  The parametric co-ordinates of qj are then used as the initial guess for an iterative procedure, used to minimise the function R(q), and so find the local thickness.

Numerous numerical methods exist for minimizing multi-dimensional functions.  Some of these, such as simplex [12], require only evaluation of the function to be minimized while others can also make use of the function's partial derivatives, often resulting in a massive decrease in the number of function evaluations necessary for convergence.  The partial derivatives of the surface with respect to u and v can easily be calculated for parametric surfaces [15] allowing the use of these faster algorithms.

Our implementation uses a bounded quasi-Newton method [17], which requires both function evaluations and the function's partial derivatives.  The partial derivatives of equation 3 can be calculated from the partial derivatives of the surface using

We note that the full Hessian matrix of surface derivatives can be easily calculated and so a standard unbounded Newton's method could be implemented.

The iteration returns the parametric coordinates of the point q* for which R(q) is minimum.  This point is considered important for a number of automated meshing procedures [4, 18].   Once the point q* has been found the thickness R(q*) can easily be evaluated using equation 3 and the centre of the maximal sphere calculated using equation 1.

 

4. Trimmed NURBS

 

Most CAD models consist of trimmed NURBS, in which some regions of the parametric domain are removed from the final surface, allowing holes and arbitrary boundary shapes to be defined.  The `trimmed' regions are demarcated by a set of curves in the parametric domain.  Winding rules or boundary crossing tests [16] can then be used to find whether a point in the parametric domain lies in one of these regions.  In order to implement our algorithm for trimmed NURBS two additional conditions must be enforced. First, no point in the search set should lie in a trimmed region of the NURBS.  This can easily be achieved by using the tests described in [16].  Second, if the iteration returns a point in a trimmed region the appropriate point on the trimming curves surrounding that region must be identified as q*. The procedure for this is identical to that described in section 3.  A set of sample points on the trimming curves [15] are searched and the sample returning the smallest value of R(q) is used as the initial guess for an iteration on the appropriate curve.

 

5. Results

 

The thickness calculated by the above algorithm can be visualised for an entire model in the following way.  Each surface in the set ∂Λ is triangulated using the methods suggested in [16, 10]. The thickness is then calculated for each vertex of the triangulation and converted into appropriate r, g, b values using a conventional colour scale.  The triangles can then be rendered using standard functions, for example those offered by OpenGL, which interpolate the colours defined at each vertex over the face of the triangle. An alternative way to visualise this information is as a coloured approximate medial axis of the object.  Equation 1 can be used to find the centres of the maximal spheres associated with each of the vertices in the surface triangulation. The sphere centres can then be used as the vertices of a new triangulation with the connectivity preserved from that of the original surface points.

Figures 2 to 4 show these representations for several trimmed NURBS models. This ‘double sided’ medial axis representation is also described in [9] where it is used for the detection of surface features.  Similar feature detection can also be performed using the medial surfaces obtained from our algorithm.

 

6. Conclusions

 

This paper describes a simple and efficient method to calculate the thickness of an object associated with any point on its surface. This is equivalent to finding the radius of the largest sphere which fits into the object and touches that surface point. The algorithms similarity to point projection allows it to be rapidly implemented in existing CAD/CAM software, often by making only small adaptations to the existing functionality.  When the thickness is calculated for a mesh of points covering the surface, an intuitive colour coded thickness map can be generated for the visualisation of thickness, as well as a medial axis type representation.

Reference

 

[1] Nina Amenta, Sunghee Choi, and Ravi Kolluri. The power crust, unions of balls, and the medial axis transform. Computational Geometry: Theory and Applications, 19(2-3):127–153, 2001.

[2] I. J. Anderson, D. A. Turner, and J. C. Mason. An efficient algorithm for solving the surface foot point problem. Information

Geometers, 1998.

[3] Pin Y. Ang and Cecil G. Armstrong. Adaptive curvature sensitive meshing of the medial axis. Engineering with computers,

18(3):253–264, 2002.

[4] C.G. Armstrong, D.J. Robinson, R.M. McKeag, T.S. Li, S.J. Bridgett, R.J. Donaghy, and C.A. McGleenan. Medials for

meshing and more. In Proceedings of 4th International Meshing Roundtable, pages 277–288, 1995.

[5] Harry Blum. Biological shape and visual science. Journal of Theoretical Biology, 38:205–287, 1973.

[6] Jonathan W. Brandt and Ralph V. Algazi. Continuous skeleton computation by voronoi diagram. Image Understanding,

55(3):329–338, 1992.

[7] Jerome H. Friedman. An algorithm for finding best matches in logarithmic expected time. ACM Transactions on Mathematical Software, 3(3):209–226, 1977.

[8] Jos Gomes and Olivier Faugeras. Reconciling distance functions and level sets. Journal of Visual Communication and

Image Representation, 11(2):209–223, 2000.

[9] M. Hisada, A. G. Belyaev, and T. L. Kunii. A skeleton-based approach for detection of perceptually salient features on

polygonal surfaces. Computer Graphics Forum, 21(4):689–700, 2002.

[10] William L. Luken. Tessellation of trimmed nurb surfaces. Computer Aided Geometric Design, 13(2):163–177, 1996.

[11] G. Malandain and S. Fernndez-Vidal. Euclidean skeletons. Image and Vision Computing, 16:317–327, 1998.

[12] J. A. Nelder and R. Mead. A simplex algorithm for function minimization. Computer Journal, 7:308–313, 1965.

[13] T.A. Osswald, L. Turng, and P.J. Gramann. Injection Molding Handbook. Carl Hanser Verlag, 2001.

[14] Sylvain Petitjean and Edmond Boyer. Regular and nonregular point sets : properties and reconstruction. Computational

Geometry, 19(2-3):101–126, Jul 2001.

[15] Les Piegl and Wayne Tiller. The NURBS book. Springer-Verlag, 1997.

[16] Les Piegl andWayne Tiller. Geometry-based triangulation of trimmed nurbs surfaces. Computer-Aided Design, 30(1):11–18, 1998.

[17] H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in C. Cambridge University Press

Geometers, 1992.

[18] Peter Sampl. Semi-structured mesh generation based on medial axis. In Proceedings, 9th International Meshing

Roundtable, pages 21–32, 2000.

[19] D. J. Sheehy, C. G. Armstrong, and D. J. Robinson. Computing the medial surface of a solid from a domain delaunay

triangulation. In Proceedings of the third ACM symposium on Solid modeling and applications, pages 201–212. ACM Press, 1995.

[20] Damian J. Sheehy, Cecil G. Armstrong, and Desmond J. Robinson. Shape description by medial surface construction.

IEEE transactions on visualisation and computer graphics, 2(1):1077–2626, 1996.

[21] E. Sherbrooke, N. Patrikalakis, and E. Brisson. Computation of the medial axis transform of 3-d polyhedra. In ACM symposium on Solid modeling and applications, pages 187–200, 1995.

[22] Evan C. Sherbrooke, Nicholas M. Patrikalakis, and F. E. Wolter. Differential and topological properties of medial

axis transforms. Graphical Models and Image Processing, 58(6):574–592, 1996.

[23] P. Soille. Morphological Image Analysis. Springer-Verlag, 2003.

[24] L. Xu, W. Xu, and Y. Chen. Plastic injection molding process optimisation using software tools. International Journal

of Vehicle Design, 25(1 / 2):53–63, 2001.

[25] P. Yushkevich, P. T. Fletcher, S. Joshi, A. Thall, and S.M. Pizer. Continuous medial representations for geometric object

modeling in 2d and 3d. Image and Vision Computing, 21(1):17–28, 2003.