{"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":36363,"databundleVersionId":4050810,"sourceType":"competition"}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install pyvista","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T17:28:51.775443Z","iopub.execute_input":"2025-12-11T17:28:51.775757Z","iopub.status.idle":"2025-12-11T17:28:58.554685Z","shell.execute_reply.started":"2025-12-11T17:28:51.775733Z","shell.execute_reply":"2025-12-11T17:28:58.553565Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nEnhanced 3D Cervical Spine Visualization Suite\nCombines multiple visualization techniques with robust processing\nOptimized for RSNA Kaggle competition data\n\"\"\"\n\nimport numpy as np\nimport pydicom\nimport os\nimport gc\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nimport matplotlib.pyplot as plt\nfrom matplotlib import cm\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nfrom skimage import measure, morphology, exposure\nfrom scipy import ndimage as ndi\n\n# Try optional imports\ntry:\n    import plotly.graph_objects as go\n    PLOTLY_AVAILABLE = True\nexcept Exception:\n    PLOTLY_AVAILABLE = False\n    print(\"⚠️ Plotly not available - interactive 3D will be skipped\")\n\n# ==================== PARAMETERS ====================\nBASE_PATH = '/kaggle/input/rsna-2022-cervical-spine-fracture-detection'\nTRAIN_IMAGES_PATH = os.path.join(BASE_PATH, 'train_images')\n\n# Processing parameters\nPATIENT_INDEX = 0\nBONE_THRESHOLD = 200          # HU threshold for bone\nSOFT_TISSUE_THRESHOLD = -300  # HU threshold for soft tissue\nRESAMPLE = True               # Resample to isotropic voxels\nTARGET_SPACING = [1, 1, 1]    # Target voxel spacing in mm\nMAX_DIM = 256                 # Max dimension for marching cubes\nLARGEST_COMPONENT = True      # Keep only largest bone structure\nAPPLY_MORPHOLOGY = True       # Apply morphological operations\n\n# Visualization parameters\nSAVE_OUTPUTS = True\nOUTPUT_DIR = '/kaggle/working'\nVERBOSE = True\n# ===================================================\n\ndef log(msg):\n    \"\"\"Print with prefix if verbose\"\"\"\n    if VERBOSE:\n        print(f\"[3D-VIZ] {msg}\")\n\ndef load_dicom_series(patient_id):\n    \"\"\"Load and sort DICOM slices for a patient\"\"\"\n    patient_path = os.path.join(TRAIN_IMAGES_PATH, patient_id)\n    \n    if not os.path.exists(patient_path):\n        raise ValueError(f\"Patient folder not found: {patient_id}\")\n    \n    dicom_files = sorted([f for f in os.listdir(patient_path) if f.endswith('.dcm')])\n    \n    if not dicom_files:\n        raise ValueError(f\"No DICOM files found for patient: {patient_id}\")\n    \n    log(f\"Loading {len(dicom_files)} DICOM slices...\")\n    \n    # Load slices\n    slices = []\n    for dicom_file in dicom_files:\n        filepath = os.path.join(patient_path, dicom_file)\n        ds = pydicom.dcmread(filepath)\n        slices.append(ds)\n    \n    # Robust multi-method sorting\n    try:\n        slices.sort(key=lambda x: int(x.InstanceNumber))\n    except Exception:\n        try:\n            slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))\n        except Exception:\n            slices.sort(key=lambda x: float(x.SliceLocation))\n    \n    return slices\n\ndef get_pixels_hu(slices):\n    \"\"\"Convert DICOM pixel data to Hounsfield Units\"\"\"\n    image = np.stack([s.pixel_array for s in slices])\n    image = image.astype(np.int16)\n    \n    # Set outside-of-scan pixels to 0\n    image[image == -2000] = 0\n    \n    # Convert to Hounsfield units\n    for slice_number in range(len(slices)):\n        intercept = float(slices[slice_number].RescaleIntercept)\n        slope = float(slices[slice_number].RescaleSlope)\n        \n        if slope != 1:\n            image[slice_number] = slope * image[slice_number].astype(np.float64)\n            image[slice_number] = image[slice_number].astype(np.int16)\n        \n        image[slice_number] += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)\n\ndef get_spacing(slices):\n    \"\"\"Extract voxel spacing from DICOM headers\"\"\"\n    try:\n        px = slices[0].PixelSpacing\n        \n        # Calculate z-spacing from consecutive slices\n        if len(slices) > 1:\n            positions = [float(ds.ImagePositionPatient[2]) for ds in slices[:5]]\n            z_diffs = [abs(positions[i+1] - positions[i]) for i in range(len(positions)-1)]\n            z_spacing = np.median(z_diffs)\n        else:\n            z_spacing = float(slices[0].SliceThickness)\n        \n        spacing = [z_spacing, float(px[0]), float(px[1])]\n        return spacing\n    except Exception as e:\n        log(f\"Warning: Could not extract spacing ({e}), using defaults\")\n        return [1.0, 1.0, 1.0]\n\ndef resample_volume(image, original_spacing, new_spacing=[1, 1, 1]):\n    \"\"\"Resample volume to isotropic voxels\"\"\"\n    spacing = np.array(original_spacing, dtype=np.float32)\n    new_spacing = np.array(new_spacing, dtype=np.float32)\n    \n    resize_factor = spacing / new_spacing\n    new_shape = np.round(image.shape * resize_factor).astype(int)\n    \n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    \n    log(f\"Resampling from {image.shape} to {tuple(new_shape)}...\")\n    image = ndi.zoom(image, real_resize_factor, order=1)\n    \n    return image, new_spacing.tolist()\n\ndef apply_bone_window(volume):\n    \"\"\"Apply bone window settings (W:1800, L:400)\"\"\"\n    lower = 400 - 1800/2\n    upper = 400 + 1800/2\n    windowed = np.clip(volume, lower, upper)\n    windowed = (windowed - lower) / (upper - lower)\n    return windowed\n\ndef apply_soft_tissue_window(volume):\n    \"\"\"Apply soft tissue window settings (W:400, L:40)\"\"\"\n    lower = 40 - 400/2\n    upper = 40 + 400/2\n    windowed = np.clip(volume, lower, upper)\n    windowed = (windowed - lower) / (upper - lower)\n    return windowed\n\ndef create_bone_mask(hu_volume, threshold, apply_morphology=True):\n    \"\"\"Create binary bone mask with optional morphological cleanup\"\"\"\n    mask = hu_volume > threshold\n    \n    if apply_morphology:\n        # Aggressive morphological operations\n        mask = ndi.binary_closing(mask, structure=np.ones((3,3,3)), iterations=2)\n        mask = ndi.binary_opening(mask, structure=np.ones((3,3,3)), iterations=1)\n        \n        # Remove small objects\n        mask = morphology.remove_small_objects(mask, min_size=1000)\n    \n    return mask\n\ndef extract_largest_component(mask):\n    \"\"\"Keep only the largest connected component\"\"\"\n    labeled = measure.label(mask, connectivity=2)\n    if labeled.max() == 0:\n        return mask\n    \n    props = measure.regionprops(labeled)\n    largest = max(props, key=lambda x: x.area)\n    \n    log(f\"Kept largest component: {largest.area:,} voxels out of {mask.sum():,}\")\n    return labeled == largest.label\n\ndef downsample_mask(mask, target_max_dim):\n    \"\"\"Downsample mask to target dimension\"\"\"\n    max_dim = max(mask.shape)\n    if max_dim <= target_max_dim:\n        return mask, 1\n    \n    factor = int(np.ceil(max_dim / target_max_dim))\n    downsampled = mask[::factor, ::factor, ::factor]\n    \n    log(f\"Downsampled from {mask.shape} to {downsampled.shape} (factor={factor})\")\n    return downsampled, factor\n\n# ==================== VISUALIZATION FUNCTIONS ====================\n\ndef render_mip_views(volume, title=\"Maximum Intensity Projection\", save_path=None):\n    \"\"\"Create Maximum Intensity Projection views\"\"\"\n    log(\"Creating MIP views...\")\n    \n    fig, axes = plt.subplots(2, 2, figsize=(14, 14))\n    fig.suptitle(title, fontsize=16, fontweight='bold')\n    \n    # Axial MIP (top-down)\n    mip_axial = np.max(volume, axis=0)\n    axes[0, 0].imshow(mip_axial, cmap='gray')\n    axes[0, 0].set_title('Axial MIP (Top View)', fontsize=12)\n    axes[0, 0].axis('off')\n    \n    # Sagittal MIP (side)\n    mip_sagittal = np.max(volume, axis=2)\n    axes[0, 1].imshow(mip_sagittal, cmap='gray')\n    axes[0, 1].set_title('Sagittal MIP (Side View)', fontsize=12)\n    axes[0, 1].axis('off')\n    \n    # Coronal MIP (front)\n    mip_coronal = np.max(volume, axis=1)\n    axes[1, 0].imshow(mip_coronal, cmap='gray')\n    axes[1, 0].set_title('Coronal MIP (Front View)', fontsize=12)\n    axes[1, 0].axis('off')\n    \n    # Volume info\n    info_text = f'Volume Shape: {volume.shape}\\n'\n    info_text += f'HU Range: [{volume.min():.0f}, {volume.max():.0f}]\\n'\n    info_text += f'Mean HU: {volume.mean():.0f}\\n'\n    info_text += f'Std HU: {volume.std():.0f}'\n    axes[1, 1].text(0.5, 0.5, info_text, ha='center', va='center', \n                    fontsize=12, family='monospace',\n                    bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))\n    axes[1, 1].axis('off')\n    \n    plt.tight_layout()\n    \n    if save_path and SAVE_OUTPUTS:\n        plt.savefig(save_path, dpi=150, bbox_inches='tight')\n        log(f\"Saved MIP views: {save_path}\")\n    \n    plt.show()\n\ndef render_slice_montage(volume, num_slices=16, window_func=None, save_path=None):\n    \"\"\"Display a montage of slices\"\"\"\n    log(\"Creating slice montage...\")\n    \n    step = max(1, volume.shape[0] // num_slices)\n    slice_indices = range(0, volume.shape[0], step)[:num_slices]\n    \n    rows = int(np.ceil(np.sqrt(num_slices)))\n    cols = int(np.ceil(num_slices / rows))\n    \n    fig, axes = plt.subplots(rows, cols, figsize=(16, 16))\n    axes = axes.flatten()\n    \n    # Apply windowing if provided\n    display_volume = window_func(volume) if window_func else volume\n    \n    for idx, slice_idx in enumerate(slice_indices):\n        axes[idx].imshow(display_volume[slice_idx], cmap='gray')\n        axes[idx].set_title(f'Slice {slice_idx}/{volume.shape[0]}', fontsize=10)\n        axes[idx].axis('off')\n    \n    # Hide unused subplots\n    for idx in range(len(slice_indices), len(axes)):\n        axes[idx].axis('off')\n    \n    plt.suptitle('CT Slices Montage', fontsize=16, fontweight='bold')\n    plt.tight_layout()\n    \n    if save_path and SAVE_OUTPUTS:\n        plt.savefig(save_path, dpi=150, bbox_inches='tight')\n        log(f\"Saved slice montage: {save_path}\")\n    \n    plt.show()\n\ndef render_multiplanar(volume, title=\"Multiplanar Reconstruction\", save_path=None):\n    \"\"\"Display orthogonal slices at volume center\"\"\"\n    log(\"Creating multiplanar reconstruction...\")\n    \n    z, y, x = volume.shape\n    mid_z, mid_y, mid_x = z//2, y//2, x//2\n    \n    fig = plt.figure(figsize=(16, 12))\n    \n    # Axial slice\n    ax1 = plt.subplot(2, 3, 1)\n    ax1.imshow(volume[mid_z, :, :], cmap='bone')\n    ax1.set_title(f'Axial (Z={mid_z})', fontsize=12)\n    ax1.axhline(mid_y, color='r', linestyle='--', alpha=0.5)\n    ax1.axvline(mid_x, color='g', linestyle='--', alpha=0.5)\n    ax1.axis('off')\n    \n    # Sagittal slice\n    ax2 = plt.subplot(2, 3, 2)\n    ax2.imshow(volume[:, :, mid_x].T, cmap='bone', aspect='auto')\n    ax2.set_title(f'Sagittal (X={mid_x})', fontsize=12)\n    ax2.axhline(mid_y, color='r', linestyle='--', alpha=0.5)\n    ax2.axvline(mid_z, color='b', linestyle='--', alpha=0.5)\n    ax2.axis('off')\n    \n    # Coronal slice\n    ax3 = plt.subplot(2, 3, 3)\n    ax3.imshow(volume[:, mid_y, :].T, cmap='bone', aspect='auto')\n    ax3.set_title(f'Coronal (Y={mid_y})', fontsize=12)\n    ax3.axhline(mid_x, color='g', linestyle='--', alpha=0.5)\n    ax3.axvline(mid_z, color='b', linestyle='--', alpha=0.5)\n    ax3.axis('off')\n    \n    # Bone window views\n    bone_windowed = apply_bone_window(volume)\n    \n    ax4 = plt.subplot(2, 3, 4)\n    ax4.imshow(bone_windowed[mid_z, :, :], cmap='gray')\n    ax4.set_title('Axial (Bone Window)', fontsize=12)\n    ax4.axis('off')\n    \n    ax5 = plt.subplot(2, 3, 5)\n    ax5.imshow(bone_windowed[:, :, mid_x].T, cmap='gray', aspect='auto')\n    ax5.set_title('Sagittal (Bone Window)', fontsize=12)\n    ax5.axis('off')\n    \n    ax6 = plt.subplot(2, 3, 6)\n    ax6.imshow(bone_windowed[:, mid_y, :].T, cmap='gray', aspect='auto')\n    ax6.set_title('Coronal (Bone Window)', fontsize=12)\n    ax6.axis('off')\n    \n    plt.suptitle(title, fontsize=16, fontweight='bold')\n    plt.tight_layout()\n    \n    if save_path and SAVE_OUTPUTS:\n        plt.savefig(save_path, dpi=150, bbox_inches='tight')\n        log(f\"Saved multiplanar view: {save_path}\")\n    \n    plt.show()\n\ndef render_3d_surface_matplotlib(verts, faces, title=\"3D Surface\", save_path=None):\n    \"\"\"Create 3D surface rendering using matplotlib\"\"\"\n    log(\"Creating 3D surface with matplotlib...\")\n    \n    fig = plt.figure(figsize=(14, 12))\n    \n    # Main view\n    ax1 = fig.add_subplot(2, 2, 1, projection='3d')\n    \n    # Subsample faces if too many\n    max_faces = 30000\n    if len(faces) > max_faces:\n        indices = np.random.choice(len(faces), max_faces, replace=False)\n        faces_plot = faces[indices]\n    else:\n        faces_plot = faces\n    \n    mesh = Poly3DCollection(verts[faces_plot], alpha=0.8, linewidths=0.1, edgecolors='darkgray')\n    mesh.set_facecolor([0.8, 0.9, 1.0])\n    ax1.add_collection3d(mesh)\n    \n    ax1.set_xlim(verts[:, 0].min(), verts[:, 0].max())\n    ax1.set_ylim(verts[:, 1].min(), verts[:, 1].max())\n    ax1.set_zlim(verts[:, 2].min(), verts[:, 2].max())\n    ax1.set_xlabel('X (mm)', fontsize=10)\n    ax1.set_ylabel('Y (mm)', fontsize=10)\n    ax1.set_zlabel('Z (mm)', fontsize=10)\n    ax1.set_title('3D View', fontsize=12)\n    ax1.view_init(elev=20, azim=45)\n    \n    # Side view\n    ax2 = fig.add_subplot(2, 2, 2, projection='3d')\n    mesh2 = Poly3DCollection(verts[faces_plot], alpha=0.8, linewidths=0.1, edgecolors='darkgray')\n    mesh2.set_facecolor([0.8, 0.9, 1.0])\n    ax2.add_collection3d(mesh2)\n    ax2.set_xlim(verts[:, 0].min(), verts[:, 0].max())\n    ax2.set_ylim(verts[:, 1].min(), verts[:, 1].max())\n    ax2.set_zlim(verts[:, 2].min(), verts[:, 2].max())\n    ax2.set_xlabel('X (mm)', fontsize=10)\n    ax2.set_ylabel('Y (mm)', fontsize=10)\n    ax2.set_zlabel('Z (mm)', fontsize=10)\n    ax2.set_title('Side View', fontsize=12)\n    ax2.view_init(elev=0, azim=90)\n    \n    # Top view\n    ax3 = fig.add_subplot(2, 2, 3, projection='3d')\n    mesh3 = Poly3DCollection(verts[faces_plot], alpha=0.8, linewidths=0.1, edgecolors='darkgray')\n    mesh3.set_facecolor([0.8, 0.9, 1.0])\n    ax3.add_collection3d(mesh3)\n    ax3.set_xlim(verts[:, 0].min(), verts[:, 0].max())\n    ax3.set_ylim(verts[:, 1].min(), verts[:, 1].max())\n    ax3.set_zlim(verts[:, 2].min(), verts[:, 2].max())\n    ax3.set_xlabel('X (mm)', fontsize=10)\n    ax3.set_ylabel('Y (mm)', fontsize=10)\n    ax3.set_zlabel('Z (mm)', fontsize=10)\n    ax3.set_title('Top View', fontsize=12)\n    ax3.view_init(elev=90, azim=0)\n    \n    # Front view\n    ax4 = fig.add_subplot(2, 2, 4, projection='3d')\n    mesh4 = Poly3DCollection(verts[faces_plot], alpha=0.8, linewidths=0.1, edgecolors='darkgray')\n    mesh4.set_facecolor([0.8, 0.9, 1.0])\n    ax4.add_collection3d(mesh4)\n    ax4.set_xlim(verts[:, 0].min(), verts[:, 0].max())\n    ax4.set_ylim(verts[:, 1].min(), verts[:, 1].max())\n    ax4.set_zlim(verts[:, 2].min(), verts[:, 2].max())\n    ax4.set_xlabel('X (mm)', fontsize=10)\n    ax4.set_ylabel('Y (mm)', fontsize=10)\n    ax4.set_zlabel('Z (mm)', fontsize=10)\n    ax4.set_title('Front View', fontsize=12)\n    ax4.view_init(elev=0, azim=0)\n    \n    plt.suptitle(title, fontsize=16, fontweight='bold')\n    plt.tight_layout()\n    \n    if save_path and SAVE_OUTPUTS:\n        plt.savefig(save_path, dpi=150, bbox_inches='tight')\n        log(f\"Saved 3D surface: {save_path}\")\n    \n    plt.show()\n\ndef render_interactive_3d_plotly(verts, faces, title=\"Interactive 3D Spine\"):\n    \"\"\"Create interactive 3D rendering using Plotly\"\"\"\n    if not PLOTLY_AVAILABLE:\n        log(\"Plotly not available - skipping interactive 3D\")\n        return\n    \n    log(\"Creating interactive 3D visualization with Plotly...\")\n    \n    x, y, z = verts.T\n    i, j, k = faces.T\n    \n    fig = go.Figure(data=[\n        go.Mesh3d(\n            x=x, y=y, z=z,\n            i=i, j=j, k=k,\n            opacity=0.7,\n            color='lightblue',\n            flatshading=False,\n            lighting=dict(\n                ambient=0.5,\n                diffuse=0.8,\n                specular=0.5,\n                roughness=0.5,\n                fresnel=0.2\n            ),\n            lightposition=dict(\n                x=100,\n                y=200,\n                z=300\n            )\n        )\n    ])\n    \n    fig.update_layout(\n        title=title,\n        scene=dict(\n            xaxis_title='X (mm)',\n            yaxis_title='Y (mm)',\n            zaxis_title='Z (mm)',\n            aspectmode='data',\n            camera=dict(\n                eye=dict(x=1.5, y=1.5, z=1.5)\n            )\n        ),\n        width=1000,\n        height=800\n    )\n    \n    if SAVE_OUTPUTS:\n        html_path = os.path.join(OUTPUT_DIR, 'interactive_3d_spine.html')\n        fig.write_html(html_path)\n        log(f\"Saved interactive 3D: {html_path}\")\n    \n    fig.show()\n\ndef save_stl(verts, faces, filename):\n    \"\"\"Save mesh to STL format\"\"\"\n    try:\n        with open(filename, 'w') as f:\n            f.write('solid mesh\\n')\n            for face in faces:\n                v0, v1, v2 = verts[face]\n                normal = np.cross(v1 - v0, v2 - v0)\n                normal = normal / (np.linalg.norm(normal) + 1e-8)\n                \n                f.write(f'  facet normal {normal[0]:.6e} {normal[1]:.6e} {normal[2]:.6e}\\n')\n                f.write('    outer loop\\n')\n                f.write(f'      vertex {v0[0]:.6e} {v0[1]:.6e} {v0[2]:.6e}\\n')\n                f.write(f'      vertex {v1[0]:.6e} {v1[1]:.6e} {v1[2]:.6e}\\n')\n                f.write(f'      vertex {v2[0]:.6e} {v2[1]:.6e} {v2[2]:.6e}\\n')\n                f.write('    endloop\\n')\n                f.write('  endfacet\\n')\n            f.write('endsolid mesh\\n')\n        log(f\"Saved STL mesh: {filename}\")\n        return True\n    except Exception as e:\n        log(f\"Failed to save STL: {e}\")\n        return False\n\n# ==================== MAIN PIPELINE ====================\n\ndef main():\n    \"\"\"Main visualization pipeline\"\"\"\n    \n    log(\"=\"*60)\n    log(\"3D SPINE VISUALIZATION SUITE\")\n    log(\"=\"*60)\n    \n    # 1. List and select patient\n    log(\"\\nAvailable patients:\")\n    patients = sorted(os.listdir(TRAIN_IMAGES_PATH))[:10]\n    for i, patient in enumerate(patients[:5]):\n        print(f\"  {i+1}. {patient}\")\n    if len(patients) > 5:\n        print(f\"  ... and {len(patients)-5} more\")\n    \n    patient_id = patients[PATIENT_INDEX]\n    log(f\"\\nSelected patient: {patient_id} (index {PATIENT_INDEX})\")\n    \n    # 2. Load DICOM series\n    slices = load_dicom_series(patient_id)\n    log(f\"Loaded {len(slices)} slices\")\n    \n    # 3. Convert to Hounsfield Units\n    volume = get_pixels_hu(slices)\n    log(f\"Volume shape: {volume.shape}\")\n    log(f\"HU range: [{volume.min()}, {volume.max()}]\")\n    \n    # 4. Get spacing\n    original_spacing = get_spacing(slices)\n    log(f\"Original spacing (Z×Y×X): {original_spacing} mm\")\n    \n    # 5. Optional resampling\n    if RESAMPLE:\n        volume, new_spacing = resample_volume(volume, original_spacing, TARGET_SPACING)\n        spacing = new_spacing\n        log(f\"Resampled to spacing: {spacing} mm\")\n    else:\n        spacing = original_spacing\n    \n    # 6. VISUALIZATION 1: MIP Views\n    log(\"\\n\" + \"=\"*60)\n    log(\"VISUALIZATION 1: Maximum Intensity Projections\")\n    log(\"=\"*60)\n    render_mip_views(\n        volume, \n        title=f\"Patient {patient_id} - MIP Views\",\n        save_path=os.path.join(OUTPUT_DIR, 'mip_views.png')\n    )\n    \n    # 7. VISUALIZATION 2: Multiplanar Reconstruction\n    log(\"\\n\" + \"=\"*60)\n    log(\"VISUALIZATION 2: Multiplanar Reconstruction\")\n    log(\"=\"*60)\n    render_multiplanar(\n        volume,\n        title=f\"Patient {patient_id} - Multiplanar Views\",\n        save_path=os.path.join(OUTPUT_DIR, 'multiplanar.png')\n    )\n    \n    # 8. VISUALIZATION 3: Slice Montage\n    log(\"\\n\" + \"=\"*60)\n    log(\"VISUALIZATION 3: Slice Montage\")\n    log(\"=\"*60)\n    render_slice_montage(\n        volume, \n        num_slices=16,\n        window_func=apply_bone_window,\n        save_path=os.path.join(OUTPUT_DIR, 'slice_montage.png')\n    )\n    \n    # 9. Create bone mask\n    log(\"\\n\" + \"=\"*60)\n    log(\"MESH GENERATION: Creating bone surface\")\n    log(\"=\"*60)\n    bone_mask = create_bone_mask(volume, BONE_THRESHOLD, APPLY_MORPHOLOGY)\n    log(f\"Initial bone voxels: {bone_mask.sum():,}\")\n    \n    # Auto-adjust threshold if needed\n    if bone_mask.sum() < 10000:\n        log(f\"Low bone count, reducing threshold from {BONE_THRESHOLD} to 100\")\n        bone_mask = create_bone_mask(volume, 100, APPLY_MORPHOLOGY)\n        log(f\"Adjusted bone voxels: {bone_mask.sum():,}\")\n    \n    # 10. Extract largest component\n    if LARGEST_COMPONENT and bone_mask.sum() > 0:\n        bone_mask = extract_largest_component(bone_mask)\n    \n    # 11. Downsample and generate mesh\n    bone_mask_small, down_factor = downsample_mask(bone_mask, MAX_DIM)\n    final_spacing = [s * down_factor for s in spacing]\n    \n    log(\"Running marching cubes algorithm...\")\n    try:\n        verts, faces, normals, vals = measure.marching_cubes(\n            bone_mask_small.astype(np.uint8),\n            level=0.5,\n            spacing=final_spacing,\n            step_size=1\n        )\n        faces = faces.astype(np.int64)\n        log(f\"Mesh created: {len(verts):,} vertices, {len(faces):,} faces\")\n        \n        # 12. Save STL\n        if SAVE_OUTPUTS:\n            stl_path = os.path.join(OUTPUT_DIR, 'spine_mesh.stl')\n            save_stl(verts, faces, stl_path)\n        \n        # 13. VISUALIZATION 4: 3D Surface (Matplotlib)\n        log(\"\\n\" + \"=\"*60)\n        log(\"VISUALIZATION 4: 3D Surface Rendering (Matplotlib)\")\n        log(\"=\"*60)\n        render_3d_surface_matplotlib(\n            verts, faces,\n            title=f\"Patient {patient_id} - 3D Bone Surface\",\n            save_path=os.path.join(OUTPUT_DIR, '3d_surface.png')\n        )\n        \n        # 14. VISUALIZATION 5: Interactive 3D (Plotly)\n        log(\"\\n\" + \"=\"*60)\n        log(\"VISUALIZATION 5: Interactive 3D (Plotly)\")\n        log(\"=\"*60)\n        render_interactive_3d_plotly(verts, faces, f\"Patient {patient_id} - Interactive 3D\")\n        \n    except Exception as e:\n        log(f\"ERROR: Mesh generation failed: {e}\")\n    \n    # Cleanup\n    del volume, bone_mask, bone_mask_small\n    if 'verts' in locals():\n        del verts, faces, normals\n    gc.collect()\n    \n    log(\"\\n\" + \"=\"*60)\n    log(\"✓ ALL VISUALIZATIONS COMPLETE!\")\n    log(\"=\"*60)\n    if SAVE_OUTPUTS:\n        log(f\"Output files saved to: {OUTPUT_DIR}\")\n        log(\"  - mip_views.png\")\n        log(\"  - multiplanar.png\")\n        log(\"  - slice_montage.png\")\n        log(\"  - 3d_surface.png\")\n        log(\"  - spine_mesh.stl\")\n        if PLOTLY_AVAILABLE:\n            log(\"  - interactive_3d_spine.html\")\n\nif __name__ == \"__main__\":\n    main()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T17:29:29.829594Z","iopub.execute_input":"2025-12-11T17:29:29.829978Z","iopub.status.idle":"2025-12-11T17:30:23.089109Z","shell.execute_reply.started":"2025-12-11T17:29:29.829944Z","shell.execute_reply":"2025-12-11T17:30:23.087585Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}