import numpy as np
from PIL import Image
from scipy.fftpack import dct

img = Image.open('/Users/user/puz/README.jpg')
Y = np.asarray(img.convert('YCbCr'))[:, :, 0].astype(float) - 128.0
h, w = Y.shape
bh, bw = h // 8, w // 8
print('image', w, 'x', h, '->', bw, 'x', bh, 'blocks =', bw * bh)

# zigzag order
zz = []
for s in range(15):
    for i in range(8):
        j = s - i
        if 0 <= j < 8:
            zz.append((i, j) if s % 2 else (j, i))

def blockdct(b):
    return dct(dct(b.T, norm='ortho').T, norm='ortho')

coeffs = np.zeros((bh, bw, 8, 8))
for by in range(bh):
    for bx in range(bw):
        coeffs[by, bx] = blockdct(Y[by*8:by*8+8, bx*8:bx*8+8])

def decode(idxs, invert, msb_first=True):
    out = bytearray()
    for by in range(bh):
        for bx in range(bw):
            bits = []
            for (i, j) in idxs:
                b = 1 if coeffs[by, bx, i, j] > 0 else 0
                bits.append(b ^ invert)
            if not msb_first:
                bits = bits[::-1]
            v = 0
            for b in bits:
                v = (v << 1) | b
            out.append(v)
    return bytes(out)

# candidate: 8 lowest-frequency AC coefficients in zigzag order (skip DC)
for name, idxs in [('zigzag AC1-8', zz[1:9]), ('zigzag inc DC', zz[0:8])]:
    for invert in (0, 1):
        for msb in (True, False):
            data = decode(idxs, invert, msb)
            head = data[:120]
            printable = sum(1 for c in head if 32 <= c < 127 or c in (10, 13))
            print(f'{name} invert={invert} msb={msb} printable={printable}/120 :: {head[:96]!r}')
