"""Tier-1: a faithful software model of the Division VelociRender PXPL5 IGC
(Pixel-Planes 5) rasteriser array, driven by geometry captured off the emulated i860.

This is the array model from sda4/DPL3/VRENDER/PXPL5SUP/IGCOPS.C + IGCTYPES.H, made
real: the screen is partitioned into 64x128 tiles (tile_x_bits=6, tile_y_bits=7); each
pixel owns a small bit-addressable memory (pxpl5_mem_chars=26 bytes) plus an enable bit.
The IGC has no per-pixel ALU program in the classic sense -- every pixel evaluates the
same linear tree in lockstep:  eval_ltree(x,y,A,B,C) = (int)(x*A + y*B + C)  (IGCOPS.C).

A triangle is drawn exactly as the hardware does (PXPL5GEO tri_zb_rgb):
  * three EDGE trees  Ei = Ai*x + Bi*y + Ci, oriented inside>=0  -> the enable register
  * a Z tree and R/G/B trees, each an Ax+By+C plane through the 3 vertices
  * per enabled pixel: znew = eval z-tree; if znew nearer than the stored z, the enable
    survives (MEM2geMEM2) and z + rgb are written into pixel memory (writes are gated by
    enable, as in hardware).
Tiles the geometry never touches are skipped (the real board bins per region first).
Finally each tile's pixel memory is read back out (read_pixmem_word) into an RGB frame.

Not a decode of the compiled bit-serial micro-code (that binary encoding is still
undecoded) -- it is the array's *computational* model, producing the pixels that
micro-code would. Validated to match shade_render.py on the same projection.
"""
import sys, struct, math, pickle, json

TILE_X_BITS, TILE_Y_BITS = 6, 7
TILE_X, TILE_Y = 1 << TILE_X_BITS, 1 << TILE_Y_BITS      # 64 x 128
PXPL5_MEM_CHARS = 26

# pixel-memory bit layout we use (fits in 26 bytes = 208 bits):
Z_BIT0, Z_BITS = 0, 24            # 24-bit depth, smaller = nearer
R_BIT0, G_BIT0, B_BIT0, C_BITS = 24, 32, 40, 8
Z_FAR = (1 << Z_BITS) - 1


class Tile:
    """One 64x128 IGC tile: per-pixel enable bit + 26-byte bit memory."""
    __slots__ = ('x0', 'y0', 'mem', 'enable', 'touched')

    def __init__(self, x0, y0):
        self.x0, self.y0 = x0, y0
        self.mem = [bytearray(PXPL5_MEM_CHARS) for _ in range(TILE_X * TILE_Y)]
        self.enable = bytearray(TILE_X * TILE_Y)
        self.touched = False
        for p in self.mem:                         # z-clear to far
            _wword(p, Z_BIT0, Z_BITS, Z_FAR)


def _rbit(pix, bit):
    return 1 & (pix[bit >> 3] >> (bit & 7))


def _wbit(pix, bit, val):
    if val:
        pix[bit >> 3] |= (1 << (bit & 7))
    else:
        pix[bit >> 3] &= ~(1 << (bit & 7))


def _rword(pix, bit0, bits):
    res = 0
    for i in range(bits):
        res |= _rbit(pix, bit0 + i) << i
    return res


def _wword(pix, bit0, bits, val):
    for i in range(bits):
        _wbit(pix, bit0 + i, val & 1)
        val >>= 1


def plane(x0, y0, v0, x1, y1, v1, x2, y2, v2):
    """A,B,C so that A*x+B*y+C == v at each of the 3 vertices (the IGC linear tree)."""
    det = (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0)
    if abs(det) < 1e-9:
        return 0.0, 0.0, v0
    A = ((v1 - v0) * (y2 - y0) - (v2 - v0) * (y1 - y0)) / det
    B = ((v2 - v0) * (x1 - x0) - (v1 - v0) * (x2 - x0)) / det
    C = v0 - A * x0 - B * y0
    return A, B, C


