Implementing integer 1/z depth buffer

I am trying to adapt the code from this answer to use 1/z depth buffer and here is my attempt

/* g++ trig.cpp -o trig -lSDL2 */

#include <algorithm>

#define SDL_MAIN_HANDLED
#include <SDL2/SDL.h>

#define SCREEN_WIDTH  600
#define SCREEN_HEIGHT 400

// TODO make an array of values to be interpolated
// TODO uv texture mapping
typedef struct
{
    int x, y, z;
    int zInv; // 1/z
    uint8_t r, g, b;
    int u, v;
} Point2D;

constexpr int SIZE = SCREEN_WIDTH * SCREEN_HEIGHT;

uint32_t pixels[SIZE];
int zDepth[SIZE];

void SetPixel(int x, int y, int zInv, uint8_t r, uint8_t g, uint8_t b)
{
    int offset = y * SCREEN_WIDTH + x;

    if (zDepth[offset] < zInv)
    {
        zDepth[offset] = zInv;
        pixels[offset] = (255 << 24) | (r << 16) | (g << 8) | b;
    }
}

int ContourX[SCREEN_HEIGHT][2]; // min X and max X for every horizontal line within the triangle
int ContourZ[SCREEN_HEIGHT][2]; // 1/z
int ContourR[SCREEN_HEIGHT][2]; // red value for every horizontal line
int ContourG[SCREEN_HEIGHT][2]; // green value for every horizontal line
int ContourB[SCREEN_HEIGHT][2]; // blue value for every horizontal line

// Scans a side of a triangle setting min X and max X in ContourX[][] (same for ContourR, ContourG, and ContourB) using the Bresenham's
void TriangleLine(const Point2D& p0, const Point2D& p1)
{
    // DDA variables
    int kx, ky, kz, kr, kg, kb; // step directions
    int dx, dy, dz, dr, dg, db; // abs delta
    kx = 0; dx = p1.x - p0.x; if (dx > 0) kx = 1; if (dx < 0) { kx = -1; dx = -dx; }
    ky = 0; dy = p1.y - p0.y; if (dy > 0) ky = 1; if (dy < 0) { ky = -1; dy = -dy; }
    kz = 0; dz = p1.zInv - p0.zInv; if (dz > 0) kz = 1; if (dz < 0) { kz = -1; dz = -dz; }
    kr = 0; dr = p1.r - p0.r; if (dr > 0) kr = 1; if (dr < 0) { kr = -1; dr = -dr; }
    kg = 0; dg = p1.g - p0.g; if (dg > 0) kg = 1; if (dg < 0) { kg = -1; dg = -dg; }
    kb = 0; db = p1.b - p0.b; if (db > 0) kb = 1; if (db < 0) { kb = -1; db = -db; }

    int n;
    n = dx;
    if (n < dy) n = dy;
    if (n < dz) n = dz;
    if (n < dr) n = dr;
    if (n < dg) n = dg;
    if (n < db) n = db;
    if (!n) n = 1;

    // target buffer according to ky direction
    int target; if (ky > 0) target = 0; else target = 1;

    // integer DDA line start point
    int x = p0.x;
    int y = p0.y;
    int z = p0.zInv;
    int r = p0.r;
    int g = p0.g;
    int b = p0.b;

    // fix endpoints just to be sure (wrong division constants by +/-1 can cause that last point is missing)
    ContourX[p0.y][target] = p0.x;
    ContourZ[p0.y][target] = p0.zInv;
    ContourR[p0.y][target] = p0.r;
    ContourG[p0.y][target] = p0.g;
    ContourB[p0.y][target] = p0.b;

    ContourX[p1.y][target] = p1.x;
    ContourZ[p1.y][target] = p1.zInv;
    ContourR[p1.y][target] = p1.r;
    ContourG[p1.y][target] = p1.g;
    ContourB[p1.y][target] = p1.b;

    int cx, cy, cz, cr, cg, cb, i;
    for (cx = cy = cz = cr = cg = cb = n, i = 0; i < n; ++i)
    {
        ContourX[y][target] = x;
        ContourZ[y][target] = z;
        ContourR[y][target] = r;
        ContourG[y][target] = g;
        ContourB[y][target] = b;

        cx -= dx; if (cx <= 0){ cx += n; x += kx; }
        cy -= dy; if (cy <= 0){ cy += n; y += ky; }
        cz -= dz; if (cz <= 0){ cz += n; z += kz; }
        cr -= dr; if (cr <= 0){ cr += n; r += kr; }
        cg -= dg; if (cg <= 0){ cg += n; g += kg; }
        cb -= db; if (cb <= 0){ cb += n; b += kb; }
    }
}

