Here is the sequence of calculations to do per frame, commonly called the "rendering pipeline":

     Clear both the pix buffer and the zbuffer
     Loop through geometry objects:
        Compute MN, the inverse transpose of animation matrix M for this geometry object
        Loop through all vertices of the geometry object:
           For each vertex (x,y,z,nx,ny,nz):
              Use matrix M to transform x,y,z to x',y',z'
              Use matrix MN to transform nx,ny,nz to nx',ny',nz'
              Use the Phong algorithm to transform nx',ny',nz' to r,g,b

        For each face of the geometry object:
           If face has more than three vertices, split into triangles

           For each triangle:
              Clip triangle to the plane z = -ε
              If clipped shape has four vertices, split into two triangles to render separately.
              Apply perspective transform to each vertex. Now you have (px,py,pz,r,g,b).
              Apply viewport transform to each vertex.  Now you have (i,j,pz,r,g,b).
              Split triangle into two trapezoids, in preparation for scan-conversion.

              For each of the two trapezoids with corners (TL, TR, BL, BR):
                 Each corner is located at pixel (i,j) and contains a value of (r, g, b, pz)
                    where pz is the perspective z coordinate.

                 Loop through each scan scan-line j of the trapezoid between TLj and BLj:
                    interpolant tj = (j - TLj) / (BLj - TLj) 
                    left edge point  L = lerp(tj, TL, BL)
                    right edge point R = lerp(tj, TR, BR)

                    Loop through pixels i of scan-line, between Li and Ri:
                       interpolant ti = (i - Li) / (Ri - Li) 
                       (r,g,b,pz) = lerp(ti, L, R)

                       If (pz > zbuffer(i,j))
                          set zbuffer(i,j) to pz
                          set pix(i,j) to r,g,b

Some helpful notes:

Here is a link to a class/method for computing the inverse of a 4×4 matrix: Invert.java

As we discussed in class, the algorithm for clipping a triangle against a plane (in this case the plane P = [0,0,-1,&epsilon] is:

    Start with an empty output polygon

    For each vertex A:

       If A is on the positive side of the polygon (ie: dot(P,A) ≥ 0),
          then add A to the output polygon.

       Consider the edge A,B, where B is the next vertex around the polygon.

       If the sign of dot(P,A) differs from the sign of dot(P,B), then do:
          t = (0 - dot(P,A)) / (dot(P,B) - dot(P,A))
          Add interpolated vertex C = lerp(t, A, b) to the output polygon.

If you have trouble with the clipping, then you can cheat by just rejecting any triangle that has at least one vertex on the positive z side3. It will give you occasional "popping" artifacts, but it's better than nothing.