We demonstrate a ray casting type technique for rendering algebraic surfaces using programmable graphics hardware (GPUs). Our approach allows for real-time exploration and manipulation of arbitrary real algebraic surfaces, with no pre-processing step, except that of a change of polynomial basis. The algorithm is based on the blossoming principle of trivariate Bernstein-Bézier functions over a tetrahedron. By computing the blossom of the function describing the surface with respect to each ray, we obtain the coefficients of a univariate Bernstein polynomial, describing the surface's value along each ray. We then use Bézier subdivision to find the first root of the curve along each ray to display the surface. These computations are performed in parallel for all rays and executed on a GPU. For a function f : R 3 → R, an implicit surface can be defined by the level set of the equation f (x, y, z) = c, where x, y, z, c ∈ R. By reordering the terms, we can, without loss of generality, just consider the zero-set of the function, e.g. f (x, y, z) − c = 0. For this work, all functions and surfaces are considered to be real-valued. If the function f is an algebraic polynomial, the resulting surface is called an algebraic surface. Such surfaces can easily describe intricate, smooth shapes with varying topology, and also allow for the interpolation and blending of such surfaces. Their polynomial nature also makes it easy to find analytic directional derivatives, normals and curvature. Our algorithm is based on the classical ray casting approach, where rays pass from the eye, and we wish to find the location of each ray's first intersection with the surface. For a given algebraic surface f (x, y, z) = 0 of total degree d, we are only interested in intersections inside a view volume V . Hence, it is natural to only consider the segment of each ray inside V . Let f and b represent the front and back coordinate of the ray as it enters and leaves V . The segment is then given by the parametrization p(t) = (1 − t)f + tb, t ∈ [0, 1]. We combine this parametrization with the equation of the algebraic surface, yielding a univariate polynomial g(t) = f (p(t)) of degree d. Now the problem of finding the closest intersection to the eye point, can be reformulated to finding the smallest By using the ray casting formulation, the process of extracting the zero set of one multivariate polynomial, has been replaced by a multitude of independent root finding problems. The GPU is well suited for such parallel computations, and we summarizes our algorithm as follows: 1. Issue spanning geometry to the GPU. The spanning geometry defines our view volume V , and initiates rendering. 2. Let the GPU rasterize the spanning geometry into a set of fragments, yielding the coordinates f and b for each ray. 3. Find the univariate polynomial, g(t), representing the value of the algebraic surface along each ray. 4. Find the smallest t such that g(t) = 0. If no zeros are found, we discard the fragment. 5. For the surviving fragments, calculate the final depth and color by using information from the algebraic surface. pixel plane For each ray segment inside the view volume V , we find the front (f) and back (b) vectors in barycentric coordinates. Then we use blossoming to find the coefficients of the univariate BernsteinBézier function along each ray. Several other works discuss ray casting like algorithms on the GPU, and current interest is high. However, direct hardware support still seems to be a long way of, and most approaches involve a considerable amount of tweaking in order to execute on a GPU. Donnelly [1] uses sphere tracing to effectively implement bounded displacement mapping using fragment shaders. Loop and Blinn [2] has demonstrated rendering of algebraic surfaces, which is very similar to our approach. Their focus has however been on rendering piecewise surfaces spanned over a large number of tetrahedrons. However, they use analytic root finding algorithms, thereby limiting the maximum degree of surfaces they can visualize to general quartic surfaces, as it is not possible to analytically solve general equations of degree 5 and higher. FIGURE 2: A quintic surface, having the maximum possible number of ordinary double points, 31. It is also known as a Dervish surface. For the space of real polynomials of total degree d, it is well known that they can be expressed in the Bernstein-Bézier (BB) basis, which is preferable for use in computations since it is numerically optimal [3] . The Bernstein basis functions B d i of degree d can be expressed as: Here, β β β = β β β(x) = (β1(x),...,βs+1(x)) denotes the barycentric coordinate of a point x ∈ R s with respect to a set of base points, (v1,...,vs+1) ∈ R s , which form a non-degenerate simplex Σs = conv(v1,...,vs+1). For instance, the 2-simplex is a triangle and the 3-simplex a tetrahedron. A polynomial then takes the form f (x) = ∑i biB d i (β β β). The de Casteljau algorithm provides a fast and numerically stable way of evaluating a BB-polynomial. Let bi denote the coefficients of an algebraic function in BB-form and compute repeated convex combinations of the base points by setting: (2) The de Casteljau algorithm also admits a subdivision formula for Bézier curves, and the control polygon of such a curve converges to the curve itself. In particular, such a curve have no more roots than there are sign changes in the control polygon. If we modify the de Casteljau algorithm by introducing a sequence of parameter points, (β β β1,...,β β βd) and use a distinct point at each level of recursion, we arrive at the blossom of the polynomial, which we denote as B. The blossom is given by: The key to finding a curve representing the value of f along the line segment xy is a surprising property of the blossom, which we summarize in the following lemma: Lemma 1 Given two points α α α = β β β(x) and γ γ γ = β β β(y), x, y ∈ Σs and a multivariate BB-polynomial f (β β β). Then the straight line segment from x to y is mapped to a univariate BB-polynomial of degree d, and its coefficients are given by repeated blossoming of the endpoints, To summarize, the blossom provides us with the algorithms needed to convert a trivariate polynomial into a set of univariate polynomials with respect to each ray.  To initiate our algorithm, we must send geometry encompassing the barycentric domain to the GPU. Theoretically any convex shape will suffice, but in practice shapes for which it is trivial to calculate the barycentric front and back coordinates for each ray are best suited. Blossoming The most computationally demanding subtask of the visualization is the calculation of the coefficients of the curve representing the value of f along each ray. We find these coefficients by blossoming and form an algorithm based on (3). If recursion had been supported by the GPU, we could have implemented (3) directly. However, current generation GPUs do not support recursion, so for a given degree d, we use the CPU to unroll the recursion and generate a long expression of shader source code, which represents the blossom. We use uniform variables to represent the coefficients of the surface, and Since the blossom is recalculated for every frame, this allows us to easily modify the surface by altering the coefficients, allowing for real-time shape manipulation. Root finding Because of its simple data structure, we decided to use a heuristic root finding method based on recursive subdivision of the coefficients, as described by Schneider [5]. The idea is to recursively subdivide the polynomial in half by using the de Casteljau algorithm (2). At each level of recursion we count the number of sign changes in the coefficients, and if both sets of coefficients contain sign changes we choose to continue to subdivide the leftmost of them,while storing the other one as we might need to unwind the recursion if the first one turns out to be a false zero. If only one set of coefficients has sign changes, we of course choose this one for further subdivision. This effectively uses the Bézier convex hull property to ensure that a zero is within an interval given by the coefficient endpoints. Since it could happen that there is a false zero we must be able to unwind the recursion and choose another interval to subdivide. On a regular CPU this would be trivial, but as the GPU lacks a stack we store the other potential coefficients and the recursion level, and when we detect a false zero we unroll to the last set of potential coefficients. The number of possible false zeros is bounded by the degree of the surface, so we know in advance how many sets of coefficients we might need to store. We terminate the recursion when we can decide that there is no root (no sign changes of the coefficients), or when we reach a fixed level of subdivisions. Our work presents a real-time rendering strategy for algebraic surfaces up to total degree five. However, as GPU technologies advance we believe it can be improved further. The two most obvious issues, and thus prime candidates for further work, are to increase rendering speed and the degree of surfaces we can visualize. Using our current root finding algorithm, we are able to stably isolate roots for polynomials up to degree 10, the limiting factor being that we run out of temporary registers on the GPU. DX10 GPUs will greatly increase the number of such registers. Our current approach is computationally intensive, and we are currently experimenting with other methods, such as B-spline knot insertion [4] or hybrid methods like Brent-Dekker. However, these methods are very dynamic, and effective GPU implementations are not trivial, and will probably be based on multipass methods. Our current blossoming algorithm is very poor at utilizing the vectorization possible within each fragment pipeline. The approach of [2] uses tensor contraction, which has a higher computational complexity than our method, but is very effective on a GPU. Efficient blocking of discrete Bspline matrices could be a way to get the best of both worlds, and be the key for effectively rendering high order surfaces.
Seland and Dokken REF rendered low-degree algebraic surfaces by using blossoms of trivariate Bernstein-Bézier functions over tetrahedra.
Real-Time Algebraic Surface Visualization
Computer Science