[SOLVED] How to use face normal for back-face culling in perspective projection scenes

I am writing low level 2D/3D rendering engine as a part of display driver for MCU platform in C++ and I hit a wall with perspective projection in 3D and face culling.

Lets assume my 3D engine is using M,V,P matrices (in the same manner as Model,View and Projection matrices in OpenGL fixed pipeline).

The idea is to convert rasterized face normal into view local coordinates (using MV) and test the sign of coordinate corresponding to view direction. The same can be done with dot product between camera view direction and the normal itself. And according to sign either rasterize or skip face. This works well for parallel rays projections ... however with perspective projections this leads to false positives (you can see faces "visually" tilted away up to some angle).

For filled surfaces this poses only performance hit as depth buffer remedies the artifacts so render looks like it should. However wire frame is a problem:

The remedy is to transform vertexes of face by MVP and do the perspective divide. And then re-compute the normal from the result and use that for face culling:

However that is a lot of more operations which on slow platforms like MCU could pose a performance problem. So my question is:

If it is possible How to use face normal for back face culling safely?

I have tried to transform the normal locally by transforming face center and its small displacement in normal direction by MVP with perspective divide and then recomputing the normal from these 2 points. Still twice the operations then using directly normal but better than 3x. However result was not correct (looked almost identical to using normal directly).

I am thinking of somehow computing the tilt angle for given projection / location and test the:

dot(normal,view_direction) >= acos(min_tilt_angle_threshold)

However not sure how to compute the tilt angle nor if it is a safe assumption, nor if the scales in x,y direction have to be applied as projection has non trivial aspect ratio as usually the screen resolution is rectangle.

My current test rendering code looks like this:

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _mesh_h
#define _mesh_h
//---------------------------------------------------------------------------
#include "vecmath.h"
//---------------------------------------------------------------------------
int view_flags;
float view_mv[16],view_p[16],view_mvp[16],view_ofs[3],view_sc[3];
static volatile float angx=0.0,angy=0.0,angz=5.0*deg;
//---------------------------------------------------------------------------
void obj2mv(float *xyz,const float *dir,float w)
    {
    // apply MV
    mat4_mul_vec3(xyz,view_mv,dir,w);
    }
//---------------------------------------------------------------------------
void obj2mvp(float *p,const float *pos)
    {
    float w;
    // apply MVP
    w=divide(1.0,mat4_mulw_vec3(p,view_mvp,pos,1.0));
    // perspective divide by w + scale x,y to pixels
    p[0]*=view_sc[0]*w;
    p[1]*=view_sc[1]*w;
    p[2]*=w;
    }
//---------------------------------------------------------------------------
void obj2scr(int &x,int &y,int &z,const float *pos)
    {
    float p[3];
    obj2mvp(p,pos);
    // center screen + usable depth range <-1,+1> -> <0,255> (8 bit depth buffer)
    p[2]++; p[2]*=view_sc[2];
    vec3_add(p,p,view_ofs);
    // float -> int
    x=p[0];
    y=p[1];
    z=p[2];
    }
//---------------------------------------------------------------------------
void mesh_view()
    {
    const float z0   =  7.5;
    const float znear=  5.0;
    const float zfar = 10.0;
//  const float FOVy = 10.0*deg;
    const float FOVy = 30.0*deg;

    const float q[3]={0.0,0.0,0.0};
    // modelview
    float m0[16]=
        {
         1.0, 0.0, 0.0,0.0,
         0.0, 1.0, 0.0,0.0,
         0.0, 0.0, 1.0,0.0,
         0.0, 0.0, -z0,1.0,
        };
    vec3_ld(view_ofs,0.5*float(lcd.xs),0.5*float(lcd.ys),0.0);
    vec3_ld(view_sc,0.5*float(lcd.xs),0.5*float(lcd.ys),0.5*255.0);
    mat4_perspective(view_p,FOVy,float(lcd.xs)/float(lcd.ys),znear,zfar);
    mat4_one(view_mv);
    mat4_rotx(view_mv,q,angx);
    mat4_roty(view_mv,q,angy);
    mat4_rotz(view_mv,q,angz);
    mat4_mul(view_mv,view_mv,m0);
    mat4_mul(view_mvp,view_mv,view_p);
    angx=fmod(angx+1.0*deg,2.0*M_PI);
    angy=fmod(angy+5.0*deg,2.0*M_PI);
    angz=fmod(angz+2.0*deg,2.0*M_PI);
    view_flags=0;
    }
