The invention relates generally to a computer apparatus and a data processing method and apparatus for the calculation of raster data from scalar fields and vector fields for the purpose of building an explicit representation of levels of the one or more fields or for visualizing the one or more fields using iso-surface and/or volume visualization methods. One specialization of the invention is the visualization of piecewise algebraic hypersurfaces. More particularly, the invention relates to a computer apparatus and data processing method which may be employed with advantage within computer aided design, computer aided manufacturing, rapid prototyping, algebraic geometry, finite element pre-processing, finite element post processing, computer animation, interpretation and visualization of medical data, and interpretation and visualization of petroleum related data where fast and topologically correct and geometrically accurate raster data is needed. The apparatus combines one or more CPUs with one or more SPUs (Stream Processing Units), with one embodiment of the SPUs being one or more Graphics Processing Unit (GPUs).
Raster data is information stored in a d-dimensional grid of cells, where d is often in the range of 1 to 4. Each cell is indexed by d discrete numbers from some finite range. All the cells contain the same number and types of information elements. A cell can contain different information element types. Raster data is popular in such fields as computer vision, image processing, and computer graphics where they are commonly used to represent images (2D raster grids), videos and MRI volumes (3D raster grids), etc. In this invention the element types used can have any interpretation. However, often they will represent point locations, gradients, intensity values, shading data, RGB values, values of the field integrated over sub-sets such as sub-volumes or line segments.
In computer graphics 3D visualization is based on making an image of an object using either perspective or parallel projection. In this process the part of space that is to be viewed is defined by a volume V, and through projections and transformations the object is mapped into standard viewing coordinates. The transformations can be found in computer graphics text books, and are not important when describing the principles of this invention.
In traditional ray tracing, the visualization is based on tracing rays starting in the eye point (centre of projection) and going through the pixels of the image to be produced. Then the first intersection of the ray with the scene visualized is calculated. Dependent of the properties of the object hit, e.g. transparency, refraction or reflection, further searching is done along the ray, or along the reflection ray. The focus of our invention is only related to what happens inside the volume V.
Within visualization of 3D data and time varying 3D data, visualization traditionally has been based on either triangulation of the geometry e.g. by the “marching cubes”-technique, or by ray tracing. Traditionally such algorithms have been run on the CPU (Central Processor Unit). However, the publication “The Ray Engine”, Graphics Hardware 2002, Sep. 1-2, 2002, Saarbrücken, Germany, The Eurographics Association, 2002, by Mathan A Carr et al, disclose a computer system for ray tracing of triangulations by using the GPU (Graphics card) as a computational resource. However, this is based on the intersection of triangles and the ray and not a segment of the ray and a scalar field as in the invention we present.
The PCT patent application PCT/NO 2005/000453 addresses a computer apparatus combining one or more CPUs with one or more SPUs (Stream Processing Units) for finding the topology and a geometric description of the intersection calculations of geometric objects and for detecting self-intersections in geometric objects. The computer apparatus proposed in the proposed invention addresses the generation of raster data, aimed at visualization or object segmentation and consequently addresses a different application area. Some of the similar mathematical approaches are used in both inventions; however, these are general in nature.
One typical example of a prior art scalar field visualization is ray tracing of level sets. Here each ray is traced from the centre of projection, through the pixel to be generated until a point in the scalar field possibly is identified with the specified scalar field level value. In the case when reflection or refraction is of interest the search is continued along the refracted ray and along the reflected ray until no more points are hit (inducing refraction and reflection) or the maximal level of recursion (number of reflections or refractions). One implementation of such an approach is described in the paper “Anders Adamson, Marc Alexa, Andrew Nealen. Adaptive Sampling of Intersectable Models Exploiting Image and Object-space Coherence” presented at the Third Eurographics Symposium on Geometry Processing, Jul. 4-6, 2005. However, in this approach the ray surface intersection is only performed on the CPU, not on the SPU (GPU) as proposed in the present invention.
Another examples is “marching cubes”, where for each cell in a volume grid, a triangulation of the iso-surfaces inside the cell is directly generated from the values in the eight corners of the cell and predefined rules on the topology of the triangulation based on the sign of the field values in the corners. The triangles generated by marching cubes can be viewed as a crude approximation to a tri-linear interpolation of the corners of the grid cell.
A related invention is disclosed in US-patent application 2004/008532 A1 as a system and method for modelling of graphical shapes, where the system comprises a central processing unit, a main processor and a graphics processor and where the calculations of the models and transformation calculations are done in the central processing unit which contains a first and a second vector processor and where the calculated data are sent to the graphics processor for rendering an later displaying on a display means.
U.S. Pat. No. 657,571 B1 describes an image system comprising a number of graphics processors, each activating their own processes based on instructions received from an external host computer and sending subsequently data to a rendering apparatus.
Briefly, and in general terms, the present invention provides improvements in creating raster data from scalar or vector fields and an apparatus whereby data processing speed and raster data quality increase substantially.
More particularly, by way of example and not necessarily by way of limitation, the present invention provides an organization of the scalar or vector field processing to be performed on at least one CPU for controlling the process, and at least on stream processor or SPU for execution of computer intensive task of the scalar field or vector field processing.
When the raster data is generated as part of a visualization process the data processing speed and visual quality increase substantially and will be within pixel accuracy. This is especially important for scalar or vector field singularities, points where the field gradient vanishes, as the level set topology change around singularities. If sub-pixel accuracy is used when generating the raster data, then improved image quality can be achieve through anti-aliasing techniques.
When the raster data is generated as part of a segmentation process the data processing speed increase substantially and the level set topology generated will be correct within the geometric distance between the elements of the raster data.
Hence, the present invention satisfies a long existing need for enhanced generation of raster data from scalar fields and vector fields.
Hence, in a first aspect of the present invention, there is provided a computer apparatus for determining high quality raster data generation of fixed and time varying scalar fields or vector fields, represented by piecewise polynomials or piecewise rational functions. The computer apparatus comprises at least one central processor (CPU) operative to do portions of the raster data generation algorithm, initializing sub-algorithms thereof, control the sub-algorithms, and possibly read back the generated raster data or transfer the raster data to other processors in the system. The computer apparatus further comprises at least on stream processing unit (SPU) operative to receive parts of the raster data algorithm from the at least one CPU and to execute sub-algorithms of the raster data algorithm, resulting in raster data that can be directly visualized, read back to the one or more CPUs or transferred to other processors.
In another aspect of the present invention is a method to create at least one set of raster data from a time varying or constant scalar field or vector field within a volume. Said scalar or vector field is described using any of polynomial, rational, piecewise polynomial and piecewise rational representation. The method comprises at least controlling execution of a raster data algorithm, initializing subalgorithms thereof, controlling said sub-algorithms and controlling the use of the results of the calculations involving said sub-algorithms. The method is characterized by controlling at least one stream processor unit (SPU) operative to receive sub-algorithms of said raster data algorithm and to execute said sub-algorithms of said raster data algorithm, resulting in raster data.
A detailed definition of the present invention is given by the attached claim set.
The invention will be described in detail below taking reference to the attached drawings where:
The invention addresses scalar fields and vector fields that are either constant in time or vary with time. However, as the invention relates to discrete time steps in time varying fields, the time component is not included in the detailed descriptions to simplify the presentation. When addressing the time varying fields the following description relates to a specific time step.
Although the invention is not limited to scalar fields or vector fields, we will describe the mathematics involved in terms of scalar fields in R3. A vector field can be described by one scalar field for each component of the vector field.
A scalar field is a function f from Rn to R, or a function f from Rn to C (the complex numbers). That is, it is a function defined on the n-dimensional Euclidean space with real values. Often it is required to be continuous, or one or more times differentiable, that is, a function of class Ck, k≧0. The scalar field can be visualized as an n-dimensional space with a real or complex number attached to each point in the space. The derivative of a scalar field results in a vector field called the gradient field.
A vector field is a function f from Rn to Rm, or a function f from Rn to Cm. A vector field can consequently be described by one scalar field for each component of the vector field. Consequently the discussion following also covers vector fields.
A trimmed scalar field or vector field is a scalar of vector field where the part of Rn addressed is defined by a number of other scalar fields ti(x), i=, . . . , M. The only part of the scalar field or vector field we address are described by points satisfying ti(x)≧0, i=, . . . , M, or alternatively ti(x)≦0, i=M.
The scalar fields we will look into are described by piecewise polynomial functions or rational functions with numerator and denominator being piecewise polynomial functions, e.g., a tensor product B-spline function or NURBS-function, or a structure of polynomial or rational functions, each described over a volume. Frequently this volume will be a simplex (a tetrahedron in R3) with a chosen order of continuity between the pieces. However, the volume is not limited to being a simplex. Other relevant volumes are finite elements used for solution of partial differential equations by the Finite Element Method.
One aspect of mathematics of the invention is related to calculations for each polynomial piece in the description of the scalar field and thus we will concentrate in the following on each polynomial piece. The part of the scalar field we are interested in is described by a volume V, e.g., a rectangular box, sphere or a tetrahedron or any volumetric shape, often the volume will have a convex shape. If the volume is non-convex the volume can be subdivided into convex sub-volumes, and the method applied to each convex sub-volume. It should be noted that the volume V is not the same as the volumes used for the description of the structure of piecewise polynomial functions.
When addressing iso-surfaces the invention combines ideas from the two approaches by:
However, the approach is not limited to only looking at the closest iso-surface within the cell or sub-volume. Other variants of the invention are to:
In the case that we just have one cell or volume and look for values of f, with f(x, y, z)=0, we have a real algebraic surface. Consequently in addition to address scalar fields the invention addresses also sampling and visualization of real algebraic surfaces.
One preferred selection for implementation is a standard PC having one or more single- or multi-core CPUs running e.g. Windows, Linux or other operating systems with one or more stream type processors, preferably of GPU-type but not restricted to GPUs. Another selected implementation is a closer integrated CPUs and stream processors such as in the CELL-processor, consequently not limiting the stream processor to be integrated into a graphics card. The approach can also be implemented on battery powered PDA-type devices such as mobile phones that combine a CPU and programmable 3D graphics acceleration.
The simplest implementation of the approach, see
The points found can both be represented as a point on the straight line (105, 312) by just remembering the parameter value of the point on the line, or be represented as a 3D point. Which information is stored depends on the later intended use of the raster data.
Visualization of level sets of scalar and vector fields is frequently used within medicine, oil & gas industry and for the interpretation of simulation results. As the above approach uses the same concepts as used in 3D computer graphics it is straight forward to use it for visualization of level sets of scalar or vector fields. We just have to make sure that the n×m points in the rectangular domain correspond to the resolution of the image we want to generate. Doing this we ensure pixel, or sub-pixel accurate visualization of the scalar field. For iso surface visualization purposes the following steps have to be added in the process.
In the case where the SPU is an integrated part of the visualization pipeline as in the case of a GPU there is a direct connection to the displays used.
With the approach proposed above the quality of the points generated can be described as within pixel resolution when viewing from the used centre of projection. If there are isolated singularities in the set, e.g. where different braches meet the sampling seldom hits exactly on the singularity but the visual image generated will reflect the behavior of the region around the singularity and thus the singularity is indirectly represented. Zooming in the points will get closer and closer to the singularity and consequently the visual impression will be correct. By changing the constant defining the level set addressed a good visual impression of the scalar field can be produced. To guarantee that singularities are not missed, the analysis of the line segments intersecting the field can be replaced by analysis of the intersection of the field by long thin swept rectangles around each line segment. The swept rectangle of each line segment should just touch the swept rectangle of adjacent straight line segments to ensure that no gaps are left.
By sampling the level set from 3 different directions,
However, if we control which points are sampled in the different projections by restricting the gradient of the scalar field to be within regular pyramids assigned to each projection as shown if
Except where the branches of the normal field penetrate the volume V, the boundary of each disjoint data set will have to be merged with the boundary of disjoint data sets from the other projections. This merging operation will typically be performed on the CPUs.
Alternatively, no removal of points is done based on gradient direction, as this removal can be performed later in the process. In this case there will be overlapping of points sets generated in the different sampling planes.
The intersection of the ray and the volume V will trim the ray to a straight line segment between a point P0 and a point P1 with P0 closer to the centre of projection than P1, see
or as a tensor product tri-variate polynomial of bi-degrees (n1, n2, n3) and total degree
The intersection of i(t) and qm(x, y, z), can be expressed as qm(l(t))=0, where qm(l(t)) is a polynomial of degree m in t. Consequently the points we are looking for are zeroes of qm(l(t)) in the interval [0,1]. The intersection of l(t) and qm
Two central tasks to be performed for calculating the points we are looking for are thus:
An efficient and numerical stable method for finding either the polynomial ƒ(t)=qm(l(t)) or , with respective degree m or degree n1+n2+n3.
It should be noted that we do not need to find all the zeroes of ƒ(t)=0 in the interval [0,1] if we are just looking for the zero closest to the centre of projection.
Dependent of the polynomial basis used for describing the scalar field different algorithms can be used for finding ƒ(t). We will now denote the scalar field q(x, y, z) and assume that it has total degree m.
One implementation is to let program components running on the one or more CPUs use Blossoming to generate the analytical expression of f(t), and let the one or more CPUs use these analytical expression for generating program segments to be run on the one or more SPUs.
Finding the Zero of ƒ(t) when Using the Bernstein Basis.
The univariate function ƒ(t) expressed in a Bernstein basis looks like
We have translated the parameter interval of the polynomial to [0,1] to simplify the description. The Bernstein representation of the polynomial has a number of nice properties:
By subdividing the function into sub-pieces represented in the Bernstein basis, all zeroes in the interval [0,1] can be efficiently determined. Following the first introduction of algorithms for finding zeros of Bernstein basis represented polynomials many different formulation of these algorithms have been published. As a central part of the proposed approach is the use of a stream processor for finding such zeroes, the actual implementation has to be adapted to the specificities of stream processors and the actual performance characteristics of the stream processor used.
When implementing the approach on a programmable graphics card the functionality of the card will be used to map the geometry into standard viewport coordinates. The n×m points in a rectangular region is mapped onto the viewport, and n×m points correspond to the viewport resolution. The intersection of the ray and the 3D box is performed before the intersection of the ray and the iso-surface. Thus we can reduce the intersection of an infinite straight line and the iso-surface, to the intersection of a parametric described straight line segment and the isosurface.
A 3D scalar field can be represented in many ways.
By a Bernstein polynomial of total degree m in barycentric coordinates over a tetrahedron with corners P1, P2, P3, P4
The part of the polynomial inside the tetrahedron satisfies β1, β2, β3, β4≧0 with β1+β2+β3+β4=1. The Cartesian coordinates of a point P=(x, y, z) represented in barycentric coordinates are calculated by P=β1+P1+β2P2+β3P3+β4P4.
of degrees n1, n2, n3.
By a tri-variate tensor product B-spline volume of degrees n1, n2, n3.
where Bi,n
The power basis represented polynomial qm(x, y, z) of total degree m can be converted to a Bernstein polynomial of total degree m and visa versa. The grid representation can be interpreted as a tri-variate B-spline, the tri-variate B-spline can be converted to a structure of tri-variate Bezier basis represented volumes. A Bezier represented volume of degrees n1, n2, n3 can be converted to a Bernstein basis represented polynomial of total degree n1+n2+n3 over a tetrahedron, or to a power basis representation.
Consequently the invention is valid for a wide variate of scalar field representations.
It also includes rational variants of the representations above.
To illustrate the principles of the approach we give some examples.
Assume that the scalar field/algebraic surface is represented by a Bernstein polynomial of total degree m in barycentric coordinates over a tetrahedron with corners P1, P2, P3, P4. Assume that the constant level we want to visualize has value c. For an algebraic surface c≡0. Consequently we can focus on the algebraic surface qm(β1, β2, β3, β4)=c=0. The tetrahedron can be used as the volume V describing the portion of the surface that we are interested in. However, the volume V can also be independent of the description of the tetrahedron e.g. a sphere.
From each pixel to be found in the image we intersect the corresponding projection ray with the view volume to find the actual line segment to be intersected with the scalar field,
In the case where the scalar field/algebraic surface is describe in a tri-variate tensor product rational Bernstein basis we address the intersection of qn
Other representations of an algebraic surface can be converted to these formats.
A typical visualization process will have multiple stages:
Each of the visible volume elements can then be treated as a scalar field represented by one rational function. However, the ray has to be clipped both by the view volume and the boundaries of the volume element.
The point grid is regarded as the control polygon of a tri-linear B-spline surface and visualized with the approach described for tri-variate rational tensor product B-splines.
In addition to the scalar field level set qm(β1, β2, β3, β4)=c, a set of other scalar fields are given ti(β1, β2, β3, β4)=0, i=1, . . . , M. The only subset of qm(β1, β2, β3, β4)=c we are interested in is described by points satisfying tl(β1, β2, β3, β4)≧0, i=1, . . . , M. The algorithmic addition to the examples above is that points found in the ray level set intersection are discarded if the point satisfies tl(β1, β2, β3, β4)<0 for any i=1, . . . , M. In this case the search for intersection points has to continue further along the ray. The trimming test can be performed as part of the algorithm running on the SPU.
If the approach of constructive solid geometry (CSG) is used for describing volume objects, algebraic or piecewise algebraic surfaces are used for describing the half-spaces defining the extent of the volume. For each piece of an implicit surface that is part of the surface of the volume object, the adjacent (immediate neighboring) implicit surfaces take on the role as trimming surfaces. Consequently this example shows how the invention can be used for generating raster data and visualization of volumes described by constructive solid geometry.
Volume visualization of scalar fields is frequently used within medicine, oil & gas industry and for the interpretation of simulation results.
In addition to the scalar field we assume that a transparency field a from Rn to R is given where the values of α all are in [0,1], e.g., all are greater or equal to zero or less or equal to 1. The value zero means that the scalar field is transparent; the value 1 means that the scalar field is opaque.
In iso-surface visualization to find each of the n×m raster data we intersect the infinite straight line with the iso-surface and select the point inside the volume V closest to the centre of projection in the volume. For scalar field volume visualization we rather calculate the raster data produced by exact integration, numeric integration or sampling by combining the scalar field values and transparency values.
Let the [a,b] be the interval on the straight line l(t) within the volume V, with a representing the end closest to the centre of projection. Let f(t) be scalar field values along the straight line and a(t) be the corresponding transparency values. Let the function h: R to Rl l>0, express how we assign a visual appearance to the scalar field. Then one alternative for calculating the value of the raster data corresponding to the straight line is to use an integrator unit—a program running on a processor using analytical formulas for calculating integrals of functions—to calculate the integral
Here the integral
expresses how visible the point l(t) is seen from the centre of projection. The value α(t)h(ƒ(t)) expresses the visual contribution of the point l(t). If exact integration is not feasible, then a numerical integrator unit—a program running on a processor using numerical integration methods for calculating the integrals of functions—can be used.
The approach is not limited to only the above integral calculation.
The approach opens up for combining the use of iso-surface and volume visualization. One approach for doing this is first to calculate the iso-surface visualization taking care to remember the location b on each straight line representing the identified iso-surface point (
Here g(b) is the raster data calculated by the iso-surface visualization, and the integral
tell show much the iso-surface visualization is obscured by the volume visualization. The raster data f(t) of the scalar field is converted to the same representation as g(b) by the function h(f(t)).
Number | Date | Country | Kind |
---|---|---|---|
20062770 | Jun 2006 | NO | national |
Filing Document | Filing Date | Country | Kind | 371c Date |
---|---|---|---|---|
PCT/NO2007/000196 | 6/7/2007 | WO | 00 | 3/31/2009 |