_3-D TEXTURE MAPPING_
by Jeremy Spiller

Listing One

// *************************************************************************
// Texture mapping, copyright (C) 02/01/1993 by Jeremy Spiller
// This file contains the implementation for the texture mapping algorithm.
// *************************************************************************

#include <stdlib.h>
#include <iostream.h>
#include <conio.h>
#include <dos.h>
#include <math.h>
unsigned _stklen = 64000U;  // --- 64K of stack space ---

// --- Global variables ---
#define max(X,Y) ( ((X) > (Y)) ? (X) : (Y) )
#define min(X,Y) ( ((X) < (Y)) ? (X) : (Y) )
#define SCR_MAX_X 319
#define SCR_MAX_Y 199
int MX = 160, MY = 100;
int MinY, MaxY, MinX [SCR_MAX_Y+1], MaxX [SCR_MAX_Y+1];
char *Map1, *Map2, *Map3, *PhysicalScreen;

// ===== Graphics functions =====
char *Pixel_Pointer, *ScreenBuffer;
// --- Set the (X, Y) position to draw pixels future pixels ---
#define Set_Pixel_Position(X,Y) \
    (Pixel_Pointer = ScreenBuffer + 320*(Y) + (X))
// --- Draw a pixel, then move the position to the right ---
#define Plot_Pixel_Accross(Color) \
    (*Pixel_Pointer++ = (Color))
// --- Call a bios function to set the screen to graphics mode ---
void ibm_graphics_mode ()
  {
    _AH =0x00;            // --- Set video mode       ---
    _AL =0x13;            // --- Mode = VGA (320x200) ---
    geninterrupt (0x10);  // --- Perform mode change  ---
  }
// --- Calls a bios function to set the screen to text mode ---
void ibm_text_mode ()
  {
    _AH =0x00;            // --- Set video mode       ---
    _AL =0x03;            // --- Mode = text          ---
    geninterrupt (0x10);  // --- Perform mode change  ---
  }
// ===== MAT3D - Matrix functions for 3D matrix operations. This is a standard
//    mathematical matrix, accessed as M [Row][Column] (or M [Y][X]) 
//    where M[0][0] is the first element.
class MAT_COLUMN
  {
    double Column [4];
    public:
    double &operator [] (int Index)
      {
    return Column [Index];
      }
  };
// --- 3D matrix class ---
class MAT3D
  {
    MAT_COLUMN Matrix [4];
    public:
    MAT_COLUMN &operator [] (int Index)
      {
    return Matrix [Index];
      }
    friend MAT3D operator * (MAT3D &A, MAT3D &B);
  };
// --- Multiply two matrices, returns A*B ---
MAT3D operator * (MAT3D &A, MAT3D &B)
  {
    int I, J, K;
    MAT3D Temp;
    for (I = 0;  I < 4;  I++)
      for (J = 0;  J < 4;  J++)
    {
      double Sum = 0;
      for (K = 0;  K < 4;  K++)
        Sum += A [I][K] * B [K][J];
      Temp [I][J] = Sum;
    }
    return Temp;
  }
// --- Produce an identity matrix ---
MAT3D Identity ()
  {
    int X, Y;
    MAT3D Temp;
    for (Y = 0;  Y < 4;  Y++)
      for (X = 0;  X < 4;  X++)
    if (X == Y)
      Temp [Y][X] = 1;
    else
      Temp [Y][X] = 0;

    return Temp;
  }
// --- Calculate pitch - rotation around the X axis ---
MAT3D RotateX (double Angle)
  {
    MAT3D Temp = Identity ();
    Temp [1][1] = Temp [2][2] = cos (Angle);
    Temp [2][1] = - ( Temp [1][2] = sin (Angle) );
    return Temp;
  }
// --- Calculate yaw - rotation around the Y axis ---
MAT3D RotateY (double Angle)
  {
    MAT3D Temp = Identity ();
    Temp [0][0] = Temp [2][2] = cos (Angle);
    Temp [0][2] = - ( Temp [2][0] = sin (Angle) );
    return Temp;
  }
