### of 3D Shapes Based on Modified Sphere Equation

##### Computer Graphics Group

Technical University of Szczecin

##### Szczecin / Poland

 1.    2.    3.    4.    5.    6.    7.

1. Introduction

The blob term is used with reference to 3D models with smooth shapes. This paper discusses a new method of obtaining models which are generalization of classical blobs. The method has been called GSE (Generalized Sphere Equation). The GSE models differ significantly, depending on parameters which have influence on their shape, size and specific surface bumps. The simplest variant of GSE method allows to obtain the classical blobs. This paper discusses the input data, mathematical interpretation and solutions which allow to render pictures of GSE models using ray tracing method.

Blob models are described as three-dimensional surfaces. There are two classifications of such surfaces: parametric and implicit.

Parametric surface is a set of points generated by three functions of given number of variables. For example, bivariate parametric surfaces are generated by functions of two variables: fx(i,j), fy(i,j), fz(i,j).

Implicit surface is a set of points which satisfy some equation in the three-dimensional space:

F(x,y,z)=0.

Implicit surface whose equation is of the second degree has been called quadric surfaces. (This term is also used with curves.) Thus, surface equation is:

Depending on values of parameters A, B, …, J equation 1 can describe plane, sphere, cone, cylinder, and the like.

Next chapter examines how the surfaces allied to quadrics have been used to making blobs. Section 3 includes the input data and the mathematical interpretation of the GSE method. In section 4 the usage of GSE method in ray tracing is discussed. There is description of implementation and results in section 5. The last chapter includes consclusions and discussions of future work.

2.1 J. F. Blinn's method

The first blobs were generated in 1982, when James F. Blinn used a method allied to quadric surfaces to solve the problem of computer aided visualization of molecular model.

In [1], a physical interpretation is given for mathematical solutions used. According to quantum mechanics, the electron in an atom can be represented as a density function of the spatial location. For a hydrogen atom, density function is:

where

,

(x1, x2, x3) is the center of an atom.

According to superposition theorem, density function can be used for many atoms. Thus the sum of density contributions of each atom should be taken into account:

where ri is the distance from (x, y, z) to the center of atom i.

A molecular surface can be defined as a set of points where density function equals some threshold value T. Using implicit surface definition a molecular surface can be described as:

In above equation bi has been specified in term of a blobbines parameter Bi.

[1] describes how to obtain the pictures of Blinn's blobs using ray tracer. Deriving ray equation begins with the various transformations into a standard viewing space. Reciprocal transformations let calculate viewing ray equation for each pixel on the screen. For each point of the ray, coordinates x and y depend on the z coordinate. Thus, calculating the z coordinate is enough to get the ray and surface intersection. J. F. Blinn has divided z value calculating into two phases:

1.      Root isolation phase: finding the range in which is a solution. In Blinn's solution for each viewing ray, a list of n values: z1, z2 ... zn-1 is made. The zi value is a z coordinate of this point on the ray in which atom i has the biggest influence Di on density function D. The list of z0, z1 ... zn-1 values is sorted in ascending order. The list is searched for the first value zi for which Di(zi) is greater than or equals T. Thus, <zi-1, zi> is a required range for numeric methods of finding zeros.

2.      Root refinement phase: finding a solution in a given range (using general numeric methods of finding zeros, like Newton method or regula falsi method). Blinn suggested using combination of two numeric methods of finding zeros: Newton method and regula falsi method.

The surface normal at a given point can be found by taking the gradient of the surface defining function, F:

Described method in previous chapter works well for small number of atoms (a few dozen). Yet it is much to slow when number of atoms is in the order of a few thousands. That is why the algorithm has been optimized.

The idea is based on the fact that only a part of atoms have influence on surface and ray intersection. Thus, the calculation can be limited to these atoms which are close enough to the viewing ray. Each atom could then be enclosed in a sphere. The enclosing spheres of all atoms are projected into screen space. Thus, it is possible to find atoms which can have influence on the color of each pixel.

2.2 Other solutions for Blinn’s method

The solutions of optimizing the algorithm of ray tracing implicit surfaces are described in [3].

While looking for intersection of viewing ray and implicit surface it is usually not possible to use an analytical method. That is because surface functions are very complicated as it was in case of the function H in J. F. Blinn’s method.

The solution is to use a numeric analysis. Usually, modifications or combinations of classical numeric methods are used, preceded by necessary transformations. Authors follow J. F. Blinn and divide calculations into root isolation phase and root refinement phase. Next two subsections describe two methods of optimizing the root isolation phase.

One of the methods of optimizing the root isolation phase is to approximate the complicated implicit function with a polynomial. In 1986, Geoff Wyvill suggested the approximation by the polynomial function C(r2):

The received models have been called soft objects. They are described in [2].

[5] describes how to use method taken from numerical analysis in root isolation phase. The method is called interval analysis and originated in 1966 [6].

To start the interval analysis, a function and a bracket of arguments must be given. The method narrows given bracket down to new bracket in which function is monotonic and has zero.

Let’s start description of the interval analysis with some definitions.

Def. An interval [a, b] is an ordered pair ab representing the range of numbers {x:axb}.

