COMPUTE AND DRAW A BICUBIC SURFACE PATCH USING FORWARD DIFFERENCES - This code implements and provides corrections to the algorithm named DrawSurfaceFwdDif presented in Fig.11.46 at page 525 of the book Computer Graphics - Principles and Practice 2.ed in C by James D.Foley et.al.

ForwardDiff_Foley_vanDam.pde

/* ============================================================================= *
 *      COMPUTE AND DRAW A BICUBIC SURFACE PATCH USING FORWARD DIFFERENCES       *
 *                                                                               *
 * This code implements and provides corrections to the algorithm named          *
 * DrawSurfaceFwdDif presented in Fig.11.46 at page 525 of the book              *
 * Computer Graphics - Principles and Practice 2.ed in C by Jamed D.Foley et.al. *
 * This algorithm was neither corrected nor repeated in the 3rd (present)        *
 * edition of Foley et.al.'s book.                                               *
 *                                                                               *
 * The algorithm, as stated in the book, does not work: There is one error and   *
 * two necessary steps are omitted. We corrected these errors (with comments)    *
 * and implemented it in the Processing Language as an easy-to-test/easy-to-show *
 * way of demonstrating the implementation. Simply copy the code into a new      *
 * Processing project and run it.                                                *
 *                                                                               *
 * All code was implemented using Structured Programming and was maintained as   *
 * explicit and syntactically near the algorithm in Foley et.al. as possible.    *
 * High-level functions from Procesing were used only where they do not relate   *
 * directly to the algorithm, such as in performing perspective projection or    *
 * in drawing the control points as spheres.                                     *
 *                                                                               *
 * All used matrix operations such as multiplication and transposition are       *
 * explicitly implemented here in the code.                                      *
 *                                                                               *
 * The implementation philosophy was to produce a code that is easily            *
 * translatable into any other programming language, given that one observes     *
 * that the matrix initializations and operations suppose a row-major memory     *
 * model for matrix organization. If you translate it into e.g. FORTRAN, you'll  *
 * have to transpose all matrix operations and initializations. It will work     *
 * without any change in Java and C++.                                           *
 *                                                                               *
 * Aldo von Wangenheim                                                           *
 * http://www.inf.ufsc.br/~awangenh                                              *
 * https://www.researchgate.net/profile/Aldo_Von_Wangenheim2                     *
 * Santa Catarina Island, October, 2014                                          *
 * ============================================================================= */

/* ----------------------------------------------------------------------------- *
   Example Constants
 * ----------------------------------------------------------------------------- */                  

// Bezier Method Matriz - you can use B-Spline, Hermite, ....    
double [][] B =  { {-1,  3, -3,  1},
                  { 3, -6,  3,  0},
                  {-3,  3,  0,  0},
                  { 1,  0,  0,  0} };

// Sample control point Matrices
double[][] Ax = {  {-1, 0, 1, 2},
                  {-1, 0, 1, 2},
                  {-1, 0, 1, 2},
                  {-1, 0, 1, 2}};
                  
double[][] Ay = {  {3, 3, 3, 3},
                  {3, -2, -2, 3},
                  {3, -2, -2, 3},
                  {3, 3, 3, 3} };

double[][] Az = {  {1, 1, 1, 1},
                  {2, 2, 2, 2},
                  {3, 3, 3, 3},
                  {4, 4, 4, 4} };

/* ----------------------------------------------------------------------------- *
   Coefficient Matrices C
 * ----------------------------------------------------------------------------- */                  
double[][] Cx = new double[4][4];      
double[][] Cy = new double[4][4];      
double[][] Cz = new double[4][4];      

/* ----------------------------------------------------------------------------- *
   Parameter Matrices
 * ----------------------------------------------------------------------------- */
double[][] Et = new double[4][4];      
double[][] Es = new double[4][4];   

/* ----------------------------------------------------------------------------- *
   Forward Differences Matrices
 * ----------------------------------------------------------------------------- */
double[][] DDx = new double[4][4];      
double[][] DDy = new double[4][4];      
double[][] DDz = new double[4][4];      

/* ============================================================================= */
void printMatrix4x4(double mat [][])
/* ============================================================================= */
{
  for (int i = 0; i<4; i++)
    {
    print("[ ");
    for (int j=0; j<4; j++)
      {
        print(mat[i][j]);
        print(" ");
      }
    println("]");
    }
  println("");
}

/* ============================================================================= */
void multMatrix4x4WithConstant(double mat [][], double mult)
/* ============================================================================= */
{
  for (int i = 0; i<4; i++)
    {
    for (int j=0; j<4; j++)
      {
        mat[i][j] = mat[i][j] * mult;
      }
    }
  
}