void DrawTriangle(const Point2D& p0, const Point2D& p1, const Point2D& p2)
{
    TriangleLine(p0, p1);
    TriangleLine(p1, p2);
    TriangleLine(p2, p0);

    int y0, y1; // min and max y
    y0 = p0.y; if (y0 > p1.y) y0 = p1.y; if (y0 > p2.y) y0 = p2.y;
    y1 = p0.y; if (y1 < p1.y) y1 = p1.y; if (y1 < p2.y) y1 = p2.y;

    int x0, z0, r0, g0, b0;
    int x1, z1, r1, g1, b1;
    int dx;
    int kz, kr, kg, kb;
    int dz, dr, dg, db;
    int z, cz;
    int r, cr;
    int g, cg;
    int b, cb;

    for (int y = y0; y <= y1; ++y)
    {
        if (ContourX[y][0] < ContourX[y][1])
        {
            x0 = ContourX[y][0];
            z0 = ContourZ[y][0];
            r0 = ContourR[y][0];
            g0 = ContourG[y][0];
            b0 = ContourB[y][0];

            x1 = ContourX[y][1];
            z1 = ContourZ[y][1];
            r1 = ContourR[y][1];
            g1 = ContourG[y][1];
            b1 = ContourB[y][1];
        }
        else
        {
            x1 = ContourX[y][0];
            z1 = ContourZ[y][0];
            r1 = ContourR[y][0];
            g1 = ContourG[y][0];
            b1 = ContourB[y][0];

            x0 = ContourX[y][1];
            z0 = ContourZ[y][1];
            r0 = ContourR[y][1];
            g0 = ContourG[y][1];
            b0 = ContourB[y][1];
        }

        dx = x1 - x0;

        kz = 0; dz = z1 - z0; if (dz > 0) kz = 1; if (dz < 0) { kz = -1; dz = -dz; }
        kr = 0; dr = r1 - r0; if (dr > 0) kr = 1; if (dr < 0) { kr = -1; dr = -dr; }
        kg = 0; dg = g1 - g0; if (dg > 0) kg = 1; if (dg < 0) { kg = -1; dg = -dg; }
        kb = 0; db = b1 - b0; if (db > 0) kb = 1; if (db < 0) { kb = -1; db = -db; }

        z = z0; cz = dx;
        r = r0; cr = dx;
        g = g0; cg = dx;
        b = b0; cb = dx;

        // x<x1 to follow top left rule (ie. don't draw bottom or right edges)
        for (int x = x0; x < x1; ++x)
        {
            SetPixel(x, y, z, r, g, b);

            cz -= dz; if (cz <= 0) { cz += dx; z += kz; }
            cr -= dr; if (cr <= 0) { cr += dx; r += kr; }
            cg -= dg; if (cg <= 0) { cg += dx; g += kg; }
            cb -= db; if (cb <= 0) { cb += dx; b += kb; }
        }
    }
}

