{"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":"import numpy as np\nimport pydicom\nimport os\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nfrom skimage import measure\nimport plotly.graph_objects as go\nfrom scipy import ndimage\n\n# Dataset path\nBASE_PATH = '/kaggle/input/rsna-2022-cervical-spine-fracture-detection'\nTRAIN_IMAGES_PATH = os.path.join(BASE_PATH, 'train_images')\n\ndef load_dicom_series(patient_id):\n    \"\"\"Load all 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    # Get all DICOM files\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    print(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    # Sort by ImagePositionPatient (Z coordinate)\n    slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))\n    \n    return slices\n\ndef get_pixels_hu(slices):\n    \"\"\"Convert DICOM pixel data to Hounsfield Units\"\"\"\n    # Stack all slices into 3D array\n    image = np.stack([s.pixel_array for s in slices])\n    \n    # Convert to int16\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 (HU)\n    for slice_number in range(len(slices)):\n        intercept = slices[slice_number].RescaleIntercept\n        slope = 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 resample_volume(image, slices, new_spacing=[1, 1, 1]):\n    \"\"\"Resample volume to isotropic voxels\"\"\"\n    # Get current spacing\n    spacing = np.array([\n        slices[0].SliceThickness,\n        slices[0].PixelSpacing[0],\n        slices[0].PixelSpacing[1]\n    ], dtype=np.float32)\n    \n    resize_factor = spacing / new_spacing\n    new_shape = np.round(image.shape * resize_factor)\n    \n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    \n    image = ndimage.zoom(image, real_resize_factor, mode='nearest')\n    \n    return image, new_spacing\n\ndef render_3d_volume_mip(volume, title=\"Maximum Intensity Projection\"):\n    \"\"\"Create Maximum Intensity Projection views\"\"\"\n    fig, axes = plt.subplots(2, 2, figsize=(12, 12))\n    fig.suptitle(title, fontsize=16)\n    \n    # Axial view (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)')\n    axes[0, 0].axis('off')\n    \n    # Sagittal view (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)')\n    axes[0, 1].axis('off')\n    \n    # Coronal view (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)')\n    axes[1, 0].axis('off')\n    \n    # 3D visualization info\n    axes[1, 1].text(0.5, 0.5, f'Volume Shape: {volume.shape}\\nHU Range: [{volume.min()}, {volume.max()}]',\n                    ha='center', va='center', fontsize=12)\n    axes[1, 1].axis('off')\n    \n    plt.tight_layout()\n    plt.show()\n\ndef render_3d_surface(volume, threshold=-300, title=\"3D Surface Rendering\"):\n    \"\"\"Create 3D surface rendering using marching cubes\"\"\"\n    print(f\"Creating 3D surface with threshold: {threshold} HU...\")\n    \n    # Use marching cubes to obtain the surface mesh\n    verts, faces, normals, values = measure.marching_cubes(volume, threshold)\n    \n    # Create 3D plot\n    fig = plt.figure(figsize=(12, 10))\n    ax = fig.add_subplot(111, projection='3d')\n    \n    # Create mesh\n    mesh = Poly3DCollection(verts[faces], alpha=0.7, linewidth=0)\n    face_color = [0.8, 0.8, 1]\n    mesh.set_facecolor(face_color)\n    ax.add_collection3d(mesh)\n    \n    # Set plot limits\n    ax.set_xlim(0, volume.shape[0])\n    ax.set_ylim(0, volume.shape[1])\n    ax.set_zlim(0, volume.shape[2])\n    \n    ax.set_xlabel('X')\n    ax.set_ylabel('Y')\n    ax.set_zlabel('Z')\n    ax.set_title(title)\n    \n    plt.tight_layout()\n    plt.show()\n\ndef render_interactive_3d(volume, threshold=-300):\n    \"\"\"Create interactive 3D rendering using Plotly\"\"\"\n    print(f\"Creating interactive 3D visualization with threshold: {threshold} HU...\")\n    \n    # Use marching cubes\n    verts, faces, normals, values = measure.marching_cubes(volume, threshold)\n    \n    # Create mesh\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.5,\n            color='lightblue',\n            flatshading=True\n        )\n    ])\n    \n    fig.update_layout(\n        title='Interactive 3D Spine CT Visualization',\n        scene=dict(\n            xaxis_title='X',\n            yaxis_title='Y',\n            zaxis_title='Z',\n            aspectmode='data'\n        ),\n        width=900,\n        height=700\n    )\n    \n    fig.show()\n\ndef render_slices_montage(volume, num_slices=16):\n    \"\"\"Display a montage of slices\"\"\"\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=(15, 15))\n    axes = axes.flatten()\n    \n    for idx, slice_idx in enumerate(slice_indices):\n        axes[idx].imshow(volume[slice_idx], cmap='gray')\n        axes[idx].set_title(f'Slice {slice_idx}')\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)\n    plt.tight_layout()\n    plt.show()\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) * 255\n    return windowed.astype(np.uint8)\n\ndef main():\n    \"\"\"Main function to run 3D volumetric rendering\"\"\"\n    \n    # Example patient IDs from the dataset\n    # You can find these by listing the train_images directory\n    print(\"Available patients:\")\n    patients = sorted(os.listdir(TRAIN_IMAGES_PATH))[:5]  # Show first 5\n    for i, patient in enumerate(patients):\n        print(f\"{i+1}. {patient}\")\n    \n    # Select a patient (change this to any patient ID in your dataset)\n    patient_id = patients[0]  # Use first patient\n    print(f\"\\nProcessing patient: {patient_id}\")\n    \n    # Load DICOM series\n    slices = load_dicom_series(patient_id)\n    print(f\"Loaded {len(slices)} slices\")\n    \n    # Get pixel data in Hounsfield Units\n    volume = get_pixels_hu(slices)\n    print(f\"Volume shape: {volume.shape}\")\n    print(f\"HU range: [{volume.min()}, {volume.max()}]\")\n    \n    # Resample to isotropic voxels (optional, for better 3D rendering)\n    print(\"\\nResampling volume to isotropic voxels...\")\n    volume_resampled, new_spacing = resample_volume(volume, slices, [1, 1, 1])\n    print(f\"Resampled volume shape: {volume_resampled.shape}\")\n    \n    # Apply bone window\n    volume_windowed = apply_bone_window(volume_resampled)\n    \n    # 1. Maximum Intensity Projection\n    print(\"\\n1. Creating Maximum Intensity Projections...\")\n    render_3d_volume_mip(volume_resampled, \"Spine CT - MIP Views\")\n    \n    # 2. Slice montage\n    print(\"\\n2. Creating slice montage...\")\n    render_slices_montage(volume_resampled, num_slices=16)\n    \n    # 3. 3D Surface rendering (bones)\n    print(\"\\n3. Creating 3D surface rendering...\")\n    # Threshold for bone: typically > 200 HU\n    render_3d_surface(volume_resampled, threshold=200, title=\"3D Bone Rendering\")\n    \n    # 4. Interactive 3D rendering\n    print(\"\\n4. Creating interactive 3D visualization...\")\n    render_interactive_3d(volume_resampled, threshold=200)\n    \n    print(\"\\n✓ All visualizations complete!\")\n\nif __name__ == \"__main__\":\n    main()","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-12-11T17:54:07.519463Z","iopub.execute_input":"2025-12-11T17:54:07.519834Z","iopub.status.idle":"2025-12-11T17:54:56.368986Z","shell.execute_reply.started":"2025-12-11T17:54:07.519802Z","shell.execute_reply":"2025-12-11T17:54:56.367394Z"}},"outputs":[],"execution_count":null}]}