/* ============================================================================= */
void multMatrix4x4(double mat1 [][], double mat2 [][], double result[][])
/* ============================================================================= */
{
  // Simple 4x4 Matrix Multiplication. Can be extended for any square matrix.
  // Makes use of the implicit parameter passing by reference/pointer employed by Java/Processing 
  // when calling functions with Arrays as arguments. Be careful when translating this into another
  // language. If you end up with an implementation where result[] is passed by value, the function 
  // won't work...
  for (int i = 0; i<4; i++)
    {
    for (int j=0; j<4; j++)
      {
         result[i][j] = 0.0;
      }
    }
  for (int i=0; i<4; i++)
    {
    for (int j=0; j<4; j++)
      {
        for (int k=0; k<4; k++)
        {
          result[i][j] = result[i][j] + (mat1[i][k] * mat2[k][j]);
        }
      }
    }
  
}

/* ============================================================================= */
void transpose(double mat [][])
/* ============================================================================= */
{
  double result[][] = new double[4][4];
  for (int i = 0; i<4; i++)
    {
    for (int j=0; j<4; j++)
      {
         result[i][j] = mat[j][i];
      }
    }
  for (int i = 0; i<4; i++)
    {
    for (int j=0; j<4; j++)
      {
         mat[i][j] = result[i][j];
      }
    }
}

