Matt Bross - March 15th, 1997

3D Programming: An Explanation of the Basics as Well as Some Advanced 
Techniques -- With Examples -- in the BASIC Programming Language

	We have all seen the effects of 3D animation, from games like DOOM to the 
movie Toy Story.  Wouldn't it be cool to be able to make your own 3D stuff!   In this little 
paper I hope to explain how it is possible to do this on your own home computer.

Converting 3D to 2D 

	This is the easy part and the first thing you need to know.  To convert 3D 
coordinates to a 2D screen use this: X, Y, and Z are 3D coordinates on a Cartesian plane, 
and SX and SY are the final screen coordinates.

SX = X / Z
SY = Y / Z

	This is very logical and easy to understand, but not something that you would just 
think of.  The farther the Z, the closer the points will be together.   With different 
variables, if the viewer is at the origin, you must give each of the variables its own 
coordinate system, relate them with variables and, then, modify the equation.  So, do it 
like this:  The screen coordinates of  programming language are not oriented around (0, 0, 
0), the origin, so you must add half the screen resolution to the end of the equation to 
center the object on the screen.  Using the same variables, VX, VY, VZ as the position of 
the viewer, and SRX and SRY as the screen's resolution.

SX = ((X - VX) / (Z - VZ)) + (SRX / 2)
SY = ((Y - VY) / (Z - VZ)) + (SRY / 2)

	Now, if you define points and then draw lines between them, this will draw a 3d 
