Graphics Notes
Graphics Notes
Graphics Notes
(2009)
COMPUTER GRAPHICS
Although this course largely treats the theoretical aspects of Computer Graphics, students are
expected to implement important algorithms in C/C++ utilising the OpenGL API.
Ref. 3 above is a classical and extensive text on the subject, but these notes have been
developed mainly from an earlier Pascal-based version of ref.1 with updates taken from ref. 2.
Page references here point to ref. 2.
Course outline
• General survey and overview of graphics systems (hardware and software)
• Fundamental algorithms for output primitives: lines, circles, ellipses, curves, area
filling, clipping
• Algorithms for controlling attributes of primitives: shading, pattern filling, colouring, anti-
aliasing
• 2D and 3D geometric transformations: translations, scalings, rotations, reflections,
shearing, inter-coordinate system transformations
• 2D and 3D viewing: World coordinate models to viewing/device coordinates
• 3D object representation/modelling: quadrics, splines and surfaces, octrees, fractal
geometry methods
• Visible surface detection, Lighting and surface rendering, Animation
Note: This document is best viewed and/or printed in colour, but students will get a b/w printed version + a pdf file
of it.
2
1.1 Introduction
Previous to the 1980s Computer Graphics (CG) was a small and specialised field. With the
advent of PCs and better graphics hardware, CG now pervades virtually all modern application
software, and is commonly deployed in interactive windows applications.
The basic task of CG is to form (synthesize) an image (on an output device e.g. CRT screen)
from a model of a real object.
The basic task of the related field of image processing, is somewhat the reverse: Analyse or
reconstruct a 2D and 3D object (i.e. an enhanced image or model) from a primitive image.
1.2 Applications of CG
i) Software design of objects (models of cars, planes etc) in science and engineering
= “computer aided design” (CAD) – usually combined with virtual reality or
animations.
ii) Presentation graphics – reports for science and business, 3D graphs
iii) Computer art – commercial advertising, “paintbrush” drawings
iv) Entertainment – movies and sound
v) Education and training
vi) Visualizations and simulations of complex processes e.g. computational fluid
dynamics (CFD) and scientific or high performance computing (HPC)
vii) Image processing – to enhance and modify existing images
viii) Graphical user interface (GUI) design – used in windows applications
Many languages provide libraries (“bindings”) implementing one or more of the above, e.g. the
early Turbo C++ for DOS uses the Borland Graphics Interface (BGI) implementing
3
So too have evolved various standards for storing graphical images as binary data on
peripheral devices, such as the Computer Graphics Metafile (CGM), Windows Bitmap (BMP),
Joint Photographic Experts Group File (JPEG), Tagged Image File Format (TIFF), Graphics
Interchange Format (GIF) and many others.
Our approach here will be to study methods for generating the fundamental building blocks on
which such libraries are constructed, namely, the “primitives” from first principles. In addition we
shall study some basic graphics algorithms employing primitives for such tasks as surface
generation, shading, geometric modelling, viewing and animation. First we review the currently
available hardware systems and software for rendering graphics.
4
screen with
inner phosphor
coating:
persistence 10-
60µs – requires
heating refresh
element
electron
emitting heated accelerating anode
(+) or mag coil electron
cathode (-)
beam
Electrons from a heated cathode are accelerated and focussed to a point beam, either with
electric fields due to +ve biased plates or with magnetic coils. When the beam hits the
phosphor coating on the screen’s inner surface, it excites its atoms to higher quantum levels,
emitting light when they fall to a lower level. Since this light generally lasts about 10-60 µs, the
screen has to be “refreshed” to maintain the displayed image. Typical refresh rates are about
60 frames per sec on standard monitors, independent of the picture complexity.
At any time the beam will hit a particular spot or point on the screen. The maximum number of
points that can be displayed without overlapping is called the resolution = the number of
points/cm that can be plotted horizontally or vertically (1-direction). The light spot intensity is
generally of Gaussian distribution, such as:-
∆
5
Thus two spots will appear separate if the separation distance d > ∆ = the width at 60% level.
Typical resolutions for high quality systems are ~ 1280 x 1024. This is larger in “high definition”
systems.
The aspect ratio:
= the ratio of the number of vertical points-to-horizontal points required to produce equal
length lines in both directions. For example,
a ratio ¾ => vertical line plotted with 3 points
= in length to horizontal line with 4 points.
....
To paint each frame takes about 1/60 – 1/80 sec = the “refresh rate” or “vertical trace rate”. A
variation involves interlacing where alternate set of lines are scanned, followed by a vertical
retrace of the other set of lines and so on. This effectively reduces the refresh rate by ½ at the
expense of some picture definition and is used with slower refresh rate devices in order to
reduce the effect of flicker.
On a b/w system with 1 bit/pixel, the frame buffer is called a bitmap. For systems with > 1
bit/pixel the FB is called a pixmap.
ending point, one line at a time, rather than scanning the entire screen. This process is called a
vector display or stroke writing or calligraphic display. The action is similar to that employed in
pen plotters and the picture (or lines composing it) can be drawn in any order required, with
refreshing done only on those lines requiring it.
The picture definition is stored as a set of line drawing commands in a section of memory
called the refresh display file or display program or display list. The system cycles through
the set of commands in the display file, drawing each line, not necessarily horizontally, one
after another. When all lines are drawn, then it cycles back to the first line command in the
buffer. It draws all component lines of the picture ± 30-60 times/sec. These displays are not
suitable for delicate shading etc but Ok for “liney” figures for which much smoother pictures
than those of raster scan displays are obtained.
electron beam
glass
Here the colour emitted depends on the depth penetrated by the beam: a slow beam
excites red only, whilst a fast beam reaches into and excites green. Intermediate
penetration gives orange and yellow – essentially 4 colours.
R R R R
G B G B G B G B
R R R R
G B G B G B G B
To achieve colour, 3 electron guns providing 3 separate beams, one to excite one
component colour of each pixel are used. Then a particular pixel colour is made up of a
combination of the 3 primary RGB colours, whose individual component intensities can
be controlled, resulting in millions of colours.
However, alignment of the composite beam so that the R-beam hits only the R-spot,
and the G-beam hits only the G-spot etc is difficult. To achieve this a metal screen is
placed before the coated surface and which contains thousands of tiny holes to allow
the electron triple beam to pass through and strike the correct component colour spots.
This metal screen is called a shadow mask. As it sits close to the screen it can used to
correctly align the 3 electron beams as shown in the figure below.
screen G B
RGB
triple dot
Better alignment can be obtained when the pixel RGB triple dots are arranged in-line
(“precision in-line delta” CRT) rather than the RGB triad format (“delta-delta”) above.
For RGB systems the colour range is huge:-
Example:
8
If 24 bits/pixel available then can use 8 bits to store intensity levels for each
colour
=> 28=256 intensity or voltage levels for each colour or component
=> ≈ 17 × 106 colours possible – called a “full” colour or “true” colour system.
Apart from the alignment problem mentioned above, the focusing of the beam by the
magnetic focusing lens (coils) presents a problem in that not all points on the screen
are in sharp focus:
screen
focusing lens
L
in focus here
L
out of focus at
L = fixed focal length edge
Making the screen too rounded to compensate gives poor images. Thus better CRTs allow
dynamic focusing where the focal length is adjusted corresponding to the position shot at on
the screen – technology used in “flat-screen” CRTs.
observer
1 2 3 4 5 6
Here, ambient or internal background light is reflected from a reflecting panel (1) to reach an
observer (O) by passing through a horizontal polarizing panel (2) with an associated horizontal
wire-grid panel which is used to align the molecules in a liquid crystal panel (4) together with a
vertical wire-grid panel (5) with an associated vertical polarizing panel (6). By setting the
current on/off on the vertical and horizontal grids, particular ‘pixel’ positions in panel 4 may
pass or block the light. Colour may be achieved by using special liquid crystal material for the
panel 4. The picture definition is usually held in a frame buffer in order to refresh the display,
since random fluctuations in the alignment of the liquid crystal molecules in panel 4 degrade
the images. For more information see the Wikipedia article:
http://en.wikipedia.org/wiki/Liquid_crystal_display
A plasma panel consists of a glass panel which encloses a gas that can be ionized, such as
neon or a combination of neon and others (for colour). The ionization is achieved by applying
firing voltages in the X- and Y-directions by conducting grids which are at right angles to each
other, thus ionizing a particular gas cell (i.e. a pixel), which then emits light. The picture
definition is held in a refresh buffer, and the refresh rate is typically 60 frames per sec.
observer
vertical grid
horizontal grid gas (x address)
(y address)
Other devices
Include thin-film electroluminescent displays, light-emitting diode (LED) panels.
PERIPHERAL
CPU DEVICES
system bus
monitor
The video controller cycles through the frame buffer (at ± 60 frames per sec), extracting pixel
data, one scan line at a time by generating memory addresses in synchronous with the scan.
Then the correct pixel intensity values are extracted via the FB port and fed via the VC to the
monitor as intensities and deflection voltages. Where there is no separate FB port for the VC,
the pixel data values are accessed via the system bus.
Raster scan
generator horizontal
+ vertical
deflection
X register Y register voltages
colour/
Frame buffer pixel register intensity
values
Operation:
Assume FB addresses range in X: 0...xmax and in Y: 0...ymax and that the scan coordinate
system is as follows:
ymax
start scan
here
(0,0) xmax
Two registers in the VC (X reg, Y reg) holding screen coordinates are initially set to
X = 0, Y= ymax
11
Then for the pixel (x,y) the memory address for the corresponding data in the FB is formed and
the data values are extracted. These are then filled into the pixel register as intensity/colour
values which are fed to CRT as control voltages for the RGB guns, thus lighting the (x,y)
position on the screen. Then X is incremented by 1 and the FB accessed again etc until
X=xmax is done, after which,
X set to 0, Y decremented by 1,
and the FB is accessed for the next scan line. When done for the scan line Y=0, the first cycle
is complete and Y is reset to ymax, and the process repeats as the next cycle.
The above process is simplistic and generally inefficient, since it assumes that we can access
memory for each pixel without delay. Since the screen is refreshed ±60 times per sec, we
require for example, for a 640 × 480 pixels with 1-bit data, an access time of 1/(640 × 480 × 60)
= 54 nanosecs. But typical RAM cycle times are much greater (~200) so we cannot access
memory every 54 nanosecs. Thus we require a VC-FB architecture that can allow access to
multiple pixel data – perhaps stored in a separate register that gets shifted out to the CRT
periodically. However, tests show that sharing memory cycles between the VC and CPU slows
down the CPU.
Hence, various other design solutions have been proposed and implemented:
i) double buffering – 2 FBs used – 1st fills, 2nd emptied by the VC, then 2nd fills, 1st
emptied out and so on
ii) use of a video look-up (LU) tabel in the VC
iii) separate graphics processor and FB
iv) advanced integrated CPU+ graphics co—processor
For i) see later, but let’s consider ii) and iii) below.
A look-up table has as many entries as there are pixel intensity values available. For example,
for a 3-colour (RGB) system with 8 bits per pixel for data there would be 28=256 entries. The
FB at any (x,y) address ~ a coordinate position on the screen holds just an index into the table
{index range 0...256] i.e. the FB content is not used to control the intensity/colour directly.
LU table intensities
FB memory at (x,y)
cells
255 1.. 0.. .. 11
1001 R
~
(x,y) 01000011 67 1001 1010 0001 1010 G
pixel 0001 B
2
1
0 0.. 0.. .. 00
This mechanism is convenient and fast when a small colour range is required, since only the
colour index need be stored in the FB for each screen coordinate.
12
system PERIPHERAL
CPU memory DEVICES
system bus
display processor
monitor
Apart from doing simple raster operations (setting pixel positions, scan converting lines, filling
etc) a DPU can be designed to do complicated operations like windowing, coordinate
transformations, 3-D imagery etc. Additional functionality such as editing FB contents before
displaying, mixing video images with FB contents etc. may also be available. The trade-off is
the expense involved as opposed to the slower but inexpensive use of the CPU for some of
these operations.
0
1
Map to device
coords
(xdc,ydc,zdc)
Model each item in
modelling coords
(xmc,ymc,zmc) monitor
Typically, when using any of the above one would to install them (if not already packaged with
the OpenGL suite) and then include the respective header files like:
#include <windows.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h> //usually omit above two if this is used
#include <cmath> // or old <math.h> etc for usual C/C++ libs
....... C/C++/OGL code
.................................
For a guide on installing and running the DevC++ compiler suite with OpenGL see the course
website.
2.6.3 GLUT display window
Setting up a somewhat minimal display window using GLUT requires:
• call to initialization routine: glutInit(&argc, argv); //args not always used
• call to create window: glutCreateWindow(“ window title here”);
• call to set up its display content, which can be specified by a graphFunc which defines the
graphics to be created: glutDisplayFunc(graphFunc); //graphFunc is yours
• call to start process, putting display contents in window, puts window program into infinite
loop, looking to process window “events”: glutMainLoop();
This should be the last function to be called.
• The window will be located in some default position and size.To use your own we call
glutInitWindowPosition(50,100); //top left corner is at (50,100)
glutInitWindowSize(400,300); //400 pixels wide, 300 high
• Other attributes may be set by calling e.g.
glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
This allows for a single refresh buffer and RGB colour mode. Mode values can be
combined with the logical and operator.
100
s
c
window title here
r
e 50
e
n 300
400
display window
//lines.cpp
//---------
//Draw lines with OpenGL + GLUT
//Compiler DevC++ 4.9.9.2 + GLUT3.7 (Win32)
#include <windows.h>
#include <GL/glut.h>
void init(void)
{
glClearColor(1.0,1.0,1.0,0.0); //display window set to white
glMatrixMode(GL_PROJECTION); //projection parameters
gluOrtho2D(0.0,200.0,0.0,150.0); //sets up WC extent
}
void lineSegments(void)
{
glClear(GL_COLOR_BUFFER_BIT); //clears display window
glColor3f(0.0, 0.0, 1.0); //line colour set to blue
glBegin(GL_LINES);
glVertex2i(180,15); glVertex2i(10,145); //line seg1
glVertex2i(145,10); glVertex2i(15,180); //line seg2
glEnd();
glFlush(); //process all OGL functions immediately
}
int main(int argc, char** argv)
{
glutInit(&argc, argv); //initialise GLUT
glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); //sets display mode
glutInitWindowPosition(50, 100); //top left display window pos
glutInitWindowSize(400, 300); //display win width and height in pixel coords
glutCreateWindow("Example01 - lines"); //now create display win
init(); //run initilization procs
glutDisplayFunc(lineSegments); //call drawing func
glutMainLoop(); //display all and wait
return 0;
}
16
A video screen coordinate system is comprised of an integer array of pixel positions (x,y),
with typical pixels shown as centred dots, starting from (0,0) at the lower left to some (xmax,
ymax) at the top right where the y-coordinate is the scan line. In practical systems, scanning
starts from the top left taken as (0,0) down to the bottom right taken as (xmax, ymax), but software
function arguments can be set to convert to these trivially from the usual coordinate system
shown below:
y
= a pixel
.
.
7
6
5
4
3
2
1
0
0 1 2 3 4 5 6 7 8 ... x
Typically, frame buffer values can be set at the (x,y) coordinates by a call to a low-level function
like,
setPixel(x,y);
Conversely, frame buffer values for the pixel at (x,y) can be retrieved by a call to a function
such as,
getPixel(x,y,color);
where color is an integer (or other data type) corresponding to the RGB colour combination set
for the pixel (x,y)..
For 3D scenes, an additional 3rd, a depth (or z) coordinate (relative to a viewing position), may
be specified.
18
We note that the coordinate values referred to above are called absolute coordinate values
(in the particular system used). Some packages employ relative coordinates, which are
offsets relative to a current position.
general 2D space
ymax
ymin
Note that when the picture constituents are defined, they must be given in absolute (world)
coordinates. Thereafter the OGL/GLUT display functions may be constructed and called. The
function call sequence above will map these to screen coordinates as on p. 15: thus xmin will be
mapped to an integral screen or pixel coordinate value and so on.
Examples:
(i) glBegin(GL_POINTS); 300
glVertex2i(50,100);
250
glVertex2i(75,150);
200
glVertex2i(100,200);
glEnd(); 150
100
(iii) glBegin(GL_POINTS);
glVertex3f(-51.25,101.56,10.34);
glVertex3f(234.61,-1213.75,170.25);
glEnd();
wcPt2D pointPos;
pointPos.x = 65.75;
pointPos.y = 32.37;
glBegin(GL_POINTS);
glVertex2f(pointPos.x,pointPos.y);
glEnd();
(v) A C++ function can be constructed to implement a setPixel() function using the
OGL point functions above (see later).
20
3.3.2 Lines
Many graphics packages provide a function for drawing a single straight line or many line
segments. In OGL, we use the constant GL_LINES in the glBegin() function segment
together with the vertex functions to define line end points. Line segments are then drawn by
taking successive pairs of vertices. Any vertex not paired with another is ignored. Other
primitive-constants are also available.
//A polyline
glBegin(GL_LINE_STRIP); p3
glVertex2iv(p1);
glVertex2iv(p2); p2
glVertex2iv(p3);
glVertex2iv(p4); p5
glVertex2iv(p5);
glEnd();
p1 p4
p1 p4
21
yend − y0
yend m= (3.2)
xend − x0
y0 b = y0 − m. x0 (3.3)
x0 xend
The algorithms are then based on equations (3.1) – (3.3) by noting that:
For given x -interval δ x along OX, the corresponding interval along OY is
δ y = m.δ x. (3.4)
Similarly, corresponding to an interval δ y along OY is, along OX, the interval
1
δ x = δ y. (3.5)
m
In analogue systems, deflection voltages are applied corresponding to equations (3.4) – (3.5),
namely,
when m < 1 , δ x is set proportional to a small horizontal deflection voltage and
δ y is calculated from (3.4).
But, when m ≥ 1 , δ y is set proportional to a small vertical deflection voltage and
δ x is calculated from (3.5).
In raster systems, lines are plotted with pixels, by stepping in units of δ x or δ y but
constrained by the pixel (grain) separation. That is, we sample the line at discrete positions and
find the nearest pixels to these positions according to the situation below:
22
sample along OY ↑
m ≤1
m >1
← calculate the y-values
along OX→ calculated pixel
calculated pixel
Quiz: Why should we not simply sample along just one coordinate (e.g. OX) and compute
along the other in both cases?
Here (3.6) – (3.7) imply L-to-R processing. For R-to-L processing we use
m ≤ 1: yk +1 = yk − m (3.8)
1
m > 1: xk +1 = xk − (3.9)
m
19 19
18 18
line path line path
17 17
16 16
15 15
14 14
13 13
12 12
11 11
10 10
10 11 12 13 14 15 16 17 18 ... x 10 11 12 13 14 15 16 17 18 ... x
y = mx + b yk +1
yk +3 dU
y
yk +2 yk dL
y k +1
yk
xk +1
10
xk xk+1 xk+2 xk+3
24
Now consider the y-coordinate value where the line xk +1 meets the mathematical line, and let
dU and d L (both > 0) be the vertical separations from y of yk +1 and yk respectively. Then
with xk +1 = xk + 1 we have
y = m ( xk + 1) + b (3.10)
=> d L = y − yk = m( xk + 1) + b − yk
and, dU = yk +1 − y = yk + 1 − m( xk + 1) − b
Hence the difference between the two separations is
d L − dU = 2m ( xk + 1) − 2 yk + 2b − 1 (3.11)
To make the right choice conveniently, we introduce a decision parameter pk for the k-th step,
involving only integer calculations:
∆y vertical endpoints separation integer
Noting that m = = =
∆x horizontal endpoints separation integer
define pk = ∆ x ( d L − d U ) .
Then pk = ∆x ( d L − dU ) = 2∆y. xk − 2∆x. yk + c , (3.12)
where c = 2 ∆y + ∆x (2b − 1) = an integer constant
with sign ( pk ) = sign ( d L − dU ) ,
since ∆x > 0 in this case.
0, if pk < 0
Thus we take pk +1 = pk + 2∆y − 2∆x × (3.14)
1, if pk > 0
Now if ( x0 . y0 ) is the starting point then,
p0 = 2∆y. x0 − 2∆x. y0 + 2∆y + 2∆x.b − ∆x
(3.15)
= 2∆y − ∆x + 2∆y. x0 − 2∆x( y0 − b) = 2∆y − ∆x
k pk ( xk +1 , yk +1 )
0 6 (21,11) y
1 2 (22,12)
2 -2 (23,12) 18
3 14 (24,13) 17
4 10 (25,14) 16
5 6 (26,15) 15
6 2 (27,16) 14
7 -2 (28,16) 13
8 14 (27,17) 12
9 10 (30,18) 11
10
10
20 21 22 23 24 25 26 27 28 29 30 31 x
The selected pixels are shown on the RHS above with the red line indicating the theoretical
(mathematical) line.
For an implementation of this algorithm, for | m |< 1 see H&B (p.98). To generalize it to lines of
arbitrary slope we note:-
i) Use symmetry between octants and quadrants to economize on calculations
ii) For slopes > 1 (+ve) step along OY and calculate successive x-values
iii) To plot from R-to-L, both x and y decrease (i.e. steps ∆ = − ve ) when slope is +ve.
To get the same pixels when processing L-to-R and R-to-L choose, during the entire
process, one of the upper or lower of the two candidate pixels whenever d L = dU .
For slopes < 0 note that one coordinate decreases whilst the other increases.
iv) For the special cases ∆y = 0 (horizontal line) or ∆x = 0 (vertical) we can load the
FB without any processing.
26
( x − xc ) + ( y − yc ) = r 2 = const
2 2
( x, y )
yc r
xc
However, this is a poor method, resulting in uneven pixel distribution and involves too many
computations in terms:
dense
sparse
( x, y )
r
yc θ
x = xc + r cosθ
,0 ≤ θ ≤ 2π (3.16)
y = yc + r sin θ
xc
calculate these
45o
use symmetry
However, employing (3.16) involves too many calculations (trigonometric functions). The
following adaptation of Bresenham’s method is faster and more efficient.
calculate these
(x,y)
45o
use symmetry
(0,0)
Thus we take the standard octant centred at (0,0) and start at the point (0,r). With this point
assumed done, we shall take unit steps along OX and by constructing decision parameters as
we did for lines, we determine the nearest of the two pixels to the circle path according to the
calculated y-values.
First, define the circle function
fc = x2 + y2 − r 2 (3.17)
Now, if ( x, y ) lies on the circle boundary, then it satisfies
f c ( x, y ) = 0.
But, if ( x, y ) lies outside then f c ( x, y ) > 0 etc. Thus for any ( x, y ) we have
0, if ( x, y ) on boundary
f c ( x, y ) = < 0,if ( x, y ) inside (3.18)
> 0,if ( x, y ) outside
Now assume that ( xk , yk ) has been done and we want the next pixel i.e. whether,
the pixel at position ( xk +1 , yk ) is closer or
the pixel at position ( xk +1 , yk −1 ) is closer:
28
yk U
yk −1 L
xk xk +1 xk + 2
midpoint line between pixel candidates U and L
To answer this question, evaluate the circle function at the midpoint y-value i.e define
2
1
pk f c xk +1 , y 1 = ( xk +1 ) + yk − − r 2
2
(3.19)
k−
2 2
at the k-the step.
Now, if
i) pk < 0 then the midpoint is inside and the pixel at scan line yk is closer
ii) pk ≥ 0 then the midpoint is outside (or on) and the pixel at scan line yk −1 is
closer.
To continue the process we obtain successive decision parameters recursively, consider the
decision parameter at the k+1st step:
2
1 1
pk +1 = f c xk +1 + 1, yk +1 − = [( xk +1 ) + 1] + yk +1 − − r 2
2
2 2
= .... (3.20)
= pk + 2( xk + 1) + ( yk2+1 − yk2 ) − ( yk +1 − yk ) + 1
where,
if pk < 0 ,then take yk +1 = yk
else if pk ≥ 0 , take yk +1 = yk − 1 .
i.e.
pk +1 = pk + ∆pk
2 xk +1 + 1, if pk < 0
(3.21)
where ∆pk = 2 xk +1 + 1 − 2 yk +1 , if pk ≥ 0
with 2 x = 2 x + 2,2 y = 2 y − 2
k +1 k k +1 k
We start the process at (0, r ) where ( x0 , y0 ) = (0, r ) with the initial decision parameter,
1
p0 = f c (1, r − )
2
2
1 5
i.e. p0 = 1 + r − − r 2 = − r ....from (3.19) (3.22)
2 4
If r = integer, we round p0 to p0 = 1 − r .
29
Remark: Since all increments are integers, and all calculations in the above are integer
calculations (i.e. no √ ‘s or trigonometric functions!), the method is fast and efficient. To recap
precisely we have the:
Example: Obtain and plot the first quadrant pixels for the circle with centre (0, 0) and radius
r = 10 .
Solution:
With the starting point ( x0 , y0 ) = (0,10) we use the initial decision parameter
p0 = 1 − r = −9 and 2 x0 = 0, 2 y0 = 20 .
The results are tabulated below, with the pixels obtained on the right. The red curve is the arc
of the mathematical circle.
k pk ( xk +1 , yk +1 ) 2 xk +1 2 yk +1 y
0 -9 (1,10) 2 20 10
1 -6 (2,10) 4 20 y=x
9
2 -1 (3,10) 6 20 8
3 6 (4,9) 8 18 7
4 -3 (5,9) 10 18 6
5 8 (6,8) 12 16 5
6 5 (7,7) 14 14 4
3
(● = calculated, ○ = by symmetry) 2
1
0
0 1 2 3 4 5 6 7 8 9 10 11 x
The selected pixels are shown on the RHS above with the red line indicating the theoretical
(mathematical) line. For a sample code see H&B p. 108 (also circle.cpp).
30
P( x, y )
d1 d2
F1
F1
d1 + d2 = constant
( x − x1 ) + ( y − y1 ) + ( x − x2 ) + ( y − y2 ) = const
2 2 2 2
only if (3.23)
i.e. we can write the equation for an ellipse as
Ax 2 + By 2 + Cxy + Dx + Ey + F = 0 (3.24)
where A, B, C , D, E , F are constants that depend on the foci and two other constants called
the semi-minor and semi-major axes.
For an ellipse in “standard form” with centre ( xc , yc ) the semi-major axis length ( rx ) and the
semi-minor length ( ry ) are measures along axes parallel to the coordinate axes, through
( xc , yc ) :
ry
yc θ
rx
xc
It can then be shown that the equation for the ellipse may be written as
2 2
x − xc y − y c
+ = 1 (3.25)
rx ry
or, in polar coordinates as:
x = xc + rx cos θ
, 0 ≤ θ ≤ 2π (3.26)
y = yc + ry sin θ
Note here that symmetry exists between quadrants, but not between octants.
31
y
region 1
ry
slope = -1 ( − x, y ) ( x, y )
region 2
x
rx
( − x, − y ) ( x, − y )
We start by considering the 1st quadrant split into two regions (region 1 and region 2),
separated by the dividing line “----“, which meets the tangent line where the slope = -1.
Now, beginning with the point (0, ry ) take unit steps
along OX whenever |slope| < 1, or
along OY whenever |slope| > 1,
For example as the above indicates, with rx < ry move clockwise until we reach the point
where slope = -1, then switch to steps along OY.
As before, at each sampling position select the next pixel along the path according to the sign
of the function f e ( x, y ) evaluated at the midpoint between the two candidate pixels.
Now, starting from (0, ry ) , we proceed in unit steps along OX, checking at each step whether
the boundary between region 1 and region 2 has been reached by testing the slope condition:
dy 2r 2 x
≡ − y2 = −1
dx 2 rx y
i.e. whether 2ry2 x = 2 rx2 y
Thus, we move out from region 1 whenever
2ry2 x ≥ 2 rx2 y (3.29)
32
We generalize the process, as we had done for a circle, by considering the midpoint values of
f e ( x, y ) as follows:
Assume that ( xk , yk ) has been done and we want the next pixel i.e. whether,
the pixel at position ( xk +1 , yk ) is closer or
the pixel at position ( xk +1 , yk −1 ) is closer:
yk U
yk −1 L
xk xk +1 xk + 2
midpoint line between pixel candidates U and L
To this end, evaluate the ellipse function at the midpoint y-value i.e define
1
2
To continue the process we obtain successive decision parameters recursively, so for the k+1st
step:
2
1 2 1
pk +1 = f e ( xk +1 + 1, yk +1 − ) = ry ( xk +1 + 1) + rx yk +1 − − rx2 ry2
2 2
2 2
(3.31)
2
1
2
1
2
= pk + 2 ry ( xk + 1) + ry + rx yk +1 − − yk −
2 2
2 2
where,
if pk < 0 ,then take yk +1 = yk
else if pk ≥ 0 , take yk +1 = yk − 1.
i.e.
pk +1 = pk + ∆pk
2 ry xk +1 + ry , if pk < 0
2 2
2 (3.32)
where ∆pk = 2 ry xk +1 + ry − 2 rx yk +1 , if pk ≥ 0
2 2
2
with 2 ry xk +1 = 2 ry xk + 2ry ,2rx yk +1 = 2rx yk − 2 rx
2 2 2 2 2
1
p0 = f e (1, ry − )
2
2
1 1
i.e. p0 = r + r ry − − rx2 ry2 = ry2 − rx2 ry + rx2
2
y x
2
....from (3.30) (3.32)
2 4
At this ( x0 , y0 ) = (0, ry ) the terms 2ry x and 2rx y are
2 2
2ry2 x = 0 (3.33)
2r y = 2 r r
x
2 2
x y (3.34)
For region 2, start at ( x0 , y0 ) = last position selected in region 1 and sample in unit steps in
the –ve y-direction and take the midpoint x-value as the line between xk and xk +1 :
yk
yk −1 L R
xk xk +1 xk + 2
midpoint line between pixel candidates L and R
Now, if
i) pk′ > 0 then the midpoint is outside and the pixel at xk is closer
ii) pk′ ≤ 0 then the midpoint is inside (or on) and the pixel at xk +1 is closer
2
1 1
2
For ellipses in non-standard positions we map the selected points by shifts (translations) and/or
rotations (see later).
Remark: Again here all calculations in the above are integer calculations so that the method is
fast and efficient.
2
5. Then for k = 0,1, 2,3,... (in region 2)
If p′k > 0 , then take next point as ( xk , yk − 1)
and set pk′ +1 = pk′ − 2rx2 yk +1 + rx2
else if p′k ≤ 0 , take next point as ( xk + 1, yk − 1)
35
Example: Obtain and plot the first quadrant pixels for the ellipse with centre (0, 0) and radii
rx = 8, ry = 6 .
Solution:
With the starting point ( x0 , y0 ) = (0, 6) we obtain
2ry2 x = 0 (with increment 2 ry2 = 72 )
2rx2 y = 2 rx2 ry (with increment −2 rx2 = −128 )
Then for region 1: The initial decision parameter is
1
p0 = ry2 − rx2 ry + rx2 = −332
4
Then we tabulate successive values as follows:
k pk ( xk +1 , yk +1 ) 2 ry2 xk +1 2 rx2 yk +1
0 -332 (1,6) 72 768
1 -224 (2,6) 144 768
2 -44 (3,6) 216 768
3 208 (4,5) 288 640
4 -108 (5,5) 360 640
5 288 (6,4) 432 512
6 244 (7,3) 504 384
2
Tabulated values for this region are then
k p′k ( xk +1 , yk +1 ) 2 ry2 xk +1 2 rx2 yk +1
0 -151 (8,2) 576 256
1 233 (8,1) 576 128
2 745 (8,0) _ _
The computed pixels nearest the ellipse boundary are shown in figure below.
36
10
9
8
7
6
5
4
3
2
1
0
0 1 2 3 4 5 6 7 8 9 10 11 x
The selected pixels are shown on the RHS above. For a sample code see H&B p. 116 (or
ellipse.cpp).
= extra point
y
18
17
16
15
14
13
12
11
10
20 21 22 23 24 25 26 27 28 29 30 31 x
37
This can be corrected by plotting only those points inside the endpoints (20,10) to (30,18), i.e.,
in this case leave out the last (top) pixel.
Similar effects are obtained for other shapes. For example, for the enclosed rectangle (0,0),
(4,0), 4,3), (0,3) we have:
y y y
5 5 5
4 4 4
3 3 3
2 2 2
1 1 1
0 0 0
0 1 2 3 4 5 x 0 1 2 3 4 5 x 0 1 2 3 4 5 x
15 15
14 14
13 13
12 12
11 (10,10) 11 (10,10)
10 10
9 9
8 8
7 7
6 6
5 5
x 4 5 6 7 8 9 10 11 12 13 14 15 x
4 5 6 7 8 9 10 11 12 13 14 15
Polygons are useful for approximating curved surfaces (‘surface tessellation’). Figures can be
constructed with them, as wire-frame models, after which it is easy to fill or render the surfaces
so that they appear realistic.
ii) Start from a point inside and paint outward until boundaries are met – used for
complicated shapes and in interactive painting software.
all
angles > 180o
< 180o
Convex polygons present no problems when designing filling algorithms, but concave ones
have to be recognised as such and then are typically split into two or more convex ones first.
Degenerate ones, however, require special treatment.
v6 v5
E5
E4
v4
E6 E3
v3
E2
E1
v1 v2
so the polygon is concave. Then for those pairs for which the z-component < 0, we split the
polygon, along the line of the first vector in the pairs. Hence here, we split along E3 .
Example:
3
E5
2 E6
E4
1
E2
E3
E1
0 1 2 3
Here
E1 = (1,0, 0), E2 = (2,1, 0) − (1,0, 0) = (1,1,0), E3 = (1, −1, 0),
E4 = (0, 2, 0), E5 = ( −3, 0,0), E6 = (0, −2, 0).
Then,
( E1 × E2 ) z = (0,0,1) z > 0,( E2 × E3 ) z = (0,0, −2) z < 0, ( E3 × E4 ) z = (0, 0, 2) z > 0,
( E4 × E5 ) z = (0,0, 6) z > 0,( E5 × E6 ) z = (0,0, 6) z > 0,( E6 × E1 ) z = (0, 0, 2) z > 0.
Thus, the polygon is concave and we split along the line E2 i.e. along the line with slope = +1
and y-intercept = -1 (confirm this). The polygon splits into 2 convex ones, but since there is no
other vector pair for which (...)z < 0, no further splits are necessary.
Other methods (e.g. the rotational method) are also used - see H&B p.126.
Count the number of edges crossed when edges around the test point P are
traversed counter-clockwise.
If any edge is crossed R-to-L set w = w+1
If any edge is crossed L-to-R set w = w-1
Then
• if w ≠ 0, take P as an interior point
• else (w = 0) P is exterior
Examples:
exterior exterior
interior interior
We see that these methods can give different results for non-standard polygons.
To determine the direction of the edge crossings (R-to-L or L-to-R) we use one of:-
i)
A
d
U
P EAB
B
ii)
A Up
d
┐
U
P
EAB
B
Here, we let Up denote the vector from B traversing the line vector U perpendicularly
i.e. it is a vector ┴U in the direction R-to-L when looking along U towards d. For
example if,
U = (ux,uy) then Up = (-uy,ux).
Now form the dot product, Up • EAB.
Then
• if Up • EAB > 0 the edge EAB is R-to-L
• if Up • EAB ≤ 0 the edge EAB is L-to-R
In general, regions or objects can be multiply connected – so we can use the above rules to
determine the interior/exterior or simply define which is which.
E3 E6
E1
S1
S2
V5
E2 V3
E4 E5
V2
V4
Other forms of the table may be used, such as one with an edge table that includes forward
pointers into the table of surface facets. This makes data checking and retrieval more efficient.
where ( x, y , z ) is a point on the plane and the A, B, C , D are fixed parameters identifying a
particular plane.
To solve for the parameters from (3.40), we must solve 4 linearly independent equations. It’s
customary to take D as an arbitrary ≠ 0 constant and then solve for A/D, B/D, C/D from
( A / D ) xk + ( B / D ) yk + (C / D ) zk = −1; k = 1, 2, 3 (3.41)
using 3 non-colinear points, such as, for example, the vertices ( xk , yk , zk ), k = 1, 2, 3 . Thus
setting
x1 y1 z1
D = − x2 y 2 z 2 (always ≠ 0 for any 3 distinct points) (3.42)
x3 y3 z3
we obtain
1 y1 z1 x1 1 z1 x1 y1 1
A = 1 y 2 z 2 , B = x2 1 z 2 , C = x2 y 2 1 (3.43)
1 y 3 z3 x3 1 z3 x3 y3 1
⇒
A = y1 ( z2 − z3 ) + y2 ( z3 − z1 ) + y3 ( z1 − z2 )
B = z1 ( x2 − x3 ) + z2 ( x3 − x1 ) + z3 ( x1 − x2 )
(3.44)
C = x1 ( y2 − y3 ) + x2 ( y3 − y1 ) + x3 ( y1 − y2 )
D = − x1 ( y2 z3 − y3 z2 ) − x2 ( y3 z1 − y1 z3 ) − x3 ( y1 z2 − y2 z1 )
These relations for the plane parameters are also valid for planes passing through (0,0,0) i.e.
for D = 0 . They get updated as vertex and other information are updated.
In addition to the above it is also required to determine for a polygon surface or plane the
orientation of a plane surface. In order to do this we employ a normal vector to it:
43
Note from (3.40), that if ( A, B, C ) is the vector of the plane’s parameters A, B, C and its 3
successive vertices are ( xk , yk , zk ), k = 1, 2, 3 then the following dot products hold
( A, B, C ) ( x1 − x2 , y1 − y2 , z1 − z2 ) = 0
( A, B, C ) ( x1 − x3 , y1 − y3 , z1 − z3 ) = 0 . (3.45)
( A, B, C ) ( x2 − x3 , y2 − y3 , z2 − z3 ) = 0
N = ( A, B, C ) , (3.46)
In many situations we need to know the “inside” and “outside” of a surface. To establish this we
apply the right-hand rule:
If the vertices bounding the surface are specified in a counter-clockwise direction, then N
points out from the top (which is defined as the outside or frontface), whilst −N points away
from the bottom (taken as the inside or backface).
N outside
V2
V1
V3 -V1
V3
V4
inside -N
Notice from the above that, another way to calculate a normal to the surface is to take the
cross-product of 2 vectors in the plane as follows (gives the correct sign for convex polygons).
Select 3 vertices along the boundary in the counter-clockwise direction e.g. V1,V2,V3 then
N = ( V2 − V1 ) × ( V3 − V1 ) (3.47)
Expanding the RHS into components then gives (check this!)
N = ( A, B, C ) (3.48)
44
namely, the plane constants! Then D can be found by substituting A, B, C from above and
one set of vertex coordinates into (3.40). Notice also, that using (3.48) the equation for the
plane (3.40) is simply,
N P = −D (3.49)
where P = ( x, y, z ) is an arbitrary point on it.
Finally, to test any point ( x, y , z ) in space, in a RH coordinate system, we observe that
Ax + By + Cz + D = 0 ⇒ point on surface
Ax + By + Cz + D < 0 ⇒ point on inside (3.50)
Ax + By + Cz + D > 0 ⇒ point on outside
scan
line
10 14 18 24
Strategy: For each scan line, find points of intersection with polygon edges, sort points L-to-R
and fill between each pair of points. For example for scan line above,
• fill between x=10 and x=14
• fill between x=18 and x=24
However, intersections with polygon vertices require special treatment. For example consider,
incr
decr
decr
1 1 y'
2
2
y
1 1 1
decr
Here,
45
• the scan line y' intersects 2+1+1 = 4 edges, passing through 1 vertex
• the scan line y intersects 2+1+1+1 = 5 edges, passing through 1 vertex.
Note that
• for the y' intersections the “interior” pixels are clear and unambiguous
• for the y intersections the “interior” pixels are not clear
and also that,
• for y', at common vertex 2 one edge is decreasing and the other increasing, i.e. 2
is a local minimum – so count this point as two vertices and fill between 2 pairs.
• for y, at common vertex 2 one edge is monotonically decreasing and so is the
other, i.e. 2 is a not a local min/max – so count this point as one vertex point
belonging to one edge only, and now can fill correctly between 2 pairs again.
To resolve polygon vertices as single or double ones, we follow the procedure below:
1. Process non-horizontal edges clockwise (or anti-clockwise)
2. For each intersection with a vertex check whether one edge has decreasing y-value at
its endpoint (wrt another before it) and the next edge has decreasing y-value at this
point (wrt another after it).
3. If true in 2, shorten the lower edge by 1 pixel y-value, ensuring only 1 point of
intersection with a scan line (~ y scan line case above).
If false in 2, leave vertex as 2 intersection points (~y' scan line case above).
In summary we have:
y increasing from e to e'
y decreasing from e' to e
shorten by 1
e' scan line e'
scan line
e y e y+1
y-1 y
shorten by 1
The calculations can be reduced by taking into account coherence properties of a scene, i.e.
the relation between one part of it relative to the next part. For example, for the case,
( xk +1 , yk +1 ) scan line yk +1
( xk , y k ) scan line yk
∆x
xk +1 = xk + .
∆y
For example:
Take an edge with slope m =7/3.
At the initial scan line set counter = 0 and here counter increment is ∆x = 3
Move up to next 3 scan lines giving counter values 3,6,9
On 3rd scan line counter > 7 ⇒
x-intersection coordinate increases by 1
counter reset to 9-7 = 2 (since now are 2 scan lines up)
Continue as above until we reach upper endpoint of the edge
(For negative slopes the process is similar)
In the following efficient algorithm implementing this scan line polygon fill we first store the
boundary as a sorted edge table. Then we proceed clockwise (or anti-clockwise) only, around
the edges, and use a bucket sort to store the edges sorted according to the smallest y-value of
each edge, in the correct scan line position. Only non-horizontal edges are stored in the table.
When the edges are processed, any intersection with a duplicate vertex point is resolved, by
shortening one edge as outlined above.
The data structure employed is shown below for a simple example. In the table, each entry for
a given scan line contains the maximum y-value for that edge and the x-intercept value (at its
lower vertex) and the inverse slope for that edge. For each scan line the edges are sorted in
the L-to-R order. The scan lines are processed from the bottom upwards, to obtain an active
edge list for each scan line intercepting the polygon boundaries. For more details see H&B p.
200.
47
scan line
number yC yB xC 1/mCB
.
B .
yD yC xD 1/mDC yE xD 1/mDE
.
.
C scan line yC
yA yE xA 1/mAE yB xA 1/mAB
C' E .
.
D
scan line yD 1
scan line yA 0
A
else
j = k + 1;
while (pts[k].y == pts[j].y)
if ((j+1) > (cnt-1))
j = 0;
else
j++;
return (pts[j].y);
}
/* Store lower-y coordinate and inverse slope for each edge. Adjust
and store upper-y coordinate for edges that are the lower member
of a monotically increasing or decreasing pair of edges */
void makeEdgeRec
(dcPt lower, dcPt upper, int yComp, Edge * edge, Edge * edges[])
{
edge->dxPerScan =
(float) (upper.x - lower.x) / (upper.y - lower.y);
edge->xIntersect = lower.x;
if (upper.y < yComp)
edge->yUpper = upper.y - 1;
else
edge->yUpper = upper.y;
insertEdge (edges[lower.y], edge);
}
p = edges[scan]->next;
while (p) {
q = p->next;
insertEdge (active, p);
p = q;
}
}
void fillScan (int scan, Edge * active)
{
Edge * p1, * p2;
int i,ix,iy;
49
p1 = active->next;
while (p1) {
p2 = p1->next;
for (i=(int)p1->xIntersect; i< (int)p2->xIntersect; i++)
setPixel(i, scan);
p1 = p2->next;
}
}
q->next = p->next;
delete p;
}
while (p)
if (scan >= p->yUpper) {
p = p->next;
deleteAfter (q);
}
else {
p->xIntersect = p->xIntersect + p->dxPerScan;
q = p;
p = p->next;
}
}
if (active->next) {
fillScan (scan, active);
updateActiveList (scan, active);
resortActiveList (active);
}
}
/* Free edge records that have been malloc'ed ... */
}
a) 4-connect fill
Test 4 neighbouring pixels (U,D,L,R) and fill.
OK for simple shapes.
b) 8 connect fill
Test the above 4 pixels (U,D,L,R) +
4 diagonal pixels and fill.
OK for more complex shapes.
Note that the 4-connect scheme may fail to correctly fill certain shapes:
fill
area
Then we can paint the interior by replacing a specified interior colour rather than searching for
a boundary colour. Again a 4-connect or 8-connect scheme can be used in an algorithm such
as:
The full coding is similar to the BFA case. Once again, it is more efficient to employ scanline
stacking rather than using a straight 4-connect pattern.
Two representation schemes for storing character fonts are those that:-
i) Use rectangular grid patterns for characters resulting in bitmapped fonts
ii) Use straight and curved line segments to form characters resulting in outline
fonts (as in Postscript)
In i) the character grid is mapped to a frame buffer segment called the font cache – requires
more storage since must allow for each variation (size, format).
In ii), less storage is required but variations have to be generated by manipulating the curve
segment definitions and then scan converting to the frame buffer – more time consuming.
i) glRecti(200,100,59,250) gives
250
200
150
100
50
the result is the same as in i). Note that the vertices are traversed in the clockwise
direction (x1,y1), (x2,y1), (x2,y2), (x1,y2) → ((x1,y1). For anti-clockwise traversal
we use the call glRectiv(vertex2, vertex1); This is required in procedures for
establishing the correct front/back face (see later).
iii) The OGL code with 6 (must be ≥ 3) vertex points, taken counter-clockwise,
will generate a closed polygon such as on RHS:
glBegin(GL_POLYGON);
glVertex2iv(p1); p6 p5
glVertex2iv(p2);
glVertex2iv(p3);
glVertex2iv(p4); p1 p4
glVertex2iv(p5);
glVertex2iv(p6);
glEnd(); p2 p3
iv) This code gives 2 unconnected triangles (note the vertices are now re-ordered!)
glBegin(GL_TRIANGLES);
p6 p5
glVertex2iv(p1);
glVertex2iv(p2);
glVertex2iv(p6); p1 p4
glVertex2iv(p3);
glVertex2iv(p4);
glVertex2iv(p5); p2 p3
glEnd();
54
vii) Here we obtain two quadrilaterals with the code (4 successive vertices form a
quad)
glBegin(GL_QUADS);
glVertex2iv(p1); p1 p4 p5 p8
glVertex2iv(p2);
glVertex2iv(p3);
glVertex2iv(p4);
glVertex2iv(p5);
glVertex2iv(p6); p2 p3 p6 p7
glVertex2iv(p7);
glVertex2iv(p8);
glEnd();
glBegin(GL_QUAD_STRIP);
glVertex2iv(p1); p1 p4 p5 p8
glVertex2iv(p2);
glVertex2iv(p4);
glVertex2iv(p3);
glVertex2iv(p5);
glVertex2iv(p6); p2 p3 p6 p7
glVertex2iv(p8);
glVertex2iv(p7);
glEnd();
55
Remark:
In all the above, the OGL filled area primitives required convex polygons. However, in the GLU
library, concave polygons (with linear boundaries) are allowed. In addition GLU provides
tessellation routines for converting such shapes into triangles, triangle meshes, fans and
straight-line segments. The latter can then be processed by the OGL functions.
z z
1 4 5
6 7
0 y 0 y
1 1
1 2
3
x x
where each pt[ ] index corresponds to a particular vertex label, indicated above.
In addition, we define the 6 faces of the cube, with either glBegin(GL_POLYGON) or
glBegin(GL_QUADS), listing the vertices in counter-clockwise direction when “looking from the
outside”:
void quad(GLint n1, GLint n2, GLint n3, GLint n4)
{
glBegin(GL_QUADS);
glVertex3iv(pt[n1]);
glVertex3iv(pt[n2]);
glVertex3iv(pt[n3]);
glVertex3iv(pt[n4]);
glEnd();
}
void cube( )
{
quad(6,2,3,7);
quad(5,1,0,4);
quad(7,3,1,5);
quad(4,0,2,6);
quad(2,0,1,3);
quad(7,5,4,6);
}
56
Judging from the above, in realistic scenes many function calls would be needed (thousands is
typical). Many shared vertices would be repeatedly defined. To take care of such issues
efficiently, and to reduce the total number of function calls, vertex arrays are employed.
The procedure to follow is:-
1. Enable the vertex-array feature in OpenGL by calling:
glEnableClientState(GL_VERTEX_ARRAY)
2. Set up the location and data format for the vertex coordinates by calling:
glVertexPointer( )
3. Render the scene, processing multiple primitives by calling a function such as:
glDrawElements( )
A sample code for the cube above is (using the previous definition of pt[ ] ):
glEnableClientState(GL_VERTEX_ARRAY);
glVertexPointer(3, GL_INT, 0, pt);
GLubyte vertIndex[ ] = {6,2,3,7,5,1,0,4,7,3,1,5,4,0,2,6,2,0,1,3,7,5,4,6};
glDrawElements(GL_QUADS, 24, GL_UNSIGNED_BYTE, vertIndex);
Note here:
• The 1st call, sets up a state, the vertex-array state, on the client machine
• The 2nd call sets up the location and format of the coordinates of the object’s vertices:
3 = three coordinates for each vertex, GL_INT = data type for each, 0 = offset between
elements (>0 allows for more info to be packed into array, here only coordinates so 0),
pt = vertex array holding coordinates.
• The 3rd call sets up the vertex index array vertIndex[ ] holding the vertex indices as
unsigned bytes, corresponding to the faces to be processed.
• The 4th call processes the cube faces, using 4 vertices at a time: the parameters
GL_QUADS => quadrilaterals, 24 = size of vertIndex or number of indices,
GL_UNSIGNED_BYTE = type of index.
Here,
• width = no. of columns, height = no. of rows of bitShape
• each element of bitShape is either 0 or 1. 1 => corresponding pixel should be lit
57
• x0, y0 = floating point nos ~ defined “origin” (relative to lower left corner of bitShape)
array
• Location in frame buffer where pattern is to be applied is called the current raster
position. The bitmap is placed with its origin (x0,y0) at this position.
• Values xOffset, yOffset (also floating point) used to update frame-buffer current raster
position after bitmap is displayed.
The current raster position may be set with a call to one of the functions
glRasterPos*(..);
similar to the glVertex*(..) functions. For more on this topic see H&B p. 144.
Similarly, a pattern defined in the colour array pixMap may be applied to a block in the frame
buffer with a function call like:
For example, the following displays a pixmap colorShape, of size 128 × 128 with elements
holding RGB colour values of type unsigned byte:
It’s possible to paste the pixmap array into a buffer, such as a depth buffer (stores object
distances from a viewing position) or a stencil buffer (stores boundary patterns for a scene).
To do this set the dataFormat parameter in the above to GL_DEPTH_COMPONENT or to
GL_STENCIL_INDEX.
Various buffers are available in OpenGl for such purposes as stereoscopic viewing (left/right
buffers), animation (double buffering) and combinations of these for use in screen refreshing.
Non-refresh buffers to hold user data can be defined too. A buffer may be selected for storing a
pixmap with a call such as:
glDrawBuffer(buffer);
• To read a block of pixels from the lower left screen coordinates (xmin,ymin) and put
into array with size width and height and with dataFormat and dataType as before
we use
glReadPixels(xmin, ymin, width, height, dataFormat, dataType, array);
• A block of pixels may be copied from a source buffer to a destination buffer by:
glCopyPixels(xmin, ymin, width, height, pixelValues);
where the lower left screen coordinates are (xmin,ymin), the size is width × height
and pixelValues = GL_COLOR, GL_DEPTH or GL_STENCIL to indicate the type of
buffer.
• In the above we must previously select buffers by glReadBuffer( ) for the source and
glDrawBuffer( ) for the destination.
Other routines are available for pixel data processing (see H&B p. 147).
colour
glutBitmapCharacter(font, character);
where font takes a symbolic GLUT constant such as for examle, GLUT_BITMAP_8_BY_13 (a
fixed width font) or GLUT_BITMAP_TIMES_ROMAN_10 (a proportionally spaced font) and
character may be specified as an ASCII code (e.g. 65) or in this instance ‘A’.
With this function each character is displayed with its lower left bitmap corner at the current
raster position. After loading the character into the frame buffer, an offset, the size of a
character width advances the raster position. Thus a text string text[36] can have its characters
consecutively displayed by code such as:
glRasterPos2i(x,y);
for (int k=0; k < 36; k++)
glutBitmapCharacter(GLUT_BITMAP_9_BY_15, text[k]);
Outline characters are generated with the polyline(GL_LINE_STRIP) boundaries and displayed
by:
glStrokeCharacter(font, character);
where font can take values such as GLUT_STROKE_ROMAN (a proportionally spaced font) or
GLUT_STROKE_MONO_ROMAN (a font with constant spacing).
wherein the list may reside on a server. Complex scenes are more easily constructed when
object descriptions are held in display lists.
Creating an OpenGL display list
The command syntax used is:
glNewList(listID, listMode);
....
....
glEndList( );
Here, a +ve integer value is assigned to listID, and listMode takes the value GL_COMPILE
(saves list for later execution) or GL_COMPILE_AND_EXECUTE (to execute the commands as
they are entered in the list, as well as keep for later execution).
Note that when the list is compiled, parameters such as coordinates and colours become set
and cannot be subsequently altered. Thus, commands such as OpenGL vertex-list pointers
should not be included in the list.
The listID is used in the call to execute the list, so to avoid accidental duplication it’s better to
generate an unused ID with the call
listID = glGenLists(1); //assigns an unused +ve integer to listID.
To query whether, a particular integer listID has been used as an ID we can use
glIsList(listID); //returns GL_TRUE or GL_FALSE
Example:
The following code segment sets up a display list to generate a regular hexagon around the
circumference of a circle and then executes the list with glCallList().
For a complete program implementing this code see DispLstHexagon.cpp. Multiple display
lists may also be processed – see H&B p. 152.
60
Any parameter affecting the way in which a primitive is to be displayed is called an “attribute”
parameter. For example, colour and size are fundamental attributes. Attributes are handled in
two ways:
i) We extend the parameter list of the primitive function to include attribute data
ii) We maintain a system list of current attributes via separate routines, which is
checked whenever any primitive function is called, before displaying the primitive.
OpenGl uses the latter approach. Such a system which maintains its own attribute list is
referred to as a state system or state machine. The attribute parameters and other related
information such as the current frame-buffer position are called state variables or parameters.
Colours codes can be stored directly in the frame buffer or put into a separate table; then the
corresponding binary pattern for each pixel is put into the frame buffer. Using for example, 3
bits per table, the number of colours is 8, with colours coded as follows in a table:
Here each bit puts a particular colour gun either on (1) or off (0) in an RGB system. Adding
more bits for each colour, of course increases the colour range. For example, for 6-bits per
pixel, we can use 4 intensity levels for each gun (00, 01, 10, 11), i.e there are 4×4×4 = 64
possible colours for a pixel. For a resolution of 1024 × 1024 a full-colour (24 bits-per-pixel)
RGB system requires 3 megabytes storage for the FB.
When FB space is limited, selected colours from the full-colour set may be held in a look-up
table where each colour is accessed by an index and where only such indices are held in the
FB. Recall section 2.2.1 (p. 11) for an example.
When a monitor is not capable of colour output we can use a table corresponding to shades of
gray. For example, intensity codes for a 4-level gray scale system would be organised as:
61
Note that any intensity near 0.33 would be mapped to int 1 (or binary 01) and so on.
Alternatively, each intensity level may be converted to a voltage level on the output device
directly.
glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
where the 1st parameter specifies a single buffer for the FB, and the 2nd specifies the RGB
mode. Using GLUT_INDEX will allow choice from an index table and GLUT_RGBA the RGBA
(enhanced) mode.
RGB and RGBA modes
The RGBA mode allows for an (optional) additional alpha coefficient to control color blending.
This is useful when primitives overlap over a pixel area, and one primitive can be set to allow
for different transparency (or opacity) levels. The current colour components are set by calling
glColor*(colorComponents);
where * is 3 (for RGB) or 4 (RGBA) followed by f or i or s or d and/or by v as in glVertex*().
Examples:
glColor3f(1.0, 0.0, 0.0); //red
glColor3fv(colorArray); //depends on vector colorArray
glColor4f(1.0, 1.0, 1.0, 1.0); //white with alpha = 1.0. OGL default
glColor3i(0, 255, 255); //for 8-bits per pixel for each color
Color-index mode
OpenGL color-index mode is set with
glIndex*(colorIndex); //colorIndex = +ve no index to table
Example: glIndexi(196); //196 held in current (x,y) FB position
OGL cannot make a LU table, but accesses window system’s. But GLUT allows user to set
values there by,
glutSetColor(index, red, green, blue);
will allow float RGB values (0.0 - 1.0) to be loaded into the table at position index.
OGL and its extension libraries have routines, comprising its Imaging Subset, for setting up
color tables. These can be enabled by the glEnable( ) function. See H&B p. 179.
pixel spans
skipped spaces
We can specify the span length and inter-span spacing by setting a pixel mask = a string of 0’s
and 1’s e.g.
1111000 => dashed line with spans of 4 pixels and spacing of 3 pixels. On a
bilevel system each bit = a bit value loaded into the frame buffer, and the computed bits are
then and-ed with these. Producing dashed lines with a fixed number of pixels can give unequal
lengths depending on the slopes:
normal larger by √2
Line width
The implementation of the line-width attribute depends on the capabilities of the output device.
For example, a heavy line on a video monitor could be adjacent parallel lines, but on a pen
plotter it would require a change of pens.
In raster systems for slopes |m| < 1 In raster systems for slopes |m| > 1
Plot pixels at (x,y) and (x,y+1) for line Plot pixels at (x,y),(x-2,y),(x-1,y),(x+1,y)
width=2. i.e have vertical extension. for line width=4 i.e. horizontal extension
This technique however gives lines of unequal thickness for different slopes.
Line endpoints
64
Can also be used to distinguish styles for thick lines. Some common endpoint styles in graphics
packages are:-
For thick line or polygon joints we can also employ styles such as:-
Another approach to line thickening, for wide lines, is to use filled rectangles for them.
Curve attributes
Methods similar to those for straight lines are employed here. See H&B p. 189.
Line width
Set with the call
glLineWidth(width); //width = 1.0, 2.0, 3.0 etc
Line style
Set with the call
glLineStipple(repeatFactor, pattern);
where,
• pattern = 16-bit integer ~ display pattern (default 0xFFFF = all bits 1 => solid line,
0x00FF => 8 bits OFF and 8 bits ON => dashed line).
• repeatFactor = number of times each bit in pattern is repeated
For example, glLineStipple(1, 0x00FF) will give dashed lines with solid segment length = 8
pixels and spacing 8 pixels. However, before this function takes effect the line-style feature
must have been enabled:
glEnable(GL_LINE_STIPPLE);
glLineStipple(repeatFactor, pattern);
......
glDisable(GL_LINE_STIPPLE); //reduce to default
Example:
The following code can be used to output 3 line graphs, each of a different style.
A fill pattern is usually defined in a rectangular array of colours, one colour for a position in the
array. Another way involves the use of a bit array which acts as a mask to be applied to the
display area. The process of filling an area with a rectangular pattern is known as tiling. For
more details see H&B p. 194.
Clearly for t < 0.5, the background colour predominates. This process is referred to as
soft filling. Other variations of (4.1) are also employed.
The polygon fill pattern will be shown up to and including the boundaries, unless the latter is
specifically set otherwise. Also, by default, the drawn polygon would be filled in solid colour.
The bit order starts from the bottom row up to the top-most row of the patern array (same order
as bitShape function). Then this pattern is enabled across the entire display window, from its
lower left corner to the top right corner. Wherever a specified polygon intersects it, the polygon
would be shown with this fill pattern.
The mask is used to set the current fill pattern with,
glPolygonStipple(fillPattern);
and then the fill routine is enabled with
glEnable(GL_POLYGON_STIPPLE);
and only then the polygon to be filled is constructed.
The pattern filling is disabled with
glDisable(GL_POLYGON_STIPPLE);
In addition, an interpolation fill, like we had for lines may be used. For example, consider the
code:
glShadeModel(GL_SMOOTH);
glBegin(GL_TRIANGLES);
glColor3f(0.0,0.0,1.0);
glVertex2i(50,50);
glColor3f(1.0,0.0,0.0);
glVertex2i(250,250);
glColor3f(0.0,1.0,0.0);
glVertex2i(350,50);
glEnd();
This will produce the shaded triangle (run TriangleShaded.cpp) where the colour is a linear
interpolation of the 3 vertex point colours:
68
Note that the option glShadeModel(GL_FLAT) will result in the single last colour (green).
To obtain different fill and edge/boundary colours we proceed as in the stub example:
For 3-D polygons, the above filling methods may produce gaps along edges. This effect is
known as stitching and is due to the different ways in scan-line filling and the edge-line
drawing algorithms meet at common pixels. One way to eliminate it is to shift the depth values
(distances from the X-Y plane) computed by the fill algorithm to avoid overlapping with the
edge depth values for the polygon. This is done by:-
glEnable(GL_POLYGON_OFFSET_FILL);
glPolygonOffset(factor1, factor2);
where, the 1st function invokes the offset routine for the scan-line filling, and, the 2nd sets the
float parameters factor1 and factor2 used to compute the amount of depth offset by the
formula:
depthOffset = factor1× maxSlope + factor2 × const (4.2)
Here maxSlope = the maximum slope in the polygon and typically factor* = 0.75 or 1.0 (some
experimentation may be needed here).
The previous code segment would then take the form:
Then the interior fill of the polygon is pushed a little back in depth, so that there is no
interference with the depth values of the edges. For other methods to handle stitching see H&B
p. 209.
69
To display a concave polygon we first split it into convex ones, in particular triangles. Then,
when the triangles are filled, but we must remove the common edges for wire-frame rendering.
For example, here,
dashed
lines
not
wanted
we can set a bit flag to GL_FALSE (to indicate a vertex is not connected to another by a
boundary edge) for a particular triangle. Thus for the triangle shown below, we obtain the
situation on its RHS,
v1 v1
v2 v2
v3 v3
glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
glBegin(GL_POLYGON);
glVertex3fv(v1);
glEdgeFlag(GL_FALSE);
glVertex3fv(v2);
glEdgeFlag(GL_TRUE);
glVertex3fv(v3);
glEnd();
Polygon edge flags can also be incorporated into vertex arrays (see H&B p. 210).
4.9 Anti-aliasing
Displayed images on raster systems are generally jagged since the final coordinates are
discrete (integers). This can cause distortion due to low-frequency or under-sampling, resulting
in the so-called aliasing error. For example consider the periodic shape or signal:
70
periodic shape
sampling positions
After sampling on a grid with wider separation (lower-frequency) we obtain the waveform:
One possible way to increase the sampling rate is to use a higher-resolution display. But, the
jaggedness is not eliminated entirely. Also, this requires a larger frame buffer – too large a
frame buffer means that the required 30-60 frames/sec refresh rate cannot be achieved.
Other ways involve varying pixel intensities along the boundaries, which is suitable for raster
systems with > 2 intensity levels (colour or grayscale). Methods used to implement this idea
are:-
• Increase the sampling rate by treating the screen as if it is covered with a finer grid
than the actual one. Then use multiple sampling points across the finer grid to
determine the actual screen pixel intensities i.e. sample at higher resolution but display
at lower resolution – called super-sampling or post-filtering
• Find the pixel intensities by calculating areas of overlap of each pixel with objects to be
displayed, by determining where object boundaries intersect individual pixel
boundaries – called area sampling or pre-filtering
• Shifting the display location of a pixel areas by micro-positioning the electron beam
in relation to the object geometry – called pixel phasing – depends on device’s
hardware capabilities.
For more technical detail on each of the above see H&B p. 215 – 221.
71
Another way involves setting up a colour ramp table, in gradations of the object colour to the
background colour (see H&B p. 222).
Combined group parameters may put onto stack with a logical or ( | ) combined parameter
e.g.,
glPushAttrib(GL_POINT_BIT | GL_POLYGON_BIT);
glPopAttrib( );
with no parameters, since all values put onto the stack are returned to the state machine.
In any graphics package we need to have routines that allow us to move objects from one
position to another or to view scenes from different angles or vantage points. These issues are
resolved by applying mathematical geometric transformations.
Translation
A 2D point ( x, y ) is translated to the point ( x′, y ′) by applying a shift (or translation) vector
(t x , t y ) to it as follows:
x′ = x + t x , y′ = y + t y (5.1)
Note that translation is a rigid body transformation i.e. objects are moved without
deformation or without altering the relative positions of points within the body.
To translate an entire straight line, we can translate the endpoints and then re-draw the line
between the new points.
To translate a polygon we
• translate each vertex
• re-draw its new line segments between the vertex pairs
old
new
P P′
To translate a circle/ellipse we translate the centre, and re-draw the figure. Similar processes
can be used for other figures.
73
class wcPt2D {
public:
GLfloat x, y;
};
void translatePolygon (wcPt2D * verts, GLint nVerts, GLfloat tx, GLfloat ty)
{
GLint k;
To delete the original polygon we can display it in the background colour before translating it.
Can use other means as well.
Rotation
Here a 2D point P( x, y ) is repositioned to the point P′( x′, y ′) by moving it along a circular
path:
new
P′ old
θ P
yr
xr
( x′, y′)
r
r ( x, y )
θ
φ
74
For rotation about an arbitrary, but fixed point ( xr , yr ) , i.e. the situation
new
P′( x ′, y ′) old
θ P( x , y )
yr
φ
xr
x′ = xr + ( x − xr )cosθ − ( y − yr )sin θ
(5.9)
y ′ = yr + ( x − xr )sin θ + ( y − yr )cosθ
Thus rotation about a general point is effectively a translation (of the pivot point) plus a rotation
about it.
Note:
• Rotations, in general, are rigid body transformations i.e. they move objects without
deformation.
• To rotate straight lines or polygons, rotate the vertices and re-draw.
• Similar processes may be applied for rotating other objects.
75
The following code shows how to implement a general pivot-point rotation of a polygon. The
polygon vertices, the pivot point and the angle of rotation must be supplied as input.
class wcPt2D {
public:
GLfloat x, y;
};
Scaling
A scaling transformation alters the size of an object. Here the new coordinates ( x′, y ′) are
given by
x ′ = s x x, y ′ = s y y (5.10)
where sx , s y > 0 are scaling factors in the X , Y directions respectively. In matrix form,
x′ sx 0 x sx 0
y ′ = 0 , or P′ = S ⋅ P with S = . (5.11)
sy y 0 s y
Note that
• for values sx , s y < 1 we have a reduction in a particular direction and for values > 1, a
magnification.
• when sx = s y the scaling is said to be uniform, otherwise (for sx ≠ s y ) the scaling is
differential
Examples:
For sx = 2, s y = 1 For sx = s y = 0.5
square → rectangle RH line → LH line
closer to axis and smaller
x′ x
76
Sometimes it is required to scale distances w.r.t. a fixed point, say, ( x f , y f ) . Then the scaled
coordinates ( x′, y ′) are computed from the original ( x, y ) by means of the relations
x′ − x f = s x ( x − x f ), y ′ − y f = s y ( y − y f ) (5.12)
which are usually written as
x′ = sx x + (1 − sx ) x f
. (5.13)
y ′ = s y y + (1 − s y ) y f
Observe that terms (1 − s x ) x f and (1 − s y ) y f are constants for all points in the object and
that the matrix form of (5.13) is :
x′ sx 0 x (1 − sx ) x f (1 − sx ) x f
= + , or P′ = S ⋅ P + . (5.14)
y ′ 0 s y y (1 − s y ) y f (1 − s y ) y f
To scale polygons, we apply (5.13) to each vertex and then redraw the line segments. For
circles with uniform scaling, we apply (5.13) to the centre and scale the radius and re-draw. For
other figures, we need to apply (5.13) to each defining point and then re-draw.
Sample code for scaling
class wcPt2D {
public:
GLfloat x, y;
};
P′ = M 1 ⋅ P + M 2 , (5.15)
where, P′, P, M 2 are 2-element column vectors and M 1 is a 2 × 2 matrix.
For a sequence of operations, we may first perform scaling, followed by a translation and then
a rotation, each time calculating newer coordinates from the old. This sequence can be
77
performed in “one go” by a composite matrix multiplication, without the additive term M 2 in
(5.15), if we employ special coordinates known as homogeneous coordinates.
Homogeneous coordinates
First, represent each ( x, y ) with the homogeneous coordinate triple ( xh , yh , h ) where
x y
x = h , y = h , with h ≠ 0 a real parameter (5.16)
h h
We can now show that all 2D geometric transformations are just matrix multiplications in
homogeneous coordinates.
Thus,
i) for translation we write (and recover (5.2) and (5.3)):
x′ 1 0 t x x
y ′ = 0 1 t y (5.17)
y
1 0 0 1 1
14243
T (t , t )
x y
i.e. we have the form
P′ = T ( t x , t y ) ⋅ P (5.18)
where, P′, P are 3-vectors and T is a 3 × 3 matrix.
Note that:
1 0 − t x
T (t x , t y ) = 0 1 −t y (Prove this!).
−1
0 0 1
ii) for rotations, about the origin, we similarly write (and recover (5.6)-(5.8)):
x′ cosθ − sin θ 0 x
y ′ = sin θ cosθ 0 y (5.19)
1 144
0 42444 0 1 1
3
R(θ )
i.e. we have the form
P′ = R(θ ) ⋅ P (5.20)
iii) for scaling about the origin as fixed point we write and recover (5.11)):
x′ sx 0 0 x
y ′ = 0 s 0 y (5.22)
y
1 0 0 1 1
14243
S (s , s )
x y
i.e. we have the form
P′ = S ( s x , s y ) ⋅ P (5.23)
Further, to rotate about a general pivot or scale w.r.t. a general fixed point, we can use a
succession of transformations about the origin. However, it’s better to use composite
transformation matrices to effect these.
1 0 t x 2 1 0 t x1 1 0 t x1 + t x 2
T ( t x 2 , t y 2 ) T ( t x1 , t y1 ) = 0 1 t y 2 0 1 t y 1 = 0 1 t y1 + t y 2
(5.26)
0 0 1 0 0 1 0 0 1
= T (t x1 + t x 2 , t y1 + t y 2 ).
Thus, 2 successive translations are additive (add their matrix arguments for
the translation parts).
P′ = R(θ 2 ) o [ R(θ1 ) ⋅ P]
(5.27)
=[R(θ 2 ) o R(θ1 )] ⋅ P
which can be verified by multiplying out the matrices giving (verify this)
R(θ 2 ) ⋅ R(θ1 ) = R(θ 2 + θ1 ) (5.28)
or, P′ = R(θ 2 + θ1 ) ⋅ P (5.29)
Thus two successive rotations are additive (add the rotation angles in the
corresponding matrix arguments).
sx 2 0 0 s x1 0 0 s x1 ⋅ s x 2 0 0
0 sy2 0 0 s y1 0 = 0 s y1 ⋅ s y 2 0 (5.30)
0 0 1 0 0 1 0 0 1
or, S ( s x 2 , s y 2 ) ⋅ S ( s x1 , s y1 ) = S ( s x1 ⋅ s x 2 , s y1 ⋅ s y 2 ) . (5.31)
Thus two successive scalings are multiplicative. For example, if we triple the
size of an object twice => final scaling is 9 × original.
Diagramatically we have:
( xr , yr ) ( xr , y r )
1 0 xr cosθ − sin θ 0 1 0 − xr
0 1 yr sin θ cosθ 0 0 1 − y
r
0 0 1 0 0 1 0 0 1
(5.32)
cosθ − sin θ xr (1 − cosθ ) + yr sin θ
= sin θ cosθ yr (1 − cosθ ) − xr sin θ
0 0 1
Diagramatically we have:
( xr , yr ) ( xr , y r )
1 0 xf sx 0 0 1 0 − x f s x 0 x f (1 − sx )
0 1 yf 0 sy 0 0 1 − y = 0 sy y f (1 − s y ) (5.34)
f
0 0 1 0 0 1 0 0 1 0 0 1
or,
T ( x f , y f ) ⋅ S ( sx , s y ) ⋅ T (− x f , − y f ) ≡ S ( x f , y f , s x , s y ) (5.35)
We can use (5.34) to write a function that accepts a general fixed point ( x f , y f ) .
81
S2
O X
θ
S1
we proceed as follows:
• Apply a rotation so that the OS1-OS2 axes coincide the OX-OY axes
• Now scale as before
• Apply reverse rotation to the original directions
(2,2)
(½,1½)
(0,1) (1,1)
(1½,½)
(0,0) (0,0)
(1,0)
Note that if the scaling above was performed w.r.t. an arbitrary fixed point rather
than the origin, then an additional translation matrix would have to be incorporated
into (5.36).
x′ rsxx rs xy trsx x
y ′ = rs rs yy trs y y (5.37)
yx
1 0 0 1 1
T (t x , t y ) ⋅ R( xc , yc ,θ ) ⋅ S ( xc , yc , sx , s y )
sx cosθ − s y sin θ xc (1 − sx cosθ ) + yc s y sin θ + t x
(5.38)
= sx sin θ s y cosθ yc (1 − s y cosθ ) − xc sx sin θ + t y
0 0 1
We note that although (5.37) ⇒ that after the matrix terms are found, 9 mults. + 6 adds are
needed for each point, after concatenation the new coordinates need only be calculated from
x′ = x ⋅ rsxx + y ⋅ rsxy + trsx
(5.39)
y ′ = x ⋅ rs yx + y ⋅ rs yy + trs y
which actually involves only 4 mults. + 4 adds. In fact, without matrix concatenation, individual
matrix operations on every point would increase the overall operations count considerably.
rxx rxy
A= (5.40)
ryx ryy
It is easily seen (by consulting (5.38)) that the following properties hold for A:
83
rxx ryx
iii) r ⋅ r = rxx ryx + rxy ryy = 0 (5.43)
xy yy
i.e. the vectors ( rxx , rxy ) and ( ryx , ryy ) map to the unit vectors (1,0) and (0,1) along the OX
and OY axes respectively.
All these properties follow from the special case of (5.38) corresponding to a rotation about
some pivot ( xr , yr ) through an angle θ followed by a translation written in the form:
cosθ − sin θ 0
R(θ ) = sin θ cosθ 0 (5.44a)
0 0 1
cosθ − sin θ
Here, is its orthogonal sub-matrix which maps the unit vectors
sin θ cosθ
(cosθ , − sin θ ) and (sin θ ,cosθ ) to the unit vectors (1,0) and (0,1) in the coordinate X and
Y directions respectively (prove this!). We can turn this argument around to determine the
rotation matrix R(θ ) from the unit vectors in the old and new coordinate axes positions as
follows:
y y
v
u
initial position
x θ x
final position
84
Suppose that after rotation thro’ θ the new position is determined by the unit vectors
u = (u x , u y ) and v = ( v x , v y ) . We now claim that the first row of the rotation sub-matrix can
be taken as u = (u x , u y ) and the second row is v = ( v x , v y ) . A proof of this can be obtained
by noting that properties (i) – (iii) above must be satisfied and u ⊥ v . Then if have found
v = ( v x , v y ) , say, we must take u = (ux , u y ) ≡ (v y , − v x ) and then the above properties are
satisfied. Moreover, check that the respective unit vectors here map to the correct (1,0) and
(0,1) of the original coordinate system.
Note that rotations involve calculations of sin θ and cosθ terms, which can be very
expensive in some applications. Many algorithms are thus designed so that a full angle θ of
rotation is achieved incrementally in steps of a small ∆θ . Then we can use one or two terms
of a power series for sin ∆θ and cos ∆θ or simply the approximations,
sin ∆θ 0
(OK for ∆θ ≤ 10 ) (5.45)
0
cos ∆θ 1
50 50
Here the triangle is first scaled, then rotated and finally translated. The composite matrix
(compMatrix) is concatenated in that order, starting with compMatrix = identity. The
complete code (from H&B p. 249) is CompositeTransformation2D.cpp and employs OpenGL
only to display the final results whilst the working part is done fully in C++.
1
original
1 0 0
0 -1 0 (5.46)
2 3
X
2′ 3′
reflected
0 0 1
Trans matrix for refl. about y=0
1′
For reflection about the Y-axis the y-values remain unchanged, but the x-values are flipped.
The path of rotation is ⊥ XY plane for all points in the body.
Y
original reflected
-1 0 0
0 0
2 2′
1 (5.47)
1 1′
3 3′ 0 0 1
X Trans matrix for refl. about x=0
For reflection about origin: Both x-values and y-values are flipped. The rotation is about axis
thro’ (0,0) ⊥ XY plane for all points in the body.
Y reflected
3′
-1 0 0
1′
2′
X 0 -1 0 (5.48)
1
2
0 0 1
Note that the above matrix is same as the rotation matrix R(θ ) with θ = 1800 , so both these
operations are equivalent.
For reflection about any other point: Same as rotation about an axis through the fixed reflection
point Pref =(xref, yref) and ⊥ XY plane for all points in the body.
86
Y reflected
3′
What is the transformation matrix
2′ for this case?
1′
Pref
1
2
original
3
X
Y original
0
3 y=x 1 0
2
1
3′ 1 0 0 (5.49)
1′
2′
0 0 1
reflected
Trans matrix for refl. about y=x
X
Prove that (5.49) is the transformation matrix by concatenating the matrices for:
• clockwise rotation thro’ 450 about (0,0), rotating line y = x to X-axis
• reflection about the X-axis
• rotating X-axis back to the line y = x
Note that this process is also equivalent to reflection about the X-axis + rotation thro’ +900.
3′
0 -1 0
2
1 1′ -1 0 0 (5.50)
X
original 3
y=-x
0 0 1
Trans matrix for refl. about y= - x
Variations
• Reflections about the coordinate axes or (0,0) can also be implemented as scaling with
negative scaling factors
• Can set elements of reflection matrix to be <, >, ±1.
o For > |±1| mirror image is shifted further away
o For < |±1| mirror image is shifted nearer axis
Shear
A shear transformation distorts the shape of an object, causing “internal layers to slide over”.
Two types here are i) a shift in x-values and ii) a shift in y-values.
x-shear:
1 shx 0
x′ = x + shx ⋅ y 0 1 0
Given by with matrix (5.51)
y′ = y
0 0 1
Here shx is any real number. For example shx = 2 changes the square below into a parm:
(0,0) (0,0)
(1,0) (1,0)
(0,0) (0,0)
(1,0) (½,0) (1½,0)
y = -1
88
1 0 0
x′ = x sh
Given by with matrix 1 − shy ⋅ xref (5.53)
y ′ = y + sh y ( x − xref ) y
0 0 1
This will shift coordinate positions vertically by an amount ∝ distance from the reference line x
= xref, e.g. with sh y = 12 relative to the line x = xref = −1 we have:-
(1,2)
(0,1½)
(1,1)
(0,1) (1,1)
(0,½)
(0,0)
(-1,0) (0,0) (1,0)
x = -1
Remark: Shears may also be expressed as compositions of the basic transformations. For
example, (5.51) may be written as a rotation + a scaling. Can you confirm this?.
y0
θ
x0 X
Y
Y′
y P
X′
y′
x′
θ X
x
Then we do a clockwise rotation via
cos θ sin θ 0
R( −θ ) = − sin θ cos θ 0 (5.55)
0 0 1
This result can also be obtained directly by deriving the relations between the respective
coordinate-distance pairs in the two systems from the figure above.
x0 X
Thus, suppose V is a point vector in the XY system and is in the same direction as the +ve Y′
coordinate axis. Then if a unit vector along the +ve Y′ coordinate axis is, say,
V
v= ≡ (vx , v y ) ) (5.57)
V
we take the unit vector along the +ve X′ coordinate axis as
u = (ux , u y ) ≡ ( v y , − v x ) (5.58)
so that the rotation matrix is (see section following equation (5.44) for the reason):
u x uy 0
R = vx vy 0 (5.59)
0 0 1
0 1 0
R = −1 0 0 . (5.60)
0 0 1
x0 X
P1 − P0 V
Then we can use the unit vector = = v ≡ (vx , v y ) (5.61)
P1 − P0 V
with u = (ux , u y ) ≡ ( v y , − v x ) (5.62)
entire
block
moved
Pmin
to new
position
P0
before after
• take the pixel array ~ the block and reverse each row
• then interchange the rows and columns
e.g.
1 2 3 3 2 1
4 5 6 6 5 4 3 6 9 12
→ → 2 5 8 11 ~ 900 rotation
7 8 9 9 8 7
1 4 7 10
10 11 12 12 11 10
For scalings and reflections similar methods (or combinations of the above) can be devised.
See H&B p.257.
ii) glReadPixels(xmin,ymin,width,height,GL_RGB,GL_UNSIGNED_BYTE,
colorArray);
Allows for the saving of a pixel block in array colorArray. If color-table indices are
stored at pixel positions, then replace GL_RGB with GL_COLOR_INDEX. To
rotate the pixel block we re-arrange the rows and columns of colorArray.
iii) glDrawPixels(width,height,GL_RGB,GL_UNSIGNED_BYTE,colorArray);
Used to put back the rotated array into the buffer, with lower-left corner at current
selected raster position.
The source buffer must have been selected with glReadBuffer(..) and the
destination buffer with glDrawBuffer(..).
iv) glPixelZoom(sx,sy);
Used to do 2D scaling with scaling factors sx and sy. These are floating-point
values: 0...< 1.0 => decrease size, > 1.0 => increase, < 0.0 => reflection of source
w.r.t. current raster position.
First call above and then follow with glCopyPixels(..) or glDrawPixels(..).
Raster ops may also be combined with logical operations (and, or, not, xor) in OpenGL.
5.11.1 Translation
92
where t x , t y , t z = translation displacements along each of X,Y,Z axes respectively. Note that
(5.63) expands to the usual component forms:
x′ = x + t x , y′ = y + t y , z′ = z + tz . (5.64)
Y
T = (t x ,t y , t z )
P′ = ( x′, y ′, z′)
P = ( x, y, z )
X
To translate a 3D object, we translate each of its defining points and then re-draw it. For
example, for polygons we translate each vertex and re-draw.
As before the inverse translation is given by the operator T = ( −t x , − t y , −t z ) with its
corresponding matrix similarly amended.
5.11.2 Rotation
Here we need to specify an axis of rotation which can be any 3D line in space and the angle of
rotation. By convention we take the counter-clockwise direction as positive when looking from
far on the X (or Y or Z or line) axis → origin O:
Y
θ = + ve observer
O X
Before considering an arbitrary 3D axis of rotation we first consider rotations about the
coordinate axes.
Y Z X
θ X θ Y θ Z
Z X Y
For rotation thro’ θ about the Z-axis recall from the 2D case, the equations:
x′ = x cosθ − y sin θ
y ′ = x sin θ + y cosθ (5.65)
z′ = z
This easily is extended to 3D homogeneous coordinates by writing
x′ cosθ − sin θ 0 0 x
y ′ sin θ cosθ 0 0 y
= or P′ = Rz (θ ) ⋅ P (5.66)
z′ 0 0 1 0 z
1 1444
0 0
424444
0 1 1
3
R (θ )
z
For the transformation equations for rotations about the X-axis and Y-axes respectively, we
can similarly obtain them by cyclic permutation of the coordinates ( x → y → z → x ) in (5.66)
and consultation of the figures above as follows:-
y ′ = y cosθ − z sin θ
z′ = y sin θ + z cosθ (5.67)
x′ = x
i.e.
x′ 1 0 0 0 x
y ′ 0 cosθ − sin θ 0 y
= . or P′ = Rx (θ ) ⋅ P (5.68)
z′ 0 sin θ cosθ 0 z
1 1444
0 0
424444
0 1 1
3
R (θ )
x
Remark:
For inverse rotations, as in 2D, we replace θ by −θ in the above. These then correspond to
clockwise rotations. Since the only changes are sin θ → − sin θ and vice versa, we find that
the inverse matrix is the same as the original but with the rows and columns interchanged, i.e.
in all the above,
R −1 = R T . (orthogonality property) (5.71)
General 3D rotations
For rotation about any other axis, we can set up the transformation matrix as a composition of
rotations about the coordinate axes.
i) As a special case for rotation about an axis || a coordinate axis we follow the
steps:
1. Translate object so rotation axis coincides the || coordinate axis
2. Do rotation about that coordinate axis
3. Translate so that rotation axis → original rotation axis position
Example:
Suppose the axis of rotation || X-axis, then we have:
95
Y Y Y Y
rotate translate
rotation axis
X X
X X
Z Z Z Z
ORIGINAL translate FINAL
ii) In the general case for rotation about an axis not || a coordinate axis we
follow the steps:
1. Translate object so rotation axis passes through (0,0,0)
2. Rotate object so rotation axis is || (i.e. coincides) one of the coordinate axes
3. Perform rotation about this coordinate axis
4. Apply inverse rotation to bring axis to original orientation but thro’ (0,0,0)
5. Apply inverse translation on object so axis returns to original position
This sequence of operations may be applied with any one of the coordinate axes. For
convenience let’s choose the Z-axis.
First, we can specify a rotation axis by 2 coordinate position vectors (P1 , P2 ) or one
position vector and direction angles (or direction cosines). Consider the case of 2 points
and a counter-clockwise angle, when looking from P2 to P1 :
Y Y Y
P2 P2
θ u u
P1 P1
X X X
Z Z Z
is –ve (clockwise) we reverse the V and u vectors (i.e. they point in the P2 to P1
direction); otherwise take θ as a negative angle.
Now:
1. We move P1 to (0,0,0) via
1 0 0 − x1
0 1 0 − y
T = 1
(5.76)
0 0 0 − z1
0 0 0 1
2. Next put u onto the Z-axis. We use the rotations
• u → XZ plane giving u′′ by rotating around the X-axis thro’ some angle α
• u′′ → Z -axis giving u z = (0,0,1) by rotating around the Y-axis thro’ some
angle β
Y Y Y Y
u u′ u
β β
α α
X X X X
u′′ u′′ = (a ,0, d )
u′′ u z = (0,0,1) u z = (0,0,1)
u z = (0,0,1)
Z Z Z Z
To find the angles α and β and hence the rotation matrices we proceed as follows:
Let u′ = (0, b, c ) be the projection of u onto the YZ plane. Then α is the angle between u′
and the Z-axis.
u′⋅ u z c
Now, cos α = = , where d = b 2 + c 2 (5.77)
u′ u z d
Also, u′ × u z = ux u′ u z sin α , where u x = (1,0, 0) (5.78)
1 0 0 0
c −b
0 0
Rx (α ) =
d d
(5.81)
b c
0 0
d d
0 0 0 1
Suppose that after this rotation the result is u′′ . Then its
97
d 0 −a 0
0 1 0 0
Ry (β ) = (5.87)
a 0 d 0
0 0 0 1
3. Now after applying (5.76), (5.81) and (5.87) we have aligned the axis of rotation along
the +ve Z-axis.
4. Next we rotate through the specified angle θ about the Z-axis via the transformation
matrix
cos θ − sin θ 0 0
sin θ cos θ 0 0
Rz (θ ) = (5.88)
0 0 1 0
0 0 0 1
5. To complete the transformation, we must transform the final rotation axis back to the
original position by applying inverses to the above transformations, giving the final
transformation matrix as
and with the following algebraic operations defined on the set of such numbers:
i) scalar multiplication –
uq = us + i (ua ) + j (ub) + k (uc ) for any real (scalar) u
ii) quaternion addition (“+”) is associative and is given by, for quats q1 , q2 –
q1 + q2 = ( s1 + s2 ) + i (a1 + a2 ) + j (b1 + b2 ) + k ( c1 + c2 )
iii) quaternion multiplication (“.”) is associative and is given by –
q1 ⋅ q2 = ( s1 + ia1 + jb1 + kc1 ) ⋅ ( s2 + ia2 + jb2 + kc2 )
= expansion with (5.91)
The operation (5.100) can be shown to be equivalent to the rotation part operation of (5.89),
with the matirx
Rx−1 (α ) ⋅ Ry−1 ( β ) ⋅ Rz (θ ) ⋅ R y ( β ) ⋅ Rx (α ) . (5.103)
As an outline of the process to be followed, write the vector part of q as v = ( a, b, c ) and
expand (5.102) to a 3 × 3 matrix equation involving p = ( x, y , z ) and p′ = ( x′, y ′, z′) , namely
p′ = M R (θ ) ⋅ p
where,
1 − 2b2 − 2c 2 2ab − 2 sc 2ac + 2 sb
M R (θ ) = 2ab + 2 sc 1 − 2a − 2c
2 2
2bc − 2 sa (5.104)
2ac − 2 sb 2bc + 2 sa 1 − 2a − 2b
2 2
99
With values for a,b,c from (5.98), u = (u x , u y , u z ) and some trigonometric identities we find
M R (θ ) =
u x2 (1 − cosθ ) + cosθ u x u y (1 − cosθ ) − uz sin θ u x uz (1 − cosθ ) + u y sin θ
(5.105)
u y u x (1 − cosθ ) + uz sin θ u (1 − cosθ ) + cosθ u y uz (1 − cosθ ) − u x sin θ
2
y
u z u x (1 − cosθ ) − u y sin θ u z u y (1 − cosθ ) + u x sin θ uz2 (1 − cosθ ) + cosθ
Finally, when the form (5.103) is fully expanded it can be shown to be just the matrix (5.105).
For more detail consult H&B p. 273.
The complete rotation about an arbitrary axis must then include translation as before to give the
form
R(θ ) = T −1 ⋅ M R (θ ) ⋅ T (5.106)
which corresponds to (5.89).
The following code shows how to construct a 3D rotation matrix (H&B p. 275).
class wcPt3D {
public:
GLfloat x, y, z;
};
typedef float Matrix4x4 [4][4];
Matrix4x4 matRot;
/* Construct the 4 by 4 identity matrix. */
void matrix4x4SetIdentity (Matrix4x4 matIdent4x4)
{
GLint row, col;
for (row = 0; row < 4; row++)
for (col = 0; col < 4 ; col++)
matIdent4x4 [row][col] = (row == col);
}
/* Premultiply matrix m1 times matrix m2, store result in m2. */
void matrix4x4PreMultiply (Matrix4x4 m1, Matrix4x4 m2)
{
GLint row, col;
Matrix4x4 matTemp;
for (row = 0; row < 4; row++)
for (col = 0; col < 4 ; col++)
matTemp [row][col] = m1 [row][0] * m2 [0][col] + m1 [row][1] *
m2 [1][col] + m1 [row][2] * m2 [2][col] +
m1 [row][3] * m2 [3][col];
for (row = 0; row < 4; row++)
for (col = 0; col < 4; col++)
m2 [row][col] = matTemp [row][col];
}
5.11.4 3D scaling
To scale p = ( x, y , z ) relative to the coordinate origin we use,
x′ sx 0 0 0 x x
y ′ 0 s 0 0 y y
P′ ≡ = y or, P′ = S ⋅ P ≡ S ⋅ (5.107)
z′ 0 0 sz 0 z z
1 1442443
0 0 0 1 1 1
S
which expands to x ′ = s x ⋅ x , y ′ = s y ⋅ y , z ′ = sz ⋅ z . (5.108)
Note that scaling changes the size of an object as well its position relative to the origin.
As before, to preserve the original shape we must have sx = s y = sz .
101
T ( x f , y f , z f ) ⋅ S ( s x , s y , sz ) ⋅ T ( − x f , − y f , − z f )
sx 0 0 (1 − sx ) x f
0 sy 0 (1 − s y ) y f (5.109)
= .
0 0 sz (1 − sz ) z f
0 0 0 1
The following code shows how to construct a 3D scaling matrix (H&B p. 277).
class wcPt3D
{
private:
GLfloat x, y, z;
public:
/* Default Constructor:
* Initialize position as (0.0, 0.0, 0.0).
*/
wcPt3D ( ) {
x = y = z = 0.0;
}
void scale3D (GLfloat sx, GLfloat sy, GLfloat sz, wcPt3D fixedPt)
{
102
Matrix4x4 matScale3D;
/* Initialize scaling matrix to identity. */
matrix4x4SetIdentity (matScale3D);
matScale3D [0][0] = sx;
matScale3D [0][3] = (1 - sx) * fixedPt.getx ( );
matScale3D [1][1] = sy;
matScale3D [1][3] = (1 - sy) * fixedPt.gety ( );
matScale3D [2][2] = sz;
matScale3D [2][3] = (1 - sz) * fixedPt.getz ( );
}
class wcPt3D {
public:
GLfloat x, y, z;
};
typedef GLfloat Matrix4x4 [4][4];
Matrix4x4 matComposite;
/* Construct the 4 by 4 identity matrix. */
void matrix4x4SetIdentity (Matrix4x4 matIdent4x4)
{
GLint row, col;
for (row = 0; row < 4; row++)
for (col = 0; col < 4 ; col++)
matIdent4x4 [row][col] = (row == col);
}
Example: The following reflection about the XY plane results in the RH → LH system (or vice
versa).
y y 1 0 0
0
z
0 1 0 0
x x M zrefl = (5.110)
0 0 −1 0
z
0 0 0 1
Shears
A 3D shear is used to modify shapes, for example a Z-axis shear relative to a reference
position zref is effected by
The effect of this transformation is to change the x and y values by an amount proportional to
the distance from the plane z = zref .
Recall that in the 2D case when we wanted to transform a scene from one Cartesian system to
another, we obtained the required transformation matrix by concatenating the matrices required
to back-transfer the new coordinate to the old position.
y y′
u′y x′
( x0 , y0 , z 0 ) u′x
(0,0,0) x
u′z
z
z′
Thus, if XYZ is one system with origin (0,0,0) and X ′Y ′Z ′ is the other where the latter has
origin ( x0 , y0 , z0 ) relative to (0,0,0) and in addition unit vectors u′x , u′y , u′z along its coordinate
axes, then the rotation matrix
u′x1 u′x 2 u′x 3 0
u ′ u ′ u ′ 0
R = y1 y2 y3 (5.112)
′ ′ ′
u z1 u z 2 u z 3 0
0 0 0 1
will transform the vectors u′x , u′y , u′z onto the X, Y and Z axes directions respectively. The final
transformation is obtained by translating ( x0 , y0 , z0 ) to (0,0,0) so that the composite matrix is
The special affine transformation, involving only translation, rotation and reflection also
• preserves angles
• preserves lengths of lines
in addition to the above properties.
106
c. glMultMatrix*(otherElements16);
//post-multiplies current matrix by matrix formed into column-major form from
//16-element array otherElements16.
Example: The code segment
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glMultMatrix(elemsM2);
glMultmatrix(elemsM1);
glGetIntegerv(GL_MODELVIEW_STACK_DEPTH, numMats);
108
//returns single integer value to array numMats, giving the current number of
matrices loaded onto the modelview stack.
b. glPushMatrix( ); //copies the current matrix to the top of the active stack. The
//second position now contains a copy of the current.
c. glPopMatrix( ); //will destroy matrix at the top of the active stack, and makes
//the second matrix the current matrix.
200
150
translated
100
original
50
200
150
100
rotated original
50
200
scale-reflected
150
100
original
50
//--------------------------------
//Version 1 - simple case
//--------------------------------
glMatrixMode (GL_MODELVIEW);
glColor3f (0.0, 0.0, 1.0);
glRecti (50, 100, 200, 150); // Display blue rectangle.
glColor3f (1.0, 0.0, 0.0);
glTranslatef (-200.0, -50.0, 0.0); // Set translation parameters.
glRecti (50, 100, 200, 150); // Display red, translated rectangle.
glLoadIdentity ( ); // Reset current matrix to identity.
glRotatef (90.0, 0.0, 0.0, 1.0); // Set 90-deg. rotation about z axis.
glRecti (50, 100, 200, 150); // Display red, rotated rectangle.
glLoadIdentity ( ); // Reset current matrix to identity.
glScalef (-0.5, 1.0, 1.0); // Set scale-reflection parameters.
glRecti (50, 100, 200, 150); // Display red, transformed rectangle.
//-------------------------------------------
//Version 2 – using matrix stacks
//------------------------------------------
glMatrixMode (GL_MODELVIEW);
glColor3f (0.0, 0.0, 1.0); // Set current color to blue.
glRecti (50, 100, 200, 150); // Display blue rectangle.
glPushMatrix ( ); // Make copy of identity (top) matrix.
glColor3f (1.0, 0.0, 0.0); // Set current color to red.
glTranslatef (-200.0, -50.0, 0.0); // Set translation parameters.
glRecti (50, 100, 200, 150); // Display red, translated rectangle.
glPopMatrix ( ); // Throw away the translation matrix.
glPushMatrix ( ); // Make copy of identity (top) matrix.
glRotatef (90.0, 0.0, 0.0, 1.0); // Set 90-deg. rotation about z axis.
glRecti (50, 100, 200, 150); // Display red, rotated rectangle.
glPopMatrix ( ); // Throw away the rotation matrix.
glScalef (-0.5, 1.0, 1.0); // Set scale-reflection parameters.
glRecti (50, 100, 200, 150); // Display red, transformed rectangle.
//---------------------------------------------------------
//Version 3 using composite transformations
//---------------------------------------------------------
class wcPt3D {
public:
GLfloat x, y, z;
};
\* Procedure for generating a matrix for rotation about an
* an axis defined with points p1 and p2.
*/
void rotate3D (wcPt3D p1, wcPt3D p2, GLfloat thetaDegrees)
{
/* Set up components for rotation-axis vector. */
float vx = (p2.x - p1.x);
float vy = (p2.y - p1.y);
float vz = (p2.z - p1.z);
/* Specfy translate-rotate-translate sequence in reverse order: */
glTranslatef (p1.x, p1.y, p1.z); // Move p1 back to original position.
/* Rotate about axis through origin: */
glRotatef (thetaDegrees, vx, vy, vz);
glTranslatef (-p1.x, -p1.y, -p1.z); // Translate p1 to origin.
}
/* Procedure for generating a matrix for a scaling
* transformation with respect to an arbitrary fixed point.
110
*/
void scale3D (GLfloat sx, GLfloat sy, GLfloat sz, wcPt3D fixedPt)
{
/* Specfy translate-scale-translate sequence in reverse order: */
/* (3) Translate fixed point back to original position: */
glTranslatef (fixedPt.x, fixedPt.y, fixedPt.z);
glScalef (sx, sy, sz); // (2) Scale with respect to origin.
/* (1) Translate fixed point to coordinate origin: */
glTranslatef (-fixedPt.x, -fixedPt.y, -fixedPt.z);
}
void displayFcn (void)
{
/* Input object description. */
/* Set up 3D viewing-transformation routines. */
/* Display object. */
glMatrixMode (GL_MODELVIEW);
/* Input translation parameters tx, ty, tz. */
/* Input the defining points, p1 and p2, for the rotation axis. */
/* Input rotation angle in degrees. */
/* Input scaling parameters: sx, sy, sz, and fixedPt. */
/* Invoke geometric transformations in reverse order: */
glTranslatef (tx, ty, tz); // Final transformation: Translate.
scale3D (sx, sy, sz, fixedPt); // Second transformation: Scale.
rotate3D (p1, p2, thetaDegrees); // First transformation: Rotate.
/* Call routines for displaying transformed objects. */
}
111
Here we study the processes involved in depicting a 2D scene or model onto a 2D screen. We
thus consider the viewing pipeline from world coordinates to device coordinates.
viewport
clipping window yvmax
ywmax
yvmin
ywmin
Recall from Chapter One, that the steps involved in the complete wc → dc 2D viewing
transformation can be stated as:
i) Construct the scene in world coordinates, using modelling coordinates for each
part.
ii) Set up a 2D viewing system with an oriented window
iii) Transform to viewing coordinates
iv) Define a viewport in normalized (0..1 or -1...1) coordinates and map to it from the
view coordinates
v) Clip all parts outside the viewport
vi) Transform to device coordinates
When all the transformations are done, clipping can be done in normalized coordinates or
device coordinates. The clipping process is fundamental in computer graphics.
yworld yview
clipping window
v
y0
R
u
xview
T
x0 xworld
World Coordinates
V
v= = (v x , v y ) = 2nd row of R
V
u = (ux , u y ) = (v y , − v x ) = 1st row of R
After this transformation, the following example illustrates the final outcome in WC:
yworld v yworld
y0 u
clipping window
x0 xworld xworld
and then the viewport transformation is applied. In the latter the viewport is specified in screen
coordinates relative to the display window position.
Here we consider a viewport defined with normalized coordinates in the range 0...1:
1
viewport
clipping window yvmax
ywmax ( xv, yv )
( xw, yw) yv
yw yvmin
ywmin
xv − xvmin xw − xwmin
=
xvmax − xvmin xwmax − xwmin
(6.2)
yv − yvmin yw − ywmin
=
yvmax − yvmin ywmax − ywmin
Re-arranging these relations give the equations in terms of scaling factors and translations
terms
xv = s x ⋅ xw + t x
(6.3)
yv = s y ⋅ yw + t y
xvmax − xvmin
sx =
xwmax − xwmin
where (6.4)
yvmax − yvmin
sy =
ywmax − ywmin
Note that when sx = s y the relative proportions (x vs. y) of individual objects are maintained,
otherwise stretching or contraction would occur in a particular direction.
114
Thus (consult Chap. 5) for step 1 and we use the scaling matrix
sx 0 xwmin (1 − sx )
S =0 sy ywmin (1 − s y ) (6.6)
0 0 1
and in step 2 we employ,
1 0 xvmin − xwmin
T = 0 1 yvmin − ywmin (6.7)
0 0 1
resulting in the composite wc → normalized viewport transformation matrix,
sx 0 tx
Μ wc.window →norm.viewport =T ⋅S = 0 sy ty . (6.8)
0 0 1
screen viewport
clipping window normalized square yvmax
ywmax 1 ( xv, yv )
( xw, yw) yv
yw -1 1
yvmin
-1
ywmin
( xnorm , ynorm )
xwmin xw xwmax xvmin xv xvmax
world coordinates normalized coordinates screen coordinates
For this process we simply have to adjust the transformation equations of section 6.3 as
follows:
115
For wc → nc we put xvmin = yvmin = −1 and xvmax = yvmax = +1 into (6.6) – (6.8) to obtain
2 xwmax + xwmin
xw − xw 0
xwmax − xwmin
max min
2 ywmax + ywmin
Μ wc.window →norm.square = 0 . (6.9)
ywmax − ywmin ywmax − ywmin
0 0 1
For nc → vc, after applying clipping, we transform the normalized square of length 2 to the
screen viewport by we putting xwmin = ywmin = −1 and xwmax = ywmax = +1 into (6.6) – (6.8)
to obtain
In choosing viewport sizes, it is important to maintain its aspect ratio as the same as the
original clipping window’s to avoid coordinate stretching and hence distortion in shapes. In the
final step we position the viewport in the display window. The usual technique is to set its
lower-left corner, at some position ( xs , ys ) relative to the lower-left position of the display
window, e.g.
100
s
c
window title here
r
e 50
e
n 300
ys
xs viewport
display window
400
Remarks:
• Character strings can be mapped by sending them through as bitmaps (for maintaining
constant character sizes), or by transforming the defining positions of outline fonts,
similar to any other primitive.
116
To make a display window the current one (on which all subsequent ops performed)
glutSetWindow(windowID); //default is last created
To redisplay a window (if it may have been damaged after a glutPopWindow(.) for example)
call
glutPostRedisplay(); //redisplays a current window
Example program:
This program will construct and display two triangles, each in its own viewport, with the second
one rotated through +900.
119
//OGLViewingProg2D.cpp
//--------------------------------
//Split-screen (dual viewport) example (H & B Chap 6)
#include <windows.h>
#include <GL/glut.h>
class wcPt2D {
public:
GLfloat x, y;
};
6.6 2D clipping
Clipping algorithms are used to remove unwanted portions of a picture – we either remove the
outside or the inside of a clipping window, which can be rectangular (most common),
polygonal or curved (complex). A study of the various algorithms for point clipping, line
clipping, polygon clipping, curved area clipping and text clipping is being deferred for later
study, if time permits, and in order to study other aspects of CG. If you are interested in this
subject you can obtain a set of (hand written) notes from me, or see H&B p. 315.
121
Viewing a scene in 3D is much more complicated than 2D viewing, where in the latter the
viewing plane on which a scene is projected from WCs is basically the screen, except for its
dimensions. In 3D, we can choose different viewing planes, directions to view from and
positions to view from. We also have a choice in how we project from the WC scene onto the
viewing plane.
In addition, depth cueing which involves colouring different surfaces or parts of them to give
the impression of different position depths, may be employed, and surface rendering,
lighting, visible surface detection and clipping algorithms come into play as well.
The view plane or projection plane is usually taken as a plane that is ⊥ zv -axis and is set at
a position zvp from the origin. Its orientation is specified by a choosing view-plane normal
vector N (w.r.t. WC system) which also specifies the direction of the positive zv direction. The
following right-handed systems are indicative of the set up typically employed.
122
yv xv
yw zv
view plane
N P0 = ( x0 , y0 , z0 )
Pref
VC system
xw
WC system
zw
The direction of viewing is usually taken as the − N (or − zv ) direction, for RH coordinate
systems (or in the opposite direction corresponding to LH coordinate systems).
Choosing the view-plane normal N :
• can take as “out from object” by taking N = P0 − OriginWC
• or, from a reference point Pref (“look at point”) in scene to P0 i.e. N = P0 − Pref
• or, define direction cosines for it using angles θ ,ϕ ,φ w.r.t. the WC X,Y,Z axes
N
n= = ( nx , n y , nz )
N
V×n
u= = (u x , u y , u z ) (7.1)
V
v = n× u = ( v x , v y , v z )
We then call this system a uvn viewing-coordinate reference frame.
scene
yv − zv
another possible view plane
v P0 = ( x0 , y0 , z0 )
xv
n u
uvn viewing reference frame
a view plane
+ zv
viewing direction...towards object
xv
yv VC
yw yw yw
P0
zv y ′ xv′ yvfinal
v
xw xvfinal
WC xw xw
zvfinal
zw zv′ zw
zw
T Rx Ry Rz
R = Rz ⋅ Ry ⋅ Rx
• Since we know the coordinate unit vectors for the VC system from (7.1), we can obtain
the final rotation matrix R by placing these vector components into its rows, by
extending the arguments of Chap. 5. Thus we can write,
ux uy uz 0
v vy vz 0
R= x (7.3)
nx ny nz 0
0 0 0 1
u x uy uz −u ⋅ P0
v vy vz − v ⋅ P0
M WC→ VC = R ⋅ T = x (7.4)
nx ny nz −n ⋅ P0
0 0 0 1
Here coordinate positions are transformed along parallel lines, which meet the view
plane, thus forming the projected image. Used in drafting and accurate engineering
drawings.
• perspective projection
Here coordinate positions are transformed along lines that converge to a point behind
the view plane, called the projection reference point or centre of projection. The
image on the view plane is formed by the intersections of the lines with it. Thus distant
objects look smaller than nearer ones of the same sizes, adding realism in the view.
P2
P1 view plane
P2
P2′
P2′ convergence point
P1
P1′
view plane P1′
where now the z coordinate may be used for depth cueing, i.e. can colour each pixel in
proportion to z.
This transformation is also known by the term orthographic parallel projection.
An orthogonal-projection view volume is used to limit the volume extent of a scene that will
be orthogonally projected onto the view plane. This is formed by defining a rectangular clipping
window in the view plane and extending this to a rectangular parallelpiped in the zv by
enclosing it between the planes zv = znear and zv = z far which are parallel to the view plane.
126
yv
view plane
far plane
( x, y )
view volume
( x, y , z ) near plane
xv
Mapping to a normalized view volume is then done in most packages, where the normalized
view volume has coordinate extents the unit cube (coordinate ranges 0...1) or a symmetric
cube (coordinate ranges -1...+1). After this, the normalized volume is mapped to screen
coordinates, which generally are LH coordinate systems.
xnorm
( −1, −1, −1)
xv xscreen
( xwmin , ywmin , znear ) screen window viewport
zv
Here, for the coordinate systems shown position ( xwmin , ywmin , znear ) → (−1, −1, −1) and
position ( xwmax , ywmax , z far ) → (1,1,1) . The procedure for obtaining the transformation
matrix from ( xv , yv , zv ) → ( xnorm , ynorm , znorm ) follows the argument for (6.9), with the
additional requirement that we transform the range ( znear , z far ) → (−1,1). The matrix that
results is
2 xwmax + xwmin
xw − xw 0 0
xwmax − xwmin
max min
2 ywmax + ywmin
0 0
M ortho→norm = ywmax − ywmin ywmax − ywmin . (7.7)
−2 znear + z far
0 0
znear − z far znear − z far
0 0 0 1
127
object
Vp
yv
view plane
( xp , y p ) xv
( x, y , z )
α
L φ
z
( x, y )
zv
Suppose that the point ( x, y , z ) → ( x p , y p ) on the view plane (at z = 0 ) and that ( x, y ) is its
orthographic projection on the plane.
Let α = angle between the line joining ( x, y ) and ( x p , y p ) on the plane and the line
( x, y , z ) → ( x p , y p ) .
Let φ = angle between the line joining ( x, y ) and ( x p , y p ) and the xv - axis. Then
x p = x + L cos φ
(7.8)
y p = y + L sin φ
128
tan α = z / L . (7.9)
Thus L = z / tan α ≡ z L1 , where L1 = 1/ tan α = cot α (7.10)
and
x p = x + L1 z cos φ
(7.11)
y p = y + L1 z sin φ
If the view plane was located at some position zvp on the zv - axis, then the above
transformation equations from VC to projected coordinates are simply affected by a translation,
so with zvp − z replacing z
x p = x + L1 ( zvp − z ) cos φ
y p = y + L1 ( zvp − z )sin φ (7.12)
z p = zvp
Remarks:
1. From (7.13) we recover the orthogonal projection (orthographic) case by setting
L1 = 0 or α = 900 .
2. The transformation (7.13) corresponds to the matrices for z-shears: planes of constant
z are sheared and projected onto the view plane. See H&B p. 365.
3. Typical values used in applications for the angle φ are φ = 300 , 450 giving a
combination of front, side, and top views (or front, side and bottom views).
4. Common values for α are α = 450 or tan α = 1 (cavalier projection) α ≈ 63.40 or
tan α = 2 (cabinet view). See H&B p.365.
Some graphics packages supply a WC→ VC transformation matrix in terms of the parallel-
projection view vector Vp . Then all WC points are projected onto the view plane with their
projection lines parallel to this vector. The matrix (7.13) can be written in terms of components
of Vp w.r.t. the view coordinate system. See H&B p.366. Finally, as in the case of orthogonal
projections we need to perform the VC → normalized → viewport transformations to obtain the
complete (composite) transformation from WC → screen viewport.
129
view plane
( x, y , z ) yv view plane
xv
( x, y , z )
( x p , y p , zvp ) ( x p , y p , zvp )
( x prp , y prp , z prp )
zv
zv zvp vc origin z prp
vc origin
Consider the general point P = ( x, y , z ) being projected onto the view plane as ( x p , y p , zvp ) ,
along the line converging to the projection reference point ( x prp , y prp , z prp ) . With the view
plane ⊥ the zv axis of the VC system and located at zvp on it, the RH figure shows the situation
when viewed from the xv axis.
Now, along the perspective projection line, the equations determining any general point
P′ = ( x′, y ′, z ′) on it may be written in parametric form as
x′ = x − ( x − x prp )u
y ′ = y − ( y − y prp )u ; 0 ≤ u ≤ 1 (7.14)
z′ = z − ( z − z prp )u
where for the parameter value u = 0 we recover the starting point P = ( x, y , z ) and for u = 1
the projection reference point ( x prp , y prp , z prp ).
Now, suppose we consider a P′ = ( x′, y ′, z ′) = ( x p , y p , zvp ) to be a point on the view plane,
i.e. z′ = zvp , then from the last of (7.14) we obtain at this position the parameter value
zvp − z
u= (7.15)
z prp − z
Substituting this into (7.14) gives the equations for the WC → VC perspective transformation
z −z z −z
x p = x prp vp + x prp vp
z −z
prp z prp − z
z prp − zvp zvp − z
yp = y + y prp (7.16)
z − z z − z
prp prp
z −z z − z
z p = zvp = z prp vp + z prp vp
z −z z − z
prp prp
Notice that if we write this equation in matrix form, then the matrix elements will contain
functions of the original z-coordinate. This does not make the transformation matrix useful for
130
the purposes of concatenation with other matrices. We resolve this problem by invoking
homogeneous coordinates by setting:
Homogeneous parameter h = z prp − z (7.17)
xh y z
Projection-homogenous coordinate relations xp = , yp = h , zp = h (7.18)
h h h
Then (7.16) becomes
Special cases:
1. When the view plane = xv − yv (i.e. u - v ) plane, then put zvp = 0 in the above.
2. When the perspective projection reference point is on the zv axis, then
x prp = y prp = 0 .
3. When the perspective projection reference point is the VC origin then
x prp = y prp = z prp = 0 .
4. A very common choice is to take the view plane as the xv − yv (i.e. u - v ) plane with
the projection reference point on the zv axis, so that x prp = y prp = zvp = 0 , giving:
z z
x p = prp x, y p = prp y (7.22)
z −z z −z
prp prp
5. If the projection reference point is on the view plane, z prp = 0 , the whole scene maps
to a single point (show this), so in general we take the view plane to lie between the
scene and the projection reference point.
6. It is also possible to obtain the transformation WC → VC matrix (7.4) in terms of other
parameters, such as spherical coordinates ( r, φ ,θ ) for positioning the viewing frame:
(See tutorial)
131
zw yv xv
viewing
direction (-N)
z0 zv
VC
frame
r P0 = ( x0 , y 0 , z0 )
WC φ
frame y0
yw
θ r sin φ
x0
xw
clipping window
projection
reference point
frustum view
volume
The frustum view volume is termed symmetric when the projection line from the projection
reference point, through the centre of the clipping window is ⊥ view plane, otherwise it is called
132
an oblique frustum. See H&B pp. 374-380 for more detail, where in particular they obtain the
transformation matrix to such a volume by concatenating M persp above with a shearing matrix
M zshear i.e. M oblique.persp = M persp ⋅ M zshear . Thereafter, the frustum volume is mapped to a
normalized view volume and a viewport in screen coordinates giving a concatenated matrix
M norm.persp for the WC
The complete transformation matrix from WC → normalized perpective-projection
coordinates then takes the form
M WC→norm.perps = M persp→norm ⋅ R ⋅ T
where R ⋅ T is the composite viewing transformation matrix (7.4).
Finally, clipping processes, visibility testing, surface rendering etc are performed followed by
the viewport transformation to screen coordinates, as in the case of oblique parallel projection.
See H&B p. 380-382. The rectangular viewport on the screen is set up as in the 2D case.
Note also that in (2D) screen coordinates, the z-coordinate value which is not used in any
position plotting, can be used for depth processing (variation of colour/intensities according to a
point’s depth in a scene etc). Typically, the screen buffer holds the (x,y) pixel values
corresponding to the 2D (x,y) position on the screen and additionally, a depth buffer holds the
depth cueing information for that position.
Viewing parameters (for the viewing transformation matrix) are supplied with:
These parameters are used in the subsequent setting up of the viewing transformation matrix.
133
The glOrtho(.) sets up a parallel projection perpendicular to the view plane (=near clipping
plane).
The call glOrtho(.,.,.,dnear = -1.0, dfar = 1.0) in a 2D situation is equivalent to calling
gluOrtho2D(.).
An oblique parallel projection function is not available in OpenGL (yet!). Can write your own by
using (7.13) – see sect. 7.7 above.
The frustum volume is symmetric about the Zview axis and the scene description is transformed
to normalized, homogeneous projection coordinates.
where
projection ref. point (Pref ) = VC origin (P0) – same as above
near clipping plane = view plane – same as above
xwmin,...ywmax = clipping window extents on near plane.
If xwmin = -xwmax, ywmin = - ywmax then have symmetric frustum with - Zview axis
as its centre line.
Default is orthogonal projection with symmetric cube as the view volume.
134
After the composite transformation (geometric + viewing), the clipping transformations routines
are applied in normalized coordinates, and then mapped to 3D screen coordnates (with Z-
coordinate is for depth information) with a viewport transformation.
To retain the same proportions of objects as in the WC scene, we set its aspect ratio as that of
the clipping window. When not invoked the default viewport is taken as the display window.
Multiple viewports may also be defined – see the 2D case.
glBegin (GL_QUADS);
glVertex3f (0.0, 0.0, 0.0);
glVertex3f (100.0, 0.0, 0.0);
glVertex3f (100.0, 100.0, 0.0);
glVertex3f (0.0, 100.0, 0.0);
glEnd ( );
glFlush ( );
}
Code outputs:
136
void init(void)
{
glViewport(0,0,width,height); // Make our viewport the whole window
glMatrixMode(GL_PROJECTION); // Select The Projection Matrix
glLoadIdentity(); // Reset The Projection Matrix
gluPerspective(45.0f,(GLfloat)width/(GLfloat)height, .5f ,150.0f);
glMatrixMode(GL_MODELVIEW); // Select The Modelview Matrix
glLoadIdentity(); / Reset The Modelview Matrix
glClearColor (1.0, 1.0, 1.0, 0.0);
}
Code outputs:
139
For the representation of solid objects, the methods used may be put into two classes
i) boundary representations (B-reps) – give a 3D object as a set of surfaces,
separating the object interior from its exterior
ii) space partitioning representations – used to describe interior properties by
partitioning the spatial region of an object with small, non-overlapping continuous
solids (cubes)
The polygon surface descriptions may be held in polygon tables, which may be a composition
of geometric data tables (containing vertices and other properties such as their orientation)
and attribute data tables (containing parameters for setting colour transparencies, surface
reflectivity, texture etc). Such a table was discussed in sect 3.9.4. Recall also, how to represent
the polygon surfaces by plane equations, in addition to which normal vectors may be used to
define plane surface orientations (sect 3.9.5).
Typically packages will provide mesh functions for tiling. Two common function types employ
triangles and quadrilaterals.
140
//Polyhedra.cpp
//-------------------
#include <windows.h>
#include <GL/glut.h>
GLsizei winWidth = 500, winHeight = 500; // Initial display-window size.
Output:
142
z z
P = ( x, y , z ) P = ( x, y , z )
r r
φ y φ y
θ r cos φ
Y θ r sin φ
x x
Ellipsoid
Z
rz
ry
Y
rx
X
This may be regarded as an extension of a sphere, with 3 radii, rx , ry , rz along each of the XYZ
axes and centred at the origin. Any point ( x, y , z ) on its surface satisfies the defining equation
2 2 2
x y z
+ + = 1. (8.4)
rx ry rz
Torus
This is a doughnut shaped object generated by rotating a circle (or other closed conic) about a
fixed axis. For example, a torus with circular cross-section centred at (0,0,0) in the XY plane
looks like:
144
(0, y , z )
rφ
Y O
O Y
θ
( x, y , z )
raxial
X
Side View (from X-axis) Top View (from Z-axis)
From the side view we see that the equation for the cross-section circle is
( y − raxial )2 + z 2 = r 2 .
When rotated about the Z-axis, the surface of a torus is generated, where a point ( x, y , z ) on it
satisfies the defining equation
( x 2 + y 2 − raxial )2 + z 2 = r 2 . (8.6)
In parametric form the equation for a torus with the above circular cross-section is
Other forms of tori are possible: for instance, we can rotate an elliptical cross-section instead of
a circular one.
8.5 Superquadrics
Surfaces may be generated by generalizations of quadratics with additional parameters,
resulting in “superquadrics”.
Superellipse
Based on the equation of an ellipse but in the generalized form:
2 2
x s y s
+ = 1 (8.8)
rx ry
where s is a real number. For s = 1 we have an ordinary ellipse. With rx = ry the following are
some shapes generated with the indicated s-values
145
Note also that the parametric form of the equation for a superellipse is
x = rx cos s θ
; − π ≤ θ ≤ π (8.9)
y = ry sin s θ
Superellipsoid
This is a generalization of the ellipsoid and is defined by the equation
s2
2
2 / s2
y
2 / s2 s1
z s1
x + + = 1; s1 , s2 = real (8.10)
r
x ry rz
s s π π
y = ry cos 1 φ sin 2 θ ; − ≤ φ ≤ , −π ≤ θ ≤ π (8.11)
s 2 2
z = rz sin 1 φ
For examples of shapes generated by the above see H&B p. 412. Superquadrics are useful
modelling components for composing more complex structures.
8.6 OpenGL, GLU and GLUT functions for quadric and cubic surfaces
1. glutWireSphere(r, nLongitudes, nLatitudes); or glutSolidSphere(...);
- constructs a sphere with radius r; nLatitudes,nLongitudes = number of latitude and
longitude lines. Surface is a quad mesh approximation. Defined in modelling coords with
centre at WC origin.
2. glutWireCone(rbase, height, nLongitudes, nLatitudes); or glutSolidCone(...);
- constructs a cone with axis along Z. Defined in modelling coords with centre at WC origin.
3. glutWireTorus(rCrossSection, rAxial, nConcentrics, nRadialSlices);
glutSolidTorus(...);
- constructs a cone by generating a circle in XY plane about Z axis. nConcentrics =
number of concentric circles with centre on Z axis for torus surface, nRadialSlices =
number of radial slices through torus surface. Centre is at WC origin and its axis is along Z.
4. glutWireTeapot(size); or glutSolidTeapot(...); ****NOT AVAILABLE IN OUR GLUT***
- constructs a teapot as a mesh of > 1000 bicubic surface patches (generated by Bezier
curve functions – see later). size = max radius for teapot bowl. Centre is at WC origin, and
vertical axis is along Y.
Historical note: The data set for the original teapot is due to Martin Newell (1975) developed
at the University of Utah. This and other well known 3D shapes are used to test many
proposed and published techniques for surface rendering in Computer Graphics, and now
form a classical set of test problems.
146
Other display modes are set with GLU_POINT (surface is a point plot), GLU_SILHOUTTE
(removes shared edges), GLU_FILL (patches are shaded filled areas).
Other shape drawing functions in place of gluSphere(...) can be:
gluCylinder(quadricName, rBase, rTop, height, nLongitudes, nLatitudes);
gluDisk(ringName, rInner, rOuter, nRadii, nRings);
gluPartialDisk(ringName, rInner, rOuter, nRadii, nRings, startAngle, sweepAngle);
In addition we have:
gluDeleteQuadric(quadricName); //frees memory when not required
gluQuadricOrientation(quadricName, normalVectorDirection);
//sets normal vector direction wirh
2nd param = GLU_OUTSIDE (sets normal out to front-face direction)
2nd param = GLU_INSIDE (sets normal in to back-face direction)
gluQuadricNormals(quadricName, generationMode); //generates surface normals
generationMode = GLU_NONE = no normals – typically no lighting applied
= GLU_FLAT = same normal for each face – same colour
= GLU_SMOOTH = normal vector for each surface vertex
gluQuadricCallback(quadricName, GLU_ERROR, function);
// a callback routine, calls function if error in creation etc
//QuadricSurfs.cpp
//-----------------------
#include <windows.h>
#include <GL/glut.h>
GLsizei winWidth = 500, winHeight = 500; // Initial display-window size.
where rk2 = xk2 + yk2 + zk2 , ak , bk are constants that adjust the amount of blobbiness (e.g.
bk < 0 gives dents instead of bumps) and T is a threshold parameter. A 3D bump, centred at 0
with standard deviation a and height b looks like:
−a a
composite shape
Other methods use density functions that fall-off to zero on a finite interval rather than decay
exponentially. For example, the metaball model is a combination of quadratic density
functions of the form:
3r 2
b 1 − 2 , 0 < r ≤ d / 3
d
3 r
2
f (r) = b 1 − , d / 3 ≤ r ≤ d (8.13)
2 d
0, r > d
A soft object model uses
22 r 2 17r 4 4 r 6
1 − + − ,0<r≤d
f (r) = 9 d 2 9 d 4 9d 6 (8.14)
0, r > d
150
A spline surface can be generated by two sets of spline curves, where each set intersects the
other orthogonally.
Cubic splines in Computer Graphics are generated by specifying a set of control points, and
choosing whether the curve passes through all the points (interpolation spline) or the curve
passes approximately through the points in a smooth sense (approximating spline):
The set of control points may be imagined to enclose a region of space which is a convex
polygon or convex hull, where every point is either on its boundary or lies inside.
q1 (u ) q2 (u )
u∗
Various degrees of smoothness can be obtained by setting continuity conditions for every pair
of curves q1 (u ) , q2 (u ) at the join u ∗ as follows:
q1 (u )
151
q1 (u ∗ ) = q2 (u ∗ ),
q1′(u ∗ ) ∝ q2′ (u ∗ )
iii) G (2) continuity (2nd order geometric continuity):
q1 (u ∗ ) = q2 (u ∗ ),
q1′(u ∗ ) ∝ q2′ (u ∗ ),
q1′′(u ∗ ) ∝ q2′′(u ∗ )
To understand the difference between these two types of continuity conditions, consider the
two-section splines on 3 control points, joined at the point P1 . The effect of imposing geometric
continuity is to pull the curve towards the slope of larger magnitude:
q2
q2
P0 P2 P0 P1 P2
P1
q1 q1
Specifying splines
A spline can be specified by
i) setting boundary conditions on each section, and/or
ii) incorporating boundary conditions into a characteristic matrix and specifying this
matrix, or
152
−1
a x 0 0 0 1 α
b 1 1 1 1 β
or, C = x= (8.20)
cx 0 0 1 0 γ
3 2 1 0 δ
d x 144 2443 {
M spline M geom
This form is also true for more general boundary conditions. We call M spline the basis matrix.
Finally, we can write the spline section (8.18) as
x (u ) = U ⋅ M spline ⋅ M geom
= α BF0 (u ) + β BF1 ( u ) + γ BF2 ( u ) + δ BF3 (u ) (by expansion) (8.21)
3
= ∑ g k ⋅ BFk ( u )
k =0
p k+1 p n −1
p1 p2 p k −1 pk pn
p0
PkR (u )
PkL (u )
Depending on the form of the boundary conditions chosen, different cubic splines can be
obtained:-
These give (4n-4) equations, but we require 4n conditions to find all the a x ,..., d x in (8.23 ).
We can fit the extreme endpoints
P1L (p 0 ) = a given value, Pn −1,R (p n ) = a given value (2 more equations) (8.25)
and set endpoint derivatives = 0
P1′L (p0 ) = 0, Pn′−1,R (p n ) = 0. (2 more equations) (8.26)
We thus have 4n-4 + 2 + 2 = 4n equations to find all the constants.
In a different approach, we can fit two more dummy points p -1 , p n+1 instead of setting
derivatives.
However, the natural spline obtained here, has a serious drawback in CG, in that, whenever a
local change is made (any control point is altered), the entire curve is affected; so other spline
types are generally preferred.
Hermite interpolation
Here, for each spline section (e.g. cubic piece) we specify
• function values at each control point
154
Cardinal splines
Quite often, it’s not easy to supply or obtain slope information at the control points, so that
Hermite interpolation becomes infeasible. Then cardinal splines are useful, with the following
defining properties:
• This is a piece-wise interpolatory cubic with tangent values specified at the endpoints
of each curve section, where these values are not given but calculated from
neighbouring position values:
P( u )
pk+1
pk
p k −1 pk +2
• The cardinal spline section P(u ) is completely specified with 4 control points by:-
P(0) = pk
P(1) = pk +1
1
P′(0) = (1 − t ) ( pk +1 − pk −1 ) (8.27)
2
1
P′(1) = (1 − t ) ( p k+2 − pk )
2
i.e. the slopes are chosen proportional to the slopes of the chords pk-1pk +1
and pk pk +2 . The parameter t is a tension parameter which controls how
loosely or tightly the curve section fits the control points:
slope pk+1
pk P( u )
slope
p k −1 pk +2 t < 0: t > 0:
looser fit tighter fit
When t = 0 the class of curves are called Catmull-Rom splines or Overhauser splines.
As before, we can convert the boundary conditions to matrix form and write the spline section
as
pk-1
p
P(u ) = u 3 u 2 u 1 ⋅ M C ⋅ k (8.28)
pk +1
pk +2
where the cardinal spline matrix is
155
− s 2 − s s − 2 s
2 s s − 3 3 − 2s − s
MC = ; s ≡ 1 (1 − t ). (8.29)
− s 0 s 0 2
0 1 0 0
Examples of what they look like graphically are depicted below for t = 0 or s = 0.5:
1 1
CAR1 (u )
CAR0 (u )
0 u 0
u
1 1
1 1
CAR2 (u )
CAR3 (u )
0 u 0 u
1 1
A Bezier curve section can be fitted to any number of control points and can be specified by
156
• boundary conditions
• a characteristic matrix, or
• blending functions
As a rule, a Bezier curve = a polynomial of degree 1 less than the number of control
points.
For example,
3 control points ⇒ Bezier curve = a parabola
4 control points ⇒ Bezier curve = a cubic
p1 p2
p3
p1
p2
p0 p0 p3
p2 p4
p1 p0
p0 p3 p2
p1
∑ BEZ
k =0
k ,n (u ) = 1 ∀n and 0 ≤ u ≤ 1. (8.40)
Then any position on the curve is a weighted sum of its control points. This
property ensures that the polynomial tends to smoothly follow the control points
without erratic oscillations.
p3 p1 = p 2 p3
p2
p0
p1 p4 p4
p0 = p5
or, a curve may pulled towards a point by choosing a double point there (2nd figure
above).
ii) When there is a large number of control points to be fitted, it’s unwise to use a
single high-degree Bezier curve, since this will introduce too many (unwanted)
oscillations, apart from the computational effort required to solve the polynomial
coefficients. Rather, we fit lower order piece-wise Bezier curves and join them
appropriately to obtain the required order of overall continuity.
iii) C 0 continuity automatically follows at the common control endpoint since both
Bezier sections pass through that point:
p1
C1
p2
p0 p′0
C2
p′3
p′1 p′2
For example, in the above one Bezier section ( C1 ) with control points p 0 ,p1 , p 2
joins another ( C2 ) with control points p′0 , p′1 ,p′2 ,p′3 at p 2 = p′0 giving continuity
at this common endpoint.
iv) To obtain C (1) continuity we select the 2nd point of the second section ( C2 ) in such
a way that the p1 − p 2 slope is the same as the p′0 − p1′ slope by setting:
p′0 = p 2
(8.41)
p1′ = p 2 + (p 2 − p1 )
159
This process must of course, be implemented on all the common control points, as
well as the composite endpoints, if required. Overall C (1) continuity also requires
that both sections must each have the same number of control points.
v) By similar means we can obtain C ( 2) continuity, but the latter is too restrictive
especially on cubic Bezier sections, so it is not usually administered.
BEZ 0,3 (u ) = (1 − u )
3
BEZ1,3 (u ) = 3u (1 − u )
2
(8.42)
BEZ 2,3 (u ) = 3u 2 (1 − u )
BEZ 3,3 (u ) = u 3
1 1
BEZ 0,3 (u ) BEZ1,3 (u )
0 u 0
u
1 1
1 1
BEZ 2,3 (u )
BEZ 3,3 (u )
0 u 0 u
1 1
Notice that
• at u = 0 only BEZ 0,3 (0) ≠ 0 , whilst all the others = 0.
• at u = 1 only BEZ 3,3 (1) ≠ 0 , whilst all the others = 0.
• consequently, from (8.37) P(0) = p 0 , P(1) = p3 i.e. Bezier curve passes thro’
endpoints
−1 3 −3 1 p0
3 −6 3 0 p1
P(u ) = u 3 u2 u 1 (8.45)
−3 3 0 0 p 2
1 0 0 0 p 3
M BEZ
//BezierCurve.cpp
//---------------
//H&B p. 435
//Computes and draws the Bezier curve thro' 4 control pts in WCs
#include <windows.h>
#include <GL/glut.h>
#include <stdlib.h>
#include <math.h>
class wcPt3D {
public:
GLfloat x, y, z;
};
GLint k, j;
for (k = 0; k <= n; k++) {
/* Compute n!/(k!(n - k)!). */
C [k] = 1;
for (j = n; j >= k + 1; j--)
C [k] *= j;
for (j = n - k; j >= 2; j--)
C [k] /= j;
}
}
glPointSize (4);
glColor3f (1.0, 0.0, 0.0); // Set point color to red.
162
Code outputs:
For a Java program that allows the user to insert the control points with mouse clicks see
Ammeraal’s Bezier.java.
m n
P(u, v ) = ∑∑ p jk BEZ j ,m ( v ) ⋅ BEZ k ,n (u ) (8.46)
j =0 k = 0
glDisable(GL_MAP1_VERTEX_3);
Here
• use suffix (*) f or d for type of values, uMin/uMax usually 0/1,
• nPts = number of control pts (> 0) in coordinate position array ctrlPts
• stride = integer offset = number of data values between beginning of one coordinate
position in ctrlPts and the beginning of the next. For just 3D coordinate positions set
stride=3. For 4D homogeneous coords set stride=4 and set symbolic const to
GL_MAP1_VERTEX_4.
2. glEvalCoord1*(uValue);
- evaluates positions along the spline path where uValue is a value in range uMin...uMax.
This function maps uValue to a u in the range 0...1.0 by means of the mapping
u − umin
u = value (8.47)
umax − umin
and then computes points on the curve using (8.32). It also generates a glVertex3 function
as u is processed. Repeated calls of glEvalCoord1*(.) are required to produce points along
the curve, which can be joined together by straight-line segments, approximating the curve.
Non-uniformly spaced points can be generated with this function as well.
The following code segment shows how to employ this process (for a full program see
OGLBezierCurve.cpp):
Then repeated calls to glEvalCoord2*(.) will display the Bezier surface with
glVertex3*(.) being generated in the process.
A sample code segment for generating a cubic Bezier surface is (from H&B):
7. GLU B-spline functions for curves and surfaces are also available. See H&B p. 467 –
473.
control points
167
8.20 Octrees
In this approach, a hiearchical tree structure (octree) is used to represent a solid object. See
H&B p. 476-479.
The amount of detail variation in the object is described by a number, its fractal dimension,
not necessarily an integer (can have fractional dimension).
Classification of fractals
a. self-similar fractals – here parts of the whole are scaled down versions of the
whole. We can use the same scale factor s for each sub-part or different ones.
If randomly different s-values are used then the fractal is said to be
statistically self-similar. For example, the latter is used to model trees,
shrubs etc.
b. self-affine fractals – here parts are formed with different scaling factors
sx , s y , sz in each of the coordinate directions. Again, random variations of
168
Fractal dimension
Detail variation in a fractal is described by a number D = the fractal dimension = a measure of
roughness or fragmentation.
Sometimes we can set a D-value and determine procedures for generating a fractal. At other
times we can only find D from properties of the constructed object.
For self-similar fractals we can obtain expressions for the fractal dimension in terms of the
scale factor s by analogy with Euclidean subdivision:
For example, choose s = ½ as a scale factor and perform repeated subdivisions of the following
basic shapes
After 1 subdivision n=2 parts
Euclidean dim DE = 1 Scale factor s = 1/ n
1 unit at start
A′ = A / n
A
V′ =V /n
V
Now, in analogy with the dimension for Euclidean objects, for a self-similar object we take the
fractal dimension D to satisfy
ns D = 1 (8.50)
169
ln n
or D= (8.51)
1
ln
s
For self-similar fractals constructed with different scaling factors D is obtained from the implicit
relation
n
∑s
k =1
D
k = 1; sk = scale factor ~ kth sub-part (8.52)
Note also, that the fractal dimension is always > the Euclidean dimension. For instance:
• Euclidean curves have DE = 1 , but a fractal curve generally has dimension
1 < D ≤ 2 , with D > 2 also possible. The closer D is to 1 the smoother it is. Other
values of D give more interesting possibilities: D = 2 gives a Peano curve that
completely fills a finite 2D region. For 2 < D < 3 the curve self-intersects, infinite
number of times – can be used to model natural boundaries e.g. shorelines.
• Fractal surfaces generally have dimension in the range 2 < D ≤ 3 with D > 3 also
possible, giving overlapping coverage – useful for terrain, clouds, water modelling.
• Fractal solids generally have dimension 3 < D ≤ 4 with D > 4 also possible, giving
self-overlapping solids – useful to depict vapour density, temperature in a region etc.
initiator generator
ln 4
The fractal dimension of the snow-flake is then D = 1.2619 . Further iterations produce
ln 3
more complicated shapes – see H&B p. 485.
In the above we can use random choices of generators or calculate coordinate positions
randomly, so that the result is an irregular and more realistic shape – used to model trees,
plants etc.
Y Y
y ( b) y ( b)
y(a ) y(a )
ymid
shape after
X X
a b a a+b b several
iterations
2
For generating terrain features we can start with a ground plane and randomly elevate the
midpoint, and continue in this way for each sub-part until several iterations are completed.
Such models can be used to depict small variations in a surface or even large variations
(mountainous terrain). See H&B p. 492-4.
Y Y
Z
e
a b
Y
ground h m f
plane
X d g c X X
Here the midpoint elevation zm may be computed by using the corner elevations plus a
random offset like:
1
zm = ( za + zb + zc + zd ) + rm (8.54)
4
Again, after several iterations a granular topography results.
Self-squaring fractals
171
Y
escaped sequence
z4
trapped sequence
z3 ′
z2 z200
COMPLEX PLANE
z3′
z1 z2′ z4′
−r X
z0 z1′ r
r =10, may be OK
z0′ some
rectangle like:
−2 ≤ x ≤ 2
−2 ≤ y ≤ 2
Choose a starting point z0 = x0 + iy0 within the rectangle and some complex constant
µ = a + ib and then generate the iterated function sequence (IFS) of (8.55) by,
z1 = z02 − µ
z1 = z02 − µ
. (8.57)
....
zn = zn −1 − µ
2
172
Now if for some n, we find zn > r we say the point z0 has escaped. The set of all starting
points z0 (e.g. z0′ shown above) which do not escape is commonly called the Julia set of the
complex polynomial F ( z ) in (8.55). Actually, the exact Julia set is the boundary of this
collection of points, but we’ll speak loosely of it.
Coding it:
In our program Julia.cpp we have used a rectangle of size [-2,2] × [2,2] and a radius of just 4,
with a=1.0 and b=0.0 in WCs. The starting points are chosen by sweeping through the window-
coordinate range in pixels values ( x p , y p ) and transforming to the WC ( x, y ) range above.
The transformation can follow from our previous theory or we can apply first principles based
on:
xmin x0 xmax
real WC X axis
-2 0 2
pixel Xp axis
0 xp maxxPix
Then a straight forward linear transformation sending xp→ x0 is given by the proportionate
lengths expression:
xp-0 x0-xmin (xmax-xmin)xp
= ⇒ x0 = + xmin (8.58)
maxxPix-0 xmax − xmin maxxPix
In our code we have also set up a 15-colour-value array of RGB values, which we use to colour
the escaping pixel.
Thus, for each pixel coordinate (xp,yp) in the pixel coordinate range, we map to the complex
value (x,y) in WCs, applying the transformation repeatedly, so that if any |(xn,yn)| > r, we light up
the escaping pixel (xp,xp). When done for the pixel coordinates in the window range, the non-lit
pixels (in the background colour) comprises the Julia set. Symmetry is used to avoid
duplicating calculations.
The output from Julia.cpp looks like:
173
Then the set of all µ = a + ib for which the origin escapes (i.e. zn > r for some n) is called
the Mandelbrot set.
This program may be coded by using real arithmetic, as in Julia.cpp. Do it as an exercise. Note
that H&B construct and employ a complex number class with overloaded operators for the
usual complex operations. You must avoid their approach and use ours in this exercise.
Self-squaring fractals
An example of generating fractals with (8.56) above is H&B’s SelfSquareFractal.cpp which
first obtains the inverse of this transformation. See H&B p. 496-499.
174
How to determine what is visible in a scene from a particular viewing position is of concern in
this chapter. For instance, which faces of an object can we see from a particular vantage
position, and what should appear obscured from this position is answered here. The techniques
involved are variously called visible surface detection methods or hidden surface removal
methods, although they may not mean exactly the same thing.
Typically, such methods work either on the object descriptions (object-space methods) or on
the pixel data (image-space methods), although a combination of both have also been
devised.
If the Vview is taken as a vector in the − zv direction with unit vector zˆ v in a right-handed
system, and object descriptions have been transformed to view-projection coordinates then the
face (9.1) is a back face if
zˆ v ⋅ N = C ≤ 0 . (9.3)
This object-space test may be used to eliminate (or cull) several but not all hidden surfaces or
parts of them, since for example some may only be partially obscured. Thus additional tests
beyond this process are usually required.
176
yv
N N
Vview
xv
camera Vview
zˆ v
zv only partly hidden
A
i.e. as z′ = z − (9.6)
C
Note that since –A/C is constant for the whole surface, succeeding z-values are found with just
one addition.
If we start with a depth value at a top vertex, we can calculate the depth values down an edge
as follows:
The beginning x-value on the next scan line is obtained from the starting x-value of the
previous scan line by (see fig. below)
1
x′ = x − ; m = edge slope
m
y adjacent scan
y-1 lines
It turns out that even with the above coherence based approach, many unnecessary
calculations are performed in this method. Thus other, better techniques have been devised.
2. Scan convert the surfaces in order starting with the surface of greatest depth (in image
space)
The process is similar to the way in which an artist will paint on a canvas – first the background
and then the foreground over it. Here the colour values for the furthest surface are entered into
the frame buffer. Then for each succeeding surface we overwrite the FB values.
With the above simplistic description, the method can fail when two or more surfaces overlap in
the viewing direction ( zv ). Thus further, more refined tests are necessary to resolve such
issues. See H&B p. 538-540. A variation of this algorithm which sub-divides the surfaces into
constitutive triangles, and which then conducts depth tests on the triangles is given in
Ammeraal p. 145. The latter can be applied to non-polygonal surfaces which would then be
broken down into polygons (triangles).
In summary,
• If both line endpoints are behind the surface then the line is entirely hidden
• If both endpoints are in front of the surface, the line is entirely visible
• If one endpoint is behind and another is in front, we calculate the intersections of the
line with the surface, find the depth values there and compare with the endpoint depths
to determine which sections of the line are hidden or visible.
4. glClear(GL_DEPTH_BUFFER_BIT);
//Initialize depth-buffer values. Clears all values to max 1.0 in norm range
0...1.0. Should be done for every new frame display. Then use
5. glEnable(GL_DEPTH_TEST);
//To enable depth-test routines
180
6. glDisable(GL_DEPTH_TEST);
//Disables them
glEnable (GL_DEPTH_TEST);
glPolygonMode (GL_FRONT_AND_BACK, GL_LINE);
glColor3f (1.0, 1.0, 1.0); //white foreground
/* Invoke the object-description routine. */
glDisable (GL_POLYGON_OFFSET_FILL);
12. OGL depth-cueing to vary brightness of object as function of distance from view
point invoked with:
glEnable(GL_FOG); //selects display mode
glFogi(GL_FOG_MODE, GL_LINEAR); //applies (9.9) with d min = 0.0, d max = 1.0
Can set own min and max values for d by
glFogf(GL_FOG_START, minDepth); //minDepth = float in range 0..1.0
glFogf(GL_FOG_END, maxDepth); //maxDepth = float in range 0..1.0
181
This point position is used in the model to find the effect of the light rays on objects in the
scene.
Intensity attenuation
From physics, it is known that the intensity (amplitude) of light drops by a factor 1/ d l2 where
d l is the distance from the light source. For a model using a point source, this can produce too
large a variation for scene parts near the source, but too little for parts that are far away. To
compensate for this effect, a linear term is added in the amplitude/intensity function to give the
form:
1
f radiation ( d l ) = (10.1)
a0 + a1d l + a2 d l2
where the values of the coefficients a0 , a1 , a2 can be suitably adjusted, even allowing for a
different set of values for each point in a multiple-point source. Further, since we cannot use
(10.1) for a source at infinity (since d l = ∞ ), knowing that all points in the scene would be
equally illuminated, we apply the adjusted form
1.0, if source is at infinity
f radiation ( d l ) = 1 (10.2)
a + a d + a d 2 , if source is near(local)
0 1 l 2 l
α Vlight
CONE OF
INFLUENCE
θl
light source
objects unlit here
Here if Vlight is a unit vector in the direction of light from a local source, the light has a cone of
influence, with angular extent θ l on either side of an axis of a cone in direction Vlight . If α is
the angle between the unit vector Vobj pointing to its location (vertex) and Vlight then
Restricting the cone’s angular extent to 00 < θ l ≤ 900 then object will be within the spotlight if
183
To allow for angular attenuation of a directional light source we can use an attenuation
function, such as for example,
f ang .atten (α ) = cosal α , 00 ≤ α ≤ θ (10.6)
where α is the angular distance from the cone axis and al is a positive attenuation coefficient
applicable when α > 0 (away from the axis where the intensity has a max value of 1.0). Then
to allow for both ordinary point sources and spotlights, we can employ the attenuation function
1.0, if source ≠ a spotlight
f l .ang .atten (α ) = 0.0, if Vobj ⋅ Vlight = cos α < cosθ l (obj. outside spotlight cone) (10.7)
( Vobj ⋅ Vlight ) , otherwise
al
Ambient light
A simple model for ambient or background lighting is to set an overall intensity level parameter
I a say, for the whole scene, resulting in the equal illumination of all surfaces in the scene, in all
directions, producing equally scattered or diffuse reflections from the surfaces.
Diffuse reflection
As in above, surfaces that reflect light equally in all directions and independent of the viewing
position are called ideal diffuse reflectors. A more realistic model will allow for intensity of the
reflected light to be dependent on the incident light (from a point-source) intensity I l as well as
the direction of incidence to the surface. Thus another useful diffuse-reflection model is based
on the equation:
where k d is a constant called the diffuse-reflection coefficient (or diffuse reflectivity) and θ the
angle of incidence.
184
N
N
θ
θ θ
incident L
light surface
Thus here, when the incident light hits the surface perpendicularly, θ = 00 , the diffuse
reflection has a maximum value of I l ,diff = k d I l . For angles away from zero, it decreases until
it reaches 0 at θ = 900 (grazes the surface with no reflection). For cosθ < 0 the source is
behind the surface.
More generally, if N is the unit normal at a position on the surface, and L is the unit vector
from this position to a point source then cosθ = N ⋅ L and the diffuse reflection equation for
single point-source illumination at a surface position can be written as
k I ( N ⋅ L), if N ⋅ L > 0
I l ,diff = d l (10.9)
0.0, if N ⋅ L ≤ 0
When the point-source is not at infinity, i.e. it is a nearby one, then for its direction vector L
above we can use (with obvious meanings for the symbols),
P -P
L = source surf (10.10)
Psource -Psurf
It is possible to combine the effects of the ambient and point-light-source intensities so that the
total diffuse reflection due to a point-source in an ambient background light is
k a I a + k d I l ( N ⋅ L), if N ⋅ L > 0
I diff = (10.11)
k a I a , if N ⋅ L ≤ 0
Typically, ka , k d are set in the range 0.0 ... 1.0 for monochromatic light and depend on the type
of surface material.
N N N
L R L L
θ θ φ V
R R
surface
shiny surface: dull surface:
ns ~ 100 ns ~ 1
On some shiny surfaces a bright spot or localized area is noticed when light shines on the
surface. When viewed from a more oblique angle it is not noticeable. This effect is called
specular reflection. To model it we introduce a specular reflection direction vector R which has
the same angle ( θ ) w.r.t. the unit normal as the vector L but on the opposite side.
185
Now if V is the vector to the viewer, and φ the specular reflection angle of the viewer w.r.t. R
then for an ideal reflector (mirror), reflected light is seen only when φ = 0 i.e. when viewed
along V = R . For a non-ideal surface, there is a range of viewing angles around R inside of
which specular reflection is seen. The model due to Phong proposes specular reflection with
intensity proportional to cosns φ for angles 00 ≤ φ ≤ 900 with ns called the specular reflection
exponent (see fig. above for its effect). Thus we can write for the reflected intensity
I l ,spec = W ⋅ I l cosns φ ; W = const.
For real materials, W the specular reflection coefficient, is a function of the incidence angle
00 ≤ θ ≤ 900 described by Fresnel’s law of reflection. Thus we write for the reflected intensity
For graphs of W (θ ) for various materials see H&B p. 569. These show that for transparent
materials (glass-like) specular reflection is appreciable around θ ≈ 900 and very small for
other angles of incidence. For opaque materials, W (θ ) ≈ k s = constant for all θ . So, in
practice it is common to set W (θ ) ≈ k s = a value in the range 0.0 ... 1.0 for each surface. Thus
a simplified specular reflection intensity expression can be obtained by noting:
• N, L and R = unit vectors
• cos φ = V ⋅ R , cos θ = N ⋅ L
• if V and L on same side of normal N ⇒ no specular reflection
• for specular reflection must have cos φ = V ⋅ R > 0 and cos θ = N ⋅ L > 0
Then,
k s I l ( V ⋅ R ) ns , if V ⋅ R > 0 and N ⋅ L > 0
I l ,spec =
0.0, if V ⋅ R < 0 and N ⋅ L ≤ 0 (10.13)
L N
N⋅L N H
L R
N⋅L L R V
it is clear that
R + L = 2( N ⋅ L) N ⇒ R = (2N ⋅ L) N − L . (10.14)
To apply the above, recall that L is calculated by taking a vector from the surface point to the
light source and then normalizing it. Similarly, V is got by taking a vector from the surface point
to the view point. Usually the latter is fixed so as an easy choice we set V = (0,0,1), the unit
vector in the +ve Z direction, as a good compromise.
A less compute-intensive amendment to the Phong model above is to use the halfway vector,
186
L+V
H= (10.15)
L+V
and replace V•R with N•H (see H&B p. 571 for an explanation).
Further, the basic models for diffuse and specular lighting may be combined to give a
composite model. In addition if there is more than one point light source we can form a
combined light source model by including additive terms for the individual ones. For example,
using (10.15) a combined diffuse, and n-source specular reflection model uses an intensity
function such as (see H&B p. 571)
n
I = I ambdiff + ∑ I l kd ( N ⋅ L) + k s ( N ⋅ L) ns . (10.16)
l =1
In the same way, it is also possible to include terms allowing for light emitters at various points
in the scene.
Colour
In the above, the intensity function is based on the assumption that we have monochromatic
lighting. To apply it to say, an RGB colour model, we employ a separate I function for each
component with correspondingly different constants kd , k s for each primary colour. Other
effects that can be incorporated into this basic model are transparency, opacity, translucency,
refraction and the effect of fog (haziness). See H&B p. 573-591 for a discussion of these
issues.
∑N k
NV = k =1
n
(10.17)
∑N
k =1
k
In 3 above, for each polygon face, the linear interpolation is done by considering the
intersections of scan lines with the polygon edges.
y
NV
3
1 p
V 4 5
2
x
For example, the intensity at the intersection point 4 is obtained by using only vertical
displacements via
y − y2 y − y4
I4 = 4 I1 + 1 I2 (10.18)
y1 − y2 y1 − y2
This applies to one of the 3 RGB components. So repeat for the other two with different I’s.
Similarly also we obtain I 5 .
Then the intensity for a point p along the scan line is computed from the equation
x − xp x − x4
Ip = 5 I4 + p I5 (10.19)
x5 − x4 x5 − x4
These calculations may be performed incrementally, as recurrence relations (see H&B p. 594)
for efficiency. For an example of how the method performs see H&B p. 594.
The normal vector at a point of intersection with a scan-line intersection (e.g. 4 above) is
obtained by vertically interpolating the normals at 1 and 2 as before by
y − y2 y − y4
N4 = 4 N1 + 1 N2 (10.20)
y1 − y2 y1 − y2
We proceed similarly to the Gouraud method for interpolation along scan lines but here for
normals. Finally we apply the illumination model for each projected pixel position.
188
This method is more accurate but more time consuming. A faster version of it is also used (see
H&B p. 596).
Computer animation is a process for simulating the time evolution of a scene and its depiction
on a viewing device. Such techniques are employed in the development of movie scenes,
computer games, scientific modelling, embedded video tracks in web pages, Java applets etc.
The methods used are broadly classified into two basic types:
• real-time animation where each stage of a sequence is displayed as it is created. The
screen refresh rate must be faster than the generation rate for compatibility. This is
generally employed for “simple” applications, except when user-interface controls are
also required (e.g. a flight simulator program).
• frame-by-frame animation where each frame of the motion is generated separately
and stored for later playback as a sequence. Generally used for large applications, e.g.
long movies and large-scale scientific simulations.
Double buffering
Here two refresh buffers are employed in alternation. Whilst one is offloaded to the screen,
the other receives the next frame construction from the processor. Then in the next cycle, the
latest updated buffer can be painted to the screen, whilst the first receives the frame
construction and so on. Which one is to be offloaded to the screen is signalled by some routine
call which swops their roles and is usually called at the end of a refresh cycle, which typically
189
takes 1/60 sec (i.e. 60 frames/sec painted), but can be called at will. If the frame being
constructed is completed before the other is offloaded, then the overall sequence is
synchronized with the refresh rate. Otherwise, if the frame construction time > 1/60 sec, then
two or 3 refresh cycles will occur from the same buffer, resulting in the same picture, before the
latest one gets offloaded. The animation sequence rates tend to be irregular when frame
construction times are multiples of the refresh time i.e. when they are near 1/60, 2/60, 3/60, ...
One way to minimize this effect is to put a delay in the construction sequence or we update
only some object motions at a time in the scene for each cycle.
Since graphics terminals are refreshed at a rate of ±60 frames/sec a 1-minute animation
sequence would require about 60×60 = 3600 frames to be constructed. Typically, about 5 in-
between frames per key frame are used, so that 720 key frames are needed.
Two important processes involved in animation sequences are changing object shapes
(“morphing”) and simulation of movement of objects (accelaration).
Examples of parts of animation sequences for morphing using linear interpolation are:
1 1'
2 3'
added
point
A simple rule or algorithm for performing the above interpolation process is given in H&B p.
741, wherein are also presented more complex pictures of morphing like those seen in motion
pictures.
The simulation of acceleration can be achieved by curve fitting techniques, employing linear or
nonlinear functions. An example is
key frame k in-between frame key frame k+1 key frame k+2
.
Here given vertex positions at the key frames, we fit a curve through them. We then advance
these positions to the next frame through a series of in-between frames. If different positions
move at different speeds (along the red dotted paths) then the curves will deform as time goes
on. For constant speed or zero acceleration the time intervals between frames is set to some
constant value ∆t which can be computed from physical considerations. For accelerating
frames we can keep adjusting ∆t .
d
F= ( mv ) or Force = mass × acceleration (for m =const)
dt
For simulating the motion of living creatures (humans, animals etc) a common technique
involves constructing frames of articulated figures like:
The articulated figure may be further enhanced by limb models as above, which are made to
bend at the joints as motion is simulated in succeeding frames.
class scrPt {
public:
GLint x, y;
};
11.5.2 A bouncing ball: This is an example of motion determined by the laws of mechanics:-
In particular, a ball bouncing on a hard surface undergoes damped harmonic motion
according to a differential equation (Newton’s 2nd law) determining its height y(x(t)),
where x(t) is its X-displacement as a function of time, and whose solution for y is:
y ( x ) = A [sin( wx + θ 0 ) ] e − kx .
194
Here A is the amplitude, w the angular frequency, θ 0 a phase constant and k the
damping constant. The ball’s trajectory then looks like:
The following is a Java program that computes and simulates this motion
//Ball.java
//------------
//Simple bouncing ball animation without double buffering
import java.awt.*;
import java.awt.event.*;
CONCLUDING APPENDIX
GOOD LUCK