Find Isosurface for Metaballs

We separate the ice particles and water particles, and use different ways to test inter-section for them. Here we introduce the method used for testing the interinter-section for water particles. The method is based on the method proposed by Kanamori et al. [KSN08], and they use a metaball to represent a water particle.

To render the shape of the metaballs smoothly, the isosurface of these metaballs can be found using ray tracing with the special principles of judgment. Each metaball has a density field, and a set of nearby metaballs represents a smooth surface by the way of combining their density fields. In the following sections, we will introduce the method to find the isosurface for the metaballs.

4.2.1 Field function of metaballs

Figure 4.3: A sketch of the field function f_i.

We use the density field function f_i to represent the metaball i based on the distance r from the central of metaball P_ito the point x, as shown in Figure 4.3. The f_i has a finite

support R_i, and R_irepresents the effective radius of metaball i. Thus, f_i(r) = 0 if r ≥ Ri. For N metaballs, the shape of the isosurface is defined by the point x∈ R³ in below:

f (x) =

N∑−1 i=0

q_if_i(∥x − Pi∥) − T = 0, (4.2)

where T is a threshold, and{qi} are the density coefficients. The normal vector at x can be derived from−∇f(x).

The Nishita and Nakamae's algorithm [NN94] converts the function f_iinto the Bézier curve form, and uses Bézier clipping to solve equation with quadratic and degree-six poly-nomial. Here we use the degree-six polynomial, and it has seven control points. The equation of Bézier curve is shown as below:

f_i(s_i) =

∑6 k=0

dⁱ_kB_k⁶(s_i), (4.3)

where s_i ∈ [0, 1] is an intersected interval, B represents the Bernstein polynomial, Bkⁿ(u) = (_n

k

)u^k(1− u)^n−k, and (k

6, dⁱ_k)(k = 0, 1, ..., 6) is the coordinates of the control points of Bézier curve. When the Bézier curve represents one metaball, the control point is set as following:

d₀ = d₁ = d₅ = d₆ = 0, d₂ = d₄ = 16

27a²_i, d₃ = (8ai+ 5)a²_i

45 , (4.4)

where a_i is the length of the intersected interval D

R²_i , and D is the discriminate from the formula of ray-metaball intersection test.

4.2.2 Framework for rendering metaballs

The isosurface of these metaballs will be found using a special ray tracing algorithm.

Figure 4.4 shows the sketch map for rendering metaballs. B₀and B₁are two nearby meta-balls, and R₀is a incident ray intersecting the surface of these two metaballs at{t0, ..., t₃}.

We detect the ray-isosurface intersection point from the first interval [t₀, t₁] to the last in-terval [t₂, t₃] using Bézier clipping. This process stops when an intersection is found or no intersection is found.

Figure 4.4: A sketch for rendering metaballs.

In each interval, the Bézier curves of the influential spheres are summed into one curve. For instance in the Figure 4.4, the interval [t₁, t₂] is effected by the B₀ and B₁ at the same time, so the Bézier curve f₀₁ in the interval [t₁, t₂] is constructed by combining two curves f0 and f1 which represents B0 and B1 respectively. It is simple to combine multiple Bézier curves, since all we need to do is add each control points dⁱ_kbelonging to each influential spheres; i.e.:

d⁰¹_k = d⁰_k+ d¹_k, (k = 0, ..., 6) (4.5)

4.2.3 Bézier clipping

Figure 4.5: A Bézier curve with five control points. (a) is the Bézier curve before clipping, and (b) is the curve after clipping.

We use Bézier clipping method [NN94] to iteratively compute the root f_i(s_i)−T = 0,

as shown in Figure 4.5. The Bézier curve is always inside the convex hull formed by its control points dⁱ_k. The two intersection points between convex hull and horizontal line s = T is s_min and s_max. Therefore, we find [s_min, s_max] iteratively and finally we can converge to the root.