Several operations on intervals were defined: addition, subtraction, multiplication, division, squaring, exponentation. They allowed to calculate the value of Blinn’s surface function when the argument is a given interval.

Interval analysis consists in the following algorithm:

1.      Calculating the interval [s0, s1]=H[t0, t1].

2.      If 0Ï[ s0, s1], it means that in the bracket <t0, t1> function doesn’t have zeros. The algorithm is finished with the adequate communicate.

3.      Calculating the interval: [r0, r1]=H[t0, t1].

4.      If 0Ï[r0, r1], it means that the function f is monotonic in bracket <t0; t1>. The algorithm is finished and the <t0t1> bracket is returned.

5.      Back to first step for intervals: [t0, t0+½(t0+t1)] and [t0+½(t0+t1), t1].

Realization of above algorithm gives monotonic segment of the function, so then the root refinement phase can start.

The idea of GSE method is based on observation that sphere is a set of points which distance from one point (center) is constant whereas ellipsoid is a set of points which sum of distances from two points is constant. The idea is to generalize this principle for three, four or more points. Thus, the generalized sphere equation is:

where

,

(x, y, z) is the point on the surface,

(xi, yi, zi) is the ith equivalent of the center.

We assume that the models made using GSE method meet the following conditions:

1)      The model surface should be determined by a set of points,

2)      The model should have soft shape,

3)      The gravity points should have only a local influence on the shape of model (to make modeling intuitive).

The models described by equation 3 meet the conditions first and second, but the condition third which concerns local influence of the gravity point on the model shape is not met. So, the GSE model equation has been expanded to:

where

Derivation of the equation 4 can be found under 3.2 . The accurate meaning of f(di), g(di) and  φi is explained in 3.1. In case when n=1, φ0=1, f(di)=1, g(di)=0 we result in a sphere equation (or circle on the plane) where the radius equals d:

3.1 Input Values

Every model is described by:

§         set of  n (n≥1) gravity points {(x0, y0, z0), ..., (xi, yi, zi), ..., (xn-1, yn-1, zn-1), } in the 3D space,

§         set of weights {φ0, ... φi, ..., φn-1}, φi Ì<-1;1>,

§         d parameter determining size of the model,

§         pair of functions f and g describing the shape and the bumps of the model, respectively.

The gravity points describe shape of the model. Every intersection of the model is a closed curve of soft shape or set of such curves (fig. 1). The weights describe influence of every single point on shape of the model. Weights range from –1 to 1, both inclusive. In fig. 2 there are model intersections with the gravity points marked. The d parameter determines size of the model. The larger the model, the smaller the influence of each gravity point (fig. 3). The functions f and g are described in section 3.1.1.

 Figure 2: The influence of the weights assigned to the gravity points on the model shape (intersections). The size of each dot is proportional to the absolute value of the weight, while the color means sign of the weight: red – positive, blue – negative

 Figure 3: The influence of the d parameter on the model size and shape (intersections)

The function f describes the overall shape of the model. It’s arguments are d0, d1, d2 ... dn-1. di stands for the distances between a given surface point and the ith gravity point. For the model to have expected shape, the function f should have smooth graph. Example equations that meet both conditions are shown below. Moreover, the function f should assign lower values to the further gravity points. (dmax stands for the distance between given surface point and the furthest gravity point).

 a) b) c)

 a) for k=2 b) for A=-3 c) Figure 4: The influence of the function f on the model shape (intersections)

The function g describes the surface bumps. In previous examples, the function g was constant and equal to 0. In fig. 5 an example of a model with non-constant function g is shown. The intersection of the same model for g=0 is marked with red.

 Figure 5: Intersection (on left) and draft 3D image of model with non-constant function

Let’s exchange the sum on the left side of equation (1) with the weighted average:

Coefficients w(d0)... w(dn-1) should meet the following conditions:

1)      The coefficients should be chosen so that the gravity points located closer to a given point of the surface have more influence on it’s location than further gravity points.

2)      The coefficients should depend on weights φi assigned to each gravity point.

3)      The coefficients should allow modifying surface bumps described by the function g.

The first condition means that coefficients should have the greatest values for the nearest gravity points and the smallest values for the furthest ones. An example solution is the inverted square of the distance:

Let’s generalize and replace the inverted square of the distance with any function f (see section 3.1) that meets the first condition:

To meet the second condition, value of each coefficient should be determined by a weight assigned to each gravity point:

To meet the third condition, we need to take into consideration the surface bumps described by the function g (see section 3.1).

This chapter discusses two calculating method: finding the intersection of a ray with a GSE model and finding normal vectors. These two methods enable the GSE models to be used in ray tracing.

4.1.1 Problem Description

Finding the intersection of a ray with a GSE model can be realized using every method used for blobs. The method used in the GSE ray tracer is similar to the method described in [1], but seems to be simpler in implementation. In future it will probably be replaced with the interval analysis method.

Let dir be the normal vector determining the ray direction in ray tracing method. Then, to find the intersection of the ray with the GSE model, we need to solve the following equations:

