#!/usr/bin/env python3
"""Forensics: is the blocky noise grid-aligned at the latent-patch scale (~16px)?

Method: temporal-difference map -> 2D FFT -> check power concentration at
patch-grid frequencies (W/76=16px, H/42=16px) vs neighboring bins,
plus high-frequency band energy. Static shot = control.
"""
import subprocess, sys, math
import numpy as np

def read_frames_gray(path, start, count, w, h):
    # decode whole stream (short clips), slice in numpy — avoids select/trim quirks
    cmd = ["ffmpeg","-v","error","-i",path,"-vf",f"format=gray","-f","rawvideo","-"]
    raw = subprocess.run(cmd, capture_output=True).stdout
    total = len(raw)//(w*h)
    fr = np.frombuffer(raw[:total*w*h], np.uint8).reshape(total,h,w).astype(np.float32)
    return fr[start:start+count]

def analyze(path, w, h, start, count, label):
    fr = read_frames_gray(path, start, count, w, h)
    if len(fr) < 3:
        print(f"{label}: NOT ENOUGH FRAMES"); return
    D = np.abs(np.diff(fr, axis=0)).mean(axis=0)          # avg temporal-diff map
    D = D - D.mean()
    win = np.outer(np.hanning(h), np.hanning(w))
    F = np.fft.fftshift(np.fft.fft2(D*win))
    P = np.abs(F)**2
    def binpow(px, py):                                   # px,py = cycles across W,H
        ix, iy = int(round(px)), int(round(py))
        cy, cx = h//2, w//2
        return P[cy+iy, cx+ix]
    # patch-grid fundamental: 76x42 patches => 76 cycles/W, 42 cycles/H
    fund = max(binpow(76,42), binpow(-76,42), binpow(76,-42), binpow(-76,-42))
    # background: mean of ring at same radius +-2 bins excluding near-axis
    r = math.hypot(76/w*2, 42/h*2)  # normalized freq radius (cycles/px *2)
    tot=0.0; cnt=0
    for dy in range(-4,5):
        for dx in range(-4,5):
            if dx==0 and dy==0: continue
            fx, fy = 76+dx*3, 42+dy*3   # sample ring around fundamental, off-grid
            v = binpow(fx,fy); tot+=v; cnt+=1
    bg = tot/max(cnt,1)
    # high-freq band energy fraction (>40 cycles across)
    cy,cx = h//2,w//2
    Y,X = np.mgrid[0:h,0:w]
    R = np.hypot((X-cx)/w,(Y-cy)/h)*2   # ~cycles/px*2 approx scale
    hf = P[R>0.05].sum(); allp = P.sum()
    print(f"{label}: frames={len(fr)}")
    print(f"  temporal-diff mean={D.mean()+abs(D.mean()):.2f} std={D.std():.2f}")
    print(f"  patch-grid peak(16px) power = {fund:.3e}")
    print(f"  local background power      = {bg:.3e}")
    print(f"  contrast ratio fund/bg      = {fund/max(bg,1e-12):.2f}x")
    print(f"  HF(>0.05cyc/px) share       = {hf/max(allp,1e-12)*100:.1f}%")

if __name__ == "__main__":
    base = "/home/kai/shared/2026-08-26/blockcache-prod/"
    analyze(base+"bcHD36_A_00001_.mp4", 1216, 672, 70, 10, "HD36-A(native+sage) hi-motion f70")
    analyze(base+"bcHD36_B_00001_.mp4", 1216, 672, 70, 10, "HD36-B(+blockcache)  hi-motion f70")
    analyze(base+"arm_A1.mp4", 768, 1344, 30, 10, "STATIC-A(native+sage) mid f30")