int main(void)
{
    // clear the screen
    std::fill(pixels, pixels + SIZE, 0);
    std::fill(zDepth, zDepth + SIZE, 0);
/*
    Point2D p0, p1, p2, p3;

    p0.x = 30;
    p0.y = 41;
    p0.r = 255;
    p0.g = 0;
    p0.b = 0;

    p1.x = 350;
    p1.y = 41;
    p1.r = 0;
    p1.g = 255;
    p1.b = 0;

    p2.x = 40;
    p2.y = 311;
    p2.r = 0;
    p2.g = 0;
    p2.b = 255;

    p3.x = 572;
    p3.y = 280;
    p3.r = 255;
    p3.g = 140;
    p3.b = 0;

    DrawTriangle(p0, p1, p2);
    DrawTriangle(p1, p2, p3);
*/

    Point2D p0, p1, p2, p3, p4, p5;

    p0.x = 10;
    p0.y = 50;
    p0.z = 10;
    p0.zInv = 0xfffff / p0.z;
    p0.r = 255;
    p0.g = 0;
    p0.b = 0;

    p1.x = 400;
    p1.y = 100;
    p1.z = 10;
    p1.zInv = 0xfffff / p1.z;
    p1.r = 255;
    p1.g = 0;
    p1.b = 0;

    p2.x = 290;
    p2.y = 380;
    p2.z = 10;
    p2.zInv = 0xfffff / p2.z;
    p2.r = 255;
    p2.g = 0;
    p2.b = 0;

    DrawTriangle(p0, p1, p2);

    p3.x = 50;
    p3.y = 350;
    p3.z = 2;
    p3.zInv = 0xfffff / p3.z;
    p3.r = 0;
    p3.g = 255;
    p3.b = 0;

    p4.x = 130;
    p4.y = 40;
    p4.z = 20;
    p4.zInv = 0xfffff / p4.z;
    p4.r = 0;
    p4.g = 255;
    p4.b = 0;

    p5.x = 380;
    p5.y = 200;
    p5.z = 5;
    p5.zInv = 0xfffff / p5.z;
    p5.r = 0;
    p5.g = 255;
    p5.b = 0;

    DrawTriangle(p3, p4, p5);

    SDL_Init(SDL_INIT_EVERYTHING);

    SDL_Window* window = SDL_CreateWindow
    (
        "Trig",
        SDL_WINDOWPOS_UNDEFINED,
        SDL_WINDOWPOS_UNDEFINED,
        SCREEN_WIDTH,
        SCREEN_HEIGHT,
        SDL_WINDOW_SHOWN
    );

    SDL_Renderer* renderer = SDL_CreateRenderer(window, -1, SDL_RENDERER_ACCELERATED);
    SDL_Texture* texture = SDL_CreateTexture(renderer, SDL_PIXELFORMAT_BGRA32, SDL_TEXTUREACCESS_STREAMING, SCREEN_WIDTH, SCREEN_HEIGHT);
    SDL_Event event;

    bool quit = false;

    while (!quit)
    {
        if (SDL_PollEvent(&event))
        {
            switch (event.type)
            {
            case SDL_QUIT:
            {
                quit = true;
                break;
            }
            }
        }

        SDL_UpdateTexture(texture, NULL, &pixels[0], SCREEN_WIDTH * 4);
        SDL_RenderCopy(renderer, texture, NULL, NULL);
        SDL_RenderPresent(renderer);
    }

    SDL_DestroyRenderer(renderer);
    SDL_DestroyWindow(window);
    SDL_Quit();

    return 0;
}

Unfortunately I didn't get the output I desired

it should look like this (I used floating point arithmetic with naive z-buffer implementation)

For reference, here is the output by OpenGL:

I am not sure what might be the problem. Maybe I need to decrease the step size for 1/z? Any pointers?

first using zInv=0xFFFF/z on 16bit int might be a problem as you forgot the sign is there too. I would use zInv=0x7FFF/z just to be safe. In case int is 32 or 64 bit on your platform/compiler then you're fine.

Other than that from a quick look at your code I do not see anything wrong with it. So it might be just problem with precision as 1/z is highly nonlinear.

To test it you can use 32bit ints for the zInv like zInv=0x7FFFFFFF/z if the behavior changes (for the better) you know its precision. Beware you have to make sure all the variables along the interpolation of zInv is 32bit too.

Another option to test is change the zInv to float if that helps a lot you know its precision (IIRC OpenGL also uses floats for this).

To help with precsion you might offset the z a bit to move from problematic part of 1/z (where |z| is small) for example:

zInv=0x7FFF/(z+z0)

so you do not divide by too small values... try for example z0=1000 if your view direction is +z ...