/* ============================================================================= */
void drawGrid()
/* ============================================================================= */
{
  // Show the control point grid
  for (int i = 0; i < 3; i++) {
    for (int j = 0; j < 4; j++) {
      pushMatrix();
      translate((float)Ax[i][j], (float)Ay[i][j], (float)Az[i][j]);  
      noStroke();
      fill(#ff3333);
      sphere(2);
      popMatrix();  
      stroke(#5050FF);
      line((float)Ax[i][j], (float)Ay[i][j], (float)Az[i][j], (float)Ax[i+1][j], (float)Ay[i+1][j], (float)Az[i+1][j]);
    }
  }
  for (int i = 0; i < 4; i++) {
    for (int j = 0; j < 3; j++) {
      pushMatrix();
      translate((float)Ax[i][j], (float)Ay[i][j], (float)Az[i][j]);  
      noStroke();
      fill(#ff3333);
      sphere(2);
      popMatrix();  
      stroke(#5050FF);
      line((float)Ax[i][j], (float)Ay[i][j], (float)Az[i][j], (float)Ax[i][j+1], (float)Ay[i][j+1], (float)Az[i][j+1]);
    }
  }
}

/* ============================================================================= */
void calculateCoefficients()
/* ============================================================================= */
{
 /* Calculates the surface coefficient matrices Cx, Cy and Cz. 
  * Given by  <method matrix> * <geometry matrix> * <transposed method matrix>
  * Bezier and B-Spline have symmetric matrices, so there's no need to transpose.
  * Uses a buffer in order not to call multMatrix4x4(Cx, B, Cx);                 */
  double buffer[][] = new double[4][4];
  println("Cx:");
  multMatrix4x4(B, Ax, buffer);  // Mult <method matrix> and <geometry matrix> and store in a buffer
  multMatrix4x4(buffer, B, Cx);  // Mult the buffer and <transposed method matrix> and store in C
  printMatrix4x4(Cx);
  println("Cy:");
  multMatrix4x4(B, Ay, buffer);
  multMatrix4x4(buffer, B, Cy);
  printMatrix4x4(Cy);
  println("Cz:");
  multMatrix4x4(B, Az, buffer);
  multMatrix4x4(buffer, B, Cz);
  printMatrix4x4(Cz);
}

void createDeltaMatrices(double delta_s, double delta_t)
{
  /* In Foley's algorithm, the deltas are explicitly initialized,                         *
   * but the initialization of the E(delta) matrices is omitted.                          *
   * This code here should be called once after the initialization                        *
   * of the delta s and delta t.                                                          */
   
  // Delta s
  Es[0][0] = 0;
  Es[0][1] = 0;
  Es[0][2] = 0;
  Es[0][3] = 1;
  
  Es[1][0] = delta_s * delta_s * delta_s;
  Es[1][1] = delta_s * delta_s;
  Es[1][2] = delta_s;
  Es[1][3] = 0;
  
  Es[2][0] = 6 * delta_s * delta_s * delta_s;
  Es[2][1] = 2 * delta_s * delta_s;
  Es[2][2] = 0;
  Es[2][3] = 0;
  
  Es[3][0] = 6 * delta_s * delta_s * delta_s;
  Es[3][1] = 0;
  Es[3][2] = 0;
  Es[3][3] = 0;

  // Delta t
  Et[0][0] = 0;
  Et[0][1] = 0;
  Et[0][2] = 0;
  Et[0][3] = 1;
  
  Et[1][0] = delta_t * delta_t * delta_t;
  Et[1][1] = delta_t * delta_t;
  Et[1][2] = delta_t;
  Et[1][3] = 0;
  
  Et[2][0] = 6 * delta_t * delta_t * delta_t;
  Et[2][1] = 2 * delta_t * delta_t;
  Et[2][2] = 0;
  Et[2][3] = 0;
  
  Et[3][0] = 6 * delta_t * delta_t * delta_t;
  Et[3][1] = 0;
  Et[3][2] = 0;
  Et[3][3] = 0;
  
  // Transpose E(delta t)
  transpose(Et);
  
} //end createDeltaMatrices

void createForwardDiffMatrices()
{
  double buffer[][] = new double[4][4];
  // DD(x)
  multMatrix4x4(Es, Cx, buffer);
  multMatrix4x4(buffer, Et, DDx);
  // DD(y)
  multMatrix4x4(Es, Cy, buffer);
  multMatrix4x4(buffer, Et, DDy);
  // DD(z)
  multMatrix4x4(Es, Cz, buffer);
  multMatrix4x4(buffer, Et, DDz);
} //end createForwardDiffMatrices

/* =================================================================
 * Draw a single curve using Forward Differences. This is verbatim 
 * of Foley's algorithm in Fig. 11.36 on page 513. No changes 
 * were made besides substituting the moveAbs(x, y, z) call by
 * a buffering of the old x,y and z values at oldx, oldy and oldz.
 * ================================================================= */
 
void DrawCurveFwdDif( int n, 
                      double x, double Dx, double D2x, double D3x,
                      double y, double Dy, double D2y, double D3y,
                      double z, double Dz, double D2z, double D3z)
{
  int i = 0;
  double oldx, oldy, oldz;
  oldx = x;
  oldy = y;
  oldz = z;
  stroke(#FFFFFF);
  for (i=1; i < n; i++) {
    x = x + Dx;  Dx = Dx + D2x;  D2x = D2x + D3x;
    y = y + Dy;  Dy = Dy + D2y;  D2y = D2y + D3y;
    z = z + Dz;  Dz = Dz + D2z;  D2z = D2z + D3z;
    line((float)oldx, (float)oldy, (float)oldz, (float)x, (float)y, (float)z); // Draw line segment 
    oldx = x;
    oldy = y;
    oldz = z;    
  } //end for
} //end DrawCurveFwdDif


/* ============================================================================= */
void UpdateForwardDiffMatrices()
/* ============================================================================= */
{
  // Implements equation 11.92 at pg. 525 of Foley & van Dam
  // This code avoids using high level matrix row by row sum functions
  
  //row1 <- row1 + row2
  DDx[0][0] =  DDx[0][0]+DDx[1][0]; DDx[0][1] = DDx[0][1]+DDx[1][1]; DDx[0][2] = DDx[0][2]+DDx[1][2]; DDx[0][3] = DDx[0][3]+DDx[1][3];
  DDy[0][0] =  DDy[0][0]+DDy[1][0]; DDy[0][1] = DDy[0][1]+DDy[1][1]; DDy[0][2] = DDy[0][2]+DDy[1][2]; DDy[0][3] = DDy[0][3]+DDy[1][3];
  DDz[0][0] =  DDz[0][0]+DDz[1][0]; DDz[0][1] = DDz[0][1]+DDz[1][1]; DDz[0][2] = DDz[0][2]+DDz[1][2]; DDz[0][3] = DDz[0][3]+DDz[1][3];
  //row2 <- row2 + row3
  DDx[1][0] =  DDx[1][0]+DDx[2][0]; DDx[1][1] = DDx[1][1]+DDx[2][1]; DDx[1][2] = DDx[1][2]+DDx[2][2]; DDx[1][3] = DDx[1][3]+DDx[2][3];
  DDy[1][0] =  DDy[1][0]+DDy[2][0]; DDy[1][1] = DDy[1][1]+DDy[2][1]; DDy[1][2] = DDy[1][2]+DDy[2][2]; DDy[1][3] = DDy[1][3]+DDy[2][3];
  DDz[1][0] =  DDz[1][0]+DDz[2][0]; DDz[1][1] = DDz[1][1]+DDz[2][1]; DDz[1][2] = DDz[1][2]+DDz[2][2]; DDz[1][3] = DDz[1][3]+DDz[2][3];
  //row3 <- row3 + row4 
  DDx[2][0] =  DDx[2][0]+DDx[3][0]; DDx[2][1] = DDx[2][1]+DDx[3][1]; DDx[2][2] = DDx[2][2]+DDx[3][2]; DDx[2][3] = DDx[2][3]+DDx[3][3];
  DDy[2][0] =  DDy[2][0]+DDy[3][0]; DDy[2][1] = DDy[2][1]+DDy[3][1]; DDy[2][2] = DDy[2][2]+DDy[3][2]; DDy[2][3] = DDy[2][3]+DDy[3][3];
  DDz[2][0] =  DDz[2][0]+DDz[3][0]; DDz[2][1] = DDz[2][1]+DDz[3][1]; DDz[2][2] = DDz[2][2]+DDz[3][2]; DDz[2][3] = DDz[2][3]+DDz[3][3];  
  
}

/* ============================================================================= */
void DrawSurfaceFwdDif( int ns, int nt)
/* ============================================================================= */
{
  /* Uses global Coeficients Cx, Cy e Cz in x, y and z  */
  double delta_s = 1.0 / (ns - 1);
  double delta_t = 1.0 / (nt - 1);
  createDeltaMatrices(delta_s, delta_t);  // Omitted in Foley's algorithm
  createForwardDiffMatrices();
  // Draw ns curves along t.
  for (int i = 0; i < ns; i++) {
    
    // Foley et.al. calls DrawCurveFwdDif using an arbitrary "n". 
    // This does not work: if you initialized the Et matrix  
    // using a delta_t calculated from nt, you
    // necessarilly will have to plot the curves in nt steps. 
    /* DrawCurveFwdDif (nt, 
                         x, Dx, D2x, D3x,
                         y, Dy, D2y, D3y,
                         z, Dz, D2z, D3z); */
    DrawCurveFwdDif (nt, 
                     DDx[0][0], DDx[0][1], DDx[0][2], DDx[0][3], 
                     DDy[0][0], DDy[0][1], DDy[0][2], DDy[0][3], 
                     DDz[0][0], DDz[0][1], DDz[0][2], DDz[0][3] );
    UpdateForwardDiffMatrices();
  } //endfor ns
  
  // Regenerate DD matrices.
  // Calling the UpdateForwardDiffMatrices() function in the 
  // loop above changed the DD matrices.
  // You MUST reinitialize them. Omitted in Foley et.al.
  createForwardDiffMatrices();
  // Transpose FwdDif DD matrices in order to be able to continue to 
  // use row by row sums to update it by calling UpdateForwardDiffMatrices(). 
  // Alternatively you could write a version of UpdateForwardDiffMatrices() 
  // that performs column by column sums and use it below instead of transposing DD.  
  transpose(DDx);
  transpose(DDy);
  transpose(DDz);  
  // Do it all again for nt: Draw nt curves along s.
  for (int i = 0; i < nt; i++) {
    
    // Foley et.al. calls DrawCurveFwdDif using an arbitrary "n". 
    // This does not work: if you initialized the Es matrix  
    // using a delta_s calculated from ns, you
    // necessarilly will have to plot the curves in ns steps. 
    /* DrawCurveFwdDif (ns, 
                         x, Dx, D2x, D3x,
                         y, Dy, D2y, D3y,
                         z, Dz, D2z, D3z); */
    DrawCurveFwdDif (ns, 
                     DDx[0][0], DDx[0][1], DDx[0][2], DDx[0][3], 
                     DDy[0][0], DDy[0][1], DDy[0][2], DDy[0][3], 
                     DDz[0][0], DDz[0][1], DDz[0][2], DDz[0][3] );
    UpdateForwardDiffMatrices();
  } //endfor nt
  
} //end DrawSurfaceFwdDif

/* ============================================================================= */
void setup() 
/* ============================================================================= */
{
  size(1200, 1050, P3D);
  
  // Scale the control points in order to better visualize
  // the generated curves.
  // IS NOT PART OF THE ALGORITHM. 
  multMatrix4x4WithConstant(Ax, 40.0);
  multMatrix4x4WithConstant(Ay, 40.0);
  multMatrix4x4WithConstant(Az, 40.0);
  printMatrix4x4(Ax);
  printMatrix4x4(Ay);
  printMatrix4x4(Az);
  
  // This is an initialization step that IS part of the algorithm.
  calculateCoefficients();
  
}


/* ============================================================================= */
void draw() 
/* ============================================================================= */

{
  // Wait until the viewport window opens and is painted black, 
  // then move the mouse around, to force a screen update and 
  // the drawing of the surface.
  // Moving the mouse vertically rotates the surface.
  // Moving the mouse horizontally zooms.
  lights();
  background(0);
  float cameraY = height/2.0;
  float fov = mouseX/float(width) * PI/4;
  float cameraZ = cameraY / tan(fov / 2.0);
  float aspect = float(width)/float(height);
  perspective(fov, aspect, cameraZ/30.0, cameraZ*30.0);
  
  translate(width/2.0, height/1.9, 0);
  rotateX(-PI/6);
  rotateY(PI/3 + mouseY/float(height) * PI);
  drawGrid();
  DrawSurfaceFwdDif(15, 15);

}