In each new interval, we want to obtain new control points for new curve, as shown in Figure 4.5(b). The de Casteljau's Algorithm is utilized here. Also, we use the Jarvis March algorithm to obtain convex hull for curve. These two algorithms will be illustrated in the following section

4.2.4 The de Casteljau's algorithm

The de Casteljau's algorithm is used to obtain new control points for the new curve.

If we have the original control points d_k(k = 0, ..., n) and the new interval [s_min, s_max], we can use this algorithm to get the new control points d^′_k(k = 0, ..., n).

The process of the de Casteljau's Algorithm as below:

1. For any n^th order Bézier curve, this algorithm starts with the point set{pk= d_k}, which contains exactly n + 1 elements.

2. It then replaces this set with an interpolated set{p^′0 = point between p₀ and p₁ at distance t, p^′₁ = point between p₁ and p₂ at distance t, ..., p^′_n−1 = point between p_n−1 and p_n at distance t}, where t is the position we want to find on the original interval. This set will be 1 element shorter than p_k.

3. Step 2 is repeated until there is only one point left, and that point will lie on the original curve, at the same t value.

In this process, the new control points are found at interval [0, t] or [t, 1]. An example is shown on Figure 4.6, assuming point set p_k(k = 3).

For the new interval [s_min, s_max], we run this algorithm two times. In the first time, we set the s_min as t, and find the new control points at interval [t, 1]. Therefore, we get a new interval [s_min, 1]. In the second time of computation, we set the s_maxas t, and we find

Figure 4.6: Finding control points for new curve at interval [0, t] or [t, 1].

the new control points at [0, t]. These new control points are the solution d^′_k at interval [s_min, s_max].

4.2.5 Jarvis march algorithm

The property of Bézier curve is that we can use a convex hull to cover the curve, and the convex hull is made by control points. Therefore, we use Jarvis March algorithm to find the convex hull for the curve. Jarvis March algorithm is a simple method that forms edges, and checks that every point is on the same side. The following is the process of algorithm:

1. We assume a virtual first point p_v to act as the initial reference.

2. We connect pv and p0 to a line L0, where p0 is the first point of point set P . p0 is regarded as a judgment point p_{f ront}.

3. We find the biggest angle for L₀and L_i, where L_iis a line constructed by p_{f ront}and the point p_ifrom the set P . The point p_i will be pushed into hull set H if L_ihas the biggest angle. Then, L_ibecomes L₀, and p_iturns into p_{f ront}.

4. Step 3 is repeated until p_{f ront} points to h₀, where h₀ is the first element in hull set H, and then we find the necessary points of the convex hull which are in the hull set H.

Here is an example with six points, shown in Figure 4.7.

Figure 4.7: Example of Jarvis March algorithm. The blue line is L₀, black lines are the necessary lines, gray lines and red line are L_i, and red line is the line with the biggest angle calculated with counterclockwise.

4.2.6 Algorithm

We illustrate the process for rendering N metaballs. We render the scene with many rays, and the algorithm is processed for each ray, as shown in Figure 4.8.

In the beginning, we shoot a ray from the view target and perform intersection tests with all the metaballs. The metaball which is intersected have even number of intersection points. We sort the intersection points of all intersected metaballs and construct k intervals.

For each interval i, we collect influential metaballs in a set and calculate the control points using de Casteljau's algorithm for these metaballs. To combine multiple Bézier curves, we add each control points dⁱ_kbelonging to the set of influential metaball and construct a new curve. If all the control points are smaller or bigger than threshold T , then there is no intersection between the convex hull and the horizontal line s = T , so we go to the interval i+1. Otherwise, we will do the Bézier clipping and find the new interval [smin, smax] until the interval size is small enough (i.e. smaller than ε, where ε is a very small number). The

Figure 4.8: Flow chart of metaballs rendering.

vector scalar t for ray direction is calculated by t = t_i + (t_i+1 − ti)∗ smin. Therefore, we can obtain the exact intersect position I_p between ray and metaballs by utilizing the following equation:

I_p = V_p+ tD, (4.6)

where V_p is the position of view target, and D is the direction vector of the ray.