import json, math, hashlib
from pathlib import Path
import numpy as np
from PIL import Image, ImageFilter
BASE=Path('.scratch/cockpit-shell-panel-layout/verify/artifacts/t35')
PILOT=BASE/'negative-estate-pilot/input-v2'
OUT=BASE/'roof-imagery-quality'
OUT.mkdir(parents=True,exist_ok=True)
configs=json.loads((PILOT/'roof-textures.json').read_text())
regions=json.loads((BASE/'roof-quality-pilot/report.json').read_text())

def pixels(path):
    return np.asarray(Image.open(path).convert('RGB'),dtype=np.float32)
def stats(arr):
    rgb=arr/255.0
    lum=.2126*arr[:,:,0]+.7152*arr[:,:,1]+.0722*arr[:,:,2]
    sat=arr.max(2)-arr.min(2)
    dx=np.abs(np.diff(lum,axis=1));dy=np.abs(np.diff(lum,axis=0))
    return {'sizePx':[int(arr.shape[1]),int(arr.shape[0]),'source channels RGB'],
      'luminancePercentiles':{str(p):round(float(np.percentile(lum,p)),2) for p in [1,5,25,50,75,95,99]},
      'darkPixelFractionLumaUnder55':round(float(np.mean(lum<55)),5),
      'brightPixelFractionLumaOver210':round(float(np.mean(lum>210)),5),
      'saturationPercentiles':{str(p):round(float(np.percentile(sat,p)),2) for p in [50,90,99]},
      'edgeGradientP95':round(float(np.percentile(np.concatenate([dx.ravel(),dy.ravel()]),95)),2),
      'rgbStddev':np.round(np.std(arr.reshape(-1,3),axis=0),2).tolist()}

textures=[]
tile_map={}
for bid,meta in configs.items():
    arr=pixels(meta['path']);rec={'id':bid,'roofCrop':stats(arr),'sha256':hashlib.sha256(Path(meta['path']).read_bytes()).hexdigest(),
      'sourceTileCount':len(meta['sourceTiles']),'sourceTiles':meta['sourceTiles']}
    textures.append(rec)
    for t in meta['sourceTiles']:tile_map[t['path']]=t['sha256']
for target in regions['results']:
    for src in target['processed']['sourceTiles']:tile_map[src['path']]=src['sha256']

# Tile-edge continuity proxy. Compare RGB values immediately beside adjacent source tile boundaries.
tiles={p:pixels(p) for p in tile_map}
seams=[]
items=sorted(tiles)
for p in items:
    parts=Path(p).parts
    try:z=int(parts[-3]);x=int(parts[-2]);tms=int(Path(p).stem)
    except (ValueError,IndexError):continue
    y=2**z-1-tms
    for axis,nx,ny in [('east',x+1,y),('south',x,y+1)]:
      q=f'D:/GIS/GIS/yuanlingPic3/{z}/{nx}/{2**z-1-ny}.png'
      if q not in tiles or p>=q:continue
      a,b=tiles[p],tiles[q]
      edgeA=a[:,-1,:] if axis=='east' else a[-1,:,:]
      edgeB=b[:,0,:] if axis=='east' else b[0,:,:]
      d=np.sqrt(np.mean((edgeA-edgeB)**2,axis=1))
      la=.2126*edgeA[:,0]+.7152*edgeA[:,1]+.0722*edgeA[:,2]
      lb=.2126*edgeB[:,0]+.7152*edgeB[:,1]+.0722*edgeB[:,2]
      seams.append({'tileA':p,'tileB':q,'orientation':axis,'meanRgbRmsDelta':round(float(np.mean(d)),2),
        'p95RgbRmsDelta':round(float(np.percentile(d,95)),2),'meanLumaDelta':round(float(np.mean(np.abs(la-lb))),2)})

report={'generatedAt':'2026-10-15','method':'Pillow/Numpy raster statistics only; not roof segmentation, radiometric correction or facade evidence',
 'roofCrops':textures,'sourceTileCount':len(tiles),'adjacentTileSeams':seams,
 'interpretationLimits':['Dark-pixel share and luminance spread combine cast shadows, roof material, capture mosaics and building occlusion.','Adjacent tile edge mismatch is a proxy and may cross different land cover; it is not a seam detector.','The 512/205px GLB atlases are roof bounding rectangles; only the roof mesh UV subset is rendered.','No source image is altered; no geometric claim follows from these statistics.'],
 'runtimeConfigChanged':False,'runtimeRendered':False}
(OUT/'roof-imagery-statistics.json').write_text(json.dumps(report,indent=2))
print(json.dumps({'passed':True,'roofCrops':[(x['id'],x['roofCrop']['darkPixelFractionLumaUnder55'],x['roofCrop']['luminancePercentiles']['50'],x['sourceTileCount']) for x in textures], 'sourceTileCount':len(tiles),'seams':[(x['orientation'],x['meanRgbRmsDelta']) for x in seams]}))