//---------------------------------------------------------------------------
void mesh_draw(const float *pnt,const int *fac) // render triangular mesh using matrix m
    {
    int mm;
    int i,i0,i1,i2;
    int x0,y0,z0,c0,x1,y1,z1,c1,x2,y2,z2,c2;
    float p[3],q[3],n[3],p0[3],p1[3],p2[3],np[3];
    for (i=0;fac[i]>=0;)
        {
        // point index for each vertexes of triangle face i
        i0=fac[i]; i++; i0+=i0+i0;
        i1=fac[i]; i++; i1+=i1+i1;
        i2=fac[i]; i++; i2+=i2+i2;
        // obtain screen space coordinates
        obj2scr(x0,y0,z0,pnt+i0);
        obj2scr(x1,y1,z1,pnt+i1);
        obj2scr(x2,y2,z2,pnt+i2);
        // computed normal
        vec3_sub(p,pnt+i1,pnt+i0);
        vec3_sub(q,pnt+i2,pnt+i1);
        vec3_mul(n,p,q);
        obj2mv(n,n,0.0);
        vec3_one(n,n);
        // normal with perspective
        obj2mvp(p0,pnt+i0);
        obj2mvp(p1,pnt+i1);
        obj2mvp(p2,pnt+i2);
        vec3_sub(p,p1,p0);
        vec3_sub(q,p2,p1);
        vec3_mul(np,p,q);
        vec3_one(np,np);

        // normal shading color
        c0=float(255.0*n[2]); c0=2+((14*c0)>>8);
        // flat shaded rendering
        c1=c2=c0;

        // back face culling
        if (np[2]<0.0) continue;

        // render face
//      mm=_mode_fill;                  // normal render
        mm=_mode_line; c1=c2=c0=15;     // wireframe render

        lcd.triangle(x0,y0,z0,c0,x1,y1,z1,c1,x2,y2,z2,c2,mm);
/*
        // render normals
        x0=(x0+x1+x2)/3;
        y0=(y0+y1+y2)/3;
        z0=(z0+z1+z2)/3;
        c0=16; c1=0;
        x1=x0+float(10.0*n[0]);
        y1=y0+float(10.0*n[1]);
        z1=z0+float(10.0*n[2]);
        lcd.line(x0,y0,z0,c0,x1,y1,z1,c1);
*/
        }
    }
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

Where mesh is hard-coded like this:

//---------------------------------------------------------------------------
float point_1[]=
    {
    -1.0,+1.0,-1.0,
    +1.0,+1.0,-1.0,
    +1.0,-1.0,-1.0,
    -1.0,-1.0,-1.0,
    -1.0,-1.0,+1.0,
    +1.0,-1.0,+1.0,
    +1.0,+1.0,+1.0,
    -1.0,+1.0,+1.0,
    };
const int points_1=sizeof(point_1)/sizeof(point_1[0]);
int face3_1[]=
    {
    0,1,2,
    0,2,3,
    4,5,6,
    4,6,7,
    3,2,5,
    3,5,4,
    2,1,6,
    2,6,5,
    1,0,7,
    1,7,6,
    0,3,4,
    0,4,7,
    -1
    };
//---------------------------------------------------------------------------

And vecmath.h is my lib for 3D vector math:

