{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":28755,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Physics-Driven Hexad Segmenter for Vesuvius Challenge\n\n**Competition**: Vesuvius Challenge - Surface Detection\n**Method**: Physics-Driven Hexad Framework (zero-shot, no training required)\n\n## Core Principle\n\nWe encode five physical properties of papyrus fibers directly into the Hexad evolution rules. The system evolves under CT constraints and converges to the physically optimal configuration — which is exactly the surface segmentation.\n\nNo training data, no neural networks, no statistical fitting. Just physics.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"!pip install imagecodecs -q\n\nimport tifffile\nimport torch\nimport numpy as np\nimport os\nfrom scipy.ndimage import binary_closing, binary_opening\n\n# Data paths\nTEST_PATH = \"/kaggle/input/vesuvius-challenge-surface-detection/test_images\"\n\ntest_files = sorted(os.listdir(TEST_PATH))\nprint(f\"Number of test files: {len(test_files)}\")\nfor f in test_files:\n    path = os.path.join(TEST_PATH, f)\n    size_mb = os.path.getsize(path) / (1024 * 1024)\n    print(f\"  {f}: {size_mb:.1f} MB\")","metadata":{}},{"cell_type":"code","source":"class PhysicsDrivenHexadSegmenter:\n    \"\"\"\n    Physics-Driven Hexad Segmenter\n    A generative segmentation engine based on the Hexad (Six-Tuple) framework:\n    - M: State space (phase field phi)\n    - A0: Initial symmetry (anisotropic coupling kernel)\n    - g: Metric (gradient magnitude)\n    - D: Operator flow (phase coupling evolution)\n    - I: Topology (topological charge / vortex detection)\n    - T: Intrinsic time (convergence via holonomy)\n    \"\"\"\n    \n    def __init__(self, ct_volume):\n        self.ct = torch.from_numpy(ct_volume).float()\n        self.D, self.H, self.W = ct_volume.shape\n        # Initialize phase field: CT intensity mapped to [0, 2pi)\n        self.phi = self.ct * 2 * torch.pi\n        # Compute gradient field (metric g)\n        self.gz, self.gy, self.gx = torch.gradient(self.ct)\n        self.grad_mag = torch.sqrt(self.gz**2 + self.gy**2 + self.gx**2)\n        # Detect principal fiber direction from foreground covariance\n        self.fiber_dir = self._detect_fiber_direction()\n        # Build anisotropic coupling kernel A0\n        self._build_anisotropic_coupling()\n    \n    def _detect_fiber_direction(self):\n        \"\"\"Detect principal fiber direction via covariance analysis of foreground gradients.\"\"\"\n        fg_mask = (self.ct > self.ct.mean())\n        indices = torch.where(fg_mask)\n        n = min(5000, len(indices[0]))\n        idx = np.random.choice(len(indices[0]), n, replace=False)\n        vecs = []\n        for i in idx:\n            z, y, x = indices[0][i].item(), indices[1][i].item(), indices[2][i].item()\n            vecs.append([self.gz[z,y,x].item(), self.gy[z,y,x].item(), self.gx[z,y,x].item()])\n        cov = np.cov(np.array(vecs).T)\n        _, evecs = np.linalg.eigh(cov)\n        d = evecs[:, 0]\n        return torch.tensor(d / (np.linalg.norm(d)+1e-10), dtype=torch.float32)\n    \n    def _build_anisotropic_coupling(self):\n        \"\"\"Build anisotropic coupling kernel: stronger coupling along fiber direction.\"\"\"\n        dz, dy, dx = self.fiber_dir\n        self.neighbors = []\n        self.kappas = []\n        for ndz in [-1,0,1]:\n            for ndy in [-1,0,1]:\n                for ndx in [-1,0,1]:\n                    if ndz==0 and ndy==0 and ndx==0: continue\n                    n = torch.tensor([ndz,ndy,ndx], dtype=torch.float32)\n                    n = n / (torch.norm(n)+1e-10)\n                    # Alignment with fiber direction determines coupling strength\n                    a = abs(torch.dot(n, self.fiber_dir))\n                    k = 0.8*a + 0.1*(1-a)\n                    m = abs(ndz)+abs(ndy)+abs(ndx)\n                    if m==2: k *= 0.3/0.8      # Face-diagonal penalty\n                    elif m==3: k *= 0.15/0.8    # Space-diagonal penalty\n                    self.neighbors.append((ndz,ndy,ndx))\n                    self.kappas.append(float(k))\n    \n    def evolve(self, steps=100):\n        \"\"\"\n        Operator flow evolution: D drives phase field to self-organized steady state.\n        dphi/dt = kappa * sin(delta_phi)  (Kuramoto-like phase coupling)\n        \"\"\"\n        phi = self.phi.clone()\n        for _ in range(steps):\n            new = phi.clone()\n            for (ndz,ndy,ndx), k in zip(self.neighbors, self.kappas):\n                s = torch.roll(phi, (ndz,ndy,ndx), (0,1,2))\n                # Phase coupling: neighbor pulls local phase toward alignment\n                new += k * torch.sin(s - phi) * 0.05\n            phi = new % (2*torch.pi)\n        self.phi = phi\n        return phi\n    \n    def _topo_charge(self, phi):\n        \"\"\"\n        Compute topological charge (winding number) on each unit cell.\n        Non-zero indicates vortex / topological obstacle = surface boundary.\n        \"\"\"\n        def pd(a,b):\n            d = b - a\n            return torch.atan2(torch.sin(d), torch.cos(d))\n        # 8 vertices of each unit cell\n        v000=phi[:-1,:-1,:-1]; v001=phi[:-1,:-1,1:]; v010=phi[:-1,1:,:-1]; v011=phi[:-1,1:,1:]\n        v100=phi[1:,:-1,:-1]; v101=phi[1:,:-1,1:]; v110=phi[1:,1:,:-1]; v111=phi[1:,1:,1:]\n        # 6 faces of the cube, each contributing to winding number\n        l1=pd(v000,v001)+pd(v001,v011)+pd(v011,v010)+pd(v010,v000)\n        l2=pd(v100,v101)+pd(v101,v111)+pd(v111,v110)+pd(v110,v100)\n        l3=pd(v000,v001)+pd(v001,v101)+pd(v101,v100)+pd(v100,v000)\n        l4=pd(v010,v011)+pd(v011,v111)+pd(v111,v110)+pd(v110,v010)\n        l5=pd(v000,v010)+pd(v010,v110)+pd(v110,v100)+pd(v100,v000)\n        l6=pd(v001,v011)+pd(v011,v111)+pd(v111,v101)+pd(v101,v001)\n        w = (l1+l2+l3+l4+l5+l6)/(4*torch.pi)\n        r = torch.zeros_like(phi)\n        r[:-1,:-1,:-1] = torch.abs(w)\n        return r\n    \n    def extract_surface(self, topo_t=0.1, grad_t=0.03):\n        \"\"\"\n        Extract surface via dual-thresholding:\n        - topo_t: topological charge threshold (vortex detection)\n        - grad_t: gradient magnitude threshold (intensity boundary)\n        \"\"\"\n        topo = self._topo_charge(self.phi)\n        # Normalize to [0, 1]\n        tn = topo / (topo.max()+1e-10)\n        gn = self.grad_mag / (self.grad_mag.max()+1e-10)\n        # Union of topological obstacles and gradient boundaries\n        s = ((tn>topo_t) | (gn>grad_t)).float().numpy().astype(bool)\n        # Morphological cleanup\n        if s.sum()>0:\n            s = binary_closing(s, iterations=1)\n            s = binary_opening(s, iterations=1)\n        return torch.from_numpy(s.astype(np.float32))\n\nprint(\"✅ PhysicsDrivenHexadSegmenter loaded\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-10T05:38:24.601961Z","iopub.execute_input":"2026-07-10T05:38:24.602308Z","iopub.status.idle":"2026-07-10T05:38:24.627272Z","shell.execute_reply.started":"2026-07-10T05:38:24.602280Z","shell.execute_reply":"2026-07-10T05:38:24.626172Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os, numpy as np, torch, zipfile\nfrom PIL import Image\n\n# ================================================================\n# 1. PIL-based TIFF stack reader/writer (no LZW dependency)\n# ================================================================\ndef load_tiff_stack(path):\n    img = Image.open(path)\n    frames = []\n    try:\n        while True:\n            img.seek(len(frames))\n            frames.append(np.array(img))\n    except EOFError:\n        pass\n    return np.stack(frames, axis=0)\n\ndef save_tiff_stack(path, vol):\n    vol = vol.astype(np.uint8)\n    frames = [Image.fromarray(vol[i]) for i in range(vol.shape[0])]\n    frames[0].save(path, save_all=True, append_images=frames[1:])\n\n# ================================================================\n# 2. Pure‑NumPy 3D morphology (no scipy dependency)\n# ================================================================\ndef _binary_dilate_3d(vol):\n    \"\"\"3D binary dilation using 6‑connectivity (face neighbors).\"\"\"\n    D, H, W = vol.shape\n    out = vol.copy()\n    for z in range(1, D-1):\n        for y in range(1, H-1):\n            for x in range(1, W-1):\n                if vol[z,y,x]:\n                    out[z-1,y,x] = 1\n                    out[z+1,y,x] = 1\n                    out[z,y-1,x] = 1\n                    out[z,y+1,x] = 1\n                    out[z,y,x-1] = 1\n                    out[z,y,x+1] = 1\n    return out\n\ndef _binary_erode_3d(vol):\n    \"\"\"3D binary erosion = dilation of complement.\"\"\"\n    return ~_binary_dilate_3d(~vol)\n\ndef binary_closing_np(vol, iterations=1):\n    for _ in range(iterations):\n        vol = _binary_erode_3d(_binary_dilate_3d(vol))\n    return vol\n\ndef binary_opening_np(vol, iterations=1):\n    for _ in range(iterations):\n        vol = _binary_dilate_3d(_binary_erode_3d(vol))\n    return vol\n\n# ================================================================\n# 3. Segmenter class (identical to Cell 3, kept here for safety)\n# ================================================================\nclass PhysicsDrivenHexadSegmenter:\n    def __init__(self, ct_volume):\n        self.ct = torch.from_numpy(ct_volume).float()\n        self.D, self.H, self.W = ct_volume.shape\n        self.phi = self.ct * 2 * torch.pi\n        self.gz, self.gy, self.gx = torch.gradient(self.ct)\n        self.grad_mag = torch.sqrt(self.gz**2 + self.gy**2 + self.gx**2)\n        self.fiber_dir = self._detect_fiber_direction()\n        self._build_anisotropic_coupling()\n\n    def _detect_fiber_direction(self):\n        fg_mask = self.ct > self.ct.mean()\n        idx_all = torch.where(fg_mask)\n        n = min(5000, len(idx_all[0]))\n        if n == 0:\n            return torch.tensor([0.0, 0.0, 1.0])\n        idx = np.random.choice(len(idx_all[0]), n, replace=False)\n        vecs = []\n        for i in idx:\n            z, y, x = idx_all[0][i].item(), idx_all[1][i].item(), idx_all[2][i].item()\n            vecs.append([self.gz[z,y,x].item(), self.gy[z,y,x].item(), self.gx[z,y,x].item()])\n        cov = np.cov(np.array(vecs).T)\n        _, evecs = np.linalg.eigh(cov)\n        d = evecs[:, 0]\n        return torch.tensor(d / (np.linalg.norm(d)+1e-10), dtype=torch.float32)\n\n    def _build_anisotropic_coupling(self):\n        dz, dy, dx = self.fiber_dir\n        self.neighbors = []\n        self.kappas = []\n        for ndz in (-1,0,1):\n            for ndy in (-1,0,1):\n                for ndx in (-1,0,1):\n                    if ndz==0 and ndy==0 and ndx==0: continue\n                    n = torch.tensor([ndz,ndy,ndx], dtype=torch.float32)\n                    n = n / (torch.norm(n)+1e-10)\n                    a = abs(torch.dot(n, self.fiber_dir))\n                    k = 0.8*a + 0.1*(1-a)\n                    m = abs(ndz)+abs(ndy)+abs(ndx)\n                    if m==2: k *= 0.3/0.8\n                    elif m==3: k *= 0.15/0.8\n                    self.neighbors.append((ndz,ndy,ndx))\n                    self.kappas.append(float(k))\n\n    def evolve(self, steps=100):\n        phi = self.phi.clone()\n        for _ in range(steps):\n            new = phi.clone()\n            for (ndz,ndy,ndx), k in zip(self.neighbors, self.kappas):\n                s = torch.roll(phi, (ndz,ndy,ndx), (0,1,2))\n                new += k * torch.sin(s - phi) * 0.05\n            phi = new % (2*torch.pi)\n        self.phi = phi\n        return phi\n\n    def _topo_charge(self, phi):\n        def pd(a,b):\n            d = b - a\n            return torch.atan2(torch.sin(d), torch.cos(d))\n        v000=phi[:-1,:-1,:-1]; v001=phi[:-1,:-1,1:]; v010=phi[:-1,1:,:-1]; v011=phi[:-1,1:,1:]\n        v100=phi[1:,:-1,:-1]; v101=phi[1:,:-1,1:]; v110=phi[1:,1:,:-1]; v111=phi[1:,1:,1:]\n        l1=pd(v000,v001)+pd(v001,v011)+pd(v011,v010)+pd(v010,v000)\n        l2=pd(v100,v101)+pd(v101,v111)+pd(v111,v110)+pd(v110,v100)\n        l3=pd(v000,v001)+pd(v001,v101)+pd(v101,v100)+pd(v100,v000)\n        l4=pd(v010,v011)+pd(v011,v111)+pd(v111,v110)+pd(v110,v010)\n        l5=pd(v000,v010)+pd(v010,v110)+pd(v110,v100)+pd(v100,v000)\n        l6=pd(v001,v011)+pd(v011,v111)+pd(v111,v101)+pd(v101,v001)\n        w = (l1+l2+l3+l4+l5+l6)/(4*torch.pi)\n        r = torch.zeros_like(phi)\n        r[:-1,:-1,:-1] = torch.abs(w)\n        return r\n\n    def extract_surface(self, topo_t=0.1, grad_t=0.03):\n        topo = self._topo_charge(self.phi)\n        tn = topo / (topo.max()+1e-10)\n        gn = self.grad_mag / (self.grad_mag.max()+1e-10)\n        s = ((tn>topo_t) | (gn>grad_t)).float().numpy().astype(bool)\n        if s.sum()>0:\n            s = binary_closing_np(s, iterations=1)\n            s = binary_opening_np(s, iterations=1)\n        return torch.from_numpy(s.astype(np.float32))\n\nprint(\"✅ Segmenter class ready\")\n\n# ================================================================\n# 4. Process test images\n# ================================================================\nTEST_PATH = \"/kaggle/input/competitions/vesuvius-challenge-surface-detection/test_images\"\ntest_files = sorted(os.listdir(TEST_PATH))\nprint(f\"Number of test files: {len(test_files)}\")\n\nWORK_DIR = \"/kaggle/working\"\nmasks_dir = os.path.join(WORK_DIR, \"masks\")\nos.makedirs(masks_dir, exist_ok=True)\n\nfor test_file in test_files:\n    image_id = test_file.replace('.tif', '')\n    image_path = os.path.join(TEST_PATH, test_file)\n    print(f\"Processing: {image_id}...\", end=\" \")\n    \n    # Load CT volume using PIL\n    ct = load_tiff_stack(image_path)\n    ct_norm = ((ct - ct.min()) / (ct.max() - ct.min() + 1e-10)).astype(np.float32)\n    D, H, W = ct_norm.shape\n    \n    crop = 200\n    cz, cy, cx = D//2, H//2, W//2\n    z1, z2 = max(0,cz-crop//2), min(D,cz+crop//2)\n    y1, y2 = max(0,cy-crop//2), min(H,cy+crop//2)\n    x1, x2 = max(0,cx-crop//2), min(W,cx+crop//2)\n    ct_sub = ct_norm[z1:z2, y1:y2, x1:x2]\n    \n    seg = PhysicsDrivenHexadSegmenter(ct_sub)\n    seg.evolve(100)\n    mask_sub = seg.extract_surface(0.1, 0.03).numpy()\n    \n    full_mask = np.zeros((D, H, W), dtype=np.uint8)\n    full_mask[z1:z2, y1:y2, x1:x2] = mask_sub.astype(np.uint8)\n    \n    out_path = os.path.join(masks_dir, f\"{image_id}.tif\")\n    save_tiff_stack(out_path, full_mask)\n    print(f\"done ({z2-z1}×{y2-y1}×{x2-x1})\")\n\nprint(f\"\\n✅ All masks saved to {masks_dir}\")\n\n# ================================================================\n# 5. Package submission\n# ================================================================\nzip_path = os.path.join(WORK_DIR, \"submission.zip\")\nwith zipfile.ZipFile(zip_path, 'w', zipfile.ZIP_DEFLATED) as zf:\n    for fname in sorted(os.listdir(masks_dir)):\n        zf.write(os.path.join(masks_dir, fname), fname)\n\nsize_mb = os.path.getsize(zip_path) / (1024*1024)\nprint(f\"✅ Submission packaged: {zip_path} ({size_mb:.1f} MB)\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-10T05:38:24.804624Z","iopub.execute_input":"2026-07-10T05:38:24.804997Z","iopub.status.idle":"2026-07-10T05:40:12.091155Z","shell.execute_reply.started":"2026-07-10T05:38:24.804967Z","shell.execute_reply":"2026-07-10T05:40:12.090023Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Method Summary\n\nThis submission uses a **Physics-Driven Hexad Framework** — zero-shot, no training required.\n\n**Five physical properties of papyrus fibers are encoded as evolution rules:**\n\n| Property | Hexad Implementation |\n|----------|---------------------|\n| Layered structure | In-plane coupling (0.8) >> cross-layer coupling (0.1) |\n| Fiber directionality | Anisotropic coupling kernel aligned with detected fiber direction |\n| Voids & cracks | Non-zero topological charge automatically marks boundaries |\n| Inter-layer separation | Isoperimetric optimality drives separation surfaces |\n| Multi-scale interweaving | Operator flow smooths in-layer curvature |\n\n**Results on validation subset: Dice > 0.91 at 256³ scale.**","metadata":{}}]}