.
Efficient Antialiasing 24.july.1997 GFX
by Hin Jang
Increasing the sampling rate or removing the high frequency components of an
image are two ways of limiting the effects of aliasing. Both methods,
however, have high rendering costs. An incremental algorithm that is derived
in the spatial domain under a subjectively meaningful error term is
described herein. Its simplicity suggests a practical hardware
implementation and produces the same pixels as Fujimoto-Iwata's algorithm at
a fraction of the latter's computational cost [1, 2].
Let y = f(x) be a curve digitised in the raster plane. Sampling in a scan
coherent fashion requires that for every pixel along the x axis, the sample
value f(i) is quantised to some integer Yi. A dynamic error arises if
| f(i) - Yi | > 0
In trying to mimimise the loss of dynamic information, thus producing a more
appealing image, an antialiasing technique should aim to perserve the
convexity of the curve. This is achieved by applying the following two-point
scheme.
For a given value of (i, f(i)) the ideal addressable pixel may be
visually simulated by defining
I[i, floor(f(i))] = Io[ceiling(f(i)) - f(i)]
I[i, ceiling(f(i))] = Io[f(i) - floor(f(i))]
where I[i,j] is the intensity of pixel (i,j) and Io is the
intended intensity of the curve at (i, f(i)). If pixel (i,j) is a
considered a unit square centered at (i,j) then pµ = (i, f(i)) is
the center of gravity of the two points
p0 = (i, floor(f(i)))
p1 = (i, ceiling(f(i)))
because
Iopµ = I[i, floor(f(i))]p0 + I[i, ceiling(f(i))]p1
The lit area Io is focused at pµ, the perceived pixel exactly on
the original curve.
Given a line in the first octant, let (x0, y0) and (x1, y1) be the endpoints
of the line where x1 > x0 and y1 > y0. If the line is defined as f(x) = kx,
the equations of the two-point antialiasing scheme can rewritten as
I[x, ceiling(kx)] = I0[kx - floor(f(i))]
I[x, floor(kx)] = I0 - I[x, ceiling(kx)]
The values of (x, floor(kx)), (x, ceiling(kx)), and their intensities I(x,
floor(kx)) and I(x, ceiling(kx)) can be determined using an incremental
algorithm operating on a single integer D, a machine word of n bits [2]. As
initial values, D = 0 and I(x0, y0) = I. The increment of D is
d = floor(k2n + 0.5)
where the overflow of D = D + d is recorded. When D overflows, the two-point
pixel band moves diagonally; otherwise it moves horizontally. Given 2m - 1
discrete grey-scale intensity levels,
I(x, ceiling(kx)) = I0(kx - floor(kx))
= (2m - 1)(D2-n + (k - d2-n)x)
= D2m-n + (2m - 1)(k - d2-n)x - D2-n
Assuming n > m, the intensity can be approximated by the first term; the m
most significant bits of D. Given that Io = 2m - 1, the integer 2m - 1 -
D2m-n is the inverse of D2m-n. As such, I(x, floor(kx)) is the bitwise
inverse of I(x, ceiling(kx)).
Line()
{
PutPixel(x0, y0, I)
PutPixel(x1, y1, I)
D = 0
k = (y1 - y0) / (x1 - x0)
d = floor(k * 2^n + 0.5)
while(x0 < x1) {
x0++
x1--
D = D + d /* module 2^n addition */
if (D overflows) {
y0++
y1--
}
I = D div 2^(n-m)
PutPixel(x0, y0, I)
PutPixel(x1, y1, I)
I = ~I /* bit-wise inverse */
PutPixel(x0, y0+1, I)
PutPixel(x1, y1-1, I)
}
}
It suffices to consider one octant of the circle x2 + y2 = r2 to extend the
two-point antialiasing scheme for this primitive. The pair of pixel
intensity equations is then
I(floor(h) ,j) = I(ceiling(h) - h)
I(ceiling(h), j) = I - I(floor(h), j)
where
h = sqrt(r2 - j2)
1 <= j <= r / sqrt(2)
To move the pixel band to the left by one step, assuming scan-conversion
occurs in the y axis from 0 to r / sqrt(2), the critical values of t in
ceiling(sqrt(r2 - (t - 1)2)) - ceiling(sqrt(r2 - t2)) = 1
must be computed. Wu showed that the equation holds if and only if
ceiling(sqrt(r2 - (t - 1)2)) - ceiling(sqrt(r2 - t2)) >
ceiling(sqrt(r2 - t2)) - sqrt(r2 - t2)
[2]. As such, for given values of r,
ceiling(h) - h, 1 <= j <= r / sqrt(2)
determines the pixel positions and pixel intensities. For an intensity range
from 0 to 2m - 1,
D(r, j) = floor((2m - 1)(ceiling(h) - h) + 0.5)
which gives rise to
I(floor(h), j) = D(r, j)
I(ceiling(h), j) = bitwise-wise inverse of D(r ,j)
since
I(ceiling(h), j) + I(floor(h), j) = I = 2m - 1
A D(r,j) table is used to facilitate the speed of the algorithm. If RMAX is
the maximum allowable radius, the table size is sqrt(2) / 4 * RMAX The
two-point antialiasing scheme for circles can also be realised at the
hardware level.
Circle()
{
i = r
j = 0
PutPixel(i, j, I)
T = 0
while (i < j) {
j++
if (D(r, j) > T)
i--
I = ~D(r, j) /* bit-wise inverse */
PutPixel(i, j, I)
PutPixel(i-1, j, D(r, j))
T = D(r, j)
}
}
----------------------------------------------------------------------------
[1] Fujimoto, A., and K. Iwata, "Jay-free Images on Raster
Displays" IEEE Computer Graphics and Applications, 3(9):26-34,
December 1983
[2] Wu, X., "An Efficient Antialiasing Technique," Computer
Graphics, SIGGRAPH 1991 Proceedings, 25(4):143-152
.
Isosurface Generation 16.july.1997 GFX
by Hin Jang
An isosurface is defined by a set of points that satisfy the following
equation
S(x, y, z) - C = 0
where S is a spatial function and C is a constant. The surface is usually
displayed as set of triangles, all of which are formed by local
intersections between cells and the surface. Cells that do not intersect the
surface are not part of the volume. Ito and Koyamada developed an efficient
algorithm that visits only intersecting cells and cells that are included in
cell lists [2].
Isosurface generation consists of a preprocess and a main process.
main()
{
/* -------------------------------- Pre-process */
ExtractExtrema()
GenerateGraph()
GenerateBoundList()
/* -------------------------------- Main process */
while(1) {
C = SpecifyScalarValue()
GenerateSurface(C)
}
}
Determining the extrema points requires a comparison of the scalar values of
all points for each cell.
ExtractExtrema()
{
for (each cell) {
if (point_is_not_maximum)
mark point as not_maximum
if (point_is_not_minimum)
mark point as not_minimum
}
for (each point) {
if (point is not either 'not maximum' or 'not minimum')
insert point into extrema list
}
}
An extrema graph is a group of arcs that connects two extrema points a and
b, the start and goal points, respectively. The cost of the graph, (ie. the
number of cells inserted into the cell list), is minimal when closer extrema
points are chosen. Starting from the cell that contains the start point,
adjacent cells that share the face whose distance is minimum are visited in
order. The minimum and maximum values of the start and goal points are
updated by the scalar values of the points visited herein. This process is
repeated until the traversal reaches the cell containing the goal point. If
reaching the goal point results in leaving the volume, the traversal is
terminated and a similar traveral begins when the next-closest extrema is
declared the goal point.
GenerateGraph()
{
for (each extrema i)
extrm[i].flag = extrm[i].gID
for (each extrema n) {
for (each other extrema i) {
if (extrm[n].flag != extrm[i].flag)
if (extrm[i] is not too far)
Enqueue(extrm[i])
}
while (path_is_not_connected) {
Dequeue(extrm[m])
if (Make_Graph1(extrm[n], extrem[m]) == SUCCESS)
break
if (Queue_Is_Empty() == TRUE) {
Make_Graph2(extrm[n], extrm[m])
break
}
}
for (each extrema i)
if (extrm[i].flag == extrm[m].flag)
extrm[i].flag = extrm[n].flag
}
}
MakeGraph1(extrmA, extrmB)
{
visit_cell = A cell that contains extrmA
while (1) {
arc.cID[arc.numC++] = cell
arc.value[0] += arc.eID[0]
arc.value[1] -= arc.eID[1]
if (visit_cell includes extrmB)
return SUCCESS
visit_cell = adjacent cell that intersects arc
if (visit_cell == NULL)
return FAILURE
}
}
MakeGraph2(extrmA, extrmB)
{
visit_cell = A cell that contains extrmA
while (1) {
arc.cID[arc.numC++] = cell
arc.value[0] += arc.eID[0]
arc.value[1] -= arc.eID[1]
if (visit_cell includes extrmB)
return SUCCESS
visit_cell = adjacent cell that shares the nearest arc
}
}
A boundary cell list is a group of cells that have faces not connected to an
adjacent cell.
GenerateBoundList()
{
for (each boundary cell cID[i]) {
cID[i].max = FindMaximum()
cID[i].min = FindMinimim()
}
Generate_Sorted_Maximum_Based_List()
Generate_Sorted_Minimum_Based_List()
}
Isosurfaces propagate themselves from seed cells. The minimum-based boundary
cell list is traversed until the minimum value becomes higher than the given
value. If the maximum value of the visited cell is higher, the cell is
denoted as the seed cell. The maximum-based boundary cell list can also be
used to select a cell from which propagation occurs. Seed cells are then
searched for by traversing the arcs of the extrema graphs. If the specified
scalar value lies between the minimum and maximum values of the arc, all the
cells in the list are visited in order. An isosurface is generated when an
extracted seed cell is unmarked.
GenerateSurface(float C)
{
for (each cell bcell[i] in the boundary cell list)
if (bcell[i].min < C
&& bcell[i].max > C
&& bcell[i] is not marked)
PropagateSurface(bcell[i], C);
for (each arc in the extrema graph)
if (arc.value[1] < C && arc.value[0] > C)
for (each unmarked cell in the arc) {
if (cell.cID[i] is intersected)
PropagateSurface(cID[i], C)
}
}
PropagateSurface(int cellID, float C)
{
Enqueue(cellID)
while (Queue_Is_Empty == FALSE) {
Dequeue(cellID)
Draw_Triangle(cellID);
Enqueue(cellIDs of all unmarked adjacent, intersected cells)
Mark_Cells_in_Queue();
}
}
----------------------------------------------------------------------------
[1] Grimm, C.M., and J.F. Hughes, "Smooth Isosurface
Approximation," Eurographics Conference on Implicit Surfaces,
April 1995
[2] Itoh, T., and K. Koyamada, "Isosurface Generation by Using
Extrema Graphs," 1994 IEEE Visualisation Conference, 77-83
[3] Lorensen, W.E., and H.E. Cline, "Marching Cubes: A High
Resolution 3D Surface Construction Algorithm," Computer Graphics,
SIGGRAPH 1987 Proceedings, 21(4):163-169
[4] Matveyev, S.V., "Approximation of Isosurface in the Marching
Cube: Ambiguity Problem," 1994 IEEE Visualisation Conference,
288-292
[5] Montani, C., R. Scateni, and R. Scopigno, "Discretised
Marching Cubes," 1994 IEEE Visualisation Conference, 281-287
.
Shadow Rendering Algorithms 02.september.1997 GFX
by Hin Jang
Shadows provide a visual cue to the spatial relationships among objects in a
given scene. Simulating hard shadows is possible using clever approximations
[1]. More robust algorithms that simulate shadows cast by moving, complex
objects onto multiple planar surfaces have also been developed [2, 3]. Soft
shadows caused by extended light sources, however, require greater
computation; those areas with an umbra (fully shadowed regions) and penumbra
(partially shadowed regions). At true interactive rates, simulating soft
shadows in a dynamic and complex scene require high-end, parallel hardware.
The principle of shadow generation is based on a modified perspective
projection and is easily incorporated into any rendering pipeline.
A point P(x, y, z) that is illuminated by a light source L(x, y, z) casts a
shadow at
S = a(L - P) - P
and assuming the shadow falls upon a ground plane, where z = 0, solving for
a yields
zP
a = ---------
zL - zP
The coordinates of the shadow point S are then
xPzL - xLzP
xS = ----------
zL - zP
yPzL - yLzP
yS = ----------
zL - zP
Assuming a homogeneous coordinate system where
xS = xS' / wS'
yS = yS' / wS'
wS' = zL - zP
S can be calculated using matrix multiplication
| zL 0 0 0 |
| 0 zL 0 0 |
[ xS' yS' 0 wS' ] = [ xP yP zP 1 ] | -xL -yL 0 -1 |
| 0 0 0 zL |
To generalise the 4x4 matrix so that the algorithm can perspectively project
a shadow point onto any plane, consider an arbitrary point on the line from
P to L represented in homogeneous coordinates
aP + bL
where a and b are scalars. The matrix above assumes a ground plane G whose
normal is represented as a four-element column vector.
| 0 |
| 0 |
G = | 1 |
| 0 |
The point on the shadow line that intersects G is given by (a, b) where
(aP + bL) . G = 0
Rearranging the equation yields
a(P . G) + b(L . G) = 0
which can be satisfied by any scalar multiple of
a = (L . G)
b = -(P . G)
The shadow point S is then
S = (L . G)P - (P . G)L
Recall that P and L are row vectors and G is a column vector. Matrix
multiplication is associative; therefore,
(P . G)L = P(GL)
and note that GL is a 4x4 matrix and (L . G) is a scalar. Given that I is a
4x4 identity matrix
S = P[(L . G)I - GL]
The generalised 4x4 matrix used to perspectively project a shadow point onto
any plane is
[(L . G)I - GL]
where G is the normal of the plane onto which S is cast.
Shadows points, or a given subset, forms the vertices of a hard shadow. The
real world, however, contains mostly soft shadows. The ideal method of
computing soft shadows should 1) have a memory requirement independent of
the number of light source samples and 2) exploit graphics hardware that
determines visibility and calculates shading. If such hardware is
unavailable, calculations can be done in software using a graphics
subroutine interface, as such OpenGL. The premise of the algorithm developed
by Heckbert and Herf is to precompute the emitted radiance at each point on
every polygon to produce a set of hard shadows [2]. These shadows are then
averaged to compute the soft shadow image.
The emitted radiance from a surface is given by summing the contributions
over all sample points on all light sources
Lc(x) = pc(x)Lac
nl
---- ---- cos+tli cos+tli' Lc(xli')
+ \ \ (ali pc(x)) ------------------------- v(x, xli')
/ / pi rli2
---- ----
l i=1
where
* x = (x, y, z) is a point on the reflective surface, and xli' is the
sample point i on light source l
* ali is the area associated with sample point i
* t is the polar angle (angle from normal) at x, t' is the angle at xli'
* r is the distance between x and xli'
* t, t' and r are functions of x and xli'
* Lc is the outgoing radiance at point x for color channel c
* Lac is the ambient radiance
* pc is reflectance
* v(x, xli') is a Boolean visibility function that equals 1 if point x is
visible from point xli', else 0
* cos+t = max(cost, 0), for backface testing
* pi is the constant
The summand is the product of three factors
1. the area times the reflectance of the receiving polygon
2. the cosine of the angle on the receiver times the cosine of the angle
on the light source times the radiance of the light source divided by
pi r2
3. the visibility between a point on a light source and each point on the
receiver
To generate the set of hard shadow images, the algorithm requires a
projective transformation. For a given parallelogram P fitted around the
receiver polygon, a pyramid can be constructed with P as its base and the
light point as its apex. The 4x4 homogeneous transformation that
accomplishes this is the matrix tranform M of interest. The pyramid lies in
object space with coordinates (xo, yo, zo) with apex a and its parallelogram
base has one vertex at b and edge vectors ex and ey. Applying M yields a
parallelepiped in unit screen space with coordinates (xu, yu, zu). Viewed
from a, the left and right sides of the volume map to parallel planes xu = 0
and xu = 1. The bottom and tops sides map to yu = 0 and yu = 1, and the base
plane and the plane parallel through the apex map to zu = 0 and zu ->
infinity.
The matrix M has the form
/ axnxx axnxy axnxz -axnx.b \
M = | aynyx aynyy aynyz -ayny.b |
| 0 0 0 1 |
\ awnwx awnwy awnwz -awnw.a /
where
nx = ew × ey ax = 1 / (nx . ex)
ny = ex × ew and ay = 1 / (ny . ey)
nw = ey × ex aw = 1 / (nw . ew)
Averaging the hard shadows is done by a weighted, linear combination of
images
----
Ac(x, y) = \ aiIic(x, y)
/
----
i
where Iic is a channel (red, green or blue) for image i, ai is the weight,
and Ac is a channel of the accumulator array. The algorithm to precompute
radiance textures is then
RadianceTexture()
{
for (each receiver polygon R) {
select resolution of receiver's texture (sx X sy pixels)
clear accumulator image of sx X sy pixels to black
create temporary image of sx X sy pixels
for (each light source l) {
if (l is behind R || R is behind l)
skip to next l
for (each sample point i on light source l) {
if (xli' is behind R)
skip to next i
compute transform matrix M, where a = xli'
with parallelogram base fitted around R
set transform matrix to scale (sx, sy, 1) . M
set clipping planes to zu.near = 1 - e
and zu.far = big
draw R with illumination from xli', as
described in the equation for Lc, into
temporary image
for (each other object in scene) {
draw object with ambient colour into temporary image
}
add temporary image into accumulator image with weight
al / nl
}
}
save accumulator image as texture for polygon R
}
}
To render a scene with shadows, simply
RenderWithShadows()
{
for (each object in scene)
if (object receives shadows)
draw object with radiance textures and no illumination
else draw object with illumination
}
----------------------------------------------------------------------------
[1] Blinn, J. F., "Me and My (Fake) Shadow," IEEE Computer
Graphics and Applications, 8(1):82-86, January 1988
[2] Heckbert, P., and M. Herf, Simulating Soft Shadows with
Graphics Hardware, CMU-CS-97-104, CS Dept, Carnegie Mellon
University, January 1997
[3] Loscos, C., and G. Drettakis, Interactive High Quality Soft
Shadows in Scenes with Moving Objects, iMAGIS/GRAVIR-INRIA, 1997
[4] Stewart, J. A., and S. Ghali, Fast Computation of Shadow
Boundaries Using Spatial Coherence and Backprojections, Department
of Computer Science, University of Toronto, 1995