class IGCArray:
    def __init__(self, W, H):
        self.W, self.H = W, H
        self.ntx = (W + TILE_X - 1) // TILE_X
        self.nty = (H + TILE_Y - 1) // TILE_Y
        self.tiles = {}
        self.tris = 0

    def _tile(self, tx, ty):
        key = (tx, ty)
        t = self.tiles.get(key)
        if t is None:
            t = Tile(tx * TILE_X, ty * TILE_Y)
            self.tiles[key] = t
        return t

    def triangle(self, v0, v1, v2):
        """Each v = (sx, sy, depth, r, g, b). depth: smaller = nearer (0..1)."""
        (x0, y0, z0, r0, g0, b0) = v0
        (x1, y1, z1, r1, g1, b1) = v1
        (x2, y2, z2, r2, g2, b2) = v2
        area = (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0)
        if abs(area) < 1e-6:
            return
        if area < 0:                                   # orient CCW so inside>=0
            x1, y1, z1, r1, g1, b1, x2, y2, z2, r2, g2, b2 = \
                x2, y2, z2, r2, g2, b2, x1, y1, z1, r1, g1, b1
        # three edge trees Ei = Ai*x + Bi*y + Ci; normalise each so the centroid is
        # inside (>=0), independent of winding / screen-Y direction.
        gcx, gcy = (x0 + x1 + x2) / 3.0, (y0 + y1 + y2) / 3.0
        E = []
        for (ax, ay), (bx, by) in (((x0, y0), (x1, y1)),
                                   ((x1, y1), (x2, y2)),
                                   ((x2, y2), (x0, y0))):
            A = (by - ay); B = -(bx - ax); C = -(A * ax + B * ay)
            if A * gcx + B * gcy + C < 0:
                A, B, C = -A, -B, -C
            E.append((A, B, C))
        # z + rgb planes
        zc = plane(x0, y0, z0, x1, y1, z1, x2, y2, z2)
        rc = plane(x0, y0, r0, x1, y1, r1, x2, y2, r2)
        gc = plane(x0, y0, g0, x1, y1, g1, x2, y2, g2)
        bc = plane(x0, y0, b0, x1, y1, b1, x2, y2, b2)
        self.tris += 1
        minx = max(0, int(min(x0, x1, x2))); maxx = min(self.W - 1, int(max(x0, x1, x2)) + 1)
        miny = max(0, int(min(y0, y1, y2))); maxy = min(self.H - 1, int(max(y0, y1, y2)) + 1)
        for ty in range(miny // TILE_Y, maxy // TILE_Y + 1):
            for tx in range(minx // TILE_X, maxx // TILE_X + 1):
                self._raster_tile(self._tile(tx, ty), E, zc, rc, gc, bc,
                                  minx, maxx, miny, maxy)

    def _raster_tile(self, t, E, zc, rc, gc, bc, minx, maxx, miny, maxy):
        (zA, zB, zC) = zc
        lx0 = max(0, minx - t.x0); lx1 = min(TILE_X - 1, maxx - t.x0)
        ly0 = max(0, miny - t.y0); ly1 = min(TILE_Y - 1, maxy - t.y0)
        (A0, B0, C0), (A1, B1, C1), (A2, B2, C2) = E
        for ly in range(ly0, ly1 + 1):
            gy = t.y0 + ly
            base = ly << TILE_X_BITS
            for lx in range(lx0, lx1 + 1):
                gx = t.x0 + lx
                # enable = inside all three edge trees
                if (A0 * gx + B0 * gy + C0) < 0: continue
                if (A1 * gx + B1 * gy + C1) < 0: continue
                if (A2 * gx + B2 * gy + C2) < 0: continue
                idx = base + lx
                pix = t.mem[idx]
                znew = int((zA * gx + zB * gy + zC) * Z_FAR)
                if znew < 0: znew = 0
                elif znew > Z_FAR: znew = Z_FAR
                if znew >= _rword(pix, Z_BIT0, Z_BITS):     # MEM2geMEM2: nearer wins
                    continue
                t.touched = True
                _wword(pix, Z_BIT0, Z_BITS, znew)
                rv = min(255, max(0, int(rc[0] * gx + rc[1] * gy + rc[2])))
                gv = min(255, max(0, int(gc[0] * gx + gc[1] * gy + gc[2])))
                bv = min(255, max(0, int(bc[0] * gx + bc[1] * gy + bc[2])))
                _wword(pix, R_BIT0, C_BITS, rv)
                _wword(pix, G_BIT0, C_BITS, gv)
                _wword(pix, B_BIT0, C_BITS, bv)

    def readout_depth(self, bg=(7, 11, 17)):
        """Read the 24-bit z stored in pixel memory back into a depth image (near=bright).
        This is the plane the SENDE z-sweep interpolates, straight out of pixel memory."""
        img = bytearray(self.W * self.H * 3)
        for y in range(self.H):
            o = y * self.W * 3
            for x in range(self.W):
                img[o] = bg[0]; img[o + 1] = bg[1]; img[o + 2] = bg[2]; o += 3
        # gather the touched z range for contrast
        zmin, zmax = Z_FAR, 0
        for t in self.tiles.values():
            if not t.touched:
                continue
            for pix in t.mem:
                z = _rword(pix, Z_BIT0, Z_BITS)
                if z != Z_FAR:
                    zmin = min(zmin, z); zmax = max(zmax, z)
        span = (zmax - zmin) or 1
        for (tx, ty), t in self.tiles.items():
            if not t.touched:
                continue
            for ly in range(TILE_Y):
                gy = t.y0 + ly
                if gy >= self.H:
                    break
                base = ly << TILE_X_BITS
                for lx in range(TILE_X):
                    gx = t.x0 + lx
                    if gx >= self.W:
                        break
                    z = _rword(t.mem[base + lx], Z_BIT0, Z_BITS)
                    if z == Z_FAR:
                        continue
                    v = 1.0 - (z - zmin) / span          # near = bright
                    o = (gy * self.W + gx) * 3
                    img[o] = int(30 + v * 120); img[o + 1] = int(90 + v * 150); img[o + 2] = int(120 + v * 130)
        return img

    def readout(self, bg=(7, 11, 17)):
        """Read pixel memory back out into an RGB framebuffer (row-major, RGB bytes)."""
        img = bytearray(self.W * self.H * 3)
        for y in range(self.H):
            o = y * self.W * 3
            for x in range(self.W):
                img[o] = bg[0]; img[o + 1] = bg[1]; img[o + 2] = bg[2]
                o += 3
        for (tx, ty), t in self.tiles.items():
            if not t.touched:
                continue
            for ly in range(TILE_Y):
                gy = t.y0 + ly
                if gy >= self.H: break
                base = ly << TILE_X_BITS
                for lx in range(TILE_X):
                    gx = t.x0 + lx
                    if gx >= self.W: break
                    pix = t.mem[base + lx]
                    if _rword(pix, Z_BIT0, Z_BITS) == Z_FAR:
                        continue
                    o = (gy * self.W + gx) * 3
                    img[o] = _rword(pix, R_BIT0, C_BITS)
                    img[o + 1] = _rword(pix, G_BIT0, C_BITS)
                    img[o + 2] = _rword(pix, B_BIT0, C_BITS)
        return img


# ------- driver: rasterise the captured 9x5 surface through the array -------
def _n(a, b, c):
    m = math.sqrt(a * a + b * b + c * c) or 1.0
    return a / m, b / m, c / m


def build_from_grid(pkl, W=620, H=560, yaw=40.0, pitch=28.0):
    objs = pickle.load(open(pkl, 'rb'))['objs']
    allv = [v for o in objs for v in o]
    xs = sorted(set(round(v['mx'], 2) for v in allv))
    zs = sorted(set(round(v['mz'], 2) for v in allv))
    grid = {(round(v['mx'], 2), round(v['mz'], 2)): v for v in allv}
    cx = sum(xs) / len(xs); cz = sum(zs) / len(zs)
    cy = sum(v['my'] for v in allv) / len(allv)
    ry, rp = math.radians(yaw), math.radians(pitch)
    cyw, syw, cp, sp = math.cos(ry), math.sin(ry), math.cos(rp), math.sin(rp)

    def rot(x, y, z):
        x, z = x * cyw + z * syw, -x * syw + z * cyw
        y, z = y * cp - z * sp, y * sp + z * cp
        return x, y, z
    L = _n(0.35, 0.55, 0.72)
    P = {}
    for (x, z), v in grid.items():
        X, Y, Z = rot(v['mx'] - cx, (v['my'] - cy) * 1.8, v['mz'] - cz)
        nx, ny, nz = rot(v['nx'], v['ny'], v['nz']); n = _n(nx, ny, nz)
        d = abs(n[0] * L[0] + n[1] * L[1] + n[2] * L[2])
        it = max(0.16, min(1.0, 0.22 + 0.85 * d))
        P[(x, z)] = (X, Y, Z, it)
    XX = [p[0] for p in P.values()]; YY = [p[1] for p in P.values()]; ZZ = [p[2] for p in P.values()]
    mnx, mxx, mny, mxy = min(XX), max(XX), min(YY), max(YY)
    mnz, mxz = min(ZZ), max(ZZ)
    pad = 0.1 * W
    s = min((W - 2 * pad) / (mxx - mnx), (H - 2 * pad) / (mxy - mny))
    ox = (W - (mxx - mnx) * s) / 2; oy = (H - (mxy - mny) * s) / 2

    def sx(p): return (p[0] - mnx) * s + ox
    def sy(p): return (mxy - p[1]) * s + oy
    def depth(p): return (p[2] - mnz) / ((mxz - mnz) or 1)      # 0 near .. 1 far

    def col(it):
        return (min(255, 28 + it * 168), min(255, 50 + it * 205), min(255, 54 + it * 122))

    arr = IGCArray(W, H)
    xsl, zsl = xs, zs
    for i in range(len(xsl) - 1):
        for j in range(len(zsl) - 1):
            a = P[(xsl[i], zsl[j])]; b = P[(xsl[i + 1], zsl[j])]
            c = P[(xsl[i], zsl[j + 1])]; d = P[(xsl[i + 1], zsl[j + 1])]

            def vtx(p):
                r, g, bb = col(p[3]); return (sx(p), sy(p), depth(p), r, g, bb)
            arr.triangle(vtx(a), vtx(b), vtx(c))
            arr.triangle(vtx(b), vtx(d), vtx(c))
    return arr


def main():
    S = r'C:\Users\cyd\AppData\Local\Temp\claude\c--VWE-TeslaRel410\4e848c76-6e89-4034-8047-d8d491cb32d8\scratchpad'
    pkl = sys.argv[1] if len(sys.argv) > 1 else S + r'\vfull.pkl'
    out = next((a for a in sys.argv[1:] if a.endswith('.ppm')), 'igc.ppm')
    W, H = 620, 560
    arr = build_from_grid(pkl, W, H)
    active = sum(1 for t in arr.tiles.values() if t.touched)
    print("IGC array: %dx%d, %d tiles (%dx%d), %d touched, %d triangles" %
          (W, H, len(arr.tiles), arr.ntx, arr.nty, active, arr.tris))
    img = arr.readout()
    with open(out, 'wb') as f:
        f.write(b'P6\n%d %d\n255\n' % (W, H)); f.write(bytes(img))
    print("wrote", out)
    # ASCII preview
    ramp = " .:-=+*#%@"
    for y in range(0, H, H // 40):
        line = "  "
        for x in range(0, W, W // 74):
            o = (y * W + x) * 3
            lum = (img[o] + img[o + 1] + img[o + 2]) / 3
            line += ramp[min(9, int(lum / 255 * 9.99))] if lum > 20 else ' '
        print(line)


if __name__ == '__main__':
    main()