Z-buffering algorithm not drawing 100% correctly

I'm programming a software renderer in Java, and am trying to use Z-buffering for the depth calculation of each pixel. However, it appears to work inconsistently. For example, with the Utah teapot example model, the handle will draw perhaps half depending on how I rotate it.

My z-buffer algorithm:

for(int i = 0; i < m_triangles.size(); i++)
{
    if(triangleIsBackfacing(m_triangles.get(i))) continue; //Backface culling
        for(int y = minY(m_triangles.get(i)); y < maxY(m_triangles.get(i)); y++)
        {
            if((y + getHeight()/2 < 0) || (y + getHeight()/2 >= getHeight())) continue; //getHeight/2 and getWidth/2 is for moving the model to the centre of the screen
            for(int x = minX(m_triangles.get(i)); x < maxX(m_triangles.get(i)); x++)
            {
                if((x + getWidth()/2 < 0) || (x + getWidth()/2 >= getWidth())) continue;
                rayOrigin = new Point2D(x, y);
                if(pointWithinTriangle(m_triangles.get(i), rayOrigin))
                {
                    zDepth = zValueOfPoint(m_triangles.get(i), rayOrigin);
                    if(zDepth > zbuffer[x + getWidth()/2][y + getHeight()/2])
                    {
                        zbuffer[x + getWidth()/2][y + getHeight()/2] = zDepth;
                        colour[x + getWidth()/2][y + getHeight()/2] = m_triangles.get(i).getColour();
                        g2.setColor(m_triangles.get(i).getColour());
                        drawDot(g2, rayOrigin);
                    }
               }
          }
     }
}

Method for calculating the z value of a point, given a triangle and the ray origin:

private double zValueOfPoint(Triangle triangle, Point2D rayOrigin) 
{
    Vector3D surfaceNormal = getNormal(triangle);
    double A = surfaceNormal.x;
    double B = surfaceNormal.y;
    double C = surfaceNormal.z;
    double d = -(A * triangle.getV1().x + B * triangle.getV1().y + C * triangle.getV1().z);
    double rayZ = -(A * rayOrigin.x + B * rayOrigin.y + d) / C;
    return rayZ;
}

Method for calculating if the ray origin is within a projected triangle:

private boolean pointWithinTriangle(Triangle triangle, Point2D rayOrigin)
{
    Vector2D v0 = new Vector2D(triangle.getV3().projectPoint(modelViewer), triangle.getV1().projectPoint(modelViewer));
    Vector2D v1 = new Vector2D(triangle.getV2().projectPoint(modelViewer), triangle.getV1().projectPoint(modelViewer));
    Vector2D v2 = new Vector2D(rayOrigin, triangle.getV1().projectPoint(modelViewer));
    double d00 = v0.dotProduct(v0);
    double d01 = v0.dotProduct(v1);
    double d02 = v0.dotProduct(v2);
    double d11 = v1.dotProduct(v1);
    double d12 = v1.dotProduct(v2);

    double invDenom = 1.0 / (d00 * d11 - d01 * d01);
    double u = (d11 * d02 - d01 * d12) * invDenom;
    double v = (d00 * d12 - d01 * d02) * invDenom;

    // Check if point is in triangle
    if((u >= 0) && (v >= 0) && ((u + v) <= 1))
    {
         return true;
    }
    return false;
}

Method for calculating surface normal of a triangle:

private Vector3D getNormal(Triangle triangle)
{
    Vector3D v1 = new Vector3D(triangle.getV1(), triangle.getV2()); 
    Vector3D v2 = new Vector3D(triangle.getV3(), triangle.getV2());
    return v1.crossProduct(v2);
}

Example of the incorrectly drawn teapot:

What am I doing wrong? I feel like it must be some small thing. Given that the triangles draw at all, I doubt it's the pointWithinTriangle method. Backface culling also appears to work correctly, so I doubt it's that. The most likely culprit to me is the zValueOfPoint method, but I don't know enough to know what's wrong with it.

My zValueOfPoint method was not working correctly. I'm unsure why :( however, I changed to a slightly different method of calculating the value of a point in a plane, found here: http://forum.devmaster.net/t/interpolation-on-a-3d-triangle-using-normals/20610/5

To make the answer here complete, we have the equation of a plane: A * x + B * y + C * z + D = 0 Where A, B, and C are the surface normal x/y/z values, and D is -(Ax0 + By0 + Cz0).

x0, y0, and z0 are taken from one of the vertices of the triangle. x, y, and z are the coordinates of the point where the ray intersects the plane. x and y are known values (rayOrigin.x, rayOrigin.y) but z is the depth which we need to calculate. From the above equation we derive: z = -A / C * x - B / C * y - D

Then, copied from the above link, we do: "Note that for every step in the x-direction, z increments by -A / C, and likewise it increments by -B / C for every step in the y-direction. So these are the gradients we're looking for to perform linear interpolation. In the plane equation (A, B, C) is the normal vector of the plane. It can easily be computed with a cross product.

Now that we have the gradients, let's call them dz/dx (which is -A / C) and dz/dy (which is -B / C), we can easily compute z everywhere on the triangle. We know the z value in all three vertex positions. Let's call the one of the first vertex z0, and it's position coordinates (x0, y0). Then a generic z value of a point (x, y) can be computed as:"

z = z0 + dz/dx * (x - x0) + dz/dy * (y - y0)

This found the Z value correctly and fixed my code. The new zValueOfPoint method is:

private double zValueOfPoint(Triangle triangle, Point2D rayOrigin) 
{
    Vector3D surfaceNormal = getNormal(triangle);
    double A = surfaceNormal.x;
    double B = surfaceNormal.y;
    double C = surfaceNormal.z;
    double dzdx = -A / C;
    double dzdy = -B / C;
    double rayZ = triangle.getV1().z * modelViewer.getModelScale() + dzdx * (rayOrigin.x - triangle.getV1().projectPoint(modelViewer).x) + dzdy * (rayOrigin.y - triangle.getV1().projectPoint(modelViewer).y);
    return rayZ;
}

We can optimize this by only calculating most of it once, and then adding dz/dx to get the z value for the next pixel, or dz/dy for the pixel below (with the y-axis going down). This means that we cut down on calculations per polygon significantly.