//---------------------------------------------------------------------------
//--- 32bit float vector math ver: 1.001 ------------------------------------
//---------------------------------------------------------------------------
#ifndef _vecmath_h
#define _vecmath_h
//---------------------------------------------------------------------------
#include <math.h>
//---------------------------------------------------------------------------
const float deg=M_PI/180.0;
const float rad=180.0/M_PI;
float vec3_tmp[4],mat4_tmp[16];
//---------------------------------------------------------------------------
float  divide       (float a,float b){ if (fabs(b)<1e-10) return 0.0; return a/b; }
float* vec3_ld      (float x,float y,float z);
float* vec3_ld      (float *p,float x,float y,float z);
void   vec3_rotx    (float *pos,const float *mid,float ang);
void   vec3_roty    (float *pos,const float *mid,float ang);
void   vec3_rotz    (float *pos,const float *mid,float ang);
void   mat4_rotx    (float *mat,const float *mid,float ang);
void   mat4_roty    (float *mat,const float *mid,float ang);
void   mat4_rotz    (float *mat,const float *mid,float ang);
void   vec3_copy    (float *c,const float *a);
void   vec3_abs     (float *c,const float *a);
void   vec3_one     (float *c,const float *a);
void   vec3_len     (float *c,const float *a,float l);
float  vec3_len     (         const float *a){ return sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])); }
float  vec3_len2    (         const float *a){ return (a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]); }
void   vec3_neg     (float *c,const float *a);
void   vec3_add     (float *c,const float *a,const float *b);
void   vec3_sub     (float *c,const float *a,const float *b);
void   vec3_mul     (float *c,const float *a,const float *b);
void   vec3_mul     (float *c,const float *a,const float  b);
void   vec3_mul     (float *c,const float  a,const float *b);
float  vec3_mul     (         const float *a,const float *b);
float* mat4_ld      (         float a0,float a1,float a2,float a3,float a4,float a5,float a6,float a7,float a8,float a9,float a10,float a11,float a12,float a13,float a14,float a15,float a16);
float* mat4_ld      (float *p,float a0,float a1,float a2,float a3,float a4,float a5,float a6,float a7,float a8,float a9,float a10,float a11,float a12,float a13,float a14,float a15,float a16);
void   mat4_copy    (float *c,const float *a);
void   mat4_set     (float *c,float  a);
void   mat4_one     (float *c);
void   mat4_add     (float *c,const float *a,const float *b);
void   mat4_sub     (float *c,const float *a,const float *b);
void   mat4_mul     (float *c,const float *a,const float *b);
void   mat4_mul_vec3(float *c,const float *a,const float *b,float w);
void   vec3_mul_mat4(float *c,const float *a,const float *b,float w);
float  mat4_mulw_vec3(float *c,const float *a,const float *b,float w);
float  vec3_mulw_mat4(float *c,const float *a,const float *b,float w);
void   mat4_subdet  (float *c,const float *a);
float  mat4_subdet  (         const float *a,int r,int s);
float  mat4_det     (         const float *a);
float  mat4_det     (         const float *a,const float *b);
void   mat4_inv     (float *c,const float *a);
void   mat4_perspective(float *per,float fovy,float aspect,float zNear,float zFar); // per[16] = perspective projection nmatrix fovy [rad] !!!
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
float* vec3_ld(float x,float y,float z)         { float *p=vec3_tmp;   p[0]=x; p[1]=y; p[2]=z; return p; }
float* vec3_ld(float *p,float x,float y,float z){                      p[0]=x; p[1]=y; p[2]=z; return p; }
void  vec3_copy(float *c,const float *a)        { for(int i=0;i<3;i++) c[i]=a[i];       }
void  vec3_abs(float *c,const float *a)         { for(int i=0;i<3;i++) c[i]=fabs(a[i]); }
void  vec3_one(float *c,const float *a)
    {
    float l=divide(1.0,sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])));
    c[0]=a[0]*l;
    c[1]=a[1]*l;
    c[2]=a[2]*l;
    }
void  vec3_len(float *c,const float *a,float l)
    {
    l=divide(l,sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])));
    c[0]=a[0]*l;
    c[1]=a[1]*l;
    c[2]=a[2]*l;
    }
void  vec3_neg(float *c,const float *a)               { for(int i=0;i<3;i++) c[i]=-a[i];     }
void  vec3_add(float *c,const float *a,const float *b){ for(int i=0;i<3;i++) c[i]=a[i]+b[i]; }
void  vec3_sub(float *c,const float *a,const float *b){ for(int i=0;i<3;i++) c[i]=a[i]-b[i]; }
void  vec3_mul(float *c,const float *a,const float *b)
    {
    float   q[3];
    q[0]=(a[1]*b[2])-(a[2]*b[1]);
    q[1]=(a[2]*b[0])-(a[0]*b[2]);
    q[2]=(a[0]*b[1])-(a[1]*b[0]);
    for(int i=0;i<3;i++) c[i]=q[i];
    }