// --- Calculate roll - rotation around the Z axis ---
MAT3D RotateZ (double Angle)
  {
    MAT3D Temp = Identity ();
    Temp [0][0] = Temp [1][1] = cos (Angle);
    Temp [1][0] = - ( Temp [0][1] = sin (Angle) );
    return Temp;
  }
// --- Calculate pitch, yaw, and roll ---
MAT3D RotateXYZ (double AngleX, double AngleY, double AngleZ)
  {
    return RotateZ (AngleZ) * RotateY (AngleY) * RotateX (AngleX);
  }
// --- Calculate Scale - enlarge or shrink model (about origin) ---
MAT3D Scale (double Scale)
  {
    MAT3D Temp = Identity ();
    Temp [0][0] = Temp [1][1] = Temp [2][2] = Scale;
    return Temp;
  }
// --- Calculate Scale for each axis (can cause distortion) ---
MAT3D ScaleXYZ (double SX, double SY, double SZ)
  {
    MAT3D Temp = Identity ();
    Temp [0][0] = SX;
    Temp [1][1] = SY;
    Temp [2][2] = SZ;
    return Temp;
  }
// --- Translate - move model ---
MAT3D Translate (double X, double Y, double Z)
  {
    MAT3D Temp = Identity ();
    Temp [0][3] = X;
    Temp [1][3] = Y;
    Temp [2][3] = Z;
    return Temp;
  }
// --- Translate a point through a matrix, returns M * [X, Y, Z, 1]t ---
void Translate_Point (MAT3D &M, double &X, double &Y, double &Z)
  {
    double TX = X*M[0][0] + Y*M[0][1] + Z*M[0][2] + M[0][3];
    double TY = X*M[1][0] + Y*M[1][1] + Z*M[1][2] + M[1][3];
    double TZ = X*M[2][0] + Y*M[2][1] + Z*M[2][2] + M[2][3];
    X = TX;
    Y = TY;
    Z = TZ;
  }
// -- Draw a line clipped on Y axis, and has only one point per scan line. --
void do_line (void Plot(long X,long Y),double X1,double Y1,double X2,double Y2)
  {
    double YL = 0, YH = 200, Temp, DeltaX;
    if ((Y1 < YL && Y2 < YL) || (Y1 > YH && Y2 > YH))
      return;
    // --- if Y1 > Y2, swap (X1, Y1) with (X2, Y2) ---
    if (Y1 > Y2)
      {
    Temp = Y1;
    Y1 = Y2;
    Y2 = Temp;
    Temp = X1;
    X1 = X2;
    X2 = Temp;
      }
    // --- Is this a horizontal line? ---
    if (Y2 - Y1 == 0)
      {
    Plot (X1, Y1);
    Plot (X2, Y2);
    return;
      }
    DeltaX = (X2 - X1) / (Y2 - Y1);
    // --- Clip points ---
    if (Y1 < YL)
      {
    X1 = X1 + (YL - Y1) * DeltaX;
    Y1 = YL;
      }
    if (Y2 > YH)
      {
    X2 = X2 + (YH - Y2) * DeltaX;
    Y2 = YH;
      }
    // --- Draw line ---
    while (Y1 <= Y2)
      {
    Plot (X1, Y1);
    X1 += DeltaX;
    Y1 += 1;
      }
  }
// --- Find the minimum and maximum bounds of the plane ---
void MinMaxPixel (long X, long Y)
  {
    if (Y >= 0 && Y < SCR_MAX_Y)
      {
    MinX [Y] = min (MinX [Y], X);
    MaxX [Y] = max (MaxX [Y], X);
      }
  }
