{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":112509,"databundleVersionId":14254895,"sourceType":"competition"}],"dockerImageVersionId":31239,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Hidden elements are difficult to detect","metadata":{}},{"cell_type":"code","source":"import os\nimport cv2\nimport torch\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom sklearn.cluster import KMeans\nfrom sklearn.preprocessing import StandardScaler\nfrom scipy.ndimage import median_filter\nimport timm\nfrom pathlib import Path\n\n# Configuration\nclass CFG:\n    device = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n    image_path = \"/kaggle/input/csiro-biomass/train/ID1012260530.jpg\"\n    n_clusters = 8\n    patch_size = 14\n    input_size = 518\n    upsample_factor = 4\n    \n    # Feature weights\n    dino_weight = 0.4\n    color_weight = 0.3\n    veg_index_weight = 0.2\n    texture_weight = 0.1\n    \ndef load_dino_backbone():\n    model = timm.create_model(\n        \"vit_base_patch14_dinov2\",\n        pretrained=True,\n        num_classes=0,\n    )\n    model.eval()\n    model.to(CFG.device)\n    print(f\"Loaded DINOv2 Base | device={CFG.device}\")\n    return model\n\ndef extract_patch_features(model, image):\n    img = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n    img = cv2.resize(img, (CFG.input_size, CFG.input_size))\n    \n    mean = np.array([0.485, 0.456, 0.406])\n    std = np.array([0.229, 0.224, 0.225])\n    img_norm = (img / 255.0 - mean) / std\n    \n    img_tensor = torch.from_numpy(img_norm).permute(2, 0, 1).float().unsqueeze(0)\n    img_tensor = img_tensor.to(CFG.device)\n    \n    with torch.no_grad():\n        features = model.forward_features(img_tensor)\n        patch_tokens = features[:, 1:, :]\n        \n    n_patches = patch_tokens.shape[1]\n    grid_h = grid_w = int(np.sqrt(n_patches))\n    \n    return patch_tokens.squeeze(0).cpu().numpy(), (grid_h, grid_w), img\n\ndef compute_vegetation_indices(image):\n    \"\"\"Compute vegetation indices with CORRECT normalization.\"\"\"\n    img = image.astype(np.float32) / 255.0\n    R, G, B = img[:, :, 0], img[:, :, 1], img[:, :, 2]\n    eps = 1e-6\n    \n    # ExG range: approximately [-2, +2] for RGB in [0,1]\n    ExG = 2*G - R - B\n    \n    # ExGR \n    ExR = 1.4*R - G\n    ExGR = ExG - ExR\n    \n    # NGRDI: already in [-1, 1]\n    NGRDI = (G - R) / (G + R + eps)\n    \n    # GLI: already in [-1, 1]\n    GLI = (2*G - R - B) / (2*G + R + B + eps)\n    \n    # VARI\n    VARI = (G - R) / (G + R - B + eps)\n    \n    # TGI range: approximately [-0.6, +0.6]\n    TGI = G - 0.39*R - 0.61*B\n    \n    return {\n        'ExG': ExG,\n        'ExGR': ExGR,\n        'NGRDI': NGRDI,\n        'GLI': GLI,\n        'VARI': VARI,\n        'TGI': TGI\n    }\n\ndef compute_texture_features(image):\n    gray = cv2.cvtColor(image, cv2.COLOR_RGB2GRAY).astype(np.float32)\n    kernel_size = 5\n    mean = cv2.blur(gray, (kernel_size, kernel_size))\n    sqr_mean = cv2.blur(gray**2, (kernel_size, kernel_size))\n    variance = sqr_mean - mean**2\n    \n    sobelx = cv2.Sobel(gray, cv2.CV_32F, 1, 0, ksize=3)\n    sobely = cv2.Sobel(gray, cv2.CV_32F, 0, 1, ksize=3)\n    edge_mag = np.sqrt(sobelx**2 + sobely**2)\n    \n    laplacian = np.abs(cv2.Laplacian(gray, cv2.CV_32F))\n    \n    return {\n        'variance': variance,\n        'edge_mag': edge_mag,\n        'laplacian': laplacian,\n    }\n\ndef extract_enhanced_features(image, grid_size):\n    h, w = image.shape[:2]\n    gh, gw = grid_size\n    \n    patch_h = h // gh\n    patch_w = w // gw\n    \n    hsv = cv2.cvtColor(image, cv2.COLOR_RGB2HSV)\n    veg_indices = compute_vegetation_indices(image)\n    texture = compute_texture_features(image)\n    \n    features = []\n    for i in range(gh):\n        for j in range(gw):\n            y1, y2 = i * patch_h, (i + 1) * patch_h\n            x1, x2 = j * patch_w, (j + 1) * patch_w\n            \n            patch_hsv = hsv[y1:y2, x1:x2]\n            \n            h_mean, s_mean, v_mean = patch_hsv.mean(axis=(0, 1))\n            h_std, s_std, v_std = patch_hsv.std(axis=(0, 1))\n            h_channel = patch_hsv[:, :, 0]\n            green_ratio = ((h_channel >= 40) & (h_channel <= 85)).mean()\n            dead_ratio = (h_channel < 40).mean()\n            \n            veg_feats = []\n            for name, idx_map in veg_indices.items():\n                patch_idx = idx_map[y1:y2, x1:x2]\n                veg_feats.extend([patch_idx.mean(), patch_idx.std()])\n            \n            tex_feats = []\n            for name, tex_map in texture.items():\n                patch_tex = tex_map[y1:y2, x1:x2]\n                tex_feats.append(patch_tex.mean())\n            \n            patch_features = [\n                h_mean, h_std, s_mean, s_std, v_mean, v_std,\n                green_ratio * 255, dead_ratio * 255,\n            ] + veg_feats + tex_feats\n            \n            features.append(patch_features)\n    \n    return np.array(features)\n\ndef upsample_features(features, grid_size, factor=4):\n    gh, gw = grid_size\n    feat_dim = features.shape[1]\n    \n    feat_map = features.reshape(gh, gw, feat_dim)\n    new_h, new_w = gh * factor, gw * factor\n    upsampled = np.zeros((new_h, new_w, feat_dim))\n    \n    for c in range(feat_dim):\n        upsampled[:, :, c] = cv2.resize(\n            feat_map[:, :, c], \n            (new_w, new_h),\n            interpolation=cv2.INTER_LINEAR\n        )\n    \n    return upsampled.reshape(-1, feat_dim), (new_h, new_w)\n\ndef cluster_patches(features, n_clusters=8):\n    scaler = StandardScaler()\n    features_norm = scaler.fit_transform(features)\n    kmeans = KMeans(n_clusters=n_clusters, random_state=42, n_init=10)\n    labels = kmeans.fit_predict(features_norm)\n    return labels, kmeans\n\ndef post_process_segmentation(label_map, grid_size):\n    labels_2d = label_map.reshape(grid_size).astype(np.float32)\n    labels_smooth = median_filter(labels_2d, size=3)\n    return labels_smooth.flatten().astype(np.int32)\n\ndef analyze_clusters_v4(image, labels, grid_size):\n    \"\"\"\n    FIXED classification with TWO-STEP approach:\n    Step 1: Is it vegetation? (ExG > 0.05 OR TGI > 0.03) AND (S > 40)\n    Step 2: If vegetation, is it Green or Dead?\n    \"\"\"\n    h, w = image.shape[:2]\n    gh, gw = grid_size\n    \n    patch_h = h // gh\n    patch_w = w // gw\n    \n    hsv = cv2.cvtColor(image, cv2.COLOR_RGB2HSV)\n    veg_indices = compute_vegetation_indices(image)\n    \n    labels_2d = labels.reshape(gh, gw)\n    \n    cluster_info = {}\n    for cluster_id in range(CFG.n_clusters):\n        mask = labels_2d == cluster_id\n        if mask.sum() == 0:\n            continue\n        \n        h_vals, s_vals, v_vals = [], [], []\n        exg_vals, tgi_vals, ngrdi_vals = [], [], []\n        \n        for i in range(gh):\n            for j in range(gw):\n                if labels_2d[i, j] == cluster_id:\n                    y1, y2 = i * patch_h, (i + 1) * patch_h\n                    x1, x2 = j * patch_w, (j + 1) * patch_w\n                    \n                    patch = hsv[y1:y2, x1:x2]\n                    h_vals.append(patch[:, :, 0].mean())\n                    s_vals.append(patch[:, :, 1].mean())\n                    v_vals.append(patch[:, :, 2].mean())\n                    \n                    exg_vals.append(veg_indices['ExG'][y1:y2, x1:x2].mean())\n                    tgi_vals.append(veg_indices['TGI'][y1:y2, x1:x2].mean())\n                    ngrdi_vals.append(veg_indices['NGRDI'][y1:y2, x1:x2].mean())\n        \n        cluster_info[cluster_id] = {\n            'H': np.mean(h_vals),\n            'S': np.mean(s_vals),\n            'V': np.mean(v_vals),\n            'ExG': np.mean(exg_vals),\n            'TGI': np.mean(tgi_vals),\n            'NGRDI': np.mean(ngrdi_vals),\n            'count': mask.sum(),\n            'percent': mask.sum() / (gh * gw) * 100\n        }\n    \n    # ADAPTIVE THRESHOLDS using percentiles\n    # Collect ExG/TGI from all non-shadow clusters\n    veg_exg = [v['ExG'] for v in cluster_info.values() if v['V'] >= 70]\n    veg_tgi = [v['TGI'] for v in cluster_info.values() if v['V'] >= 70]\n    \n    if len(veg_exg) > 0:\n        exg_green_thr = np.percentile(veg_exg, 60)\n        exg_dead_thr = np.percentile(veg_exg, 35)\n        tgi_green_thr = np.percentile(veg_tgi, 60)\n        tgi_dead_thr = np.percentile(veg_tgi, 35)\n    else:\n        exg_green_thr, exg_dead_thr = 0.15, 0.10\n        tgi_green_thr, tgi_dead_thr = 0.08, 0.05\n    \n    print(\"\\n=== Enhanced Cluster Analysis v5 (Adaptive) ===\")\n    print(f\"Thresholds: ExG[green>{exg_green_thr:.3f}, dead<{exg_dead_thr:.3f}]\")\n    print(f\"            TGI[green>{tgi_green_thr:.3f}, dead<{tgi_dead_thr:.3f}]\")\n    print(\"-\" * 70)\n    \n    classifications = {}\n    for cid, vals in sorted(cluster_info.items(), key=lambda x: -x[1]['percent']):\n        h_val, s_val, v_val = vals['H'], vals['S'], vals['V']\n        exg, tgi, ngrdi = vals['ExG'], vals['TGI'], vals['NGRDI']\n        \n        # STEP 1: Shadow check\n        if v_val < 70:\n            label = \"SHADOW\"\n        # STEP 2: Is it vegetation? (ExG > 0.05 OR TGI > 0.03) AND (S > 40)\n        elif not (((exg > 0.05) or (tgi > 0.03)) and (s_val > 40)):\n            label = \"SOIL\"\n        else:\n            # STEP 3: GREEN vs DEAD using adaptive thresholds\n            # GREEN: High ExG AND high TGI (above 60th percentile)\n            if (exg >= exg_green_thr) and (tgi >= tgi_green_thr):\n                label = \"GREEN\"\n            # DEAD: Low ExG OR low TGI, plus yellowish hue\n            elif (20 <= h_val <= 55) and ((exg < exg_dead_thr) or (tgi < tgi_dead_thr) or (ngrdi < 0.03)):\n                label = \"DEAD\"\n            # CLOVER: Pink/purple\n            elif (h_val < 20 or h_val > 150) and s_val > 60:\n                label = \"CLOVER\"\n            # Borderline: use ExG to decide\n            elif exg >= 0.15:\n                label = \"GREEN\"\n            elif exg < 0.12 and h_val < 45:\n                label = \"DEAD\"\n            else:\n                label = \"GREEN\"  # Default for ambiguous\n        \n        classifications[cid] = label\n        cluster_info[cid]['label'] = label\n        \n        print(f\"  Cluster {cid}: [{label:6s}]\")\n        print(f\"    HSV: H={h_val:.1f}, S={s_val:.1f}, V={v_val:.1f}\")\n        print(f\"    VegIdx: ExG={exg:.3f}, TGI={tgi:.3f}, NGRDI={ngrdi:.3f}\")\n        print(f\"    Area: {vals['percent']:.1f}%\")\n    \n    return cluster_info, classifications\n\ndef visualize_segmentation_v4(image, labels, grid_size, classifications, veg_indices, save_path):\n    labels_2d = labels.reshape(grid_size)\n    \n    label_map = cv2.resize(\n        labels_2d.astype(np.float32), \n        (image.shape[1], image.shape[0]),\n        interpolation=cv2.INTER_NEAREST\n    ).astype(np.int32)\n    \n    # Color map with SIMPLE strings (no emojis)\n    color_map = {\n        \"GREEN\": [0, 200, 0],\n        \"DEAD\": [255, 200, 0],      # Yellow-orange for better visibility\n        \"CLOVER\": [255, 105, 180],\n        \"SHADOW\": [50, 50, 50],\n        \"SOIL\": [139, 90, 43],      # Brown for soil\n    }\n    \n    seg_colored = np.zeros_like(image)\n    for cluster_id, label in classifications.items():\n        mask = label_map == cluster_id\n        color = color_map.get(label, [128, 128, 128])\n        seg_colored[mask] = color\n    \n    fig, axes = plt.subplots(2, 3, figsize=(15, 10))\n    \n    axes[0, 0].imshow(image)\n    axes[0, 0].set_title(\"Original Image\", fontsize=12)\n    axes[0, 0].axis(\"off\")\n    \n    axes[0, 1].imshow(seg_colored)\n    axes[0, 1].set_title(f\"Segmentation ({CFG.n_clusters} clusters)\", fontsize=12)\n    axes[0, 1].axis(\"off\")\n    \n    overlay = cv2.addWeighted(image, 0.6, seg_colored, 0.4, 0)\n    axes[0, 2].imshow(overlay)\n    axes[0, 2].set_title(\"Overlay\", fontsize=12)\n    axes[0, 2].axis(\"off\")\n    \n    exg_resized = cv2.resize(veg_indices['ExG'], (image.shape[1], image.shape[0]))\n    axes[1, 0].imshow(exg_resized, cmap='RdYlGn', vmin=-0.5, vmax=0.5)\n    axes[1, 0].set_title(\"ExG (Excess Green)\", fontsize=12)\n    axes[1, 0].axis(\"off\")\n    \n    tgi_resized = cv2.resize(veg_indices['TGI'], (image.shape[1], image.shape[0]))\n    axes[1, 1].imshow(tgi_resized, cmap='RdYlGn', vmin=-0.3, vmax=0.3)\n    axes[1, 1].set_title(\"TGI (Triangular Greenness)\", fontsize=12)\n    axes[1, 1].axis(\"off\")\n    \n    # Create legend\n    legend_img = np.ones((image.shape[0], image.shape[1], 3), dtype=np.uint8) * 255\n    y_offset = 50\n    for name, color in color_map.items():\n        cv2.rectangle(legend_img, (20, y_offset-20), (60, y_offset+10), color, -1)\n        cv2.putText(legend_img, name, (70, y_offset), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0,0,0), 2)\n        y_offset += 50\n    axes[1, 2].imshow(legend_img)\n    axes[1, 2].set_title(\"Legend\", fontsize=12)\n    axes[1, 2].axis(\"off\")\n    \n    plt.tight_layout()\n    plt.savefig(save_path, dpi=150, bbox_inches='tight')\n    plt.close()\n    print(f\"\\nSaved visualization to {save_path}\")\n    \n    return label_map\n\ndef estimate_biomass_v4(cluster_info):\n    print(\"\\n\" + \"=\"*70)\n    print(\"BIOMASS ESTIMATION v4:\")\n    print(\"=\"*70)\n    \n    green_pct = 0\n    dead_pct = 0\n    clover_pct = 0\n    soil_pct = 0\n    shadow_pct = 0\n    \n    for cid, info in cluster_info.items():\n        pct = info['percent']\n        label = info.get('label', '')\n        \n        if label == \"GREEN\":\n            green_pct += pct\n        elif label == \"DEAD\":\n            dead_pct += pct\n        elif label == \"CLOVER\":\n            clover_pct += pct\n        elif label == \"SHADOW\":\n            shadow_pct += pct\n        else:\n            soil_pct += pct\n    \n    print(f\"\\n  RAW AREA:\")\n    print(f\"    GREEN: {green_pct:.1f}%\")\n    print(f\"    DEAD: {dead_pct:.1f}%\")\n    print(f\"    CLOVER: {clover_pct:.1f}%\")\n    print(f\"    SOIL: {soil_pct:.1f}%\")\n    print(f\"    SHADOW: {shadow_pct:.1f}%\")\n    \n    biomass_total = green_pct + dead_pct + clover_pct\n    if biomass_total > 0:\n        print(f\"\\n  NORMALIZED (biomass only):\")\n        print(f\"    GREEN: {green_pct/biomass_total*100:.1f}%\")\n        print(f\"    DEAD: {dead_pct/biomass_total*100:.1f}%\")\n        print(f\"    CLOVER: {clover_pct/biomass_total*100:.1f}%\")\n    else:\n        print(\"\\n  No biomass detected (all soil/shadow)\")\n    \n    return {\n        'green': green_pct,\n        'dead': dead_pct,\n        'clover': clover_pct,\n        'soil': soil_pct,\n        'shadow': shadow_pct,\n    }\n\ndef main():\n    print(\"=\"*70)\n    print(\"DINOv2 Pseudo-Segmentation v4 - FIXED\")\n    print(\"=\"*70)\n    \n    if not os.path.exists(CFG.image_path):\n        train_dir = Path(\"/kaggle/input/csiro-biomass/train\")\n        if train_dir.exists():\n            images = list(train_dir.glob(\"*.jpg\"))\n            if images:\n                CFG.image_path = str(images[0])\n        else:\n            print(\"No images found\")\n            return\n    \n    print(f\"\\n1. Loading image: {CFG.image_path}\")\n    image = cv2.imread(CFG.image_path)\n    if image is None:\n        print(\"Could not load image\")\n        return\n    \n    h, w = image.shape[:2]\n    left_half = image[:, :w//2, :]\n    print(f\"   Image shape: {left_half.shape}\")\n    \n    print(\"\\n2. Loading DINOv2...\")\n    model = load_dino_backbone()\n    \n    print(\"\\n3. Extracting DINOv2 features...\")\n    dino_features, grid_size, img_rgb = extract_patch_features(model, left_half)\n    print(f\"   DINOv2: {dino_features.shape}, grid: {grid_size}\")\n    \n    print(\"\\n4. Computing vegetation indices...\")\n    veg_indices = compute_vegetation_indices(img_rgb)\n    \n    print(\"\\n5. Extracting enhanced features...\")\n    enhanced_features = extract_enhanced_features(img_rgb, grid_size)\n    print(f\"   Enhanced: {enhanced_features.shape}\")\n    \n    print(\"\\n6. Combining features...\")\n    scaler = StandardScaler()\n    dino_norm = scaler.fit_transform(dino_features)\n    enhanced_norm = StandardScaler().fit_transform(enhanced_features)\n    \n    combined = np.concatenate([\n        dino_norm * CFG.dino_weight,\n        enhanced_norm * (1 - CFG.dino_weight)\n    ], axis=1)\n    \n    print(f\"\\n7. Upsampling features ({CFG.upsample_factor}x)...\")\n    upsampled, new_grid = upsample_features(combined, grid_size, CFG.upsample_factor)\n    \n    print(f\"\\n8. Clustering (K={CFG.n_clusters})...\")\n    labels, kmeans = cluster_patches(upsampled, n_clusters=CFG.n_clusters)\n    \n    print(\"\\n9. Post-processing...\")\n    labels = post_process_segmentation(labels, new_grid)\n    \n    print(\"\\n10. Analyzing clusters (Two-Step)...\")\n    labels_small = cv2.resize(\n        labels.reshape(new_grid).astype(np.float32), \n        (grid_size[1], grid_size[0]),\n        interpolation=cv2.INTER_NEAREST\n    ).astype(np.int32).flatten()\n    \n    cluster_info, classifications = analyze_clusters_v4(img_rgb, labels_small, grid_size)\n    \n    print(\"\\n11. Visualizing...\")\n    save_path = \"/kaggle/working/segmentation_v4.png\"\n    visualize_segmentation_v4(img_rgb, labels, new_grid, classifications, veg_indices, save_path)\n    \n    print(\"\\n12. Estimating biomass...\")\n    biomass = estimate_biomass_v4(cluster_info)\n    import matplotlib.pyplot as plt\n    import matplotlib.image as mpimg\n\n    img_result = mpimg.imread('/kaggle/working/segmentation_v4.png')\n    plt.figure(figsize=(12, 8))\n    plt.imshow(img_result)\n    plt.axis('off')\n    plt.title(\"Segmentación de Biomasa V4\")\n    plt.show()\n    \n    return biomass, cluster_info\n\nif __name__ == \"__main__\":\n    main()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-01T00:38:19.952489Z","iopub.execute_input":"2026-01-01T00:38:19.954156Z","iopub.status.idle":"2026-01-01T00:38:47.296841Z","shell.execute_reply.started":"2026-01-01T00:38:19.954109Z","shell.execute_reply":"2026-01-01T00:38:47.295420Z"}},"outputs":[],"execution_count":null}]}