""" skiss_fourier_codec.py ================================================================================ 2D Skiss-Fourier Image Manifold Encoder / Decoder Theory: - Explicit pixel grids are skiss-violations (massive, inert, unstoreable). - Images are retro-compressed into a finite parametric Fourier surface: I(x,y; θ) = Σ_p Σ_q θ_pq · u_p(y) · v_q(x) - The pixel-coordinate path is erased; only the spectral residue θ survives. - Generation is performed on demand at any real-valued (x,y). ================================================================================ """ import numpy as np from pathlib import Path from typing import Union, Tuple, Dict, Optional # ============================================================================== # Stationary Basis Builders (the "Skiss Kernels") # ============================================================================== def _build_real_fourier_basis_1d(coords: np.ndarray, K: int, period: float) -> np.ndarray: """ Build the real 1D Fourier basis matrix on coordinate vector `coords`. Columns: [1, cos(2π·1·c/p), ..., cos(2π·K·c/p), sin(2π·1·c/p), ..., sin(2π·K·c/p)] Returns: B: ndarray of shape (len(coords), 2*K+1) """ coords = np.asarray(coords, dtype=np.float64) k = np.arange(1, K + 1, dtype=np.float64) ones = np.ones((coords.shape[0], 1), dtype=np.float64) if K == 0: return ones cos_part = np.cos(2.0 * np.pi * k * coords[:, None] / period) # (N, K) sin_part = np.sin(2.0 * np.pi * k * coords[:, None] / period) # (N, K) return np.concatenate([ones, cos_part, sin_part], axis=1) # (N, 2K+1) # ============================================================================== # Encoder: Explicit Image → Skiss-Complete Residue θ # ============================================================================== class SkissEncoder: """ Retro-compresses an explicit H×W image into a finite parameter tensor θ. The fitting exploits the separable Kronecker structure: min_Θ || I - U Θ Vᵀ ||_F² solved via the normal equations: A Θ Bᵀ = Uᵀ I V, where A = UᵀU, B = VᵀV. """ def __init__(self, K: Optional[int] = None, epsilon: Optional[float] = None, reg_lambda: float = 1e-6): """ Args: K: Manual critical rank (number of Fourier modes per axis). If None, epsilon must be provided. epsilon: Fraction of spectral energy to retain (0 < epsilon < 1). Used to auto-estimate K if K is not given. reg_lambda: Tiny ridge for numerical stability near degenerate modes. """ if K is None and epsilon is None: raise ValueError("Specify either K or epsilon.") self.K_manual = K self.epsilon = epsilon self.reg_lambda = reg_lambda self.theta: Optional[Dict] = None @staticmethod def _estimate_K(image: np.ndarray, epsilon: float) -> int: """Auto-estimate critical rank from FFT energy retention.""" if image.ndim == 3: gray = image.mean(axis=2) else: gray = image H, W = gray.shape F = np.fft.fft2(gray.astype(np.float64)) mag = np.abs(F).flatten() vals = np.sort(mag)[::-1] cumsum = np.cumsum(vals) total = cumsum[-1] + 1e-12 # Find how many coefficients to keep for (1-epsilon) energy n_keep = int(np.searchsorted(cumsum, (1.0 - epsilon) * total)) + 1 # Heuristic: n_keep ≈ (2K+1)² for the real Fourier basis K_est = int(np.ceil((np.sqrt(n_keep) - 1.0) / 2.0)) return max(K_est, 1) def encode(self, image: np.ndarray, verbose: bool = True) -> Dict: """ Compress an image into the skiss residue θ. Args: image: ndarray, shape (H, W) or (H, W, C). Values typically in [0, 255]. verbose: Print compression statistics. Returns: theta: Serializable dict containing the spectral residue. """ if image.ndim == 2: image = image[:, :, None] H, W, C = image.shape I = image.astype(np.float64) K = self.K_manual if self.K_manual is not None else self._estimate_K(I, self.epsilon) # Guard against Nyquist degeneracy (sin(kπ) = 0 at k = N/2) K = min(K, min(H, W) // 2 - 1) K = max(K, 0) # Sample coordinates on the original lattice y_coords = np.arange(H, dtype=np.float64) x_coords = np.arange(W, dtype=np.float64) # Stationary basis matrices (the unchanging "laws") U = _build_real_fourier_basis_1d(y_coords, K, float(H)) # (H, 2K+1) V = _build_real_fourier_basis_1d(x_coords, K, float(W)) # (W, 2K+1) # Normal-equation matrices A = UᵀU, B = VᵀV # For a uniform integer grid 0..N-1 with period N, these are diagonal: # diag(A) ≈ [N, N/2, ..., N/2] (cos & sin blocks) A = U.T @ U B = V.T @ V # Ridge regularization protects against singularities at Nyquist or small images if self.reg_lambda > 0: A += self.reg_lambda * np.eye(A.shape[0]) B += self.reg_lambda * np.eye(B.shape[0]) # Solve A Θ Bᵀ = RHS for each channel # Step 1: Z = A⁻¹ RHS # Step 2: Θ = Z (B⁻¹)ᵀ → Θᵀ = B⁻¹ Zᵀ coeffs = np.zeros((C, 2 * K + 1, 2 * K + 1), dtype=np.float64) for c in range(C): RHS = U.T @ I[:, :, c] @ V # (2K+1, 2K+1) Z = np.linalg.solve(A, RHS) # A Z = RHS Theta_c = np.linalg.solve(B, Z.T).T # B Θᵀ = Zᵀ coeffs[c] = Theta_c # Assemble the skiss-complete residue self.theta = { 'coefficients': coeffs, # (C, 2K+1, 2K+1) 'K': K, 'shape': (H, W), 'channels': C, 'period_y': float(H), 'period_x': float(W), 'basis': 'real_fourier_separable', 'reg_lambda': float(self.reg_lambda), } if verbose: orig_size = H * W * C param_size = coeffs.size ratio = orig_size / param_size print(f"[SkissEncoder] Image: {H}×{W}×{C}") print(f"[SkissEncoder] K={K} | Params: {param_size} | " f"Pixels: {orig_size} | Compression: {ratio:.1f}×") return self.theta def save(self, path: Union[str, Path]): """Save theta as a compressed .npz archive.""" if self.theta is None: raise RuntimeError("Nothing to save. Call encode() first.") np.savez_compressed(str(path), **self.theta) @staticmethod def load_theta(path: Union[str, Path]) -> Dict: """Load a theta dict from disk.""" data = np.load(str(path), allow_pickle=False) # Convert 0-d arrays back to scalars where convenient theta = {k: (v.item() if v.ndim == 0 else v) for k, v in data.items()} return theta # ============================================================================== # Decoder: Skiss Residue θ → On-Demand Image Generation # ============================================================================== class SkissDecoder: """ Generates pixels from the skiss residue θ at arbitrary coordinates. The explicit pixel path is never reconstructed unless explicitly requested. """ def __init__(self, theta: Dict): self.theta = theta self.K = int(theta['K']) self.coeffs = theta['coefficients'] # (C, 2K+1, 2K+1) self.C = int(theta['channels']) self.period_y = float(theta['period_y']) self.period_x = float(theta['period_x']) self.H_orig, self.W_orig = theta['shape'] def generate(self, x: np.ndarray, y: np.ndarray, channel: Optional[int] = None) -> np.ndarray: """ Evaluate the parametric Fourier surface at arbitrary real coordinates. Args: x: ndarray of arbitrary shape (float, in units of original pixel grid). y: ndarray, same shape as x. channel: None → return all channels (..., C). int → return single channel (same shape as x). Returns: values: ndarray of generated pixel intensities. """ x = np.asarray(x, dtype=np.float64) y = np.asarray(y, dtype=np.float64) if x.shape != y.shape: raise ValueError("x and y must have identical shape.") flat_x = x.ravel() flat_y = y.ravel() n_pts = flat_x.size # Evaluate stationary basis at the query coordinates (the "skiss kernels") V = _build_real_fourier_basis_1d(flat_x, self.K, self.period_x) # (n, 2K+1) U = _build_real_fourier_basis_1d(flat_y, self.K, self.period_y) # (n, 2K+1) if channel is not None: Theta = self.coeffs[channel] # (2K+1, 2K+1) # Batch evaluation: I_n = Σ_{p,q} U[n,p] * Theta[p,q] * V[n,q] # = np.sum( (U @ Theta) * V, axis=1 ) vals = np.sum((U @ Theta) * V, axis=1) return vals.reshape(x.shape) else: out = np.zeros((n_pts, self.C), dtype=np.float64) for c in range(self.C): Theta = self.coeffs[c] out[:, c] = np.sum((U @ Theta) * V, axis=1) return out.reshape(x.shape + (self.C,)) def reconstruct(self, shape: Optional[Tuple[int, int]] = None) -> np.ndarray: """ Reconstruct a full pixel grid. If shape is None, uses original dimensions. """ if shape is None: shape = (self.H_orig, self.W_orig) H, W = shape y_grid, x_grid = np.mgrid[0:H, 0:W] return self.generate(x_grid, y_grid) def reconstruct_patch(self, x_min: float, x_max: float, y_min: float, y_max: float, res_x: int, res_y: int) -> np.ndarray: """ Generate a continuous viewport patch at arbitrary resolution. Useful for zooming beyond native resolution without interpolation artifacts. """ x = np.linspace(x_min, x_max, res_x) y = np.linspace(y_min, y_max, res_y) y_grid, x_grid = np.meshgrid(y, x, indexing='ij') return self.generate(x_grid, y_grid) # ============================================================================== # Diagnostics # ============================================================================== def skiss_metrics(original: np.ndarray, reconstructed: np.ndarray) -> Dict[str, float]: """Return MSE and PSNR between original and skiss-reconstructed images.""" original = original.astype(np.float64) reconstructed = reconstructed.astype(np.float64) mse = float(np.mean((original - reconstructed) ** 2)) if mse == 0: psnr = float('inf') else: max_val = 255.0 if original.max() > 1.0 else 1.0 psnr = 20.0 * np.log10(max_val / np.sqrt(mse)) return {'mse': mse, 'psnr': psnr, 'skiss_complete': psnr > 30.0} # ============================================================================== # Example / Self-Test # ============================================================================== if __name__ == "__main__": print("=" * 70) print("Skiss-Fourier Image Codec — Self-Test") print("=" * 70) # 1. Generate a synthetic test image (grayscale + periodic structure) H, W = 256, 256 y, x = np.mgrid[0:H, 0:W] # A "recognizable" pattern: smooth blobs + sharp cross (high-freq path to erase) test_img = ( 128.0 + 64.0 * np.sin(2 * np.pi * 3 * x / W) * np.cos(2 * np.pi * 2 * y / H) + 32.0 * np.exp(-((x - W/2)**2 + (y - H/2)**2) / (2 * (30.0)**2)) + 16.0 * (np.abs(x - W/2) < 5) # vertical line + 16.0 * (np.abs(y - H/2) < 5) # horizontal line ) test_img = np.clip(test_img, 0, 255).astype(np.uint8) print(f"\nSynthetic test image: {H}×{W}") # 2. ENCODE: retro-compress the explicit pixel grid into θ encoder = SkissEncoder(K=16) # Critical rank 16 → 33×33 = 1089 coeffs theta = encoder.encode(test_img) encoder.save("skiss_residue.npz") # 3. DECODE: generate from skiss residue on demand theta_loaded = SkissEncoder.load_theta("skiss_residue.npz") decoder = SkissDecoder(theta_loaded) # a) Full grid reconstruction (the explicit path temporarily regenerated) rec = decoder.reconstruct() metrics = skiss_metrics(test_img, rec) print(f"\nFull reconstruction metrics: {metrics}") # b) On-demand continuous sampling (not aligned to original pixel lattice) # Query 5× sub-sampled coordinates at HALF-INTEGER offsets fine_x = np.linspace(0, W - 1, W * 2) + 0.5 # 2× supersample, shifted fine_y = np.linspace(0, H - 1, H * 2) + 0.5 fine_y_grid, fine_x_grid = np.meshgrid(fine_y, fine_x, indexing='ij') super_res = decoder.generate(fine_x_grid, fine_y_grid) print(f"\nSuper-resolution sample shape: {super_res.shape} " f"(generated at 2× native resolution from θ alone)") # c) Partial viewport generation (only pay work for visible region) patch = decoder.reconstruct_patch(x_min=50, x_max=80, y_min=50, y_max=80, res_x=300, res_y=300) print(f"Viewport patch shape: {patch.shape} " f"(arbitrary resolution from skiss manifold)") # 4. Demonstrate progressive recognition (sending low-K first) print("\n--- Progressive Skiss Transmission ---") for K_test in [0, 2, 4, 8, 16]: enc_tmp = SkissEncoder(K=K_test) th_tmp = enc_tmp.encode(test_img, verbose=False) dec_tmp = SkissDecoder(th_tmp) rec_tmp = dec_tmp.reconstruct() m = skiss_metrics(test_img, rec_tmp) status = "✅ Solid" if m['psnr'] > 35 else ("⚠️ Liquid" if m['psnr'] > 20 else "❌ Gas") print(f" K={K_test:2d} | Params={th_tmp['coefficients'].size:5d} " f"| PSNR={m['psnr']:6.2f} dB | {status}") print("\nDone. Residue saved to: skiss_residue.npz")