// --- Fill up the min/max outline ---
void Find_Outline (MAT3D &M, double Xp, double Yp)
  {
    double P1x, P1y, P2x, P2y, P3x, P3y, P4x, P4y;
    double Z1 = 0, Z2=0, Z3=0, Z4=0;
    int CountY;
    MaxY = 0;
    MinY = SCR_MAX_Y;
    // --- Initialize array values ---
    for (CountY = 0;  CountY < SCR_MAX_Y;  CountY++)
      {
    MinX [CountY] =  32000;
    MaxX [CountY] = -32000;
      }
    // --- Project four corners ---
    P1x = 0;  P1y = 0;
    P2x = Xp;  P2y = 0;
    P3x = Xp;  P3y = Yp;
    P4x = 0;  P4y = Yp;
    Translate_Point (M, P1x, P1y, Z1);
    Translate_Point (M, P2x, P2y, Z2);
    Translate_Point (M, P3x, P3y, Z3);
    Translate_Point (M, P4x, P4y, Z4);
    // --- Projection formula ---
    P1x = P1x/fabs(Z1) + MX;
    P1y = P1y/fabs(Z1) + MY;
    P2x = P2x/fabs(Z2) + MX;
    P2y = P2y/fabs(Z2) + MY;
    P3x = P3x/fabs(Z3) + MX;
    P3y = P3y/fabs(Z3) + MY;
    P4x = P4x/fabs(Z4) + MX;
    P4y = P4y/fabs(Z4) + MY;
    // --- Calculate min and max values of outline ---
    do_line (MinMaxPixel, P1x, P1y, P2x, P2y);
    do_line (MinMaxPixel, P2x, P2y, P3x, P3y);
    do_line (MinMaxPixel, P3x, P3y, P4x, P4y);
    do_line (MinMaxPixel, P4x, P4y, P1x, P1y);
  }
// ==== Project a texture map onto the screen.  The size of the map is
// [0..Xp] [0..Yp]. It will be transformed through the matrix M,
// translated by the variables MX and MY, and drawn on the screen.
void Project_Plane (MAT3D &M, char *Map, int Xp, int Yp)
  {
    int GridX, GridY;
    long X, Y, Z, DeltaX, DeltaY, DeltaZ;
    int Ypos, Xpos, MaxXpos, LineLength;
    Find_Outline (M, Xp, Yp);
    // --- Calculate inverse of M * [X, Y, 0, 1]t ---
    double Sx1 = M [0][0], Sy1 = M [0][1], T1 = M [0][3];
    double Sx2 = M [1][0], Sy2 = M [1][1], T2 = M [1][3];
    double Sx3 = M [2][0], Sy3 = M [2][1], T3 = M [2][3];
    const long FIXED_POINT = 64;
    // --- Calculate X axis (Scale X) ---
    long SXx = (T2*Sy3 - T3*Sy2) * FIXED_POINT;   // Scale X
    long SXy = (T3*Sy1 - T1*Sy3) * FIXED_POINT;   // Scale Y
    long SXt = (T1*Sy2 - T2*Sy1) * FIXED_POINT;   // Translate
    // --- Calculate Y axis (Scale Y) ---
    long SYx =  (T3*Sx2 - T2*Sx3) * FIXED_POINT;  // Scale X
    long SYy =  (T1*Sx3 - T3*Sx1) * FIXED_POINT;  // Scale Y
    long SYt =  (T2*Sx1 - T1*Sx2) * FIXED_POINT;  // Translate
    // --- Calculate Z axis (Scale Z) ---
    long SZx = (Sx3*Sy2 - Sx2*Sy3) * FIXED_POINT;   // Scale X
    long SZy = (Sx1*Sy3 - Sx3*Sy1) * FIXED_POINT;   // Scale Y
    long SZt = (Sx2*Sy1 - Sx1*Sy2) * FIXED_POINT;   // Translate
    for (Ypos = 0;  Ypos < SCR_MAX_Y;  Ypos += 1)
      {
    Xpos = max (0, MinX [Ypos]);
    MaxXpos = min (SCR_MAX_X, MaxX [Ypos]);
    LineLength = MaxXpos - Xpos + 1;
    if (Xpos <= MaxXpos)
      {
        X = ((Xpos-MX)*SXx + (Ypos-MY)*SXy + SXt);
        DeltaX = SXx;
        Y = ((Xpos-MX)*SYx + (Ypos-MY)*SYy + SYt);
        DeltaY = SYx;
        Z = ((Xpos-MX)*SZx + (Ypos-MY)*SZy + SZt) | 1;
        DeltaZ = SZx & ~1;    // --- Force Z to stay odd ---
        Set_Pixel_Position (Xpos, Ypos);
        while (--LineLength >= 0)
          {
        GridX = X / Z;
        GridY = Y / Z;
        X += DeltaX;
        Y += DeltaY;
        Z += DeltaZ;
        Plot_Pixel_Accross (Map [256*GridY + GridX]);
          }
      }
      }
  }