void  vec3_mul(float *c,const float *a,const float  b){ for(int i=0;i<3;i++) c[i]=a[i]*b; }
void  vec3_mul(float *c,const float  a,const float *b){ for(int i=0;i<3;i++) c[i]=a*b[i]; }
float vec3_mul(         const float *a,const float *b){ float c=0; for(int i=0;i<3;i++) c+=a[i]*b[i]; return c; }
//---------------------------------------------------------------------------
void vec3_rotx(float *pos,const float *mid,float ang)
    {
    float y,z,c,s;
    y=pos[1]-mid[1];
    z=pos[2]-mid[2];
    c=cos(ang);
    s=sin(ang);
    pos[1]=+(y*c)+(z*s)+mid[1];
    pos[2]=-(y*s)+(z*c)+mid[2];
    }
//---------------------------------------------------------------------------
void vec3_roty(float *pos,const float *mid,float ang)
    {
    float x,z,c,s;
    x=pos[0]-mid[0];
    z=pos[2]-mid[2];
    c=cos(ang);
    s=sin(ang);
    pos[0]=+(x*c)+(z*s)+mid[0];
    pos[2]=-(x*s)+(z*c)+mid[2];
    }
//---------------------------------------------------------------------------
void vec3_rotz(float *pos,const float *mid,float ang)
    {
    float x,y,c,s;
    x=pos[0]-mid[0];
    y=pos[1]-mid[1];
    c=cos(ang);
    s=sin(ang);
    pos[0]=+(x*c)+(y*s)+mid[0];
    pos[1]=-(x*s)+(y*c)+mid[1];
    }
//---------------------------------------------------------------------------
void mat4_rotx(float *mat,const float *mid,float ang)
    {
    const float c=cos(ang),s=sin(ang);
    const float rx[16]=
        { 1, 0, 0, 0,
          0, c, s, 0,
          0,-s, c, 0,
          0, 0, 0, 1 };
    vec3_sub(mat+12,mat+12,mid);
    mat4_mul(mat,mat,rx);
    vec3_add(mat+12,mat+12,mid);
    }
void mat4_roty(float *mat,const float *mid,float ang)
    {
    const float c=cos(ang),s=sin(ang);
    const float ry[16]=
        { c, 0, s, 0,
          0, 1, 0, 0,
         -s, 0, c, 0,
          0, 0, 0, 1 };
    vec3_sub(mat+12,mat+12,mid);
    mat4_mul(mat,mat,ry);
    vec3_add(mat+12,mat+12,mid);
    }
void mat4_rotz(float *mat,const float *mid,float ang)
    {
    const float c=cos(ang),s=sin(ang);
    const float rz[16]=
        { c, s, 0, 0,
         -s, c, 0, 0,
          0, 0, 1, 0,
          0, 0, 0, 1 };
    vec3_sub(mat+12,mat+12,mid);
    mat4_mul(mat,mat,rz);
    vec3_add(mat+12,mat+12,mid);
    }
