# -*- coding: utf-8 -*- """Comprehensive Model Divergence & Composition Drift Evaluation Suite. Computes 6 gold-standard quantitative metrics between ground-truth reference images (unquantized BF16) and test candidate images (quantized variants): 1. PSNR (Peak Signal-to-Noise Ratio) 2. SSIM (Structural Similarity Index) 3. SigLIP Global Cosine Similarity (Semantic Alignment / Instruction Following) 4. SigLIP Spatial Patch Similarity (Composition & Layout Drift) 5. Sobel Edge Acuity (Edge Sharpness vs Blur) 6. Chroma Spill Index (SECAM-like Cb/Cr Color Bleed across Luminance Boundaries) 7. Radial 2D FFT High-Frequency Spectral Power (Microscopic Texture Preservation) """ import os import sys import json import math import numpy as np import cv2 import torch import torch.nn.functional as F from PIL import Image from transformers import AutoProcessor, AutoModel _SIGLIP_MODEL = None _SIGLIP_PROCESSOR = None def get_siglip(device: str = "cuda:0"): global _SIGLIP_MODEL, _SIGLIP_PROCESSOR if _SIGLIP_MODEL is None: model_name = "google/siglip-base-patch16-224" _SIGLIP_PROCESSOR = AutoProcessor.from_pretrained(model_name, local_files_only=True) _SIGLIP_MODEL = AutoModel.from_pretrained(model_name, local_files_only=True).to(device).eval() return _SIGLIP_MODEL, _SIGLIP_PROCESSOR def compute_psnr_ssim(img_ref_np: np.ndarray, img_cand_np: np.ndarray) -> tuple[float, float]: """Compute PSNR and SSIM between two RGB images in uint8 [0, 255].""" if img_ref_np.shape != img_cand_np.shape: img_cand_np = cv2.resize(img_cand_np, (img_ref_np.shape[1], img_ref_np.shape[0])) # PSNR mse = np.mean((img_ref_np.astype(np.float64) - img_cand_np.astype(np.float64)) ** 2) if mse == 0: psnr = 100.0 else: psnr = float(20 * math.log10(255.0 / math.sqrt(mse))) # SSIM on grayscale ref_gray = cv2.cvtColor(img_ref_np, cv2.COLOR_RGB2GRAY).astype(np.float64) cand_gray = cv2.cvtColor(img_cand_np, cv2.COLOR_RGB2GRAY).astype(np.float64) C1 = (0.01 * 255) ** 2 C2 = (0.03 * 255) ** 2 mu1 = cv2.GaussianBlur(ref_gray, (11, 11), 1.5) mu2 = cv2.GaussianBlur(cand_gray, (11, 11), 1.5) mu1_sq = mu1 * mu1 mu2_sq = mu2 * mu2 mu1_mu2 = mu1 * mu2 sigma1_sq = cv2.GaussianBlur(ref_gray * ref_gray, (11, 11), 1.5) - mu1_sq sigma2_sq = cv2.GaussianBlur(cand_gray * cand_gray, (11, 11), 1.5) - mu2_sq sigma12 = cv2.GaussianBlur(ref_gray * cand_gray, (11, 11), 1.5) - mu1_mu2 ssim_map = ((2 * mu1_mu2 + C1) * (2 * sigma12 + C2)) / ((mu1_sq + mu2_sq + C1) * (sigma1_sq + sigma2_sq + C2)) ssim = float(np.mean(ssim_map)) return psnr, ssim def compute_composition_alignment( img_ref: Image.Image, img_cand: Image.Image, device: str = "cuda:0" ) -> tuple[float, float]: """Compute SigLIP global semantic cosine similarity and spatial patch layout similarity.""" model, processor = get_siglip(device) inputs_ref = processor(images=img_ref, return_tensors="pt").to(device) inputs_cand = processor(images=img_cand, return_tensors="pt").to(device) with torch.no_grad(): out_ref = model.vision_model(pixel_values=inputs_ref.pixel_values) out_cand = model.vision_model(pixel_values=inputs_cand.pixel_values) # Pooled output for global semantic instruction alignment pooled_ref = F.normalize(out_ref.pooler_output, dim=-1) pooled_cand = F.normalize(out_cand.pooler_output, dim=-1) global_sim = float(torch.sum(pooled_ref * pooled_cand).item()) # Patch tokens for spatial layout / composition alignment (excluding CLS/pooling if any) # last_hidden_state: [1, 196, 768] patches_ref = F.normalize(out_ref.last_hidden_state, dim=-1) patches_cand = F.normalize(out_cand.last_hidden_state, dim=-1) # Cosine similarity across corresponding spatial patches patch_sim = float(torch.mean(torch.sum(patches_ref * patches_cand, dim=-1)).item()) return global_sim, patch_sim def compute_edge_acuity(img_np: np.ndarray) -> float: """Compute Sobel edge acuity (mean gradient magnitude on luminance).""" gray = cv2.cvtColor(img_np, cv2.COLOR_RGB2GRAY).astype(np.float32) gx = cv2.Sobel(gray, cv2.CV_32F, 1, 0, ksize=3) gy = cv2.Sobel(gray, cv2.CV_32F, 0, 1, ksize=3) mag = np.sqrt(gx**2 + gy**2) return float(np.mean(mag)) def compute_chroma_spill(img_np: np.ndarray) -> float: """Compute Chroma Spill: standard deviation of Cb/Cr in flat (low luminance gradient) regions.""" ycrcb = cv2.cvtColor(img_np, cv2.COLOR_RGB2YCrCb).astype(np.float32) y, cr, cb = ycrcb[:, :, 0], ycrcb[:, :, 1], ycrcb[:, :, 2] gx = cv2.Sobel(y, cv2.CV_32F, 1, 0, ksize=3) gy = cv2.Sobel(y, cv2.CV_32F, 0, 1, ksize=3) grad_mag = np.sqrt(gx**2 + gy**2) # Flat region threshold (bottom 30% of gradient magnitude) threshold = np.percentile(grad_mag, 30) flat_mask = grad_mag < threshold cb_std = float(np.std(cb[flat_mask])) cr_std = float(np.std(cr[flat_mask])) return (cb_std + cr_std) / 2.0 def compute_hf_fourier_power(img_np: np.ndarray) -> float: """Compute 2D FFT radial energy in the outer 50% frequency band.""" gray = cv2.cvtColor(img_np, cv2.COLOR_RGB2GRAY).astype(np.float32) h, w = gray.shape f = np.fft.fft2(gray) fshift = np.fft.fftshift(f) magnitude = np.abs(fshift) # Radial mask for high frequencies cy, cx = h // 2, w // 2 y_idx, x_idx = np.ogrid[:h, :w] r = np.sqrt((x_idx - cx) ** 2 + (y_idx - cy) ** 2) r_max = np.sqrt(cx**2 + cy**2) hf_mask = r > (0.5 * r_max) total_energy = np.sum(magnitude) hf_energy = np.sum(magnitude[hf_mask]) return float((hf_energy / total_energy) if total_energy > 0 else 0.0) def evaluate_pair(ref_path: str, cand_path: str, device: str = "cuda:0") -> dict: """Run full evaluation comparing a candidate image against the unquantized reference.""" ref_pil = Image.open(ref_path).convert("RGB") cand_pil = Image.open(cand_path).convert("RGB") ref_np = np.array(ref_pil) cand_np = np.array(cand_pil) psnr, ssim = compute_psnr_ssim(ref_np, cand_np) global_sim, patch_sim = compute_composition_alignment(ref_pil, cand_pil, device=device) ref_acuity = compute_edge_acuity(ref_np) cand_acuity = compute_edge_acuity(cand_np) acuity_ratio = cand_acuity / (ref_acuity + 1e-6) ref_spill = compute_chroma_spill(ref_np) cand_spill = compute_chroma_spill(cand_np) spill_ratio = cand_spill / (ref_spill + 1e-6) ref_hf = compute_hf_fourier_power(ref_np) cand_hf = compute_hf_fourier_power(cand_np) hf_ratio = cand_hf / (ref_hf + 1e-6) return { "psnr_db": round(psnr, 2), "ssim": round(ssim, 4), "global_semantic_sim": round(global_sim, 4), "spatial_composition_sim": round(patch_sim, 4), "edge_acuity": { "ref": round(ref_acuity, 2), "cand": round(cand_acuity, 2), "ratio": round(acuity_ratio, 4), }, "chroma_spill": { "ref": round(ref_spill, 2), "cand": round(cand_spill, 2), "ratio": round(spill_ratio, 4), }, "hf_power": { "ref": round(ref_hf * 1000, 2), "cand": round(cand_hf * 1000, 2), "ratio": round(hf_ratio, 4), } } if __name__ == "__main__": if len(sys.argv) < 3: print("Usage: python eval_divergence_metrics.py ") sys.exit(1) res = evaluate_pair(sys.argv[1], sys.argv[2]) print(json.dumps(res, indent=2))