// -- If you are compiling this for a 286 machine, or with a 286 compiler, 
//     replace the two lines (GridX = X/Z) and (GridY = Y/Z) with the following
//     assembly code. This will speed up the process by using a machine 
//     language divide instruction instead of a subroutine.  
//      asm   mov    al, byte ptr Y+3
//      asm   cbw
//      asm   mov    dx, ax
//      asm   mov    cx, word ptr Z+1
//      asm   mov    ax, word ptr Y+1
//      asm   idiv   cx
//      asm   mov    word ptr GridY, ax;
//
//      asm   mov    al, byte ptr X+3
//      asm   cbw
//      asm   mov    dx, ax
//      asm   mov    cx, word ptr Z+1
//      asm   mov    ax, word ptr X+1
//      asm   idiv   cx
//      asm   mov    word ptr GridX, ax
//
// --- Is the viewer on the visible (or invisible) side of the plane? ---
int Visible (MAT3D &M)
  {
    // --- Point1 = M * [0, 0, 0, 1]t ---
    double X1 = M [0][3];
    double Y1 = M [1][3];
    double Z1 = M [2][3];
    // --- Point2 = M * [0, 1, 0, 1]t ---
    double X2 = M [0][1] + M [0][3];
    double Y2 = M [1][1] + M [1][3];
    double Z2 = M [2][1] + M [2][3];
    // --- Point3 = M * [1, 0, 0, 1]t ---
    double X3 = M [0][0] + M [0][3];
    double Y3 = M [1][0] + M [1][3];
    double Z3 = M [2][0] + M [2][3];
    double D = -X1*(Y2*Z3-Y3*Z2) - X2*(Y3*Z1-Y1*Z3) - X3*(Y1*Z2-Y2*Z1);
    // --- If the viewer is on the positive side, the plane is visible ---
    return D > 0;
  }
// --- Return a matrix that distorts the model ---
MAT3D Distort_Function ()
  {
    MAT3D Distort = Identity ();
    // Distort [0][0] = 2;              // Stretch on X axis
    // Distort [1][0] = Distort [0][1] = .4;    // Parallelogram
    // Distort [3][1] = .005;           // Pyramid
    return Distort;
  }
