{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"},{"sourceId":9187072,"sourceType":"datasetVersion","datasetId":5504483}],"dockerImageVersionId":30761,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"The provided code demonstrates how to achieve this mapping by leveraging the positional metadata embedded in DICOM files, such as ImagePositionPatient and PixelSpacing. These values represent the physical location of each slice in a 3D coordinate system and allow us to relate the z-coordinate of an axial slice to the y-coordinate in a sagittal image. The code calculates the corresponding y-coordinate on the sagittal slice and then visualizes this mapping by drawing a red line at the appropriate position on the sagittal image.","metadata":{}},{"cell_type":"markdown","source":"predict point, show the point instead of drawing the line","metadata":{}},{"cell_type":"code","source":"import glob\nimport os\nimport numpy as np\nimport pydicom\nimport matplotlib.pyplot as plt\nfrom matplotlib.patches import Polygon\n\ndef map_axial_to_sagittal(axial_folder_path, sagittal_image_path, predict_point=False):\n    def convert_to_8bit(x):\n        lower, upper = np.percentile(x, (1, 99))\n        x = np.clip(x, lower, upper)\n        x = x - np.min(x)\n        x = x / np.max(x) \n        return (x * 255).astype(\"uint8\")\n\n    def load_dicom_stack(dicom_folder, plane, reverse_sort=False):\n        dicom_files = glob.glob(os.path.join(dicom_folder, \"*.dcm\"))\n        dicoms = [pydicom.dcmread(f) for f in dicom_files]\n        plane = {\"sagittal\": 0, \"coronal\": 1, \"axial\": 2}[plane.lower()]\n        positions = np.asarray([float(d.ImagePositionPatient[plane]) for d in dicoms])\n        idx = np.argsort(-positions if reverse_sort else positions)\n        ipp = np.asarray([d.ImagePositionPatient for d in dicoms]).astype(\"float\")[idx]\n        array = np.stack([d.pixel_array.astype(\"float32\") for d in dicoms])\n        array = array[idx]\n        iop = np.asarray(dicoms[0].ImageOrientationPatient).reshape(2, 3)\n        return {\"array\": convert_to_8bit(array), \"positions\": ipp, \"pixel_spacing\": np.asarray(dicoms[0].PixelSpacing).astype(\"float\"), \"orientation\": iop}\n\n    # Load axial and sagittal images\n    ax_t2 = load_dicom_stack(axial_folder_path, plane=\"axial\", reverse_sort=True)\n    sag_t2 = pydicom.dcmread(sagittal_image_path)\n    sag_array = convert_to_8bit(sag_t2.pixel_array)\n\n    # Get sagittal image position and pixel spacing\n    sag_position = np.array(sag_t2.ImagePositionPatient)\n    sag_pixel_spacing = np.array(sag_t2.PixelSpacing)\n    sag_orientation = np.array(sag_t2.ImageOrientationPatient).reshape(2, 3)\n\n    # Calculate sagittal image corners in patient space\n    sag_height, sag_width = sag_array.shape\n    sag_corners = np.array([\n        sag_position,\n        sag_position + sag_orientation[0] * sag_width * sag_pixel_spacing[0],\n        sag_position + sag_orientation[0] * sag_width * sag_pixel_spacing[0] + sag_orientation[1] * sag_height * sag_pixel_spacing[1],\n        sag_position + sag_orientation[1] * sag_height * sag_pixel_spacing[1]\n    ])\n\n    # Calculate axial slice corners and intersections with sagittal plane\n    intersections = []\n    for ax_pos in ax_t2[\"positions\"]:\n        ax_corners = np.array([\n            ax_pos,\n            ax_pos + ax_t2[\"orientation\"][0] * ax_t2[\"array\"].shape[2] * ax_t2[\"pixel_spacing\"][0],\n            ax_pos + ax_t2[\"orientation\"][0] * ax_t2[\"array\"].shape[2] * ax_t2[\"pixel_spacing\"][0] + \n                     ax_t2[\"orientation\"][1] * ax_t2[\"array\"].shape[1] * ax_t2[\"pixel_spacing\"][1],\n            ax_pos + ax_t2[\"orientation\"][1] * ax_t2[\"array\"].shape[1] * ax_t2[\"pixel_spacing\"][1]\n        ])\n        \n        # Calculate intersection of axial slice with sagittal plane\n        normal = np.cross(sag_orientation[0], sag_orientation[1])\n        d = -np.dot(normal, sag_position)\n        \n        intersection = []\n        for i in range(4):\n            p1, p2 = ax_corners[i], ax_corners[(i+1)%4]\n            t = -(d + np.dot(normal, p1)) / np.dot(normal, p2 - p1)\n            if 0 <= t <= 1:\n                point = p1 + t * (p2 - p1)\n                intersection.append(point)\n        \n        if len(intersection) == 2:\n            intersections.append(intersection)\n\n    # Convert intersection points to sagittal image coordinates\n    sag_intersections = []\n    for inter in intersections:\n        sag_inter = []\n        for point in inter:\n            v = point - sag_position\n            x = np.dot(v, sag_orientation[0]) / sag_pixel_spacing[0]\n            y = np.dot(v, sag_orientation[1]) / sag_pixel_spacing[1]\n            sag_inter.append([x, y])\n        sag_intersections.append(sag_inter)\n\n    # Visualize the mapping\n    plt.figure(figsize=(12, 6))\n    plt.subplot(1, 2, 1)\n    plt.imshow(sag_array, cmap=\"gray\")\n\n    if predict_point:\n        # Predict the center point of each intersection\n        for inter in sag_intersections:\n            x = sum(p[0] for p in inter) / len(inter)\n            y = sum(p[1] for p in inter) / len(inter)\n            plt.plot(x, y, 'ro', markersize=3)\n        plt.title(\"Sagittal View with Predicted Axial Slice Centers\")\n    else:\n        for inter in sag_intersections:\n            x = [p[0] for p in inter]\n            y = [p[1] for p in inter]\n            plt.plot(x, y, color=\"red\", linestyle=\"--\", alpha=0.5)\n        plt.title(\"Sagittal View with Axial Slice Intersections\")\n\n    plt.xlim(0, sag_width)\n    plt.ylim(sag_height, 0)\n\n    plt.subplot(1, 2, 2)\n    plt.imshow(ax_t2[\"array\"][len(ax_t2[\"array\"]) // 2], cmap=\"gray\")\n    plt.title(\"Middle Axial Slice\")\n\n    plt.tight_layout()\n    plt.show()\n\n    return sag_intersections\n\naxial_folder_path = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/1057173941/3608799236'\nsagittal_image_path = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/1057173941/2918184860/9.dcm'\nintersections = map_axial_to_sagittal(axial_folder_path, sagittal_image_path, predict_point=False)","metadata":{"execution":{"iopub.status.busy":"2024-08-31T20:04:04.072221Z","iopub.execute_input":"2024-08-31T20:04:04.073424Z","iopub.status.idle":"2024-08-31T20:04:05.904935Z","shell.execute_reply.started":"2024-08-31T20:04:04.073359Z","shell.execute_reply":"2024-08-31T20:04:05.902997Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For mapping the sagittal to axial","metadata":{}},{"cell_type":"code","source":"import glob\nimport os\nimport numpy as np\nimport pydicom\nimport matplotlib.pyplot as plt\nfrom matplotlib.patches import Polygon\n\ndef map_sagittal_to_axial(sagittal_folder_path, axial_image_path, predict_point=False):\n    def convert_to_8bit(x):\n        lower, upper = np.percentile(x, (1, 99))\n        x = np.clip(x, lower, upper)\n        x = x - np.min(x)\n        x = x / np.max(x) \n        return (x * 255).astype(\"uint8\")\n\n    def load_dicom_stack(dicom_folder, plane):\n        dicom_files = glob.glob(os.path.join(dicom_folder, \"*.dcm\"))\n        dicoms = [pydicom.dcmread(f) for f in dicom_files]\n        plane = {\"sagittal\": 0, \"coronal\": 1, \"axial\": 2}[plane.lower()]\n        positions = np.asarray([float(d.ImagePositionPatient[plane]) for d in dicoms])\n        idx = np.argsort(positions)\n        ipp = np.asarray([d.ImagePositionPatient for d in dicoms]).astype(\"float\")[idx]\n        array = np.stack([d.pixel_array.astype(\"float32\") for d in dicoms])\n        array = array[idx]\n        iop = np.asarray(dicoms[0].ImageOrientationPatient).reshape(2, 3)\n        return {\"array\": convert_to_8bit(array), \"positions\": ipp, \"pixel_spacing\": np.asarray(dicoms[0].PixelSpacing).astype(\"float\"), \"orientation\": iop}\n\n    # Load sagittal and axial images\n    sag_t2 = load_dicom_stack(sagittal_folder_path, plane=\"sagittal\")\n    ax_t2 = pydicom.dcmread(axial_image_path)\n    ax_array = convert_to_8bit(ax_t2.pixel_array)\n\n    # Get axial image position and pixel spacing\n    ax_position = np.array(ax_t2.ImagePositionPatient)\n    ax_pixel_spacing = np.array(ax_t2.PixelSpacing)\n    ax_orientation = np.array(ax_t2.ImageOrientationPatient).reshape(2, 3)\n\n    # Calculate axial image corners in patient space\n    ax_height, ax_width = ax_array.shape\n    ax_corners = np.array([\n        ax_position,\n        ax_position + ax_orientation[0] * ax_width * ax_pixel_spacing[0],\n        ax_position + ax_orientation[0] * ax_width * ax_pixel_spacing[0] + ax_orientation[1] * ax_height * ax_pixel_spacing[1],\n        ax_position + ax_orientation[1] * ax_height * ax_pixel_spacing[1]\n    ])\n\n    # Calculate sagittal slice corners and intersections with axial plane\n    intersections = []\n    for sag_pos in sag_t2[\"positions\"]:\n        sag_corners = np.array([\n            sag_pos,\n            sag_pos + sag_t2[\"orientation\"][0] * sag_t2[\"array\"].shape[2] * sag_t2[\"pixel_spacing\"][0],\n            sag_pos + sag_t2[\"orientation\"][0] * sag_t2[\"array\"].shape[2] * sag_t2[\"pixel_spacing\"][0] + \n                      sag_t2[\"orientation\"][1] * sag_t2[\"array\"].shape[1] * sag_t2[\"pixel_spacing\"][1],\n            sag_pos + sag_t2[\"orientation\"][1] * sag_t2[\"array\"].shape[1] * sag_t2[\"pixel_spacing\"][1]\n        ])\n        \n        # Calculate intersection of sagittal slice with axial plane\n        normal = np.cross(ax_orientation[0], ax_orientation[1])\n        d = -np.dot(normal, ax_position)\n        \n        intersection = []\n        for i in range(4):\n            p1, p2 = sag_corners[i], sag_corners[(i+1)%4]\n            t = -(d + np.dot(normal, p1)) / np.dot(normal, p2 - p1)\n            if 0 <= t <= 1:\n                point = p1 + t * (p2 - p1)\n                intersection.append(point)\n        \n        if len(intersection) == 2:\n            intersections.append(intersection)\n\n    # Convert intersection points to axial image coordinates\n    ax_intersections = []\n    for inter in intersections:\n        ax_inter = []\n        for point in inter:\n            v = point - ax_position\n            x = np.dot(v, ax_orientation[0]) / ax_pixel_spacing[0]\n            y = np.dot(v, ax_orientation[1]) / ax_pixel_spacing[1]\n            ax_inter.append([x, y])\n        ax_intersections.append(ax_inter)\n\n    # Visualize the mapping\n    plt.figure(figsize=(12, 6))\n    plt.subplot(1, 2, 1)\n    plt.imshow(ax_array, cmap=\"gray\")\n\n    if predict_point:\n        # Predict the center point of each intersection\n        for inter in ax_intersections:\n            x = sum(p[0] for p in inter) / len(inter)\n            y = sum(p[1] for p in inter) / len(inter)\n            plt.plot(x, y, 'ro', markersize=3)\n        plt.title(\"Axial View with Predicted Sagittal Slice Centers\")\n    else:\n        for inter in ax_intersections:\n            x = [p[0] for p in inter]\n            y = [p[1] for p in inter]\n            plt.plot(x, y, color=\"red\", linestyle=\"--\", alpha=0.5)\n        plt.title(\"Axial View with Sagittal Slice Intersections\")\n\n    plt.xlim(0, ax_width)\n    plt.ylim(ax_height, 0)\n\n    plt.subplot(1, 2, 2)\n    plt.imshow(sag_t2[\"array\"][len(sag_t2[\"array\"]) // 2], cmap=\"gray\")\n    plt.title(\"Middle Sagittal Slice\")\n\n    plt.tight_layout()\n    plt.show()\n\n    return ax_intersections\n\nsagittal_folder_path = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/4003253/1054713880'\naxial_image_path = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/4003253/2448190387/11.dcm'\nintersections = map_sagittal_to_axial(sagittal_folder_path, axial_image_path, predict_point=False)","metadata":{"execution":{"iopub.status.busy":"2024-08-31T20:03:46.386307Z","iopub.execute_input":"2024-08-31T20:03:46.386813Z","iopub.status.idle":"2024-08-31T20:03:47.686997Z","shell.execute_reply.started":"2024-08-31T20:03:46.386766Z","shell.execute_reply":"2024-08-31T20:03:47.685757Z"},"trusted":true},"execution_count":null,"outputs":[]}]}