shape; not very impressive but definitely a start.  (There are other ways of converting 3D 
to 2D (ortho, where coordinates (X, Y, Z) become (X,Y), and parallel projection, where 
(X, Y, Z) becomes (X + Z, Y + Z), are two) but the way I showed you, perspective 
projection, is the best.  The more impressive things are rotation and translation.  

Rotation and Translation

	Now, for rotation,  Everything up until now has been pretty easy if you think about 
it, unless math really isn't your thing.  Now comes the first, relatively speaking, hard part.  
A rotation rotates all the points, and a translation moves all points over in the same 
direction and the same amount.  Translations are simpler, so I'll describe them first.  
Simply add a translation value to every point's respective X, Y, or Z coordinate, 
depending on what direction you are translating.  When programming, you will want to 
store this in a separate variable to make the object's own coordinate system.  If you don't 
understand how this gives the object its own coordinate system, then think about this:  if 
you don't use the translation values, then don't add the amounts of translation and you are 
left with an object surrounding  the origin.  In order to save space, SX equals screen x; SY 
equals screen y; X, Y, Z means virtual coordinates x, y, z; TX, TY, TZ are translations 
along the x, y, and z axes; and SRX, SRY are the screen's resolution; and VX, VY, VZ 
are the coordinates of the viewer.

SX = ((TX + X -VX) / (TZ + Z - VZ)) + (SRX / 2)
SY = ((TY + Y - VY) / (TZ + Z - VZ)) + (SRY / 2)

Translations can also be done with matrices following, but it is more complicated.   This  
also gives an introduction to matrix mathematics (see appendix D), which is used for 
rotation.

(    1,     0,      0, 0) = T1  
(    0,     1,      0, 0)
(    0,     0,      1, 0)
(-TX, -TY, -TZ, 1)

The same matrix can also be expressed in spherical coordinates (which are used for 
rotation and are explained next) to keep everything in the same form.

(                                     1,                                      0,                      0, 0) = T1  
(                                     0,                                      1,                      0, 0)
(                                     0,                                      0,                      1, 0)
(-roe(COS theta)(SIN phi), -roe(SIN theta)(SIN phi), -roe(COS phi), 1)

For rotation you must use spherical coordinates.  Spherical coordinates are coordinates 
that are expressed as distances from the origin and in angles, rather than in distances from 
the X, Y, and Z axes.  Spherical coordinates are noted (theta, phi, roe) as roe(a Greek 
letter) describes the distance from the origin, theta (another Greek letter) is the angle from 
the XY plane and phi (yet another Greek letter) is the angle from the Z axis. The 
conversions between the standard Cartesian plane and the spherical coordinates are:

	X = roe (SIN theta)(COS phi) 		theta = arcTAN (X / Y)
	Y = roe (SIN theta)(SIN phi)		phi = arcCOS (Z / Roe)
	Z = roe (COS theta)			roe = SQR(X ^ 2 + Y ^ 2 + Z ^ 2)

For 3D animation, the distance from the origin never changes, so the only variables are 
theta and phi.  So, for rotation around just one axis you must first find each point's 
distance from the origin.  Using the equation above, which is derived from the 
Pythagorean theorem (in the case that a and b are legs of a right triangle and c is the 
hypotenuse of the triangle, a ^ 2 + b ^ 2 = c ^ 2, but since there is another variable - depth 
- you must add another variable to the equation that I will call "e" (this is not official).  So, 
the equation becomes (a ^ 2 + b ^ 2 + e ^ 2) = c2.   If this is on a 2D plane e = 0 and e ^ 2 
= 0. So, it has no use, but in 3D, e  0, all the time.)  Using this theorem, convert the point 
to spherical coordinates, once again using the equations above, add the amount of rotation 
to phi and then convert it back to normal Cartesian coordinates and plot  the point on the 
screen.  This can be done with the following 4 * 4 matrices.

Z axis clockwise rotation

(   SIN theta,  COS phi, 0, 0) = T2
(-COS theta, SIN theta, 0, 0)
(                0,             0, 1, 0)
(                0,             0, 0, 1)

X axis counter-clockwise rotation

(1,              0,            0, 0)
(0, -COS phi,  SIN phi, 0)
(0,    SIN phi, COS phi, 0)
(0,              0,             0, 1)

To get the final matrix that will be used, and converted into 2D rectangular and plotted on 
the screen, multiply all the preceding matrices together.  F = T1 * T2 * T3.

F = ( -SIN theta, -(COS  theta)(COS phi), -(COS theta)(SIN phi), 0)
      ( COS theta,    -(SIN theta)(COS phi),  -(SIN theta)(SIN phi), 0)
      (               0,                           SIN phi,                    -COS phi, 0)
      (               0,                                     0,                              roe, 1)

So  the original (X, Y, Z, 1) translated and rotated = F * (X, Y, Z,1) = (Xf, Yf, Zf, 1).  

Xf = -X(SIN theta) + COS theta
Yf = -X(COS theta)(COS phi) - Y(SIN theta)(COS phi) + Z(SIN phi)
Zf = -X(COS theta)(SIN phi) - Y(SIN theta)(SIN phi) - Z(COS phi) + roe

A simpler matrix might  be

Rotation around the Y axis where Yan is the angle of rotation around the Y axis
(COS Yan, 0,   SIN Yan) = T1
(             0, 1,               0)
(-SIN Yan, 0, COS Yan)

Rotation around the Z axis where An is the angle of rotation around the Z axis
(1,              0,             0) = T2
(0, COS An, -SIN An)
(0,   SIN An, COS An)

Rotation around the X axis where Xan is the angle of rotation around the X axis
(COS Xan, -SIN Xan, 0) = T3
(-SIN Yan, COS Yan, 0)
(             0,              0, 1)

The final Matrix (T1 * T2 * T3) where S1 = SIN Yan, C2 = COS An, etc. is  
(  C1C3 + S1S2S3, -C1S3 + C3S1S2,  C2S1) = F
(                   C2S3,                    C2C3,     -S2)
(-C3S1 + C1S2S3,   S1S3 + C1C3S2, C1C2)

So 

Xf = x(s1s2s3 + c1c3) + y(c2s3) + z(c1s2s3 - c3s1)
Yf = x(c3s1s2 - c1s3) + y(c2c3) + z(c1c3s2 + s1s3)
Zf = x(c1s2s3 - c3s1) + y(-s2)  + z(c1c2)

or simpler and faster (some operations are repeated) where R1 = the angle of rotation 
around the X axis, R2 is around the Y, and R3 is the Z.

Rotation around the Y axis
TEMPX = X * COS R2 - Z * SIN R2 
TEMPZ = X * SIN R2 + Z * COS R2 

Rotation around the X axis
Zf = TEMPZ * COS R1 - Y * SIN R1
TEMPY = TEMPZ * SIN R1 + Y * COS R1

Rotation around the Z axis
Xf = TEMPX * COS R3 + TEMPY * SIN R3
Yf = TEMPY * COS R3 - TEMPX * SIN R3

Now, if you understand this so far, then you'll have no problem programming your own 
3D program in QBASIC, a programming language that comes free with DOS.  If you are 
having trouble, don't worry; maybe reading this example will help you.  Otherwise, just 
copy this code and run it with QBASIC.

'WIRE3DUO.BAS by Matt Bross, 1997
DEFINT A-Z 'defines all unmarked variables as integers
TYPE PointType 'defines user defined type PointType
  x AS INTEGER 'X coordinate
  Y AS INTEGER 'Y coordinate
  z AS INTEGER 'Z coordinate
END TYPE 'end definition of type PointType
TYPE LineType 'defines user defined type LineType
  C1 AS INTEGER  'number of the first point of a line
  C2 AS INTEGER  'number of the second point of a line
  CLR AS INTEGER 'color of the line
END TYPE 'end definition of type LineType
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%INFO%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SCREEN 0, 0, 0, 0 'set to screen mode 0
WIDTH 80, 25'set height to 25 and width to 80 c
CLS 'clear the screen
COLOR 15 'print in color 15
PRINT "WIRE3DUO.BAS by Matt Bross, 1997" 'print the text in quotes
PRINT
COLOR 7
PRINT "3D WIREFRAME ANIMATION FOR THE QBASIC LANGUAGE"
PRINT "USE FREELY AS LONG AS I RECIEVE CREDIT DUE"
PRINT
COLOR 15
PRINT "--------CONTROLS--------"
COLOR 7
PRINT "   0 - reset rotation"
PRINT "   5 - stop rotation"
PRINT "   S - reset location"
PRINT "   A - stop translation"
PRINT "2, 8 - rotate around x axis"
PRINT "4, 6 - rotate around y axis"
PRINT "-, + - rotate around z axis"
PRINT CHR$(24); ", "; CHR$(25); " - translate vertically"
PRINT CHR$(27); ", "; CHR$(26); " - translate horizontally"
PRINT "Z, X - translate depthwise"
PRINT " Esc - exit"
COLOR 15
PRINT "----CASE INSENSITIVE----"
PRINT
PRINT "LOADING  0%"
'------------------------------>BEGIN INIT<-----------------------------
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%LOAD VARIABLES%%%%%%%%%%%%%%%%%%%%%%%%%%%
'make fast lookup variables for calling values (e.g. if x = 1 is used 
'often,x = 1, is slower than one = 1: x = one)
ZERO = 0: ONE = 1: THREESIXD = 360
intMAX = 32767: intMIN = -32768
bitsX10 = 1024
A$ = "A": S$ = "S": z$ = "Z": x$ = "X"
ZERO$ = "0": TWO$ = "2": FOUR$ = "4": FIVE$ = "5": SIX$ = "6"
EIGHT$ = "8"
PLUS$ = "+": MINUS$ = "-": EQUALS$ = "="
UP$ = CHR$(0) + "H": down$ = CHR$(0) + "P"
LEFT2$ = CHR$(0) + "K": RIGHT2$ = CHR$(0) + "M"
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%SCREEN VARIABLES%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = 300: SD = 60           'distance and perspective
MaxSpin = 20: MaxSpeed = 2 'max. translate and rotate
ScnMode = 7: m = 1: b = 0  'screen mode and pages
IF ScnMode <> 7 AND ScnMode <> 9 THEN ScnMode = 7
IF ScnMode = 7 THEN DX = 160: DY = 100 'center of screen mode 7
IF ScnMode = 9 THEN DX = 320: DY = 175 'center of screen mode 9
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%MAKE TABLES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DIM SINx(360) AS LONG 'initialize array for sines
DIM COSx(360) AS LONG 'initialize array for cosines
FOR i = 0 TO THREESIXD 'loop 360 times
'Create sine and cosine tables of all degrees 0 to 360 in radians and 
'scale by 1024 or 10 bits and store as a 4-byte integer.  The numbers 
'will be scaled down again in the main loop.  This is called fixed point 
'math and is faster that floating point, or math with a decimal.
  SINx(i) = SIN(i * (22 / 7) / 180) * bitsX10
  COSx(i) = COS(i * (22 / 7) / 180) * bitsX10
IF i MOD 40 = 0 THEN LOCATE 22, 9: PRINT STR$(i \ 4) + "%"
NEXT
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%LOAD POINTS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RESTORE PointData 'Set pointer to read from label PointData.
READ MaxPoints 'Reads MaxPoints from appended data.
DIM points(1 TO MaxPoints) AS PointType    'points at start
DIM ScnPnts(1 TO MaxPoints) AS PointType   'points after rotation
DIM SX(1 TO MaxPoints), SY(1 TO MaxPoints) 'points drawn to screen
FOR i = 1 TO MaxPoints: READ points(i).x, points(i).Y, points(i).z: NEXT
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%LOAD LINES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RESTORE LineData 'Set pointer to read from label LineData.
READ MaxLines 'Reads MaxLines from appended data.
DIM l(1 TO MaxLines) AS LineType 'line data
DIM OLX1(1 TO MaxLines), OLY1(1 TO MaxLines) 'old line data
DIM OLX2(1 TO MaxLines), OLY2(1 TO MaxLines)
FOR i = 1 TO MaxLines: READ l(i).C1, l(i).C2, l(i).CLR: NEXT
'-------------------------------->END INIT<-----------------------------
LOCATE 22, 9: PRINT "100%"
PRINT "Press a Key"
DO: LOOP UNTIL INKEY$ <> "" 'Do a loop until a key is pressed.
SCREEN ScnMode, 0, 0, 0 'set screen mode ScnMode
'----------------------------->BEGIN MAIN LOOP<-------------------------
DO
'*********************************GET KEY*******************************
K$ = UCASE$(INKEY$) 'get the upper case equivalent of input from the 
keyboard
SELECT CASE K$ 'select which case the string variable k$ is
  CASE ZERO$
    R1 = ZERO: R2 = ZERO: R3 = ZERO 'stop rotation and reset it
    D1 = ZERO: D2 = ZERO: D3 = ZERO
  CASE FIVE$
    D1 = ZERO: D2 = ZERO: D3 = ZERO 'stop rotation
  CASE A$
    MX = ZERO: MY = ZERO: MZ = ZERO 'stop translation
  CASE S$
    MX = ZERO: MY = ZERO: MZ = ZERO 'stop translation and reset it
    MMX = ZERO: MMY = ZERO: MMZ = ZERO
  CASE TWO$
    D1 = D1 - ONE 'rotate counter-clockwise around the x axis
  CASE EIGHT$
    D1 = D1 + ONE 'rotate clockwise around the x axis
  CASE FOUR$
    D2 = D2 - ONE 'rotate counter-clockwise around the y axis
  CASE SIX$
    D2 = D2 + ONE 'rotate clockwise around the y axis
  CASE PLUS$, EQUALS$
    D3 = D3 - ONE 'rotate clockwise around z axis
  CASE MINUS$
    D3 = D3 + ONE 'rotate counter-clockwise around the z axis
  CASE UP$
    MY = MY + ONE 'translate positively along the y axis
  CASE down$
    MY = MY - ONE 'translate negatively along the y axis
  CASE LEFT2$
    MX = MX + ONE 'translate positively along the x axis
  CASE RIGHT2$
    MX = MX - ONE 'translate negatively along the x axis
  CASE z$
    MZ = MZ + ONE 'translate positively along the z axis
  CASE x$
    MZ = MZ - ONE 'translate negatively along the z axis
  CASE CHR$(27)
    GOTO ShutDown 'end program
END SELECT
'*********************************ROTATION******************************
'keep sines and cosines in limits of the arrays
R1 = (R1 + D1) MOD THREESIXD
R2 = (R2 + D2) MOD THREESIXD
R3 = (R3 + D3) MOD THREESIXD
IF R1 < ZERO THEN R1 = THREESIXD + R1
IF R2 < ZERO THEN R2 = THREESIXD + R2
IF R3 < ZERO THEN R3 = THREESIXD + R3
'********************************TRANSLATION****************************
MMX = MMX + MX: MMY = MMY + MY: MMZ = MMZ + MZ
'keep variables within limits of integers
IF MMX > intMAX THEN MMX = intMAX
IF MMY > intMAX THEN MMX = intMAX
IF MMZ > intMAX THEN MMZ = intMAX
IF MMX < intMIN THEN MMX = intMIN
IF MMY < intMIN THEN MMX = intMIN
IF MMZ < intMIN THEN MMZ = intMIN
'******************************ROTATE POINTS****************************
S1& = SINx(R1): S2& = SINx(R2): S3& = SINx(R3)
C1& = COSx(R1): C2& = COSx(R2): C3& = COSx(R3)
FOR i = 1 TO MaxPoints
'Rotate points around the y axis.
  TEMPX = (points(i).x * C2& - points(i).z * S2&) \ bitsX10
  TEMPZ = (points(i).x * S2& + points(i).z * C2&) \ bitsX10
'Rotate points around the x axis.
  ScnPnts(i).z = (TEMPZ * C1& - points(i).Y * S1&) \ bitsX10
  TEMPY = (TEMPZ * S1& + points(i).Y * C1&) \ bitsX10
'Rotate points around the z axis.
  ScnPnts(i).x = (TEMPX * C3& + TEMPY * S3&) \ bitsX10
  ScnPnts(i).Y = (TEMPY * C3& - TEMPX * S3&) \ bitsX10
'*****************************CONVERT 3D TO 2D**************************
TEMPZ = ScnPnts(i).z + MMZ - SD
IF TEMPZ < ZERO THEN  'calculate points visible to the viewer
  SX(i) = (ScnPnts(i).x + MMX) * D \ TEMPZ + DX
  SY(i) = (ScnPnts(i).Y + MMY) * D \ TEMPZ + DY
END IF
NEXT
'******************************DRAW POLYGONS****************************
FOR i = 1 TO MaxLines
  coord1 = l(i).C1: coord2 = l(i).C2
  'erase old line
  LINE (OLX1(i), OLY1(i))-(OLX2(i), OLY2(i)), 0
  'get new points
  OLX1(i) = SX(coord1): OLY1(i) = SY(coord1)
  OLX2(i) = SX(coord2): OLY2(i) = SY(coord2)
  'Draw line from (SX1, SY1) to (SX2, SY2) in color CLR
  LINE (SX(coord1), SY(coord1))-(SX(coord2), SY(coord2)), l(i).CLR
NEXT
'******************************FRAME COUNTER****************************
F = F + 1
IF TIMER >= T# THEN 
	T# = TIMER + 1
	FPS = F	
	LOCATE 1, 1
	PRINT FPS: F = 0
END IF
'wait for vertical retrace, the electron beam to finish scanning the 
'monitor
WAIT &H3DA, 8
LOOP
'----------------------------->END MAIN LOOP<---------------------------
ShutDown:
SCREEN 0, 0, 0, 0: WIDTH 80, 25: CLS
PRINT "Final Location of Points"
PRINT " X", " Y", " Z": PRINT STRING$(31, "-")
FOR i = 1 TO MaxPoints: PRINT ScnPnts(i).x, ScnPnts(i).Y, ScnPnts(i).z: 
NEXT
PRINT : PRINT "Free space": PRINT " String Array   Stack": PRINT 
STRING$(21, "-")
PRINT FRE(""); FRE(-1); FRE(-2): END
PointData:
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%CUBE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'number of points
DATA 8
'Location of points (x, y, z)
DATA -10, 10,-10,  -10,-10,-10,  -10, 10, 10,  -10,-10, 10
DATA  10, 10, 10,   10,-10, 10,   10, 10,-10,   10,-10,-10
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%PYRAMID%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DATA 5
DATA   0, 10,  0,    0,  0,-10,   10,  0,  0,    0,  0, 10
DATA -10,  0,  0
LineData:
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%CUBE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'number of lines
DATA 12
'The points above can be numbered, the first data statement is 1.  The 
'points listed make lines from point 1 to point 2 in color 3.
DATA 1,2, 2,    1,3, 2,   3,4, 2,   2,4, 2,   1,7, 1,   3,5, 1,   5,7, 1
DATA 5,6, 1,    7,8, 1,   6,8, 1,   4,6, 1,   2,8, 1
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%PYRAMID%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DATA 8
DATA 1,2, 1,    1,3, 1,   1,4, 1,   1,5, 1,   2,3, 1,   3,4, 1,   4,5, 1
DATA 2,5, 1

Look over this until it makes sense before you proceed,  because compared to what's next, 
this was easy.  

Introduction to Shading

	The next thing to know is what a normal vector is.  First of a vector, in 3d space, 
is a line segment with a direction and a length, called a magnitude.  A "normal" vector is a 
vector perpendicular to the a plane.  For all practical purposes, always use triangles as 
planes because you don't have to worry about noncoplanar (not on the same plane) points, 
because any three points make a plane.

	To find the dot product, do the following to the three points of a triangle: subtract 
x, y, z from the one point to the other two points and cross multiply them and then set 
them to a variable to be stored in memory or if you want to set the equations
	
This finds the vectors
X coordinate of vector a = point 1 x - point 2 x
Y coordinate of vector a = point 1 y -  point 2 y
Z coordinate of vector a = point 1 z - point 2 z	
X coordinate of vector a = point 1 x - point 3 x
Y coordinate of vector a = point 1 y -  point 3y
Z coordinate of vector a = point 1 z - point 3 z

This finds their cross product
Normal vector x = (Y vector a)(Z vector b) - (Z vector b)(Y vector a)
Normal vector y = (Z vector a)(X vector b) - (X vector a)(Z vector b)
Normal vector z = (X vector a)(Y vector b) - (Y vector a)(X vector b) 

	Now, you might wonder: why do you need to know the normal vector?  You use 
the normal vector to detect if a plane is visible or not.  If a plane is not visible, then it is 
culled and does not have to be drawn.  Backface-culling, or searching for culled polygons, 
saves a lot of time when you use advanced shading systems that will be demonstrated 
later.  The normal vector also determines the shade of the plane if you use Lambert's Law 
of Shading.  Lambert's Law states  that the shade of the plane is related to the angle at 
which the light passes through the normal vector of a plane.  To find the angle of the light 
and the normal vector use this formula, with Vectors A and B, Ax is the X coordinate of 
the vector A, etc., and LA is the length of Vector A and LB is the length of B:
	
	COS(angle) = ((Ax)(Bx) + (Ay)(By) + 	(Az)(Bz)) /  (LA * LB)

If you forgot the 3D distance formula, it is d = SQR(x ^  2 + y ^ 2 + z ^ 2).  If the angle is 
greater than 0 (degrees) and less than 90 or greater than 270 and less than 360, then it is 
visible.   Then, for even less math, reduce the magnitude of both vectors to 1 to eliminate 
a multiply and a divide.  To do this, find the length of the normal  vector and divide the 
normal vector by its magnitude.  For more on vector mathematics see appendix C

After this, you must sort the planes by their Z coordinate, unless you use a Z buffer.  To 
do this, you must find the average of the Z coordinates.  the first Z coordinate plus the 
next Z coordinate plus the last Z coordinate, and then divide that quantity by 3.   Then use 
bubble sort to sort the coordinates.  Now, you should begin to understand the following 
QBASIC Example by Rich Geldreich, a 3D program with flat, has one, solid color per 
plane, shading.

'Shaded 3-D animation with shadows [solid5.bas] for QB4.5/PDS
'By Rich Geldreich 1992
'Notes...
'   This version uses some floating  point math in the initialization
'code for shading, but after initialization floating point math is not
'used at all.
'   The  shading  employs Lambert's Law to determine the intensity of
'each visible polygon.  Three simple  lookup tables are calculated at
'initialization time  which  are  used  to  eliminate  multiples  and
'divides in the main animation code.
'   The hidden face  detection  algorithm  was  made  by Dave Cooper.
'It's fast, and does not require any multiples and divides under most
'cases.  The "standard" way of detecting hidden faces, by finding the
'dot product of the normal of each polygon and  the  viewing  vector,
'was not just good (or fast) enough for me!
'   The PolyFill routine is the major  bottleneck  of  this  program.
'QB's  LINE  command isn't as fast as I would like it to be...  On my
'286-10, the speed isn't that bad (after all, this is all-QB!).  On a
'386 or 486, this thing should fly...  [hopefully]
'   The  shadows  are  calculated by projecting a line with the light
'source's vector through each of the points on the solid.  Where this
'line hits the ground  plane(which  has  a  constant Y coordinate) is
'where the new shadow point is plotted.
'   This program is 100% public domain- but  of  course  please  give
'some credit if you use anything from this program.  Thanks!
DEFINT A-Z
DECLARE SUB CullPolygons ()
DECLARE SUB DrawLine (xs%, ys%, xe%, ye%, EdgeList() AS ANY)
DECLARE SUB DrawObject ()
DECLARE SUB DrawShadows ()
DECLARE SUB EdgeFill (EdgeList() AS ANY, YLow%, YHigh%, C%)
DECLARE SUB FindNormals ()
DECLARE SUB PolyFill (x1%, y1%, x2%, y2%, x3%, y3%, C%)
DECLARE SUB RotatePoints ()
DECLARE SUB ShadePolygons ()
 
CONST True = -1, False = 0
 
TYPE EdgeType              'for fast polygon rasterization
    Low         AS INTEGER
    High        AS INTEGER
END TYPE
TYPE PointType
    XObject     AS INTEGER 'original coordinate
    YObject     AS INTEGER
    ZObject     AS INTEGER 'rotated coordinate
    XWorld      AS INTEGER
    YWorld      AS INTEGER
    ZWorld      AS INTEGER
    XView       AS INTEGER 'rotated & translated coordinate
    YView       AS INTEGER
    XShadow     AS INTEGER 'coordinates projected onto the ground plane
    YShadow     AS INTEGER
END TYPE
TYPE PolyType
    P1          AS INTEGER '3 points which make up the polygon (points
    P2          AS INTEGER 'to the point list array)
    P3          AS INTEGER
    Culled      AS INTEGER 'True if plane not visible
    ZCenter     AS INTEGER 'Z center of polygon
    ZOrder      AS INTEGER 'Used in the shell sort of the ZCenters
    Intensity   AS INTEGER 'Intensity of polygon
    WorldXN     AS INTEGER 'Contains the coordinates of the point
    WorldYN     AS INTEGER ' which is both perpendicular and a constant
    WorldZN     AS INTEGER ' distance from the polygon
    NormalX     AS INTEGER 'Normal of polygon -128 to 128
    NormalY     AS INTEGER ' (used for fast Lambert shading)
    NormalZ     AS INTEGER
END TYPE
TYPE LineType
    P1          AS INTEGER 'Used for shadow projection
    P2          AS INTEGER
END TYPE
 
DIM SHARED EdgeList(199) AS EdgeType
DIM SHARED SineTable(359 + 90) AS LONG 'cos(x)=sin(x+90)
DIM SHARED R1, R2, R3, ox, oy, oz
DIM SHARED MaxPoints, MaxPolys, MaxLines
 
DIM SHARED lines(100) AS LineType
DIM SHARED Polys(100) AS PolyType
DIM SHARED Points(100) AS PointType
DIM SHARED lx(256), ly(256), lz(256)  'lookup tables for Lambert shading

DIM SHARED s, XLow(1), XHigh(1), YLow(1), YHigh(1)
DIM SHARED ShadowXLow(1), ShadowXHigh(1), ShadowYLow(1), ShadowYHigh(1)
DIM SHARED lx, ly, lz
 
PRINT "QuickBASIC/PDS 3-D Solid Animation": PRINT "By Rich Geldreich 
1992"
PRINT : PRINT "Keys: [Turn NUMLOCK on]"
PRINT "Left.....................Angle 1 -"

PRINT "Right....................Angle 1 +"
PRINT "Up.......................Angle 2 -"
PRINT "Down.....................Angle 2 +"
PRINT "-........................Angle 3 -"
PRINT "+........................Angle 3 +"
PRINT "5........................Rotation Stop"
PRINT "0........................Rotation Reset"
PRINT "Up Arrow.................Forward"
PRINT "Down Arrow...............Backward"
PRINT "Left Arrow...............Left"
PRINT "Right Arrow..............Right"
PRINT : PRINT "Initializing..."
 
MaxPoints = 4  'Pyramid.
'Points follow...
DATA -100,0,100, -100,0,-100, 100,0,-100, 100,0,100, 0,-290,0
MaxPolys = 5
'Polygons follow (they must be specified in counterclockwise
'order for correct hidden face removal and shading)
DATA 4,0,3, 4,3,2, 4,1,0, 4,2,1, 3,0,1, 3,1,2
MaxLines = 7
'Lines follow for shadow computation...
DATA 4,0, 4,1, 4,2, 4,3, 0,1, 1,2, 2,3, 3,0
 
'MaxPoints = 7 'Cube.
'DATA -100,100,100
'DATA 100,100,100
'DATA 100,100,-100
'DATA -100,100,-100
'DATA -100,-100,100
'DATA 100,-100,100
'DATA 100,-100,-100
'DATA -100,-100,-100
'MaxPolys = 11
'DATA 5,4,0, 5,0,1
'DATA 6,2,3, 3,7,6
'DATA 6,5,1, 6,1,2
'DATA 7,0,4, 7,3,0
'DATA 6,7,4, 6,4,5
'DATA 0,3,2, 1,0,2
'MaxLines = 11
'DATA 0,1, 1,2, 2,3, 3,0
'DATA 4,5, 5,6, 6,7, 7,4
'DATA 4,0, 5,1, 6,2, 7,3
 
'MaxPoints = 15 'Weird pencil-like shape...
'DATA 0,0,0, 250,0,0, 400,40,0, 400,60,0, 250,100,0, 0,100,0,-20,90,0
'DATA  -20,10,0
'DATA 0,0,-50, 250,0,-50, 400,40,-50, 400,60,-50, 250,100,-50, 0,10
'DATA -50, -20,90,-50, -20,10,-50
'MaxPolys = 27
'DATA 1,0,7, 1,7,2, 2,7,6, 2,6,3, 3,6,4, 4,6,5
'DATA 9,15,8, 9,10,15, 10,14,15, 10,11,14, 11,13,14, 11,12,13
'DATA 8,7,0, 8,15,7, 8,0,1, 9,8,1, 9,1,2, 10,9,2, 10,2,3, 11,10,3
'DATA 12,11,4, 11,3,4, 4,5,13, 4,13,12
'DATA 5,6,14, 5,14,13, 14,6,7, 14,7,15
'MaxLines = 23
'DATA 0,1, 1,2, 2,3, 3,4, 4,5, 5,6, 6,7, 7,0
'DATA 8,9, 9,10, 10,11, 11,12, 12,13, 13,14, 14,15, 15,0
'DATA 0,8, 1,9, 2,10, 3,11, 4,12, 5,13, 6,14, 7,15
 
FOR a = 0 TO MaxPoints
    READ Points(a).XObject, Points(a).YObject, Points(a).ZObject
    X = X + Points(a).XObject: Y = Y + Points(a).YObject: Z = Z + 
Points(a).ZObject
NEXT
'Center the object
X = X \ (MaxPoints + 1): Y = Y \ (MaxPoints + 1): Z = Z \ (MaxPoints +1)
FOR a = 0 TO MaxPoints
    Points(a).XObject = Points(a).XObject - X
    Points(a).YObject = Points(a).YObject - Y
    Points(a).ZObject = Points(a).ZObject - Z
NEXT
FOR a = 0 TO MaxPolys
    READ Polys(a).P1, Polys(a).P2, Polys(a).P3
NEXT
FOR a = 0 TO MaxLines
    READ lines(a).P1, lines(a).P2
NEXT
 
'Precalculate the normal point of each polygon for fast Lambert shading

FindNormals
 
'Precalculate the sine table
a = 0

FOR a! = 0 TO (359 + 90) / 57.29 STEP 1 / 57.29
    SineTable(a) = SIN(a!) * 1024: a = a + 1
NEXT
 
'Some light source configurations won't work that great!
l1 = 70: l2 = 40           'light source's spherical coordinates
a1! = l1 / 57.29: a2! = l2 / 57.29

s1! = SIN(a1!): c1! = COS(a1!)
s2! = SIN(a2!): c2! = COS(a2!)
lx = 128 * s1! * c2!        'convert spherical coordinates to a vector
ly = 128 * s1! * s2!        'scale up by 128 for integer math
lz = 128 * c1!
 
FOR a = -128 TO 128         'precalculate the three light source tables
    lx(a + 128) = lx * a    'for fast Lambert shading
    ly(a + 128) = ly * a
    lz(a + 128) = lz * a
NEXT
 
PRINT "Strike a key...": DO: LOOP WHILE INKEY$ = ""
 
R1 = 0: R2 = 0: R3 = 0      'three angles of rotation
ox = 0: oy = -50: oz = 1100 'object's origin (this program cannot 	
			    'currently handle the object when it goes 		
		    'behind the viewer!)
s = 1: t = 0
 
SCREEN 7, , 0, 0
OUT &H3C8, 0                'set 16 shades
FOR a = 0 TO 15
    OUT &H3C9, (a * 34) \ 10
    OUT &H3C9, (a * 212) \ 100
    OUT &H3C9, (a * 4) \ 10
    IF a = 7 THEN OUT &H3C7, 16: OUT &H3C8, 16
NEXT
LINE (0, 100)-(319, 199), 9, BF
LINE (0, 0)-(319, 99), 3, BF
SCREEN 7, , 1, 0
LINE (0, 100)-(319, 199), 9, BF
LINE (0, 0)-(319, 99), 3, BF
 
YHigh(0) = -32768: ShadowYHigh(0) = -32768
YHigh(1) = -32768: ShadowYHigh(1) = -32768
DO
    'Flip active and work pages so user doesn't see our messy drawing
    SCREEN 7, , s, t: SWAP s, t
 
    'Wait for vertical retrace to reduce flicker
    WAIT &H3DA, 8
 
    'Erase the old image from the screen
    IF YHigh(s) <> -32768 THEN
        IF YHigh(s) < 100 THEN
            LINE (XLow(s), YLow(s))-(XHigh(s), YHigh(s)), 3, BF
        ELSEIF YLow(s) < 100 THEN
            LINE (XLow(s), YLow(s))-(XHigh(s), 99), 3, BF
            LINE (XLow(s), 100)-(XHigh(s), YHigh(s)), 9, BF
        ELSE
            LINE (XLow(s), YLow(s))-(XHigh(s), YHigh(s)), 9, BF
        END IF
    END IF
    IF ShadowYHigh(s) <> -32768 THEN
        LINE (ShadowXLow(s), ShadowYLow(s))-(ShadowXHigh(s), 
ShadowYHigh(s)), 9, BF
    END IF
    RotatePoints
    CullPolygons
    ShadePolygons
 
    XLow(s) = 32767: XHigh(s) = -32768
    YLow(s) = 32767: YHigh(s) = -32768
    DrawShadows
    DrawObject
 
    R1 = (R1 + D1) MOD 360: IF R1 < 0 THEN R1 = R1 + 360
    R2 = (R2 + D2) MOD 360: IF R2 < 0 THEN R2 = R2 + 360
    R3 = (R3 + D3) MOD 360: IF R3 < 0 THEN R3 = R3 + 360
    oz = oz + dz: ox = ox + dx
    IF oz < 600 THEN
        oz = 600: dz = 0
    ELSEIF oz > 8000 THEN
        oz = 8000: dz = 0
    END IF
    IF ox < -4000 THEN
        ox = -4000: dx = 0
    ELSEIF ox > 4000 THEN
        ox = 4000: dx = 0
    END IF
    a$ = INKEY$
    SELECT CASE a$
    CASE "4"
        D1 = D1 - 2
    CASE "6"
        D1 = D1 + 2
    CASE "8"
        D2 = D2 - 2
    CASE "2"
        D2 = D2 + 2
    CASE "5"
        D1 = 0: D2 = 0: D3 = 0
    CASE "0"
        R1 = 0: R2 = 0: R3 = 0
        D1 = 0: D2 = 0: D3 = 0
    CASE "+"
        D3 = D3 + 2
    CASE "-"
        D3 = D3 - 2
    CASE CHR$(27)
        END
    CASE CHR$(0) + CHR$(72)
        dz = dz - 20
    CASE CHR$(0) + CHR$(80)
        dz = dz + 20
    CASE CHR$(0) + CHR$(77)
        dx = dx - 20
    CASE CHR$(0) + CHR$(75)
        dx = dx + 20
    END SELECT
LOOP
 
'Shades the polygons using Lambert shading. Notice the total lack of
'floating point math- only 1 divide is required for each polygon shaded.
'(This divide can be eliminated, but the shading won't be as accurate.)

'"Culls" the polygons which aren't visible to the viewer. Also shades
'each polygon using Lambert's law.
SUB CullPolygons
    'This algorithm for removing hidden faces was developed by Dave 
'Cooper.
    'There is another method, by finding the dot product of the
    'plane's normal and the viewing vector, but this algorithm is
    'much faster because of its simplicity(and lack of floating point
    'calculations).
    FOR a = 0 TO MaxPolys
        P1 = Polys(a).P1
        P2 = Polys(a).P2
        P3 = Polys(a).P3
 
        IF Points(P1).YView <= Points(P2).YView THEN
            IF Points(P3).YView < Points(P1).YView THEN
                PTop = P3
                PNext = P1
                PLast = P2
            ELSE
                PTop = P1
                PNext = P2
                PLast = P3
            END IF
        ELSE
            IF Points(P3).YView < Points(P2).YView THEN
                PTop = P3
                PNext = P1
                PLast = P2
            ELSE
                PTop = P2
                PNext = P3
                PLast = P1
            END IF
        END IF
 
        XLow = Points(PTop).XView
        YLow = Points(PTop).YView
 
        XNext = Points(PNext).XView
        XLast = Points(PLast).XView
 
        IF XNext <= XLow AND XLast >= XLow THEN
            Polys(a).Culled = True
        ELSEIF XNext >= XLow AND XLast <= XLow THEN
            Polys(a).Culled = False
        ELSE
            YNext = Points(PNext).YView
            YLast = Points(PLast).YView
IF((YNext-YLow)*256&)\(XNext-XLow)<((YLast-YLow)*256&)\(XLast-XLow) THEN
                Polys(a).Culled = False
            ELSE
                Polys(a).Culled = True
            END IF
        END IF
 
    NEXT
END SUB

'Enters a line into the edge list. For each scan line, the line's
'X coordinate is found. Notice the lack of floating point math in this
'subroutine.
SUB DrawLine (xs, ys, xe, ye, EdgeList() AS EdgeType)
 
    IF ys > ye THEN SWAP xs, xe: SWAP ys, ye
 
    IF ye < 0 OR ys > 199 THEN EXIT SUB

    IF ys < 0 THEN
        xs = xs + ((xe - xs) * -ys) \ (ye - ys)
        ys = 0

    END IF
 
    xd = xe - xs
    yd = ye - ys
 
    IF yd <> 0 THEN xi = xd \ yd: xrs = ABS(xd MOD yd)
 
    xr = -yd \ 2
 
    IF ye > 199 THEN ye = 199
 
    xdirect = SGN(xd) + xi
 
    FOR Y = ys TO ye
        IF xs < EdgeList(Y).Low THEN EdgeList(Y).Low = xs
        IF xs > EdgeList(Y).High THEN EdgeList(Y).High = xs
 
        xr = xr + xrs
        IF xr > 0 THEN
            xr = xr - yd
            xs = xs + xdirect
        ELSE
            xs = xs + xi
        END IF
    NEXT
 
END SUB

SUB DrawObject
 
    'Find the center of each visible polygon, and prepare the order 
list.
    NumPolys = 0
    FOR a = 0 TO MaxPolys
        IF Polys(a).Culled = False THEN 'is this polygon visible?
            Polys(NumPolys).ZOrder = a
            NumPolys = NumPolys + 1
            Polys(a).ZCenter = Points(Polys(a).P1).ZWorld + 
Points(Polys(a).P2).ZWorld + Points(Polys(a).P3).ZWorld
        END IF
    NEXT
    'Sort the visible polygons by their Z center using a shell sort.
    NumPolys = NumPolys - 1
    Mid = (NumPolys + 1) \ 2
    DO
        FOR a = 0 TO NumPolys - Mid
            CompareLow = a
            CompareHigh = a + Mid
		DO WHILE Polys(Polys(CompareLow).ZOrder).ZCenter < 
Polys(Polys(CompareHigh).ZOrder).ZCenter
                SWAP Polys(CompareLow).ZOrder, Polys(CompareHigh).ZOrder
                CompareHigh = CompareLow
                CompareLow = CompareLow - Mid
                IF CompareLow < 0 THEN EXIT DO
            LOOP
        NEXT
        Mid = Mid \ 2
    LOOP WHILE Mid > 0
    'Plot the visible polygons.
    FOR Z = 0 TO NumPolys
        a = Polys(Z).ZOrder 'which polygon do we plot?
        P1 = Polys(a).P1: P2 = Polys(a).P2: P3 = Polys(a).P3
        PolyFill (Points(P1).XView), (Points(P1).YView), 
(Points(P2).XView), (Points(P2).YView), (Points(P3).XView), 
(Points(P3).YView), (Polys(a).Intensity)
    NEXT
END SUB

SUB DrawShadows
    YLow = 32767: YHigh = -32768
    XLow = 32767: XHigh = -32768
    FOR a = 0 TO MaxPoints
        'Project the 3-D point onto the ground plane...
        temp& = (Points(a).YWorld - 200)
        X = Points(a).XWorld - (temp& * lx) \ ly
        Y = 200 'ground plane has a constant Y coordinate
        Z = Points(a).ZWorld - (temp& * lz) \ ly
        'Put the point into perspective
        xTemp = 160 + (X * 400&) \ Z
        yTemp = 100 + (Y * 300&) \ Z
 
        Points(a).XShadow = xTemp
        Points(a).YShadow = yTemp
 
        'Find the lowest & highest X Y coordinates
        IF yTemp < YLow THEN YLow = yTemp

        IF yTemp > YHigh THEN YHigh = yTemp
        IF xTemp < XLow THEN XLow = xTemp
        IF xTemp > XHigh THEN XHigh = xTemp
    NEXT
 
    'Store lowest & highest coordinates for later erasing...
    ShadowXLow(s) = XLow: ShadowYLow(s) = YLow

    ShadowXHigh(s) = XHigh: ShadowYHigh(s) = YHigh
    IF XHigh < 0 OR XLow > 319 OR YLow > 199 OR YHigh < 0 THEN EXIT SUB
    IF YHigh > 199 THEN YHigh = 199
    IF YLow < 0 THEN YLow = 0
 
    'Initialize the edge list
    FOR a = YLow TO YHigh
        EdgeList(a).Low = 32767
        EdgeList(a).High = -32768
    NEXT
 
    'Enter the lines into the edge list
    FOR a = 0 TO MaxLines
        P1 = lines(a).P1
        P2 = lines(a).P2
        DrawLine (Points(P1).XShadow), (Points(P1).YShadow), 
(Points(P2).XShadow), (Points(P2).YShadow), EdgeList()
        'LINE ((Points(P1).XShadow),(Points(P1).YShadow))-
((Points(P2).XShadow), (Points(P2).YShadow)), 0
    NEXT
 
    'Fill the polygon
    EdgeFill EdgeList(), YLow, YHigh, 3
 
END SUB

SUB EdgeFill (EdgeList() AS EdgeType, YLow, YHigh, C)
    FOR a = YLow TO YHigh
        LINE (EdgeList(a).Low, a)-(EdgeList(a).High, a), C
    NEXT
END SUB

'This routine initializes the data required by the fast Lambert shading
'algorithm. It calculates the point which is both perpendicular
'and a constant distance from each polygon and stores it. This point
'is then rotated with the rest of the points. When it comes time for
'shading, the normal to the polygon is looked up in a simple lookup
'table for maximum speed.
SUB FindNormals
    FOR a = 0 TO MaxPolys
        P1 = Polys(a).P1: P2 = Polys(a).P2: P3 = Polys(a).P3
 
        'find the vectors of 2 lines inside the polygon
        ax! = Points(P2).XObject - Points(P1).XObject
        ay! = Points(P2).YObject - Points(P1).YObject
        az! = Points(P2).ZObject - Points(P1).ZObject
 
        bx! = Points(P3).XObject - Points(P2).XObject
        by! = Points(P3).YObject - Points(P2).YObject
        bz! = Points(P3).ZObject - Points(P2).ZObject
 
        'find the cross product of the 2 vectors
        nx! = ay! * bz! - az! * by!
        ny! = az! * bx! - ax! * bz!
        nz! = ax! * by! - ay! * bx!
 
        'normalize the vector so it ranges from -1 to 1
        M! = SQR(nx! * nx! + ny! * ny! + nz! * nz!)
        IF M! <> 0 THEN nx! = nx! / M!: ny! = ny! / M!: nz! = nz! / M!
        'store the vector for later rotation(notice that it is scaled
        'up by 128 so it can be stored as an integer variable)
        Polys(a).WorldXN = nx! * 128 + Points(P1).XObject
        Polys(a).WorldYN = ny! * 128 + Points(P1).YObject
        Polys(a).WorldZN = nz! * 128 + Points(P1).ZObject
    NEXT
END SUB

'Draws a polygon to the screen. Simply finds the start and stop X
'coordinates for each scan line within the polygon and uses the
'LINE command for filling.
SUB PolyFill (x1, y1, x2, y2, x3, y3, C) 'for QB 4.5 guys
 
    'find lowest and high X & Y coordinates
    IF y1 < y2 THEN YLow = y1 ELSE YLow = y2
    IF y3 < YLow THEN YLow = y3
    IF y1 > y2 THEN YHigh = y1 ELSE YHigh = y2
    IF y3 > YHigh THEN YHigh = y3
 
    IF x1 < x2 THEN XLow = x1 ELSE XLow = x2
    IF x3 < XLow THEN XLow = x3
    IF x1 > x2 THEN XHigh = x1 ELSE XHigh = x2
    IF x3 > XHigh THEN XHigh = x3
 
    IF YLow < 0 THEN YLow = 0
    IF YHigh > 199 THEN YHigh = 199
 
    IF XLow < XLow(s) THEN XLow(s) = XLow
    IF XHigh > XHigh(s) THEN XHigh(s) = XHigh
 
    IF YLow < YLow(s) THEN YLow(s) = YLow
    IF YHigh > YHigh(s) THEN YHigh(s) = YHigh
 

    'check for polygons which cannot be visible
    IF YHigh < 0 OR YLow > 199 OR XLow > 319 OR XHigh < 0 THEN EXIT SUB
 
    'initialize the edge list
    FOR a = YLow TO YHigh
        EdgeList(a).Low = 32767
        EdgeList(a).High = -32768
    NEXT
 
    'Remember the lowest & highest X and Y coordinates drawn to the
    'screen for later erasing
 
    'Find the start and stop X coordinates for each scan line
    DrawLine (x1), (y1), (x2), (y2), EdgeList()
    DrawLine (x2), (y2), (x3), (y3), EdgeList()
    DrawLine (x3), (y3), (x1), (y1), EdgeList()
    EdgeFill EdgeList(), YLow, YHigh, C
 
END SUB

'Rotates the points of the object and the object's normals.
'Avoids floating point math for speed.
SUB RotatePoints
 
    'lookup the sine and cosine of each angle...
    s1& = SineTable(R1): c1& = SineTable(R1 + 90)
    s2& = SineTable(R2): c2& = SineTable(R2 + 90)
    s3& = SineTable(R3): c3& = SineTable(R3 + 90)
 
    'rotate the points of the object
    FOR a = 0 TO MaxPoints
        xo = Points(a).XObject
        yo = Points(a).YObject
        zo = Points(a).ZObject
        GOSUB Rotate3D
 
        Points(a).XView = 160 + (x2 * 400&) \ z3
        Points(a).YView = 100 + (y3 * 300&) \ z3
        'IF y3 > 300 THEN STOP
 
        Points(a).XWorld = x2
        Points(a).YWorld = y3
        Points(a).ZWorld = z3
    NEXT
    'rotate the normals of each polygon...
    FOR a = 0 TO MaxPolys
        xo = Polys(a).WorldXN
        yo = Polys(a).WorldYN
        zo = Polys(a).WorldZN
        GOSUB Rotate3D
        P1 = Polys(a).P1
        'unorigin the point
        x2 = x2 - Points(P1).XWorld
        y3 = y3 - Points(P1).YWorld
        z3 = z3 - Points(P1).ZWorld
        'check the bounds just in case of a round off error
        IF x2 < -128 THEN x2 = -128 ELSE IF x2 > 128 THEN x2 = 128
        IF y3 < -128 THEN y3 = -128 ELSE IF y3 > 128 THEN y3 = 128
        IF z3 < -128 THEN z3 = -128 ELSE IF z3 > 128 THEN z3 = 128
        'store the normal back; it's now ready for the shading
        'calculations (which are simplistic now)
        Polys(a).NormalX = x2 + 128
        Polys(a).NormalY = y3 + 128
        Polys(a).NormalZ = z3 + 128
    NEXT
    EXIT SUB
 
Rotate3D:
    x1 = (xo * c1& - zo * s1&) \ 1024 'yaw
    z1 = (xo * s1& + zo * c1&) \ 1024
 
    z3 = (z1 * c3& - yo * s3&) \ 1024 + oz 'pitch
    y2 = (z1 * s3& + yo * c3&) \ 1024
 
    x2 = (x1 * c2& + y2 * s2&) \ 1024 + ox 'roll
    y3 = (y2 * c2& - x1 * s2&) \ 1024 + oy
 
RETURN
END SUB

SUB ShadePolygons
    FOR a = 0 TO MaxPolys
        IF Polys(a).Culled = False THEN

         'lookup the polygon's normal for shading
         '(128*128)\15 = 1092
         Intensity = (lx(Polys(a).NormalX) + ly(Polys(a).NormalY) + 
lz(Polys(a).NormalZ)) \ 1092
         IF Intensity < 0 THEN Intensity = 0
         Intensity = Intensity + 5
         IF Intensity > 15 THEN Intensity = 15

         Polys(a).Intensity = Intensity
        END IF
    NEXT
END SUB

Gouraud and Phong Shading 

Gouraud and phong are two different methods of shading.  Gouraud shading is linear and 
less realistic than phong, but it is faster.  Phong shading is based on a cosinual wave.  
Linear interpolation is a big long word for something simple once you understand it.  
Hopefully, I will succeed in describing it to you.  Linear interpolation takes two values and 
finds the slope between them (e.g. the linear interpolation of  points (0, 0) and (10, 5) is 
the y delta (a Greek letter that means change, in this case 0 - 5) over the X delta (0 - 10) 
the increment between the two points is -5 / - 10 or 1 / 2, for every 1 x over the y goes up 
1 / 2.  This is almost like slope intercept form that you learned in 7th or 8th grade (Like 
me!).  This slope intercept will be modified in Algebra II to y - Y1 = m(x - X1) which is 
much easier to work with.)

The use Gouraud shading you must find the slope of the points and the slope of the colors.  
To find the color, do the same as you did for the shade of a polygon in the last QBASIC 
example, but this time find the intensity of the point instead of the normal.  Then linearly 
interpolate the color and the y coordinates and plot your pixel, then increase the point and 
the color by the slope and continue until your polygon is filled.  Here is an ASCII 
representation of how gouraud shading works.  The numbers represent palette values in 
hexadecimal , to save space.

		10                           10                    10
	          / |                           d/ |e                  d/e|e
                     /  |                          a/  |c                 a/ab|c
                    /   |                         7/   |a                7/89|c
                   /    |                        4/    |8               4/56 |8
                  /     |                       1/     |6              1/245|6
              0 /___|4                    0/___|4            0/1234|4
                                                1223                 1223

Phong shading finds the normal to each pixel and the dot product of the normal and the 
light vector, and then plots the point.  This is very slow because you must use one square 
root per pixel to find its distance from the origin but it is 100% accurate.  

Setting up your palette

There are two different methods for setting up your palette in a program.  Linearly, just 
shading a color by a constant increment, or the phong illumination model that is:

	color =  ambient + (diffuse)(COS a) + (specular)(COS a) ^ n

Ambient is the color of the low lights.  Diffuse controls the actual color.  Specular 
determines the highlight.  N hold special properties for the texture of the surface, and a is 
an angle between 0 and 90, the angles light can hit a plane.  

Polygon Clipping

This basically clips off the ends of the polygons that aren't on the screen.  If you try to plot 
a point off the screen, your computer might crash.  This doesn't happen often.  Usually 
either nothing happens or the program crashes.  The worst case scenario is it will 
permanently damage your monitor.  But with clipping this will never happen.  All you have 
to use to avoid this disaster are a few conditional statements.  If the x coordinate of the 
point to be plotted is greater than the screen width or less than zero, don't plot the point.  
If the y coordinate of the point to be plotted is greater than the screen height or less than 
zero, then don't plot the point.  Pretty simple. 

Here is my example showing all of these techniques.  The GIF routine is by Rich 
Geldreich;  the Sorting algorithm is based on one by Ryan Wellman, and the Gouraud 
shading is based on one by Luke Molnar.

Painter's Algorithm

	An Algorithm is a problem that is solved by a computer.  All of this stuff has been 
algorithms.  The Painter's Algorithm states that the viewer can only see the polygons 
closest to him/her.  So if you sort the polygons by their average Z coordinate ((z1+ z2 + 
z3 + ... + zn) / n) back to front, then the picture will look right, and then draw it in that 
order.  When you look at a person, you don't usually see their internal organs because 
their skin is in the way;  in exactly the same way that you don't see a polygon that is 
behind another polygon.  The main drawback to this is overlapping polygons with equal or 
close Z coordinates cannot be drawn correctly.  That is why you need...

Hidden Face Removal

	One of the best ways to optimize, make faster, 3D programs is hidden face 
removal.  A lot of time is spent rendering - or drawing - the shape after it has been 
calculated.  This is very wasteful because it takes a lot of time to write individual pixels to 
the screen especially when polygons aren't visible to the viewer because they are covered 
up by other polygons.

Backface-Culling

	Backface-Culling is a way to determine if a plane is visible of not.  To find this you 
use the generally you use the dot product of the normal vector of the plane and vector of 
light, but there are other ways of doing this (list the points in counterclockwise order and 
then see if they are still in that order).  This only works for enclosed shapes though, and 
does not work even then all the time, because when you look at a plane in 3D space you 
see two sides of it, the front and the back.  This determines if the front or back side of the 
polygon is showing or not, the dot product is negative.  The dot product is the following, 
where (X1, Y1, Z1);  (X2, Y2, Z2); (X3, Y3, Z3) are the first three points of a polygon:

	(X3((Z1Y2) - (Y1Z2)) + Y3((X1Z2)-(Z1X2)) + ((Z3(Y1X2)-(X1Y2))

Z-Buffering and S-Buffering

	Z-buffering and S-buffering are extremely advanced systems that replace Z sorts 
and backface culling and more accurate.  

The Z-buffer uses a copy of the screen to remember if a close pixel has been placed. The 
major drawback to this technique is that it use an incredible amount of array space.  The 
must have an element for ever pixel on the screen times the number of bytes you use to 
store each element.  One byte has 256 different combinations, so your world can 256 
points deep.  This isn't very much.  Two bytes give you 65536 (256) different 
combinations but for just screen mode 13h (hexadecimal) this is 128k (320 by 200 by 2 
bytes) which is a huge array.  The way Z buffering works is it finds where to plot a pixel 
and it interpolates the Z coordinates, if the pixel is the closest at that point to the view 
then it plots the point and changes the element for the pixel point to the new Z value else it 
moves onto the next point.  This technique can draw intersecting planes with ease and 
saves time with slow sorting algorithms.  

An S-buffer is basically a compressed Z buffer, so it uses less space, and therefore is 
better.  S- buffers finds the z of the parts of the plane its looking at and adds it to a linked 
list of the planes and sorts them in memory before drawing them to the screen so no time 
is wasted with graphics.  

The main disadvantages of s-buffering compared to Z-buffering are that curved surfaces 
aren't rendered well, and that s-buffering is very disorganized and harder to manage.

BSP Trees

	A BSP (of Binary Space Partitioning) tree is a very advanced and versatile 
technique.  A  BSP tree can be used to remove hidden surfaces and raytrace along with 
many other things.  What a BSP tree does is it stores all space and if any part of space has 
a "hyperplane" in it, the space is divided in half.  A hyperplane in 3D space is a plane (in 
2D a line).  If the plane is split there are two ways to go (this is where binary in the name 
comes from) and this branches off like a tree.  To first build a tree, select any plane in 
space and set polygons on it.  If there is a plane in any part of space then this creates a 
new area to make a tree on,  repeat this.  The most common use of this technique is 
hidden face removal.  To use this for hidden face removal, just follow the tree backwards 

DEFINT A-Z
DECLARE FUNCTION GetByte% ()
DECLARE SUB BufferWrite (a%)
DECLARE SUB MakeGif (a$, ScreenX%, ScreenY%, Xstart%, YStart%, Xend%, 
Yend%, NumColors%, AdaptorType%)
DECLARE SUB PutByte (a%)
DECLARE SUB PutCode (a%)
DECLARE SUB pal (c%, R%, G%, B%)
CONST TRUE = -1, FALSE = NOT TRUE

'GS3DO.BAS by Matt Bross, 1997
'The sorting algorithm was originally written by Ryan Wellman, which I
'modified for my own purposes.  I made the 3D program with help from '3D 
tutorials by Lithium /VLA, Shade3D.BAS by Rich Geldreich; and 'Gouraud 
fill with Luke Molnar's (of M/K Productions) gorau.bas. The GIF 'support 
is from Rich Geldreich's MakeGif.BAS.
'
'completely RANDOMIZE
RANDOMIZE TIMER: DO: RANDOMIZE TIMER: LOOP UNTIL RND > .5
'ON ERROR GOTO ErrorHandler
TYPE PointType
  x AS SINGLE     'X coordinate
  y AS SINGLE     'Y coordinate
  z AS SINGLE     'Z coordinate
  shade AS INTEGER  'shade of points
  dis AS SINGLE     'distance from the origin (0, 0, 0)
END TYPE
TYPE PolyType
  C1 AS INTEGER     'number of the first point of a polygon
  C2 AS INTEGER     'number of the second point of a polygon
  C3 AS INTEGER     'number of the third point of a polygon
  culled AS INTEGER 'TRUE if the polygon isn't visible
  AvgZ AS INTEGER   'used to sort Z coordinates of polygons
END TYPE
TYPE FillType
  Y1 AS INTEGER     'starting Y coordinate
  Y2 AS INTEGER     'ending Y coordinate
  clr1 AS INTEGER   'starting color
  clr2 AS INTEGER   'ending color
END TYPE
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%INFO%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SCREEN 0, 0, 0, 0: WIDTH 80, 25: CLS
PRINT "GS3DO.BAS by Matt Bross, 1997"
PRINT
PRINT "3D ANIMATION FOR THE QBASIC LANGUAGE"
PRINT "COPYRIGHT MATT BROSS.  USE FREELY AS"
PRINT "LONG AS CREDIT IS GIVEN."
PRINT
PRINT "--------CONTROLS--------"
PRINT "   0 - reset rotation"
PRINT "   5 - stop rotation"
PRINT "   S - reset location"
PRINT "   A - stop translation"
PRINT "2, 8 - rotate around x axis"
PRINT "4, 6 - rotate around y axis"
PRINT "-, + - rotate around z axis"
PRINT CHR$(24); ", "; CHR$(25); " - translate vertically"
PRINT CHR$(27); ", "; CHR$(26); " - translate horizontally"
PRINT "Z, X - translate depthwise"
PRINT " Esc - exit"
PRINT "----CASE INSENSITIVE----"
PRINT
INPUT "OBJECT TO LOAD", file$
IF file$ = "" THEN file$ = "pyramid.txt"
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%VARIABLES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'SRX = the screen's x resolution
'SRY = the screen's y resolution
SRX = 320: SRY = 200
'DX = the X coordinate of the center of the screen
'DY = the Y coordinate of the center of the screen
DX = SRX \ 2: DY = SRY \ 2
'D = the viewer's distance then object: SD = controls perspective
D = 350: SD = 140
'MaxSpin = controls the maximum rotation speed
'MaxSpeed = controls the maximum translation speed
MaxSpin = 20: MaxSpeed = 10
'NumClr = the number of palette values to assign to shading each color
NumClr = 63
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%GIF STUFF%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DIM SHARED OutBuffer$, OStartAddress, OAddress, OEndAddress, Oseg
DIM SHARED CodeSize, CurrentBit, Char&, BlockLength
DIM SHARED Shift(7) AS LONG
DIM SHARED x, y, Minx, MinY, MaxX, MaxY, Done, GIFFile, LastLoc&
ShiftTable:
DATA 1, 2, 4, 8, 16, 32, 64, 128
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%SIN TABLES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'create SINe and COSine tables for 360 degrees in radians, and then 
'scale 1024 times for faster math.
'$STATIC
DIM SINx(360) AS LONG, COSx(360) AS LONG
FOR i = 0 TO 360
  SINx(i) = SIN(i * (22 / 7) / 180) * 1024 'use 1024 to shift binary 
digits
  COSx(i) = COS(i * (22 / 7) / 180) * 1024 'over 6 bits.
NEXT i
'%%%%%%%%%%%%%%%%%%%%%%%%%%GOURAUD SHADE ARRAYS%%%%%%%%%%%%%%%%%%%%%%%%%
DIM scan(320) AS FillType 'DIM gouraud shading array
DIM coord(1 TO 3)
'%%%%%%%%%%%%%%%%%%%%%%%%DOUBLE BUFFERING ARRAYS%%%%%%%%%%%%%%%%%%%%%%%%
DIM SHARED aofs&
DIM SHARED ScnBuf(32001) 'DIM array to serve as page in SCREEN 13
ScnBuf(0) = 320 * 8 'set length (x)
ScnBuf(1) = 200     'set height (y)
DEF SEG = VARSEG(ScnBuf(2)) 'get segment of beginning of array data
aofs& = VARPTR(ScnBuf(2))   'get offset of beginning of array data
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%LIGHT TABLES%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DIM LX(256), LY(256), LZ(256)
'Location of light source in spherical coordinates
l1 = 70: l2 = 40: a1! = l1 / 57.29: a2! = l2 / 57.29
s1! = SIN(a1!): C1! = COS(a1!): s2! = SIN(a2!): C2! = COS(a2!)
LX = 128 * s1! * C2!: LY = 128 * s1! * s2!: LZ = 128 * C1!
'find length of segment from light source to origin (0, 0, 0)
ldis! = SQR(LX * LX + LY * LY + LZ * LZ) / 128
FOR a = -128 TO 128
  LX(a + 128) = LX * a 'make light source lookup tables for shading
  LY(a + 128) = LY * a
  LZ(a + 128) = LZ * a
NEXT a
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%LOAD OBJECT%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
OPEN file$ FOR INPUT AS #1
'Load Points Data
INPUT #1, MaxPoints, MaxPolys
DIM POINTS(MaxPoints) AS PointType     'at start
DIM ScnPnts(MaxPoints) AS PointType    'after rotation
DIM SX(MaxPoints), SY(MaxPoints)       'points drawn to screen
FOR i = 1 TO MaxPoints
INPUT #1, x!, y!, z!: POINTS(i).x = x!: POINTS(i).y = y!: POINTS(i).z=z!
  'find distance from point to the origin (0, 0, 0)
  dis! = SQR(x! * x! + y! * y! + z! * z!)
  POINTS(i).dis = dis! * ldis!: ScnPnts(i).dis = dis! * ldis!
NEXT i
'Load Polygon Data
DIM SHARED P(MaxPolys) AS PolyType 'stores all polygon data
FOR i = 1 TO MaxPolys
INPUT #1, P(i).C1, P(i).C2, P(i).C3
NEXT i: CLOSE
PRINT "Press a Key"
DO: LOOP UNTIL INKEY$ <> ""
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%SET PALETTE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SCREEN 13: CLS
s! = 0: ClrStep! = 63 / NumClr
FOR a = 0 TO NumClr
  pal a, c, c, c
  s! = s! + ClrStep!: c = INT(s!)
NEXT a
'%%%%%%%%%%%%%%%%%%%%%%%%%%%LOOK UP VARIABLES%%%%%%%%%%%%%%%%%%%%%%%%%%%
ZERO = 0: ONE = 1: THREE6D = 360
'------------------------------>BEGIN MAIN LOOP<------------------------
DO
'*********************************GET KEY*******************************
k$ = UCASE$(INKEY$)
SELECT CASE k$
  CASE "0"
    R1 = ZERO: R2 = ZERO: R3 = ZERO
    D1 = ZERO: D2 = ZERO: D3 = ZERO
  CASE "5"
    D1 = ZERO: D2 = ZERO: D3 = ZERO
  CASE "A"
    MX = ZERO: MY = ZERO: MZ = ZERO
  CASE "S"
    MX = ZERO: MY = ZERO: MZ = ZERO
    MMX = ZERO: MMY = ZERO: MMZ = ZERO
  CASE "2"
    D1 = D1 - ONE
  CASE "8"
    D1 = D1 + ONE
  CASE "4"
    D2 = D2 - ONE
  CASE "6"
    D2 = D2 + ONE
  CASE "+", "="
    D3 = D3 - ONE
  CASE "-"
    D3 = D3 + ONE
  CASE CHR$(0) + "H"
    MY = MY + ONE
  CASE CHR$(0) + "P"
    MY = MY - ONE
  CASE CHR$(0) + "K"
    MX = MX + ONE
  CASE CHR$(0) + "M"
    MX = MX - ONE
  CASE "Z"
    MZ = MZ + ONE
  CASE "X"
    MZ = MZ - ONE
  CASE CHR$(27)
    GOTO ShutDown
  CASE "G"
    INPUT "SAVE SCREEN as .GIF"; s$
    IF LEFT$(UCASE$(s$), 1) = "Y" THEN
      INPUT "Input Filename:", file$
      IF RIGHT$(UCASE$(file$), 4) <> ".GIF" THEN file$ = file$ + ".GIF"
      PUT (0, 0), ScnBuf, PSET
      MakeGif file$, 320, 200, 0, 0, 319, 199, 256, 2
    END IF
END SELECT
'*********************************ROTATION******************************
'keep rotation speed under control
IF D1 > MaxSpin THEN D1 = MaxSpin
IF D2 > MaxSpin THEN D2 = MaxSpin
IF D2 > MaxSpin THEN D2 = MaxSpin
IF D1 < -MaxSpin THEN D1 = -MaxSpin
IF D2 < -MaxSpin THEN D2 = -MaxSpin
IF D2 < -MaxSpin THEN D2 = -MaxSpin
'keep SINes and COSines in array limits
R1 = (R1 + D1) MOD THREE6D: IF R1 < ZERO THEN R1 = THREE6D + R1
R2 = (R2 + D2) MOD THREE6D: IF R2 < ZERO THEN R2 = THREE6D + R2
R3 = (R3 + D3) MOD THREE6D: IF R3 < ZERO THEN R3 = THREE6D + R3
'********************************TRANSLATION****************************
'Keep translation speed from becoming uncontrollable
IF MX > MaxSpeed THEN MX = MaxSpeed
IF MY > MaxSpeed THEN MY = MaxSpeed
IF MZ > MaxSpeed THEN MZ = MaxSpeed
IF MX < -MaxSpeed THEN MX = -MaxSpeed
IF MY < -MaxSpeed THEN MY = -MaxSpeed
IF MZ < -MaxSpeed THEN MZ = -MaxSpeed
MMX = MMX + MX: MMY = MMY + MY: MMZ = MMZ + MZ
'Keeps variables within limits of integers
IF MMX > 32767 THEN MMX = 32767
IF MMY > 250 THEN MMY = 250
IF MMZ > 120 THEN MMZ = 120
IF MMX < -32767 THEN MMX = -32767
IF MMY < -120 THEN MMY = -120
IF MMZ < -327 THEN MMZ = -327
'*******************************MOVE OBJECT*****************************
FOR i = 1 TO MaxPoints
'rotate points around the Y axis
  TEMPX = (POINTS(i).x * COSx(R2) - POINTS(i).z * SINx(R2)) \ 1024
  TEMPZ = (POINTS(i).x * SINx(R2) + POINTS(i).z * COSx(R2)) \ 1024
'rotate points around the X axis
  ScnPnts(i).z = (TEMPZ * COSx(R1) - POINTS(i).y * SINx(R1)) \ 1024
  TEMPY = (TEMPZ * SINx(R1) + POINTS(i).y * COSx(R1)) \ 1024
'rotate points around the Z axis
  ScnPnts(i).x = (TEMPX * COSx(R3) + TEMPY * SINx(R3)) \ 1024
  ScnPnts(i).y = (TEMPY * COSx(R3) - TEMPX * SINx(R3)) \ 1024
'******************************CONVERT 3D TO 2D*************************
TEMPZ = ScnPnts(i).z + MMZ - SD
IF TEMPZ < ZERO THEN  'only calculate points visible by viewer
  SX(i) = (D * ((ScnPnts(i).x + MMX) / TEMPZ)) + DX
  SY(i) = (D * ((ScnPnts(i).y + MMY) / TEMPZ)) + DY
END IF
'*******************************SHADE POINTS****************************
X1 = ScnPnts(i).x: Y1 = ScnPnts(i).y: Z1 = ScnPnts(i).z
s = CINT((X1 * LX + Y1 * LY + Z1 * LZ) \ ScnPnts(i).dis) + 128
  IF s < ZERO THEN s = ZERO
  IF s > 256 THEN s = 256
shade = (LX(s) + LY(s) + LZ(s)) \ 3
  IF shade < ZERO THEN shade = ZERO
  IF shade > NumClr THEN shade = NumClr
ScnPnts(i).shade = shade
NEXT
FOR i = 1 TO MaxPolys
'*************************CULL NON-VISIABLE POLYGONS********************
'this isn't perfect yet so I REMmed it, so for more speed unrem the 
following
coord(1) = P(i).C1: coord(2) = P(i).C2: coord(3) = P(i).C3
X1 = ScnPnts(coord(1)).x: X2 = ScnPnts(coord(2)).x: X3 = 
ScnPnts(coord(3)).x
Y1 = ScnPnts(coord(1)).y: Y2 = ScnPnts(coord(2)).y: Y3 = 
ScnPnts(coord(3)).y
Z1 = ScnPnts(coord(1)).z: Z2 = ScnPnts(coord(2)).z: Z3 = 
ScnPnts(coord(3)).z
cull1 = X3 * ((Y1 * Z2) - (Z1 * Y2)): cull2 = Y3 * ((X1 * Z2) - (Z1 * 
X2))
cull3 = Z3 * ((X1 * Y2) - (Y1 * X2))
IF cull1 + cull2 + cull3 = ZERO THEN P(i).culled = TRUE ELSE P(i).culled 
= FALSE
'******************FIND AVERAGE Z COORDINATE OF EACH POLYGON************
P(i).AvgZ = (Z1 + Z2 + Z3) \ 3
NEXT i
'******************SORT POLGONS BY THEIR AVERAGE Z COORDINATE***********
increment = MaxPolys + 1
DO
increment = increment \ 2
FOR index = 1 TO MaxPolys - increment
  IF P(index).AvgZ > P(index + increment).AvgZ THEN
  SWAP P(index), P(index + increment)
    IF index > increment THEN
    cutpoint = index
    DO
    index = (index - increment): IF index < 1 THEN index = 1
      IF P(index).AvgZ > P(index + increment).AvgZ THEN
        SWAP P(index), P(index + increment)
      ELSE
        index = cutpoint: EXIT DO
      END IF
    LOOP
    END IF
  END IF
NEXT index
LOOP UNTIL increment <= 1
'******************************DRAW POLYGONS****************************
'clear screen buffer.  Use a 320 by 200 BLOADable graphic for a 
background.
ERASE ScnBuf: ScnBuf(0) = 2560: ScnBuf(1) = SRY

FOR i = 1 TO MaxPolys
  IF P(i).culled = FALSE THEN
    'load points
    coord(1) = P(i).C1: coord(2) = P(i).C2: coord(3) = P(i).C3
    'find highest and lowest Y coordinates
    xmin = SRX: xmax = ZERO
    IF SX(coord(1)) > xmax THEN xmax = SX(coord(1))
    IF SX(coord(2)) > xmax THEN xmax = SX(coord(2))
    IF SX(coord(3)) > xmax THEN xmax = SX(coord(3))
    IF SX(coord(1)) < xmin THEN xmin = SX(coord(1))
    IF SX(coord(2)) < xmin THEN xmin = SX(coord(2))
    IF SX(coord(3)) < xmin THEN xmin = SX(coord(3))
    'keep min's and max's in the limits of the screen
    IF xmin < ZERO THEN xmin = ZERO
    IF xmax > SRX THEN xmax = SRX
    IF xmin > SRX THEN EXIT FOR
    IF xmax < ZERO THEN EXIT FOR
    IF SY(coord(1)) AND SY(coord(2)) AND SY(coord(3)) < ZERO THEN EXIT 
FOR
    IF SY(coord(1)) AND SY(coord(2)) AND SY(coord(3)) > SRY THEN EXIT 
FOR

    ERASE scan

    FOR j = 1 TO 3
      k = j + 1: IF k > 3 THEN k = 1
      VAL1 = coord(j): VAL2 = coord(k)
      IF SX(VAL1) > SX(VAL2) THEN SWAP VAL1, VAL2
      Y1 = SY(VAL1): X1 = SX(VAL1): Y2 = SY(VAL2): X2 = SX(VAL2)
      col1 = ScnPnts(VAL1).shade: Col2 = ScnPnts(VAL2).shade
      XDelta = X2 - X1: YDelta = Y2 - Y1: CDelta = Col2 - col1
      IF XDelta <> ZERO THEN
        YSlope = (YDelta / XDelta) * 128
        CSlope = (CDelta / XDelta) * 128
      ELSE
        YSlope = ZERO
        CSlope = ZERO
      END IF

      YVal& = Y1 * 128: CVal& = col1 * 128
      IF X1 < ZERO THEN X1 = ZERO
      IF X2 > SRX THEN X2 = SRX

      FOR f = X1 TO X2
        IF scan(f).Y1 = ZERO THEN
          scan(f).Y1 = YVal& \ 128
          scan(f).clr1 = CVal& \ 128
        ELSE
          scan(f).Y2 = YVal& \ 128
          scan(f).clr2 = CVal& \ 128
        END IF
        YVal& = YVal& + YSlope
        CVal& = CVal& + CSlope
      NEXT f
  NEXT j

  FOR f = xmin TO xmax
 
    IF scan(f).Y1 > scan(f).Y2 THEN
      Y1 = scan(f).Y2: Y2 = scan(f).Y1
      col1 = scan(f).clr2: Col2 = scan(f).clr1
    ELSE
      Y1 = scan(f).Y1: Y2 = scan(f).Y2
      col1 = scan(f).clr1: Col2 = scan(f).clr2
    END IF

    YDelta = Y2 - Y1: CDelta = Col2 - col1
    IF YDelta = ZERO THEN YDelta = 1
    CSlope = (CDelta / YDelta) * 128: CVal& = col1 * 128

    FOR j = scan(f).Y1 TO scan(f).Y2
    'clip polygon to screen (set boundaries)
    IF f < SRX AND f > ZERO AND j > ZERO AND j < SRY THEN
      pixel = CVal& \ 128: IF pixel > NumClr THEN pixel = NumClr
      'write pixel to screen buffer
      POKE aofs& + f + j * 320&, pixel
    END IF
      CVal& = CVal& + CSlope
    NEXT j
  NEXT f
END IF
NEXT i

PUT (ZERO, ZERO), ScnBuf, PSET 'dump array to screen, like PCOPY
'******************************FRAME COUNTER****************************
'LOCATE 1, 1: PRINT fps: frame = frame + 1
'LOCATE 2, 1: PRINT TIMER - D#: D# = TIMER
'IF TIMER > t# THEN t# = TIMER + 1: fps = frame: frame = zero
LOOP
'------------------------------>END MAIN LOOP<--------------------------
ShutDown:
DEF SEG
SCREEN 0, 0, 0, 0: WIDTH 80, 25: CLS
PRINT "GS3DO.BAS by Matt Bross, 1997"
PRINT : PRINT "THERE WERE"; MaxPoints; "POINTS AND"; MaxPolys; 
"POLYGONS"
PRINT : PRINT "Free space"
PRINT " String Array   Stack"
PRINT STRING$(21, "-")
PRINT FRE(""); FRE(-1); FRE(-2): END

ErrorHandler:
RESUME NEXT

'Puts a byte into the disk buffer... when the disk buffer is full it is
'dumped to disk.
SUB BufferWrite (a) STATIC

    IF OAddress = OEndAddress THEN  'are we at the end of the buffer?
        PUT GIFFile, , OutBuffer$   ' yup, write it out and
        OAddress = OStartAddress    '  start all over
    END IF
    POKE OAddress, a                'put byte in buffer
    OAddress = OAddress + 1         'increment position
END SUB

'This routine gets one pixel from the display.
FUNCTION GetByte STATIC

    GetByte = POINT(x, y)           'get the "byte"
    x = x + 1                       'increment X coordinate
    IF x > MaxX THEN                'are we too far?
        LINE (Minx, y)-(MaxX, y), 0 'a pacifier for impatient users
        x = Minx                    'go back to start
        y = y + 1                   'increment Y coordinate
        IF y > MaxY THEN            'are we too far down?
            Done = TRUE             ' yup, flag it then
        END IF
    END IF
END FUNCTION

'
'-----------------------------------------------------------------------
'   PDS 7.1 & QB4.5 GIF Compression Routine v1.00 By Rich Geldreich 1992
'-----------------------------------------------------------------------
'
'A$          = output filename
'ScreenX     = X resolution of screen(320, 640, etc.)
'ScreenY     = Y resolution of screen(200, 350, 480, etc.)
'XStart      = <-upper left hand corner of area to encode
'YStart      = < "                                      "
'Xend        = <-lower right hand corner of area to encode
'Yend        = < "                                       "
'NumColors   = # of colors on screen(2, 16, 256)
'AdaptorType = 1 for EGA 2 for VGA
'NOTE: EGA palettes are not supported in this version of MakeGIF.
'
SUB MakeGif (a$, ScreenX, ScreenY, Xstart, YStart, Xend, Yend, 
NumColors, AdaptorType)
    'hash table's size - must be a prime number!
    CONST Table.Size = 7177

    DIM Prefix(Table.Size - 1), Suffix(Table.Size -1),code(Table.Size-1)
   
    'The shift table contains the powers of 2 needed by the
    'PutCode routine. This is done for speed. (much faster to
    'look up an integer than to perform calculations...)
    RESTORE ShiftTable
    FOR a = 0 TO 7: READ Shift(a): NEXT
   
    'MinX, MinY, MaxX, MaxY have the encoding window
    Minx = Xstart: MinY = YStart
    MaxX = Xend: MaxY = Yend
   
    'Open GIF output file
    GIFFile = FREEFILE                  'use next free file
    OPEN a$ FOR BINARY AS GIFFile
   
    'Put GIF87a header at beginning of file
    B$ = "GIF87a"
    PUT GIFFile, , B$
   
    'See how many colors are in this image...
    SELECT CASE NumColors
        CASE 2              'monochrome image
            BitsPixel = 1   '1 bit per pixel
            StartSize = 3   'first LZW code is 3 bits
            StartCode = 4   'first free code
            StartMax = 8    'maximum code in 3 bits

        CASE 16             '16 colors images
            BitsPixel = 4   '4 bits per pixel
            StartSize = 5   'first LZW code is 5 bits
            StartCode = 16  'first free code
            StartMax = 32   'maximum code in 5 bits
       
        CASE 256            '256 color images
            BitsPixel = 8   '8 bits per pixel
            StartSize = 9   'first LZW code is 9 bits
            StartCode = 256 'first free code
            StartMax = 512  'maximum code in 9 bits
    END SELECT
    'This following routine probably isn't needed- I've never
    'had to use the "ColorBits" variable... With the EGA, you
    'have 2 bits for Red, Green, & Blue. With VGA, you have 6 bits.
    SELECT CASE AdaptorType
        CASE 1
            ColorBits = 2               'EGA
        CASE 2
            ColorBits = 6               'VGA
    END SELECT
    
    PUT GIFFile, , ScreenX  'put screen's dimensions
    PUT GIFFile, , ScreenY
    'pack colorbits and bits per pixel
    a = 128 + (ColorBits - 1) * 16 + (BitsPixel - 1)
    PUT GIFFile, , a
    'throw a zero into the GIF file
    a$ = CHR$(0)
    PUT GIFFile, , a$
    'Get the RGB palette from the screen and put it into the file...
    SELECT CASE AdaptorType
    CASE 1
        STOP
        'EGA palette routine not implemented yet
    CASE 2
        OUT &H3C7, 0
        FOR a = 0 TO NumColors - 1
          'Note: a BIOS call could be used here, but then we have to use
            'the messy CALL INTERRUPT subs...
            R = (INP(&H3C9) * 65280) \ 16128 'C=R * 4.0476190(for 0-255)
            G = (INP(&H3C9) * 65280) \ 16128
            B = (INP(&H3C9) * 65280) \ 16128
            a$ = CHR$(R): PUT GIFFile, , a$
            a$ = CHR$(G): PUT GIFFile, , a$
            a$ = CHR$(B): PUT GIFFile, , a$
        NEXT
    END SELECT
   
    'write out an image descriptor...
    a$ = ","                        '"," is image seperator
    PUT GIFFile, , a$               'write it
    PUT GIFFile, , Minx             'write out the image's location
    PUT GIFFile, , MinY
    ImageWidth = (MaxX - Minx + 1)  'find length & width of image
    ImageHeight = (MaxY - MinY + 1)
    PUT GIFFile, , ImageWidth       'store them into the file
    PUT GIFFile, , ImageHeight
    a$ = CHR$(BitsPixel - 1)        '# bits per pixel in the image
    PUT GIFFile, , a$
   
    a$ = CHR$(StartSize - 1)        'store the LZW minimum code size
    PUT GIFFile, , a$
   
    'Initialize the vars needed by PutCode
    CurrentBit = 0: Char& = 0
    
    MaxCode = StartMax          'the current maximum code size
    CodeSize = StartSize        'the current code size
    ClearCode = StartCode       'ClearCode & EOF code are the
    EOFCode = StartCode + 1     '  first two entries
    StartCode = StartCode + 2   'first free code that can be used
    NextCode = StartCode        'the current code
   
  OutBuffer$ = STRING$(5000, 32)  'output buffer; for speedy disk writes
    a& = SADD(OutBuffer$)           'find address of buffer
    a& = a& - 65536 * (a& < 0)
    Oseg = VARSEG(OutBuffer$) + (a& \ 16)   'get segment + offset >> 4
    OAddress = a& AND 15                  'get address into segment
    OEndAddress = OAddress + 5000         'end of disk buffer
    OStartAddress = OAddress              'current location in disk 
buffer
    DEF SEG = Oseg

    GOSUB ClearTree                       'clear the tree & output a
    PutCode ClearCode                     '  clear code
   
    x = Xstart: y = YStart          'X & Y have the current pixel
    Prefix = GetByte                'the first pixel is a special case
    Done = FALSE                    'True when image is complete
    
    DO 'while there are more pixels to encode

        DO 'until we have a new string to put into the table

           IF Done THEN 'write out the last pixel, clear the disk buffer
                      'and fix up the last block so its count is correct
               
                PutCode Prefix  'write last pixel
                PutCode EOFCode 'send EOF code
               
                IF CurrentBit <> 0 THEN
                    PutCode 0       'flush out the last code...
                END IF
                PutByte 0

                OutBuffer$ = LEFT$(OutBuffer$, OAddress - OStartAddress)
                PUT GIFFile, , OutBuffer$
           a$ = ";" + STRING$(8, &H1A) 'the 8 EOF chars is not standard,
                                       'but many GIF's have them, so how
                                       'much could it hurt?
                PUT GIFFile, , a$
               
            a$ = CHR$(255 - BlockLength) 'correct the last block's count
                PUT GIFFile, LastLoc&, a$
               
                CLOSE GIFFile
                EXIT SUB
            ELSE 'get a pixel from the screen and see if we can find
                 'the new string in the table
                Suffix = GetByte
                GOSUB Hash        'is it there?
             IF Found = TRUE THEN Prefix = code(index) 'yup, replace the
                                    'prefix:suffix string with whatever
                                    'code represents it in the table
            END IF
        LOOP WHILE Found  'don't stop unless we find a new string
        
        PutCode Prefix    'output the prefix to the file
       
        Prefix(index) = Prefix  'put the new string in the table
        Suffix(index) = Suffix
  code(index) = NextCode  'we've got to keep track if what code this is!
       
   Prefix = Suffix         'Prefix=the last pixel pulled from the screen
       
        NextCode = NextCode + 1   'get ready for the next code
        IF NextCode = MaxCode + 1 THEN  'can an output code ever exceed
                                        'the current code size?
            'yup, increase the code size

            MaxCode = MaxCode * 2
            
        'Note: The GIF89a spec mentions something about a deferred clear
        'code. When the clear code is deferred, codes are not entered
        'into the hash table anymore. When the compression of the image
        'starts to fall below a certain threshold, the clear code is
        'sent and the hash table is cleared. The overall result is
        'greater compression, because the table is cleared less often.
        'This version of MakeGIF doesn't support this, because some GIF
        'decoders crash when they attempt to enter too many codes
        'into the string table.

            IF CodeSize = 12 THEN       'is the code size too big?
                PutCode ClearCode       'yup; clear the table and
                GOSUB ClearTree         'start over
                NextCode = StartCode
                CodeSize = StartSize
                MaxCode = StartMax


            ELSE
                CodeSize = CodeSize + 1 'just increase the code size if
            END IF                      'it's not too high( not > 12)
        END IF
        
    LOOP 'while we have more pixels
ClearTree:
    FOR a = 0 TO Table.Size - 1       'clears the hashing table
        Prefix(a) = -1                '-1 = invalid entry
        Suffix(a) = -1
        code(a) = -1
    NEXT
RETURN
'this is only one of a plethora of ways to search the table for
'a match! I used a binary tree first, but I switched to hashing
'cause it's quicker(perhaps the way I implemented the tree wasn't
'optimal... who knows!)
Hash:
    'hash the prefix & suffix(there are also many ways to do this...)
     '?? is there a better formula?
    index = ((Prefix * 256&) XOR Suffix) MOD Table.Size
    '
    '(Note: the table size(7177 in this case) must be a prime number, or
    'else there's a chance that the routine will hang up... hate when
    'that happens!)
    '
    'Calculate an offset just in case we don't find what we want on the
    'first try...
    IF index = 0 THEN    'can't have Table.Size-0 !
        Offset = 1
    ELSE
        Offset = Table.Size - index
    END IF
    
 DO 'until we (1) find an empty entry or (2) find what we're lookin for
                                  
       
        IF code(index) = -1 THEN  'is this entry blank?
            Found = FALSE         'yup- we didn't find the string
            RETURN
        'is this entry the one we're looking for?
        ELSEIF Prefix(index) = Prefix AND Suffix(index) = Suffix THEN
            'yup, congrats you now understand hashing!!!
    
            Found = TRUE
            RETURN
        ELSE
            'shoot! we didn't find anything interesting, so we must
            'retry- this is what slows hashing down. I could of used
            'a bigger table, that would of speeded things up a little
            'because this retrying would not happen as often...
            index = index - Offset
            IF index < 0 THEN   'too far down the table?
                'wrap back the index to the end of the table
                index = index + Table.Size
            END IF
        END IF
    LOOP
END SUB

SUB pal (c, R, G, B)
OUT &H3C8, c
OUT &H3C9, R
OUT &H3C9, G
OUT &H3C9, B
END SUB

'Puts a byte into the GIF file & also takes care of each block.
SUB PutByte (a) STATIC
    BlockLength = BlockLength - 1   'are we at the end of a block?
    IF BlockLength <= 0 THEN        ' yup,
        BlockLength = 255           'block length is now 255
   LastLoc& = LOC(1) + 1 + (OAddress - OStartAddress) 'remember the pos.
    BufferWrite 255                                    'for later fixing
    END IF
    BufferWrite a   'put a byte into the buffer
END SUB

'Puts an LZW variable-bit code into the output file...
SUB PutCode (a) STATIC
    Char& = Char& + a * Shift(CurrentBit) 'put the char were it belongs;
 CurrentBit = CurrentBit + CodeSize  ' shifting it to its proper place
DO WHILE CurrentBit > 7              'do we have a least one full byte?
   PutByte Char& AND 255             ' yup! mask it off and write it out
   Char& = Char& \ 256               'shift the bit buffer right 8 bits
   CurrentBit = CurrentBit - 8       'now we have 8 less bits
    LOOP 'until we don't have a full byte
END SUB
'--------------------------------->CUT HERE PYRAMID.TXT<----------------------------------
5 36
10 0 0
0 0 10
0 0 -10
-10 0 0
0 10 0
5 1 3
5 3 1
3 1 5
3 5 1
1 3 5
1 5 3
5 1 2
5 2 1
2 1 5
2 5 1
1 2 5
1 5 2
5 2 4
5 4 2
2 4 5
2 5 4
4 5 2
4 2 5
5 3 4
5 4 3
3 4 5 
3 5 4
4 5 3
4 3 5
4 3 1
4 1 3
3 1 4
3 4 1
1 3 4
1 4 3
4 1 2
4 2 1
2 1 4
2 4 1
1 2 4
1 4 2
'----------------------------------------->END SNIP<--------------------------------------------- 

Texture Mapping

Texture mapping is adding a picture and curves it to a 3D surface.  No mapping algorithm 
is 100% correct but this adds many incredible effects to 3D programs and for advanced 
3D raytracers used for art and what not.  The easiest way to do this is to use a 4 point 
plane, or a 2D quadrilateral.  This technique is used in many DOOM-like games (Although 
many DOOM-like games, including DOOM, are 2D raytracers) for texturing walls and 
floors.  All you have to do is scan convert all the sides and fill it in.  Scan conversion 
divides each side up into equal parts, e.g. divide the distance of every side of the polygon 
by any same number, and then fill the points in each section with the pixel color of the 
texture.  But this only textures a quadrilateral, so for other shades than squares or 
rectangles you must create virtual points to use as the points of the quadrilateral and this 
screws up the perspective.  This is called affine texture mapping.    
ex.
'<HTML>
' <HEAD>
'  <TITLE>
'   Texturemapping
'  </TITLE>
' </HEAD>
' <BODY BGCOLOR=FFFFFF>
'  <B>
'   Texturemapping
'  </B>
'  <BR>
'  Here is some source code I found on a BBS by my house and I guess it 
'came from rec.games.programmer way back in 1994!  (Geeze thats old!)  
'This code looks like it will work fine with the 3D code that I have 
'on-line as long as you define the faces with 4 points.<BR>
'  This is not my code.<BR>
'  <OL>
'   <FONT COLOR=00AC00>
'    <CODE>
'     <PRE>
'From rec.games.programmer Fri. Mar 18 13:17:50 1994
'Path:govonca!torn!howland.reston.ans.net!wupost!gumby!yale!yale.edu!nig
'el.msen.com!ilium!rcsuna.gmr.com!mloeff01.elec.mid.gmeds.com!sun85.elec
'.mid.gmeds.com!not-for-mail
'From: holz@cpchq.cpc.gmeds.com (White Flame)
'Newsgroups: rec.games.programmer
'Subject: Fast texture mapping algo!
'Date: 15 Mar 1994 09:31:10 -0500
'Organization: North American Operations, G.M. Corp.
'Lines: 354
'Sender: holz@sun85.elec.mid.gmeds.com
'Distribution: world
'Message-ID: <2m4gre$bak@sun85.elec.mid.gmeds.com>
'Reply-To: holz@cpchq.cpc.gmeds.com
'NNTP -Posting - Host: sun85
'
'
'Okay, here's that texture-mapping algorithm that
'Alan Parker posted in comp.graphics.algorithms.
'it's pretty fast, but I think that I might have
'copied a line wrong or something, as you can see
'by the top and bottom faces of the square.  They're
'a pixel off.  Oh, well, here it is:
'
'--------------8<-----------------------
DEFSNG A-Z
' Texture Mapping example.
' Maps 4 sided picture onto a 4 sided polygon.

' Written in GFA basic version 2

' By Alan Parker (E-Mail: selap@dmu.ac.uk)
' On 1/2/94

' Ported to Microsoft QBASIC v1.0 by White Flame on 2/9/94

' This texture mapping method was 'borrowed' from Ben of Chaos.
' I've no idea who originally thought of this method.

' There's not many detailed comments, I've already dome my best when I
' explained the algorithm before and nobody understood it, now it's up 
to you!

' Note: All variables with % after them are Integers, all other are 
'Reals.
'   Y co-ords go from top(0) to bottom(199) of the screen.
'   'picture' refers to the original bitmapped picture.
'   'polygon' refers to the polygon drawn on screen.
'   The picture to be mapped must be in the top, left (0,0) corner of 
the
'   screen.
'   This is not a perspective map, but you would only have to modify the
'   scan converter to generate the picture x,y points depending on the
'   z co-ords of the polygon. The texture mapping routine would stay the
'   same.
'   I've half managed to do a perspective map, but it's still bugged!
'   It will be posted when (if?) I manage to get it working...

DECLARE SUB GetPolygonPoints ()
DECLARE SUB FindSmallLargeY ()
DECLARE SUB ScanConvert (X1%, Y1%, X2%, Y2%, Pside$)
DECLARE SUB TextureMap ()
DECLARE SUB ScanLeftSide (X1%, X2%, Ytop%, Lineheight%, Pside$)
DECLARE SUB ScanRightSide (X1%, X2%, Ytop%, Lineheight%, Pside$)

' *********************************************
' Machine specific code
' Find screen address
Screenbase% = 0

' Set screen address and set to low resolution

SCREEN 13
DEF SEG = &HA000
'Set Palette
OUT &H3C8, 0
FOR x = 0 TO 63
  OUT &H3C9, x
  OUT &H3C9, 0
  OUT &H3C9, 0
NEXT
FOR x = 0 TO 63
  OUT &H3C9, 63 - x
  OUT &H3C9, x
  OUT &H3C9, 0
NEXT
FOR x = 0 TO 63
  OUT &H3C9, 0
  OUT &H3C9, 63 - x
  OUT &H3C9, x
NEXT
FOR x = 0 TO 63
  OUT &H3C9, x
  OUT &H3C9, x
  OUT &H3C9, 63
NEXT

'Load in the picture to be texture mapped (must be in top,left of 
'screen)

0
CLS
FOR Y% = 0 TO 15
  FOR x% = 0 TO 15
    PSET (x%, Y%), x% + 16 * Y%
  NEXT
NEXT

' *********************************************
' General Texture Mapping code.
' Machine independent, probably :-)

' Variables and arrays.

DIM SHARED Lefttable(400, 2), Righttable(400, 2)  'Scan converter tables 
for polygon
DIM SHARED Miny%, Maxy%

' p.s - replace both the 400 values above with you max. screen height in 
pixels

DIM SHARED Polypoints%(3, 1) ' Array for polygon co-ords, 4 pairs(x,y) 
co-ords
DIM SHARED Pwidth%, Pheight%

Pwidth% = 15   'original picture width in pixels
Pheight% = 15  'original picture height in pixels

Miny% = 32767  'set initial smallest y co-ord of polygon after rotation
Maxy% = 0      'set initial largest y co-ord of polygon after rotation

' Main program

GetPolygonPoints


FindSmallLargeY

FOR x% = 0 TO 3
  PSET (Polypoints%(x%, 0), Polypoints%(x%, 1)), 63
NEXT

' Send polygon points to the scan converter
X1% = Polypoints%(0, 0)
Y1% = Polypoints%(0, 1)
X2% = Polypoints%(1, 0)
Y2% = Polypoints%(1, 1)
ScanConvert X1%, Y1%, X2%, Y2%, "top"     'scan top of picture
X1% = Polypoints%(1, 0)
Y1% = Polypoints%(1, 1)
X2% = Polypoints%(2, 0)
Y2% = Polypoints%(2, 1)
ScanConvert X1%, Y1%, X2%, Y2%, "right"   'scan right of picture
X1% = Polypoints%(2, 0)
Y1% = Polypoints%(2, 1)
X2% = Polypoints%(3, 0)
Y2% = Polypoints%(3, 1)
ScanConvert X1%, Y1%, X2%, Y2%, "bottom"  'scan bottom of picture
X1% = Polypoints%(3, 0)
Y1% = Polypoints%(3, 1)
X2% = Polypoints%(0, 0)
Y2% = Polypoints%(0, 1)
ScanConvert X1%, Y1%, X2%, Y2%, "left"    'scan left of picture

' Do the actual texture mapping
TextureMap

' The points for polygon.
' These points in the form x,y define the shape of a four sided polygon
' rotated 45 degrees. These are 2d points that would normally come from 
' a 3d routine after perspective had been applied.
' The points must be defined clockwise....

Polygonpoints:
DATA 100,0
DATA 190,10
DATA 180,100
DATA 90,90


' an alternative polygon to try

'Polygonpoints:
'DATA 50,50
'DATA 150,50
'DATA 150,150
'DATA 50,150



DO: LOOP WHILE INKEY$ = ""
GOTO 0

'--------------8<-----------------------
'
'Stupid "^M"'s are showing up at the end of every line here on xvnews
'under OpenWindows.  It always does that on DOS text file, I hope
'it doesn't screw anything up...



'White Flame
'</PRE>
'    </CODE>
'    </FONT>
'   </OL>
'   <BR>
'   <IMG SRC="images/bastex.gif"><BR>
'  | <A HREF="index.html">Back</A> |<BR>
' </BODY>
'</HTML>

SUB FindSmallLargeY
  FOR Count% = 0 TO 3
    Ycoord% = Polypoints%(Count%, 1)

    IF Ycoord% < Miny% THEN       ' is this the new lowest y co-ord?
      Miny% = Ycoord%             ' Yes...
    END IF

    IF Ycoord% > Maxy% THEN       ' is this the new highest y co-ord?
      Maxy% = Ycoord%             ' Yes...
    END IF
  NEXT
END SUB

SUB GetPolygonPoints
  RESTORE Polygonpoints:        ' start of un-rotated polygon co-ords

  FOR Count% = 0 TO 3
    READ Polypoints%(Count%, 0)  ' read x co-ord
    READ Polypoints%(Count%, 1)  ' read y co-ord
   
    RANDOMIZE TIMER
    Polypoints%(Count%, 0) = RND * 320
    Polypoints%(Count%, 1) = RND * 200
  NEXT
END SUB

SUB ScanConvert (X1%, Y1%, X2%, Y2%, Pside$)

' This procedure takes as defined by X1%,Y1%,x2,Y2%.
' It also takes a string telling it which side of the picture we are 
'mapping.  The strings are top,right,bottom,left. This routine decides 
'which 'side' of the polygon the line is on, and then calls the 
'appropriate routine.

IF Y2% < Y1% THEN
  SWAP X1%, X2% 'swap these variable round so we always scan from top           
SWAP Y1%, Y2% 'to bottom
  Lineheight% = Y2% - Y1%
  ScanLeftSide X1%, X2%, Y1%, Lineheight%, Pside$
ELSE
  Lineheight% = Y2% - Y1%
  ScanRightSide X1%, X2%, Y1%, Lineheight%, Pside$
END IF
END SUB

SUB ScanLeftSide (X1%, X2%, Ytop%, Lineheight%, Pside$)

' This procedure calculates the x points for the left side of the 
'polygon. It also calculates the x,y co-ords of the picture for the left 
'side of the polygon.

Lineheight% = Lineheight% + 1       ' No divide by zero
Xadd = (X2% - X1%) / Lineheight%

IF Pside$ <> "top" AND Pside$ <> "bottom" AND Pside$ <> "left" AND 
Pside$ <> "right" THEN PRINT "I'M MESSED UP!"

IF Pside$ = "top" THEN
  Px = Pwidth% - 1
  Py = 0
  Pxadd = -Pwidth% / Lineheight%
  Pyadd = 0
END IF
IF Pside$ = "right" THEN
  Px = Pwidth%
  Py = Pheight%
  Pxadd = 0
  Pyadd = -Pheight% / Lineheight%
END IF
IF Pside$ = "bottom" THEN
  Px = 0
  Py = Pheight%
  Pxadd = Pwidth% / Lineheight%
  Pyadd = 0
END IF
IF Pside$ = "left" THEN
  Px = 0
  Py = 0
  Pxadd = 0
  Pyadd = Pheight% / Lineheight%
END IF

x = X1%
FOR Y% = 0 TO Lineheight%
  PSET (x, Ytop% + Y%), 63
  Lefttable(Ytop% + Y%, 0) = x     'polygon x
  Lefttable(Ytop% + Y%, 1) = Px    'picture x
  Lefttable(Ytop% + Y%, 2) = Py    'picture y
  x = x + Xadd                      'Next polygon x
  Px = Px + Pxadd                   'Next picture x
  Py = Py + Pyadd                   'Next picture y
NEXT
END SUB

SUB ScanRightSide (X1%, X2%, Ytop%, Lineheight%, Pside$)

' This procedure calculates the x points for the right side of the ' 
polygon. It also calculates the x,y co-ords of the picture for the
' right side of the polygon.

Lineheight% = Lineheight% + 1    ' No divide by zero
Xadd = (X2% - X1%) / Lineheight%

IF Pside$ <> "top" AND Pside$ <> "bottom" AND Pside$ <> "left" AND 
Pside$ <> "right" THEN PRINT "I'M MESSED UP!"

IF Pside$ = "top" THEN
  Px = 0
  Py = 0
  Pxadd = Pwidth% / Lineheight%
  Pyadd = 0
END IF
IF Pside$ = "right" THEN
  Px = Pwidth%
  Py = 0
  Pxadd = 0
  Pyadd = Pheight% / Lineheight%
END IF
IF Pside$ = "bottom" THEN
  Px = Pwidth%
  Py = Pheight%
  Pxadd = -Pwidth% / Lineheight%
  Pyadd = 0
END IF
IF Pside$ = "left" THEN
  Px = 0
  Py = Pheight%
  Pxadd = 0
  Pyadd = -Pheight% / Lineheight%
END IF

x = X1%
FOR Y% = 0 TO Lineheight%
  PSET (x, Ytop% + Y%), 63
  Righttable(Ytop% + Y%, 0) = x    'polygon x
  Righttable(Ytop% + Y%, 1) = Px   'picture x
  Righttable(Ytop% + Y%, 2) = Py   'picture y
  x = x + Xadd                      'Next polygon x
  Px = Px + Pxadd                   'Next picture x
  Py = Py + Pyadd                   'Next picture y
NEXT Y%
END SUB

SUB TextureMap

' This is the actual mapping routine
' It takes the co-ords that have been calculated by the scan converter
' and 'traces' across the original picture in between them looking at 
' the pixel color and then plotting a pixel in that color in the 
' current position within the polygon.

FOR Y% = Miny% TO Maxy%
  Polyx1 = Lefttable(Y%, 0)  'get left polygon x
  Px1 = Lefttable(Y%, 1)      'get left picture x
  Py1 = Lefttable(Y%, 2)      'get left picture y

  Polyx2 = Righttable(Y%, 0) 'get right polygon x
  Px2 = Righttable(Y%, 1)     'get right picture x
  Py2 = Righttable(Y%, 2)     'get right picture y
  Linewidth% = Polyx2 - Polyx1  'what is the width of this polygon line
  IF Linewidth% = 0 THEN
  Linewidth% = Linewidth% + 1 'QUICK fix so it doesn't do divide by zero
  END IF
  Pxadd = (Px2 - Px1) / Linewidth% 'squash picture xdist into poly xdist
  Pyadd = (Py2 - Py1) / Linewidth% 'squash picture ydist into poly ydist

  FOR x% = Polyx1 TO Polyx2
    Col% = POINT(Px1, Py1)   'get color of pixel at current pos in pic
    PSET (x%, Y%), Col%     ' plot the pixel in the correct color
    Px1 = Px1 + Pxadd        ' move x picture co-ord
    Py1 = Py1 + Pyadd        ' move y picture co-ord
  NEXT

NEXT
END SUB

To make a texture for a sphere, first us spherical coordinates to first define a sphere.    For 
all values of theta and phi with r as a constant radius of the sphere.  

X = r * SIN(theta) * COS(phi)
Y = r * SIN(theta) * SIN (phi)
Z = r * COS(theta)  

Now in terms of u and v, generally used for variables in texture mapping.

X = r * SIN(v * PI) * COS(u * 2 * PI)
Y = r * SIN(v * PI) * SIN(u * 2 * PI)
Z = r * COS(v * PI)

Now solving for u, v

U = arcSIN(z / r) / PI
V = (arcCOS(x / (r * SIN(v * PI)))) / (2 * PI)

Other types of texture mapping are perspective correct texture mapping and bump 
mapping.

Appendix A: Binary, Decimal, and Hexadecimal

Binary, as you might know, is used by computers.  Decimal is what most people use and 
hexadecimal is a base system used by programmers to in easier to remember than binary.  
Binary is base two; decimal is base ten; hexadecimal is base sixteen.  What the base means 
is what the highest number can be in was column.  As you know in decimal you have the 
1's, 10's, 100's, and so on, columns, what you might not have realized is the these numbers 
are powers of 10 (this where the base 10 comes form).  

100  = 1
101 = 10 
102 = 100

The number in each column is multiplied by the base of the number system to the columns 
position to the left of zero'th power.  Binary has a 1's then 2's then 4's, 8's, and so on, 
columns.  Hexadecimal has a 1's, 16's, 256's, and so on, column's.  Hexadecimal is noted 
with and h after the number, or in BASIC with a &H before it. 

 Here are some examples in all three number systems.

Decimal	Binary	Hexadecimal	
1	1	1	
2	10	2	
3	11	3	
4	100	4	
5	101	5	
6	110	6	
7 	111	7	
8	1000	8	
9	1001	9	
10	1010	10	
11	1011	a	
12	1100	b	
13	1101	c	
14	1110	d	
15	1111	e	
16	10000	f	

Appendix B: Bitwise Operators

A Bitwise operator is comparing of each bit of a number's binary equivalence.  A bit is 
either a 1 or a 0, the numbers used in binary.  There are six operators that I know of, they 
are listed as the following:  AND, OR, XOR, NOT, EQV, IMP.  

The AND operator sees if both values being compared have a 1 in the same bit position.  
ex.  9 AND 5 = x
9 = 1001 in binary
5 = 0101 in binary
x = 0001 in binary because in bit position one both numbers were 1
x = 1 in decimal

The OR, inclusive, operator sees if one or both of the values is/are 1.
ex. 9 OR 5 = x
9 = 1001 in binary
5 = 0101 in binary
x = 1101 in binary because there are 1's in bit positions 1, 3, and 4
x = 13 in decimal

The XOR, exclusive, operator sees if one and only one of the bit positions is 1
ex. 9 XOR 5 = x
9 = 1001 in binary
5 = 0101 in binary
x = 1100 in binary because in positions 3 and 4 1 appeared once
x = 12 in decimal

The NOT operator looks if the first value has a 0 in any bit position.
ex. 9 NOT 5 = x
9 = 1001 in binary
5 = 0101 in binary
x = 0110 in binary because 0's appeared in bit positions 2 and 3
x= 6 in decimal

The EQV operator sees if both bit positions are the same, either 1 or 0.
ex. 9 EQV 5 = x
9 = 1001 in binary
5 = 0101 in binary
x = 0011 in binary because in bit positions 1 and 2 the bits were the same
x = 3 in decimal

The IMP operator sees if the first value's bit isn't 1 and the second values bit isn't 0


ex. 9 IMP 5 = x
9 = 1001 in binary 
5 = 0101 in binary
x = 0111 in binary
x = 7 in decimal

Appendix C: Vector Mathematics

	Although this is usually taught in Calculus and Linear Algebra college courses, it is 
very easy to understand.  A vector is really a quantity, and it is represented by a line 
segment with a direction and a size(or length, or module, or norm, or magnitude).  In this 
paper vectors have been used as rays from the origin, point (0, 0, 0)  through point p.  
Vectors are usually named a, b, c, ....  All the basic (add, subtract, multiply, divide) 
operations can be done to a vector.  

	Vector addition is easy.  Just perform each operation to each point of the same 
type, x to x, y to y, z to z.  ex.  given vectors A and B, A + B means the ordered triplet 
(vector A's X coordinate + vector B's X coordinate, vector A's Y coordinate + vector B's 
Y coordinate, vector A's Z coordinate + vector B's Z coordinate).  To make it even easier, 
I'll substitute number's for variables.  Given vector A = (10, 11, 12) and B = (3, 4, 5) and 
A + B = C, C = (10 + 3, 11 + 4, 12 + 5), or more simply, C = (13, 15, 17).  The same 
thing applies for subtraction,  since subtraction is the same as adding the additive inverse, 
or opposite.  

	There are more than one different vector multiplication, and divisions obviously.  
The first is scaling a vector, or increasing or decreasing, the magnitude of the vector.  
Given Vector A and a scale factor S, S * A means (A's X * S, A's Y * S, A's Z * S).  If S 
= -1 then this is written   -A instead of -1 * A. Next there is the dot product.  The dot 
product is a multiplication of two vectors. The result is a quantity (not a vector) and is 
written, with vectors A and B, A dot B.  Once again with vectors A and B, the dot 
product would be Ax *  Bx + Ay * By + Az * Bz.  The dot product can be used for many 
things.  The magnitude, written |A| for vector A. has a value of (A dot A) ^ (1 / 2), this is 
3D distance formula in a different form.  Dividing any vector by its length is called 
normalizing the vector. Normalizing a vector reduces it's length to 1, this eventually leads 
to faster and less math.  The dot product can also be used to find the angle between any 
two vectors.  For vectors A and B, and theta being the angle between them, A dot B = |A| 
* |B| * COS theta.  The other type of vector multiplication is the cross product.  The cross 
product of two vectors, once again I will use A and B, is (Ay * Bz - Az * Ay, Az * Bx - 
Bz * Ax, Ax * By - Ay * Bx).  The cross product of two vectors is perpendicular or 
orthonormal (forms a right angle) to the plane of the two vectors, and has a length of 
|A||B|SIN theta, when theta is the angle between the two vectors.

Appendix D: Matrices

	A matrix is way to express systems of linear equations.  The actual matrix looks 
something like a grid.  

        (m11, m12, m13)
M = (m21, m22, m23)
        (m31, m32, m33)

The above is an example of a 3 x 3 matrix.  It has 3 rows and 3 columns.  It can also be 
written M = (mij), where i is the row and j is the column.  If this looks hard, don't worry.  
You have all used matrices if you knew it of not.  The multiplication tables from second 
grade are matrices.  If a matrix has the same number of rows and columns then it is a 
square matrix.  An identity matrix is a square matrix that (mij) = 1 if i = j.  It looks like this

(1, 0, 0)
(0, 1, 0)
(0, 0, 1)

Matrix addition, along with subtraction, is almost like adding vectors.  For two matrices 
with then same dimensions, add the same type.  Matrix addition is written A + B = C, 
which means (aij) + (bij) = (cij).  ex.

(1, 0, 1)    (2, 3, 4)    (1 + 2, 0 + 3, 1 + 4)     (3, 3, 5)
(1, 0, 2) + (4, 5, 6) = (1 + 4, 0 + 5, 2 + 6) = (5, 5, 8)
(3, 4, 5)    (6, 3, 1)    (3 + 6, 4 + 3, 5 + 1)     (9, 7, 6)

With Matrices there are, like in vector mathematics, two ways of multiplying matrices.  
Matrices can be multiplied by a scalar or another matrix.  Multiplication by a scalar is 
written, where s is the scalar, M * s, which means to multiply each element, or part, of the 
matrix by k.  ex.

(11, 12, 13)            (110, 120, 130)
(21, 22, 23) * 10 = (210, 220, 230)
(31, 32, 33)            (310, 320, 330)

Matrix multiplication is a little hard to understand at first, but is just as easy as the last two 
operations.  For matrices A, B, and C, given that they have all the same dimensions, A * B 
= C, means that (cij) equals the sum of all (aik) * (bkj), for k = 1 to the number of 
columns.  This is rather vague, so I hope this example helps.  (notice that the first matrix's 
first position is the first position of the final matrix and the second matrix's first position is 
the final matrix's second position.)

(1, 2)(10, 11)    (1 * 10 + 2 * 13, 1 * 11 + 2 * 14)     (36, 39)
(4, 5)(13, 14) = (4 * 10 + 5 * 13, 4 * 11 + 5 * 14) = (105, 114) 

An important note is that matrix multiplication is not commutative, A * B <> B * A.  Also 
any matrix multiplied by an identity matrix is the original, non-identity, matrix. This is why 
it is called an identity matrix because it returns the same identity of the original matrix.  

Appendix E: QBASIC

	QBASIC is a high level programming language that is very easy to learn, and I 
make many references to, especially the math set.  The mathematics symbols are 

+ add 
- subtract
* multiply
/ divide
\ integer divide (rounds to closes integer, estimates)
^ n to the nth power
SQR(n) square root of n
SIN(n) sine of n
COS(n) cosine of n
TAN(n) tangent of n
MOD modulo, divides by number and returns remainder

	The above and following are some things someone new to QBASIC will need to 
understand the programming examples.

	There are different types of data formats for variables.  There is the INTEGER, 
called short in some other programming languages.  The INTEGER uses to bytes to save 
a signed number.  A signed number is a number that can be positive or negative.  An 
unsigned number is positive.  With two bytes (16 bits) you can save 65536 numbers.  This 
is because a bit can be either 1 or 0, or two combinations.  16 bits have 2 ^ 16 many 
combinations.  But since an INTEGER is a signed number, then it can have a range from
-32,768 to 32,767.  A LONG is an integer that uses four bytes (32 bits) and has 2 ^ 32 or 
4294967296 different combinations.  A LONG is also a signed number so it has a range of 
-2,147,483,648 to 2,147,483,647.  A SINGLE is a floating point number.  This means that 
it has a decimal.  This makes math slow because the interpreter needs to remember where 
the decimal point is.  A SINGLE use four bytes (32 bits).  A DOUBLE is the same as a 
single but it uses eight bytes (64 bits) to store it's number.  A STRING is a variable that 
stores text.  QBASIC has some shortcuts to defining the variables' data types.  A variable 
can have a suffix of % for INTEGER, & for LONG, ! for SINGLE,  # for DOUBLE, or $ 
for STRING.

	An array is a variable that has more than one value stored in it.  A value is 
accessed by way of elements, or parts of, of the array's dimensions that are in parenthesis.  
ex.

array(0) = 5 
This looks at the value of array at 0 at dimension 1

An array can have more than one dimension like the following: Array (1, 4). But first an 
array must be DIMensioned by the DIM command.  DIM array( lower bound of 
dimension 1 TO upper bound of dimension 1).  If the lower bound is omitted then it is  0 
by default. The DIM command can also be used to define a variable's, or array's, data type. 

Appendix F: Futher Reading

Lithuim\VLA 3drotate.doc
---. 3dshade.doc
Garms, Christian.  3d Graphics in BASIC
Bourke, Paul.  Parametric Equation of a sphere and texture mapping.
Barrett,Sean.  Texture Mapping
Helminen, Hannu .  Free Direction Texture Mapping
Loisel, Sebastien .  ZED3D.doc.
FAQ: 3-D Information for the Programmer
Denthor of Asphyxia.  VGA trainer 14
---.  VGA Trainer 17
---.  VGA Trainer 20
---.  VGA Trainer 21
---.  VGA Trainer 8
---.  VGA Trainer 9
Voltaire/OTM.  Polygons and Shading.
---.  Real-time Phong Shading
S-Buffers.  Paul Nettle
BSP Tree Frequently Asked Questions
comp.graphics.algorithms FAQ
