{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.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":71549,"databundleVersionId":8561470,"sourceType":"competition"}],"dockerImageVersionId":30732,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import glob\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport numpy as np\nimport os\nimport pandas as pd\nimport pydicom","metadata":{"execution":{"iopub.status.busy":"2024-06-10T11:12:27.109945Z","iopub.execute_input":"2024-06-10T11:12:27.110347Z","iopub.status.idle":"2024-06-10T11:12:28.680630Z","shell.execute_reply.started":"2024-06-10T11:12:27.110314Z","shell.execute_reply":"2024-06-10T11:12:28.678992Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"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\ndef 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    # if reverse_sort=False, then increasing array index will be from RIGHT->LEFT and CAUDAL->CRANIAL\n    # thus we do reverse_sort=True for axial so increasing array index is craniocaudal\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    return {\"array\": convert_to_8bit(array), \"positions\": ipp, \"pixel_spacing\": np.asarray(dicoms[0].PixelSpacing).astype(\"float\")}","metadata":{"execution":{"iopub.status.busy":"2024-06-10T11:12:28.683662Z","iopub.execute_input":"2024-06-10T11:12:28.687261Z","iopub.status.idle":"2024-06-10T11:12:28.721698Z","shell.execute_reply.started":"2024-06-10T11:12:28.687195Z","shell.execute_reply":"2024-06-10T11:12:28.719998Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv(\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_series_descriptions.csv\")\nstudy = df.loc[df.study_id == df.study_id.iloc[0]]\nstudy","metadata":{"execution":{"iopub.status.busy":"2024-06-10T11:12:28.732022Z","iopub.execute_input":"2024-06-10T11:12:28.732959Z","iopub.status.idle":"2024-06-10T11:12:28.779185Z","shell.execute_reply.started":"2024-06-10T11:12:28.732914Z","shell.execute_reply":"2024-06-10T11:12:28.777963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_dir = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/\"\n\nfor row in study.itertuples():\n    if row.series_description == \"Sagittal T2/STIR\":\n        sag_t2 = load_dicom_stack(os.path.join(image_dir, str(row.study_id), str(row.series_id)), plane=\"sagittal\")\n    elif row.series_description == \"Sagittal T1\":\n        sag_t1 = load_dicom_stack(os.path.join(image_dir, str(row.study_id), str(row.series_id)), plane=\"sagittal\")\n    elif row.series_description == \"Axial T2\":\n        ax_t2 = load_dicom_stack(os.path.join(image_dir, str(row.study_id), str(row.series_id)), plane=\"axial\", reverse_sort=True)","metadata":{"execution":{"iopub.status.busy":"2024-06-10T11:12:28.780748Z","iopub.execute_input":"2024-06-10T11:12:28.781145Z","iopub.status.idle":"2024-06-10T11:12:30.563541Z","shell.execute_reply.started":"2024-06-10T11:12:28.781107Z","shell.execute_reply":"2024-06-10T11:12:30.562265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.subplot(1, 3, 1)\nplt.imshow(sag_t2[\"array\"][len(sag_t2[\"array\"]) // 2], cmap=\"gray\")\nplt.subplot(1, 3, 2)\nplt.imshow(sag_t1[\"array\"][len(sag_t1[\"array\"]) // 2], cmap=\"gray\")\nplt.subplot(1, 3, 3)\nplt.imshow(ax_t2[\"array\"][len(ax_t2[\"array\"]) // 2], cmap=\"gray\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-10T11:12:30.565038Z","iopub.execute_input":"2024-06-10T11:12:30.565468Z","iopub.status.idle":"2024-06-10T11:12:31.121923Z","shell.execute_reply.started":"2024-06-10T11:12:30.565423Z","shell.execute_reply":"2024-06-10T11:12:31.120595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We make use of the following information:\n\n1. ImagePositionPatient specifies the xyz-coordinate (in mm) of the TOP LEFTHAND corner of the image.\n2. PixelSpacing specifies the size (in mm) of individual pixels in the image. \n3. By convention, coordinate values INCREASE (i.e., become more positive) from RIGHT to LEFT and TOE to HEAD. \n\nIn the \"world space\" (3D space measured in real-world millimeters) for lumbar spine studies, each axial slice is along the z-axis (ImagePositionPatient[2]), each coronal (front to back) slice is on the y-axis (ImagePositionPatient[1]), and each sagittal slice is along the x-axis (ImagePositionPatient[0]). \n\nThus the world space z-axis corresponds to the y-axis on a sagittal IMAGE. If we want to map the axial slices to the sagittal slices, we are essentially looking for the y-axis coordinate that each axial slice occupies in the sagittal image.","metadata":{}},{"cell_type":"code","source":"top_left_hand_corner_sag_t2 = sag_t2[\"positions\"][len(sag_t2[\"array\"]) // 2]\ntop_left_hand_corner_sag_t2","metadata":{"execution":{"iopub.status.busy":"2024-06-10T11:12:31.123457Z","iopub.execute_input":"2024-06-10T11:12:31.123863Z","iopub.status.idle":"2024-06-10T11:12:31.133582Z","shell.execute_reply.started":"2024-06-10T11:12:31.123827Z","shell.execute_reply":"2024-06-10T11:12:31.132255Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is the coordinate of the top left hand corner of a sagittal T2 image. \n\nRemember, ImagePositionPatient2 (i.e., the 3rd coordinate) is the z-axis in world space which is the y-axis in a sagittal image.\n\nThus y=0 in pixel space in the sagittal image is z=-282.82 in world space. \n\nNow, we essentially want to map each y-coordinate in pixel space to its z-coordinate in world space.\n\nEach pixel has its own dimensions in world space, using the PixelSpacing attribute in the DICOM metadata. So for each increment of 1 on the y-axis in pixel space we can use an increment using the PixelSpacing. ","metadata":{}},{"cell_type":"code","source":"sag_y_axis_to_pixel_space = [top_left_hand_corner_sag_t2[2]]\nwhile len(sag_y_axis_to_pixel_space) < sag_t2[\"array\"].shape[1]: \n    sag_y_axis_to_pixel_space.append(sag_y_axis_to_pixel_space[-1] - sag_t2[\"pixel_spacing\"][1])\n    # Coordinates DECREASE (i.e. become more negative) from HEAD to TOE","metadata":{"execution":{"iopub.status.busy":"2024-06-10T11:12:31.135374Z","iopub.execute_input":"2024-06-10T11:12:31.135828Z","iopub.status.idle":"2024-06-10T11:12:31.150431Z","shell.execute_reply.started":"2024-06-10T11:12:31.135789Z","shell.execute_reply":"2024-06-10T11:12:31.149334Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we have a corresponding z-coordinate in world space for each y-coordinate in the sagittal plane, so we can map the z-coordinates in world space from the axial images (ImagePositionPatient[2]) to the sagittal image y-coordinates. ","metadata":{}},{"cell_type":"code","source":"sag_y_coord_to_axial_slice = {}\nfor ax_t2_slice, ax_t2_pos in zip(ax_t2[\"array\"], ax_t2[\"positions\"]):\n    diffs = np.abs(np.asarray(sag_y_axis_to_pixel_space) - ax_t2_pos[2])\n    sag_y_coord = np.argmin(diffs)\n    sag_y_coord_to_axial_slice[sag_y_coord] = ax_t2_slice","metadata":{"execution":{"iopub.status.busy":"2024-06-10T11:12:31.152034Z","iopub.execute_input":"2024-06-10T11:12:31.152507Z","iopub.status.idle":"2024-06-10T11:12:31.167060Z","shell.execute_reply.started":"2024-06-10T11:12:31.152469Z","shell.execute_reply":"2024-06-10T11:12:31.165827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sag_midline_slice = sag_t2[\"array\"][len(sag_t2[\"array\"]) // 2]\nplt.imshow(sag_midline_slice, cmap=\"gray\")\nfor k in [*sag_y_coord_to_axial_slice]: \n    plt.axhline(y=k, color=\"red\", linestyle=\"--\")\n\nplt.show()\n\nfor k, v in sag_y_coord_to_axial_slice.items():\n    plt.subplot(1, 2, 1)\n    plt.imshow(sag_midline_slice, cmap=\"gray\")\n    plt.axhline(y=k, color=\"red\", linestyle=\"--\")\n    plt.subplot(1, 2, 2)\n    plt.imshow(v, cmap=\"gray\")\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-10T11:12:31.170826Z","iopub.execute_input":"2024-06-10T11:12:31.171290Z","iopub.status.idle":"2024-06-10T11:12:50.674602Z","shell.execute_reply.started":"2024-06-10T11:12:31.171250Z","shell.execute_reply":"2024-06-10T11:12:50.673524Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we can do something similar for mapping sagittal images to axial images.","metadata":{}},{"cell_type":"code","source":"top_left_hand_corner_ax_t2 = ax_t2[\"positions\"][len(ax_t2[\"array\"]) // 2]\ntop_left_hand_corner_ax_t2\n\nax_x_axis_to_pixel_space = [top_left_hand_corner_ax_t2[0]]\nwhile len(ax_x_axis_to_pixel_space) < ax_t2[\"array\"].shape[2]: \n    ax_x_axis_to_pixel_space.append(ax_x_axis_to_pixel_space[-1] + ax_t2[\"pixel_spacing\"][0])\n    # **NOTE: Coordinates INCREASE (i.e. become more positive) from RIGHT to LEFT\n\nax_x_coord_to_sag_slice = {}\nfor sag_t2_slice, sag_t2_pos in zip(sag_t2[\"array\"], sag_t2[\"positions\"]):\n    diffs = np.abs(np.asarray(ax_x_axis_to_pixel_space) - sag_t2_pos[0])\n    ax_x_coord = np.argmin(diffs)\n    ax_x_coord_to_sag_slice[ax_x_coord] = sag_t2_slice\n    \nmid_ax_slice = ax_t2[\"array\"][len(ax_t2[\"array\"]) // 2]\nplt.imshow(mid_ax_slice, cmap=\"gray\")\nfor k in [*ax_x_coord_to_sag_slice]: \n    plt.axvline(x=k, color=\"red\", linestyle=\"--\")\n\nplt.show()\n\nfor k, v in ax_x_coord_to_sag_slice.items():\n    plt.subplot(1, 2, 1)\n    plt.imshow(mid_ax_slice, cmap=\"gray\")\n    plt.axvline(x=k, color=\"red\", linestyle=\"--\")\n    plt.subplot(1, 2, 2)\n    plt.imshow(v, cmap=\"gray\")\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-10T11:12:50.676117Z","iopub.execute_input":"2024-06-10T11:12:50.676549Z","iopub.status.idle":"2024-06-10T11:12:57.786831Z","shell.execute_reply.started":"2024-06-10T11:12:50.676509Z","shell.execute_reply":"2024-06-10T11:12:57.785672Z"},"trusted":true},"execution_count":null,"outputs":[]}]}