// --- rotate a cube ---
main ()
  {
    int I;
    long Fx, Fy;
    // --- Set up memory ---
    PhysicalScreen = (char *)MK_FP (0xA000, 0x0000);
    ScreenBuffer = new char [320U*200];
    Map1 = new char [256U*200];
    if (!( Map2 = new char [256U*200] )) Map2 = Map1;
    if (!( Map3 = new char [256U*200] )) Map3 = Map1;
    if (ScreenBuffer == 0 || Map1 == 0)
      {
    cout << "Sorry, not enough memory to run this program!\n";
    exit (1);
      }
    // --- Draw a picture of crosses (X*Y) ---
    for (Fx = 0;  Fx < 256;  Fx++)
      for (Fy = 0;  Fy < 199;  Fy++)
    Map1 [256*Fy + Fx] = (Fx - 128) * (Fy - 100) / 256 + 128;
    // --- Draw a picture of circles (X*X + Y*Y) ---
    for (Fx = 0;  Fx < 256;  Fx++)
      for (Fy = 0;  Fy < 199;  Fy++)
    {
      long Dx = (Fx - 128);
      long Dy = (Fy - 100);
      Map2 [256*Fy + Fx] = (Dx*Dx + Dy*Dy) / 128 + 16;
    }
    // --- Draw some lines (X) ---
    for (Fx = 0;  Fx < 256;  Fx++)
      for (Fy = 0;  Fy < 199;  Fy++)
    Map3 [256*Fy + Fx] = Fx + 17;
     // --- Build 6 sides of a cube relative to the origin ---
     MAT3D Flip, S1, S2, S3, S4, S5, S6;
     Flip = RotateY (180 * 3.141592654 / 180);
     Flip = Translate (256, 0, 0) * Flip;
     S1 = Translate (0, 0, 0);
     S2 = Flip;
     S2 = RotateX (-90 * 3.141592654 / 180) * S2;
     S3 = Flip;
     S3 = Translate (-56, 0, 0) * S3;
     S3 = RotateY (90 * 3.141592654 / 180) * S3;
     S4 = Flip;
     S4 = Translate (0, 0, 200) * S4;
     S5 = RotateX (-90 * 3.141592654 / 180);
     S5 = Translate (0, 200, 0) * S5;
     S6 = RotateY (90 * 3.141592654 / 180);
     S6 = Translate (256, 0, 0) * S6;
     ibm_graphics_mode ();
     for (I = 0;  I < 20000;  I++)
       {
     // --- Move center of cube to origin, rotate cube, move   ---
     // --- cube back (away from viewer), and scale for screen ---
     MAT3D Cube = Translate (-128, -100, -100);
     Cube = Distort_Function () * Cube;
     Cube = RotateXYZ (I*0.051, I*0.01, I*0.037) * Cube;
     Cube = Translate (0, 0, sin (-I*0.0121) * 350 + 350 + 400) * Cube;
     Cube = ScaleXYZ (230, 200, 1) * Cube;  // --- scale for screen ---
     // --- Transform each side of the cube ---
     MAT3D Side1 = Cube*S1;
     MAT3D Side2 = Cube*S2;
     MAT3D Side3 = Cube*S3;
     MAT3D Side4 = Cube*S4;
     MAT3D Side5 = Cube*S5;
     MAT3D Side6 = Cube*S6;
     // --- Draw each side of the cube (if visible) ---
     memset (ScreenBuffer, 0, 320U*200);  // Clear screen buffer
     if (Visible (Side1))
       Project_Plane (Side1, Map1, 255, 199);
     if (Visible (Side2))
       Project_Plane (Side2, Map2, 255, 199);
     if (Visible (Side3))
       Project_Plane (Side3, Map3, 199, 199);
     if (Visible (Side4))
       Project_Plane (Side4, Map1, 255, 199);
     if (Visible (Side5))
       Project_Plane (Side5, Map2, 255, 199);
     if (Visible (Side6))
       Project_Plane (Side6, Map3, 199, 199);
     memcpy (PhysicalScreen, ScreenBuffer, 320U*200);  // Copy to screen
     if (kbhit ()) if (getch () == 27) break;  // Break at users request
       }
    ibm_text_mode ();
    cout << "Texture mapped cube copyright (C) 1993 by Jeremy Spiller.\n";
  }


Figure 5. 

(a)
NewPoint = Translate(0,0,100)*Scale(2)*RotateX(30 Deg)*Translate(DCX,DCY,DCZ)*OldPoint


(b)

T = Translate(0,0,100)*Scale(2)*RotateX(30 Degrees)*Translate(DCX,DCY,DCZ)


Figure 6. 

Xscreen = (AX+BY+C)/(GX+HY+I) 
Yscreen = (DX+EY+F)/(GX+HY+I) 


Figure 7 
	
A = Y1(Z2-Z3)+Y2(Z3-Z1)+Y3(Z1-Z2)
B = Z1(X2-X3)+Z2(X3-X1)+Z3(X1-X2)
C = X1(Y2-Y3)+X2(Y3-Y1)+X3(Y1-Y2)
D = -X1(Y2*Z3-Y3*Z2)-X2(Y3*Z1-Y1*Z3)-X3(Y1*Z2-Y2*Z1)