The last three of the above equations are the parametrical description of a line in space. After inserting x, y and z, we result in:

The solution is to find the minimal positive value of the t parameter. Due to the complexity of equations 5, instead of searching an analytical solution method, a numerical method was used. As in [1] it has been divided in to two phases: root isolation phase and root refinement phase.

4.1.2 Root Isolation Phase

For a given t, let’s denote the equations 5 solution error with ε(t) (di as in equation 5):

The precise solution is the smallest of the error function’s zero points. Thus, we have to numerically find the smallest t where ε(t)≈0.

 Figure 6: Isolation of the first zero

Let’s consider n planes π0, π 1,... π n-1 crossing the appropriate gravity points and perpendicular to the direction of a given viewing ray. The ray intersects the π 0, π 1,... π n-1 planes in n points determined by the n values of the t parameter: t0, t1,... tn-1. We can calculate the intersections of ray with planes by combining the ray equation with the appropriate plane equation. The mathematical description can be found in [4].

From the geometrical interpretation we know that in intervals where the ray runs inside the model, the error function has positive values, while in intervals where the ray runs outside the model the error function has negative values.

Thus, we can deduce that the smallest zero point of the function can be located (see fig. 6 ):

§         in the <0; t0> bracket if ε(0)>0 and ε(t0)≥0,

§         in the <ti-1; ti> bracket if ε(t0), ε(t1), ... , ε(ti-1)<0 and ε(ti) ≥0.

To recap, the algorithm for root isolation phase can be divided in to the following steps:

1)        Finding the ti values for the points where the ray intersects the planes that include the successive gravity points.

2)        Sorting the ti values.

3)        Checking if the following condition is met for successive ti:

If yes, we can stop checking and use the secant method in the following interval:

§           <T1;T2>=<ti-1; t i> if i>0, or

§           <T1;T2>=< t 0-r; t 0> if i=0.

4.1.3 Root Refinement Phase

After finishing the above algorithm, we result in an interval <T1;T2> where ε(T1)<0 and ε(T2)>0.

Now, we can find the function’s zero point using one of classical numeric method. The iterative secant method, which is especially well suited, can be found in [7].

This section discusses two methods of finding the normal vectors for a GSE models. The first method is quite quick but it works only if g(di)=0. The second method is slow but it works also if g(di)≠0.

4.2.1 Method Based on the Dummy Centers

If assume that g=0 the method of finding the normal vectors can be based on the idea of "dummy centers" defined in the following way:

 Def. For a given surface point P(x, y, z) of a solid, a dummy center is a C point such so C-P is a normal vector beginning in P.

For the solids described here, a dummy center for any surface point P is a weighted average of the gravity points:

For a sphere, after inserting n=1, φ0=1, f(di)=1, we result in the geometrical center coordinates:

Fig. 7 a) shows results of the described method. Normal vectors are marked with blue lines and the dummy centers with green dots.

 Figure 7: GSE normal vectors calculated using the dummy centers (on left) and the gradient (on right)

4.2.2 Method Based on the Gradient

In case when the g≠0 the normal vector can not be found using dummy centers method. However every implicit surface normal can be found by taking the gradient of the surface defining function. Thus, according to equation 3 we have:

For example the x component will be

,

where

The result for the GSE model intersection is shown in fig. 7 b).

Two calculating methods described in section 4 has been implemented in ray tracer supporting GSE models. The GSE ray tracer  has been written in the C++ Builder and exist as a standalone  application (see screenshot of the ray tracer in fig.8). User interface allows changing the following input data:

§         position and direction of the camera,

§         light position,

§         activation/deactivation of  antialiasing,

§         selection of example scene,

§         position of solids (GSE models and spheres) can be seen in 2D projections.

On the output the program generates picture of scene composed of spheres and GSE models. Sample results are showed in fig. 9.

 Figure 8: Editing with the GSE models ray tracer

 Figure 9: Pictures of GSE models generated by raytracer

To test the ray tracer implementation, an example GSE (see fig. 10) shape was rendered in different light condition and the results were evaluated. We did not notice any deformation or other faults.

 Figure 10: The same test model rendered for six different light positions

Moreover the correctness of GSE shapes were conformed by comparison of rendered pictures and intersections of the same GSE models (fig. 11).

 Figure 11: GSE shape intersection and picture generated by the raytracer

Tests showed that finding the intersection of a ray with a GSE model and finding normal vectors give expected results.

Execution time has been measured for pictures in resolution 400x400. The results are presented in tab.1. Achieved results let treat calculating methods as quite efficient.

 number of gravity points time for AMD K6 300 MHz time for AMD Duron 600 MHz 2 9 s 3 s 3 11 s 4 s 6 24 s 9 s 12 53 s 22 s

Presented results show that GSE method allows getting models which resemble Blinn's blobs. According to expectations, the method lets control the shape of model and modify surface bumps.

The main future work is to ray trace the GSE models with surface bumps (g≠0). In the figure 12 draft pictures of such models are presented.

Another future goal is to evaluate the possibilities of efficient calculation of texture coordinates. It would allow to map textures on the GSE shapes.  Performance optimization are also taken into consideration.