127 lines
4.2 KiB
Python
127 lines
4.2 KiB
Python
import struct
|
|
|
|
d = r"C:\Project\MAG160C\analysis\pairs_verify_20260813"
|
|
npix = 19200
|
|
|
|
def load(path, fmt, n):
|
|
b = open(path, "rb").read()
|
|
return struct.unpack_from("<%d%s" % (n, fmt), b, 0)
|
|
|
|
# 1) verify interpolation: thrw vs interp(thr0, thr1, t)
|
|
thr0 = load(d + r"\pair_000.thr0", "h", 19200*2 + 512) # (nsegs-1)*npix + 1024/2
|
|
thr1 = load(d + r"\pair_000.thr1", "h", 19200*2 + 512)
|
|
thrw = load(d + r"\pair_000.thrw", "h", 19200*2)
|
|
T = [8304, 18390, 28495, 33637, 38730, 47664]
|
|
dev54 = 29289
|
|
sel = 2
|
|
t = ((dev54 - T[sel]) << 12) // (T[sel+1] - T[sel])
|
|
print("t =", t)
|
|
diff = 0
|
|
for i in range(19200*2):
|
|
a = thr0[i]
|
|
b = thr1[i]
|
|
w = thrw[i]
|
|
calc = a + ((b - a) * t >> 12)
|
|
# clamp to int16 (out-of-range branch)
|
|
if calc > 32767: calc = 32767
|
|
if calc < -32768: calc = -32768
|
|
calc16 = struct.unpack("<h", struct.pack("<H", calc & 0xFFFF))[0]
|
|
if calc16 != w:
|
|
diff += 1
|
|
if diff <= 5:
|
|
print(f" thr mismatch i={i} a={a} b={b} calc={calc16} w={w}")
|
|
print(f"thr interp: diff={diff}/38400")
|
|
|
|
g0 = load(d + r"\pair_000.g0", "H", 19200*3*2)
|
|
g1 = load(d + r"\pair_000.g1", "H", 19200*3*2)
|
|
gw = load(d + r"\pair_000.gw", "H", 19200*3*2)
|
|
diff = 0
|
|
for i in range(19200*3*2):
|
|
a = g0[i]; b = g1[i]; w = gw[i]
|
|
calc = a + ((b - a) * t >> 12)
|
|
if calc > 65535: calc = 65535
|
|
if calc < 0: calc = 0
|
|
if calc != w:
|
|
diff += 1
|
|
if diff <= 5:
|
|
print(f" gain mismatch i={i} a={a} b={b} calc={calc} w={w}")
|
|
print(f"gain interp: diff={diff}/115200")
|
|
|
|
# 2) NUC verify with working tables + blind
|
|
f20 = load(d + r"\pair_000.f20", "H", npix)
|
|
ref = load(d + r"\pair_000.ref", "H", npix)
|
|
raw = load(d + r"\pair_000.raw", "H", npix)
|
|
blind = open(d + r"\pair_000.blind", "rb").read()
|
|
recs = []
|
|
for i in range(len(blind)//40):
|
|
target, typ = struct.unpack_from("<II", blind, i*40)
|
|
neigh = struct.unpack_from("<8I", blind, i*40+8)
|
|
recs.append((target, typ, neigh))
|
|
nuc = [0]*npix
|
|
for i in range(npix):
|
|
d2 = (int(f20[i]) - int(ref[i])) >> 1
|
|
seg = 0
|
|
if d2 > thrw[i*2]: seg = 1
|
|
if seg == 1 and d2 > thrw[i*2+1]: seg = 2
|
|
p = (seg*npix + i)*2
|
|
v = gw[p+1] + ((gw[p] * d2) >> 12)
|
|
if v < 0: v = 0
|
|
if v > 65535: v = 65535
|
|
nuc[i] = v
|
|
out = list(nuc)
|
|
for (tgt, typ, ng) in recs:
|
|
if typ == 8:
|
|
out[tgt] = sum(out[ng[j]] for j in range(8)) >> 3
|
|
elif typ == 7:
|
|
out[tgt] = sum(out[ng[j]] for j in range(7)) // 7
|
|
elif typ == 6:
|
|
out[tgt] = sum(out[ng[j]] for j in range(6)) // 6
|
|
elif typ == 5:
|
|
out[tgt] = sum(out[ng[j]] for j in range(5)) // 5
|
|
elif typ == 4:
|
|
out[tgt] = sum(out[ng[j]] for j in range(4)) >> 2
|
|
elif typ == 3:
|
|
out[tgt] = sum(out[ng[j]] for j in range(3)) // 3
|
|
diff = sum(1 for i in range(npix) if out[i] != raw[i])
|
|
sad = sum(abs(out[i]-raw[i]) for i in range(npix))
|
|
print(f"NUC+blind vs official raw: diff={diff}/19200 MAE={sad/npix:.4f}")
|
|
|
|
# 3) window + gray: LUT1024 + idx
|
|
win = open(d + r"\pair_000.win", "rb").read()
|
|
hi, lo, fmax, fmin, mean, std, axmax, axmin = struct.unpack_from("<8I", win, 0)
|
|
print(f"win: hi={hi} lo={lo} fmax={fmax} fmin={fmin} mean={mean} std={std}")
|
|
lut = open(d + r"\pair_000.lut", "rb").read()
|
|
S = 0xFFC00000 // (hi - lo)
|
|
gray = []
|
|
for i in range(npix):
|
|
v = out[i]
|
|
if v <= lo: idx = 0
|
|
elif v >= hi: idx = 1023
|
|
else: idx = (v - lo) * S >> 22
|
|
if idx > 1023: idx = 1023
|
|
gray.append(lut[idx])
|
|
# compare with official gray (320x240) - downsample center? official gray is 2x
|
|
g320 = load(d + r"\pair_000.gray", "B", 320*240)
|
|
# check the 2x relation: gray320[2y][2x] should equal gray160[y][x] (or the interp)
|
|
match = 0
|
|
tot = 0
|
|
for y in range(0, 120):
|
|
for x in range(0, 160):
|
|
g2 = g320[y*2*320 + x*2]
|
|
if g2 == gray[y*160 + x]: match += 1
|
|
tot += 1
|
|
print(f"gray 2x top-left parity: {match}/{tot}")
|
|
# what mapping? try LUT1024 idx directly vs gray320[0..] offset scan
|
|
best = None
|
|
for dy in range(4):
|
|
for dx in range(4):
|
|
m = 0
|
|
for y in range(0, 120):
|
|
for x in range(0, 160):
|
|
yy = y*2+dy; xx = x*2+dx
|
|
if yy < 240 and xx < 320:
|
|
if g320[yy*320+xx] == gray[y*160+x]: m += 1
|
|
if best is None or m > best[0]:
|
|
best = (m, dy, dx)
|
|
print("best offset:", best)
|