//---------------------------------------------------------------------------
float* mat4_ld (         float a0,float a1,float a2,float a3,float a4,float a5,float a6,float a7,float a8,float a9,float a10,float a11,float a12,float a13,float a14,float a15,float a16) { float *p=mat4_tmp; p[0]=a0; p[1]=a1; p[2]=a2; p[3]=a3; p[4]=a4; p[5]=a5; p[6]=a6; p[7]=a7; p[8]=a8; p[9]=a9; p[10]=a10; p[11]=a11; p[12]=a12; p[13]=a13; p[14]=a14; p[15]=a15; return p; }
float* mat4_ld (float *p,float a0,float a1,float a2,float a3,float a4,float a5,float a6,float a7,float a8,float a9,float a10,float a11,float a12,float a13,float a14,float a15,float a16) {                      p[0]=a0; p[1]=a1; p[2]=a2; p[3]=a3; p[4]=a4; p[5]=a5; p[6]=a6; p[7]=a7; p[8]=a8; p[9]=a9; p[10]=a10; p[11]=a11; p[12]=a12; p[13]=a13; p[14]=a14; p[15]=a15; return p; }
void mat4_copy(float *c,const float *a){ for(int i=0;i<16;i++) c[i]=a[i]; }
void mat4_set (float *c,const float  a){ for(int i=0;i<16;i++) c[i]=a; };
void mat4_one (float *c){ for(int i=0;i<16;i++) c[i]=0; c[0]=1; c[5]=1; c[10]=1; c[15]=1; };
void mat4_add (float *c,const float *a,const float *b){ for(int i=0;i<16;i++) c[i]=a[i]+b[i]; }
void mat4_sub (float *c,const float *a,const float *b){ for(int i=0;i<16;i++) c[i]=a[i]-b[i]; }
void mat4_mul (float *c,const float *a,const float *b)
    {
    float q[16];
    q[ 0]=(a[ 0]*b[ 0])+(a[ 1]*b[ 4])+(a[ 2]*b[ 8])+(a[ 3]*b[12]);
    q[ 1]=(a[ 0]*b[ 1])+(a[ 1]*b[ 5])+(a[ 2]*b[ 9])+(a[ 3]*b[13]);
    q[ 2]=(a[ 0]*b[ 2])+(a[ 1]*b[ 6])+(a[ 2]*b[10])+(a[ 3]*b[14]);
    q[ 3]=(a[ 0]*b[ 3])+(a[ 1]*b[ 7])+(a[ 2]*b[11])+(a[ 3]*b[15]);
    q[ 4]=(a[ 4]*b[ 0])+(a[ 5]*b[ 4])+(a[ 6]*b[ 8])+(a[ 7]*b[12]);
    q[ 5]=(a[ 4]*b[ 1])+(a[ 5]*b[ 5])+(a[ 6]*b[ 9])+(a[ 7]*b[13]);
    q[ 6]=(a[ 4]*b[ 2])+(a[ 5]*b[ 6])+(a[ 6]*b[10])+(a[ 7]*b[14]);
    q[ 7]=(a[ 4]*b[ 3])+(a[ 5]*b[ 7])+(a[ 6]*b[11])+(a[ 7]*b[15]);
    q[ 8]=(a[ 8]*b[ 0])+(a[ 9]*b[ 4])+(a[10]*b[ 8])+(a[11]*b[12]);
    q[ 9]=(a[ 8]*b[ 1])+(a[ 9]*b[ 5])+(a[10]*b[ 9])+(a[11]*b[13]);
    q[10]=(a[ 8]*b[ 2])+(a[ 9]*b[ 6])+(a[10]*b[10])+(a[11]*b[14]);
    q[11]=(a[ 8]*b[ 3])+(a[ 9]*b[ 7])+(a[10]*b[11])+(a[11]*b[15]);
    q[12]=(a[12]*b[ 0])+(a[13]*b[ 4])+(a[14]*b[ 8])+(a[15]*b[12]);
    q[13]=(a[12]*b[ 1])+(a[13]*b[ 5])+(a[14]*b[ 9])+(a[15]*b[13]);
    q[14]=(a[12]*b[ 2])+(a[13]*b[ 6])+(a[14]*b[10])+(a[15]*b[14]);
    q[15]=(a[12]*b[ 3])+(a[13]*b[ 7])+(a[14]*b[11])+(a[15]*b[15]);
    for(int i=0;i<16;i++) c[i]=q[i];
    }
void mat4_mul_vec3(float *c,const float *a,const float *b,float w)
    {
    float q[3];
    q[0]=(a[ 0]*b[0])+(a[ 4]*b[1])+(a[ 8]*b[2])+(a[12]*w);
    q[1]=(a[ 1]*b[0])+(a[ 5]*b[1])+(a[ 9]*b[2])+(a[13]*w);
    q[2]=(a[ 2]*b[0])+(a[ 6]*b[1])+(a[10]*b[2])+(a[14]*w);
    for(int i=0;i<3;i++) c[i]=q[i];
    }
void vec3_mul_mat4(float *c,const float *a,const float *b,float w)
    {
    float q[3];
    q[0]=(a[0]*b[ 0])+(a[1]*b[ 1])+(a[2]*b[ 2])+(w*b[ 3]);
    q[1]=(a[0]*b[ 4])+(a[1]*b[ 5])+(a[2]*b[ 6])+(w*b[ 7]);
    q[2]=(a[0]*b[ 8])+(a[1]*b[ 9])+(a[2]*b[10])+(w*b[11]);
    for(int i=0;i<3;i++) c[i]=q[i];
    }
float mat4_mulw_vec3(float *c,const float *a,const float *b,float w)
    {
    float q[4];
    q[0]=(a[ 0]*b[0])+(a[ 4]*b[1])+(a[ 8]*b[2])+(a[12]*w);
    q[1]=(a[ 1]*b[0])+(a[ 5]*b[1])+(a[ 9]*b[2])+(a[13]*w);
    q[2]=(a[ 2]*b[0])+(a[ 6]*b[1])+(a[10]*b[2])+(a[14]*w);
    q[3]=(a[ 3]*b[0])+(a[ 7]*b[1])+(a[11]*b[2])+(a[15]*w);
    for(int i=0;i<3;i++) c[i]=q[i];
    return q[3];
    }
float vec3_mulw_mat4(float *c,const float *a,const float *b,float w)
    {
    float q[4];
    q[0]=(a[0]*b[ 0])+(a[1]*b[ 1])+(a[2]*b[ 2])+(w*b[ 3]);
    q[1]=(a[0]*b[ 4])+(a[1]*b[ 5])+(a[2]*b[ 6])+(w*b[ 7]);
    q[2]=(a[0]*b[ 8])+(a[1]*b[ 9])+(a[2]*b[10])+(w*b[11]);
    q[3]=(a[0]*b[12])+(a[1]*b[13])+(a[2]*b[14])+(w*b[15]);
    for(int i=0;i<3;i++) c[i]=q[i];
    return q[3];
    }
void mat4_subdet(float *c,const float *a)
    {
    float   q[16];
    int     i,j;
    for (i=0;i<4;i++)
     for (j=0;j<4;j++)
      q[j+(i<<2)]=mat4_subdet(a,i,j);
    for (i=0;i<16;i++) c[i]=q[i];
    }
float mat4_subdet(const float *a,int r,int s)
    {
    float   c,q[9];
    int     i,j,k;
    k=0;                            // q = sub matrix
    for (j=0;j<4;j++)
     if (j!=s)
      for (i=0;i<4;i++)
       if (i!=r)
        {
        q[k]=a[i+(j<<2)];
        k++;
        }
    c=0;
    c+=q[0]*q[4]*q[8];
    c+=q[1]*q[5]*q[6];
    c+=q[2]*q[3]*q[7];
    c-=q[0]*q[5]*q[7];
    c-=q[1]*q[3]*q[8];
    c-=q[2]*q[4]*q[6];
    if (int((r+s)&1)) c=-c;       // add signum
    return c;
    }
float mat4_det(const float *a)
    {
    float c=0;
    c+=a[ 0]*mat4_subdet(a,0,0);
    c+=a[ 4]*mat4_subdet(a,0,1);
    c+=a[ 8]*mat4_subdet(a,0,2);
    c+=a[12]*mat4_subdet(a,0,3);
    return c;
    }
float mat4_det(const float *a,const float *b)
    {
    float c=0;
    c+=a[ 0]*b[ 0];
    c+=a[ 4]*b[ 1];
    c+=a[ 8]*b[ 2];
    c+=a[12]*b[ 3];
    return c;
    }
void mat4_inv(float *c,const float *a)
    {
    float   d[16],D;
    mat4_subdet(d,a);
    D=mat4_det(a,d);
    if (D) D=1.0/D;
    for (int i=0;i<16;i++) c[i]=d[i]*D;
    }
void mat4_perspective(float *per,float fovy,float aspect,float zNear,float zFar)
    {
    float f;
    for (int i=0;i<16;i++) per[i]=0.0;
    // original gluProjection
//  f=divide(1.0,tan(0.5*fovy))
//  per[ 0]=f/aspect;
//  per[ 5]=f;
    // corrected gluProjection
    f=divide(1.0,tan(0.5*fovy*aspect));
    per[ 0]=f;
    per[ 5]=f*aspect;
    // z range
    per[10]=divide(zFar+zNear,zNear-zFar);
    per[11]=-1.0;
    per[14]=divide(2.0*zFar*zNear,zNear-zFar);
//  zNear=divide(-per[11],per[10]);     // get znear from perspective projection matrix
    }
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

The lcd class is not included (it would not fit 30K limit and it is not important anyway as its just the low level rasterization) The only important thing here is its xs,ys members which are resolution of the screen and depth buffer is based on 8 bit unsigned integers.

Usage is simple just in timer I call

mesh_view();
mesh_draw(point_1,face3_1);