{"metadata":{"kaggle":{"accelerator":"none","dataSources":[{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"}],"dockerImageVersionId":30761,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.14"},"papermill":{"default_parameters":{},"duration":35.69904,"end_time":"2024-06-10T11:14:46.501526","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-06-10T11:14:10.802486","version":"2.5.0"}},"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":{"papermill":{"duration":1.258583,"end_time":"2024-06-10T11:14:15.437301","exception":false,"start_time":"2024-06-10T11:14:14.178718","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-09-25T15:24:09.552998Z","iopub.execute_input":"2024-09-25T15:24:09.553438Z","iopub.status.idle":"2024-09-25T15:24:10.292805Z","shell.execute_reply.started":"2024-09-25T15:24:09.553393Z","shell.execute_reply":"2024-09-25T15:24:10.291524Z"},"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    iop = np.asarray([d.ImageOrientationPatient 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, \"orientations\": iop, \n            \"pixel_spacing\": np.asarray(dicoms[0].PixelSpacing).astype(\"float\")}\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    idx = np.argsort(-positions if reverse_sort else positions)\n    return [dicoms[i] for i in idx]","metadata":{"papermill":{"duration":0.021474,"end_time":"2024-06-10T11:14:15.464487","exception":false,"start_time":"2024-06-10T11:14:15.443013","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-09-25T15:24:10.295401Z","iopub.execute_input":"2024-09-25T15:24:10.296111Z","iopub.status.idle":"2024-09-25T15:24:10.313527Z","shell.execute_reply.started":"2024-09-25T15:24:10.296051Z","shell.execute_reply":"2024-09-25T15:24:10.312311Z"},"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\")\n#study = df.loc[df.study_id == df.study_id.iloc[0]]\nstudy = df.loc[df.study_id == 4003253]\n#study = df.loc[df.study_id == 1140848367]\nstudy","metadata":{"papermill":{"duration":0.056945,"end_time":"2024-06-10T11:14:15.526771","exception":false,"start_time":"2024-06-10T11:14:15.469826","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-09-25T15:24:10.315265Z","iopub.execute_input":"2024-09-25T15:24:10.315735Z","iopub.status.idle":"2024-09-25T15:24:10.391824Z","shell.execute_reply.started":"2024-09-25T15:24:10.315692Z","shell.execute_reply":"2024-09-25T15:24:10.390407Z"},"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        print (len(sag_t2))\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        print (len(sag_t1))\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)\n        print (len(ax_t2))","metadata":{"papermill":{"duration":1.821294,"end_time":"2024-06-10T11:14:17.353902","exception":false,"start_time":"2024-06-10T11:14:15.532608","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-09-25T15:24:10.395507Z","iopub.execute_input":"2024-09-25T15:24:10.396415Z","iopub.status.idle":"2024-09-25T15:24:11.061117Z","shell.execute_reply.started":"2024-09-25T15:24:10.396349Z","shell.execute_reply":"2024-09-25T15:24:11.059744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.subplot(1, 3, 1)\nplt.imshow(sag_t2[len(sag_t2) // 2].pixel_array, cmap=\"gray\")\nplt.subplot(1, 3, 2)\nplt.imshow(sag_t1[len(sag_t1) // 2].pixel_array, cmap=\"gray\")\nplt.subplot(1, 3, 3)\nplt.imshow(ax_t2[len(ax_t2) // 2].pixel_array, cmap=\"gray\")\nplt.show()","metadata":{"papermill":{"duration":0.563875,"end_time":"2024-06-10T11:14:17.923371","exception":false,"start_time":"2024-06-10T11:14:17.359496","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-09-25T15:24:11.062692Z","iopub.execute_input":"2024-09-25T15:24:11.063169Z","iopub.status.idle":"2024-09-25T15:24:11.762491Z","shell.execute_reply.started":"2024-09-25T15:24:11.063118Z","shell.execute_reply":"2024-09-25T15:24:11.761229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"parent_dir = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification\"\n\ntrain_imgs_dir = \"train_images\"\ntrain_descr_name = \"train_series_descriptions.csv\"\ntrain_name = \"train.csv\"\ntrain_xy_name = \"train_label_coordinates.csv\"\n\ndf_cords = pd.read_csv(os.path.join(parent_dir, train_xy_name))\ndf_descr = pd.read_csv(os.path.join(parent_dir, train_descr_name))\n\nstudy_cords = df_cords[df_cords.study_id == study.study_id.iloc[0]]\n\nstudy_cords = study_cords.merge(df_descr, on=[\"study_id\", \"series_id\"], how=\"left\")\n\nprint (study_cords.series_description.value_counts())\n\nstudy_cords.head()","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:11.764182Z","iopub.execute_input":"2024-09-25T15:24:11.764566Z","iopub.status.idle":"2024-09-25T15:24:11.951938Z","shell.execute_reply.started":"2024-09-25T15:24:11.764523Z","shell.execute_reply":"2024-09-25T15:24:11.950485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"study_cords[study_cords.series_description == 'Axial T2']","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:11.954036Z","iopub.execute_input":"2024-09-25T15:24:11.95445Z","iopub.status.idle":"2024-09-25T15:24:11.972849Z","shell.execute_reply.started":"2024-09-25T15:24:11.954406Z","shell.execute_reply":"2024-09-25T15:24:11.971594Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def point_to_plane_distance(point, plane_point, plane_normal):\n    # Calculate the distance from the point to the plane\n    point = np.array(point)\n    plane_point = np.array(plane_point)\n    plane_normal = np.array(plane_normal)\n    \n    d_plane = np.abs(np.dot(plane_normal, point - plane_point)) / np.linalg.norm(plane_normal)\n    \n    return d_plane\n\ndef project_point_onto_plane(point, plane_point, plane_normal):\n    # Project the point onto the plane\n    point = np.array(point)\n    plane_point = np.array(plane_point)\n    plane_normal = np.array(plane_normal)\n    \n    d_plane = np.dot(plane_normal, point - plane_point) / np.linalg.norm(plane_normal)**2\n    projected_point = point - d_plane * plane_normal\n    \n    return projected_point\n\ndef point_in_square(projected_point, square_vertices):\n    # Check if the projected point lies within the square\n    v1 = square_vertices[1] - square_vertices[0]\n    v2 = square_vertices[3] - square_vertices[0]\n    \n    vp = projected_point - square_vertices[0]\n    \n    dot11 = np.dot(v1, v1)\n    dot12 = np.dot(v1, v2)\n    dot22 = np.dot(v2, v2)\n    dot1p = np.dot(v1, vp)\n    dot2p = np.dot(v2, vp)\n    \n    inv_denom = 1 / (dot11 * dot22 - dot12 * dot12)\n    u = (dot22 * dot1p - dot12 * dot2p) * inv_denom\n    v = (dot11 * dot2p - dot12 * dot1p) * inv_denom\n    \n    return (u >= 0) and (v >= 0) and (u <= 1) and (v <= 1)\n\ndef distance_point_to_segment(point, v1, v2):\n    # Compute the distance from a point to a line segment\n    v1, v2, point = np.array(v1), np.array(v2), np.array(point)\n    segment_vector = v2 - v1\n    point_vector = point - v1\n    segment_length = np.linalg.norm(segment_vector)\n    \n    if segment_length == 0:\n        return np.linalg.norm(point_vector)\n    \n    t = max(0, min(1, np.dot(point_vector, segment_vector) / (segment_length**2)))\n    projection = v1 + t * segment_vector\n    \n    return np.linalg.norm(projection - point)\n\ndef point_to_square_distance(point, square_vertices):\n    # Ensure the vertices are numpy arrays\n    square_vertices = [np.array(v) for v in square_vertices]\n    \n    # Compute the plane normal\n    v1 = square_vertices[1] - square_vertices[0]\n    v2 = square_vertices[3] - square_vertices[0]\n    plane_normal = np.cross(v1, v2)\n    \n    # Calculate the distance to the plane\n    d_plane = point_to_plane_distance(point, square_vertices[0], plane_normal)\n    \n    # Project the point onto the plane\n    projected_point = project_point_onto_plane(point, square_vertices[0], plane_normal)\n    \n    # Check if the projected point is within the square\n    if point_in_square(projected_point, square_vertices):\n        return d_plane, projected_point\n    \n    return np.inf, None\n\n# Example usage\npoint = [0, 0, 10]\nsquare_vertices = [\n    [ 0,  0,  0],\n    [10,  0,  0],\n    [10, 10,  0],\n    [ 0, 10,  0]\n]\n\ndistance, _ = point_to_square_distance(point, square_vertices)\nprint(f\"Distance from point to square: {distance}\")","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:11.974925Z","iopub.execute_input":"2024-09-25T15:24:11.975432Z","iopub.status.idle":"2024-09-25T15:24:12.001368Z","shell.execute_reply.started":"2024-09-25T15:24:11.975362Z","shell.execute_reply":"2024-09-25T15:24:11.999893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ax_t2[1].InstanceNumber","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:12.003827Z","iopub.execute_input":"2024-09-25T15:24:12.004344Z","iopub.status.idle":"2024-09-25T15:24:12.027891Z","shell.execute_reply.started":"2024-09-25T15:24:12.004287Z","shell.execute_reply":"2024-09-25T15:24:12.026507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_corners_world_cords(dicom_file):\n    Oxyz = dicom_file.ImagePositionPatient\n    Xxyz = dicom_file.ImageOrientationPatient[:3]\n    Yxyz = dicom_file.ImageOrientationPatient[3:]\n    height = dicom_file.height if hasattr(dicom_file, 'height') else dicom_file.pixel_array.shape[0]\n    width = dicom_file.width if hasattr(dicom_file, 'width') else dicom_file.pixel_array.shape[1]\n    spX, spY = dicom_file.PixelSpacing\n    \n    Oxyz = np.asarray(Oxyz).astype(\"float\")\n    Xxyz = np.asarray(Xxyz).astype(\"float\")\n    Yxyz = np.asarray(Yxyz).astype(\"float\")\n    \n    # Return the 4 corners in world coordinates\n    return [Oxyz, \n            Oxyz + width * Xxyz * spX, \n            Oxyz + width * Xxyz * spX + height * Yxyz * spY, \n            Oxyz + height * Yxyz * spY]\n\ndef get_point_world_cords(dicom_file, x, y):\n    Oxyz = dicom_file.ImagePositionPatient\n    Xxyz = dicom_file.ImageOrientationPatient[:3]\n    Yxyz = dicom_file.ImageOrientationPatient[3:]\n    spX, spY = dicom_file.PixelSpacing\n    \n    Oxyz = np.asarray(Oxyz).astype(\"float\")\n    Xxyz = np.asarray(Xxyz).astype(\"float\")\n    Yxyz = np.asarray(Yxyz).astype(\"float\")\n    \n    # Calculate the world coordinates of the pixel (x, y)\n    return Oxyz + (x * Xxyz * spX) + (y * Yxyz * spY)\n\ndef wordl_cords2pixel_cords(dicom_file, point_xyz):\n    Oxyz = dicom_file.ImagePositionPatient\n    Xxyz = dicom_file.ImageOrientationPatient[:3]\n    Yxyz = dicom_file.ImageOrientationPatient[3:]\n    spX, spY = dicom_file.PixelSpacing\n    \n    Oxyz = np.asarray(Oxyz).astype(\"float\")\n    Xxyz = np.asarray(Xxyz).astype(\"float\")\n    Yxyz = np.asarray(Yxyz).astype(\"float\")\n    \n    diff = point_xyz - Oxyz\n    \n    # Convert world coordinates back to pixel coordinates\n    x = np.dot(diff, Xxyz) / spX\n    y = np.dot(diff, Yxyz) / spY\n    \n    return np.array([x, y])\n\ndef project_ax2sag(sagit, axial, x, y):\n    # Get the 3D world coordinates of the point on the axial plane\n    point3D = get_point_world_cords(axial, x, y)\n    \n    # Get the corners of the sagittal plane in 3D world coordinates\n    sagit_corners3D = get_corners_world_cords(sagit)\n    \n    # Project the 3D point from the axial plane onto the sagittal plane\n    distance, proj_point3D = point_to_square_distance(point3D, sagit_corners3D)\n    \n    if proj_point3D is not None:\n        # Convert the projected point's world coordinates back to pixel coordinates in the sagittal plane\n        return distance, wordl_cords2pixel_cords(sagit, proj_point3D)\n    else:\n        return np.inf, None\n\nx, y = study_cords.query(\"series_description == 'Axial T2' and instance_number == 3\")[[\"x\", \"y\"]].iloc[0]\n\nmy_sag = sag_t2[sag_t2.__len__() // 2]\nmy_ax = ax_t2[1]\n\ndistance, (new_x, new_y) = project_ax2sag(my_sag, my_ax, x, y)\n\ndistance, (new_x, new_y)","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:48.166736Z","iopub.execute_input":"2024-09-25T15:24:48.167631Z","iopub.status.idle":"2024-09-25T15:24:48.202708Z","shell.execute_reply.started":"2024-09-25T15:24:48.16758Z","shell.execute_reply":"2024-09-25T15:24:48.201242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(1, 2)\n\nfig.set_size_inches((20, 10))\n\nax[0].imshow(my_ax.pixel_array, cmap=plt.cm.gray)\nax[0].plot(x, y, 'x', markersize=15, markeredgewidth=3, color=\"green\")  \n\nax[1].imshow(my_sag.pixel_array, cmap=plt.cm.gray)\nax[1].plot(int(new_x), int(new_y), \"x\", markersize=15, markeredgewidth=3, color=\"green\") ","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:48.205294Z","iopub.execute_input":"2024-09-25T15:24:48.205966Z","iopub.status.idle":"2024-09-25T15:24:49.325407Z","shell.execute_reply.started":"2024-09-25T15:24:48.205896Z","shell.execute_reply":"2024-09-25T15:24:49.324006Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"study","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:49.327085Z","iopub.execute_input":"2024-09-25T15:24:49.327482Z","iopub.status.idle":"2024-09-25T15:24:49.339372Z","shell.execute_reply.started":"2024-09-25T15:24:49.327439Z","shell.execute_reply":"2024-09-25T15:24:49.33814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_dicom_dict(dicom_folder):\n    dicom_files = glob.glob(os.path.join(dicom_folder, \"*.dcm\"))\n    dicoms = [pydicom.dcmread(f) for f in dicom_files]\n    return {int(dicom.InstanceNumber):dicom for dicom in dicoms}\n\nimage_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_dict(os.path.join(image_dir, str(row.study_id), str(row.series_id)))\n        print (len(sag_t2))\n    elif row.series_description == \"Sagittal T1\":\n        sag_t1 = load_dicom_dict(os.path.join(image_dir, str(row.study_id), str(row.series_id)))\n        print (len(sag_t1))\n    elif row.series_description == \"Axial T2\":\n        ax_t2 = load_dicom_dict(os.path.join(image_dir, str(row.study_id), str(row.series_id)))\n        print (len(ax_t2))","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:49.340922Z","iopub.execute_input":"2024-09-25T15:24:49.341317Z","iopub.status.idle":"2024-09-25T15:24:49.494909Z","shell.execute_reply.started":"2024-09-25T15:24:49.341274Z","shell.execute_reply":"2024-09-25T15:24:49.493538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"study_cords","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:49.497733Z","iopub.execute_input":"2024-09-25T15:24:49.498296Z","iopub.status.idle":"2024-09-25T15:24:49.524398Z","shell.execute_reply.started":"2024-09-25T15:24:49.498237Z","shell.execute_reply":"2024-09-25T15:24:49.523222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def connect_meta(df, dirname=\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/\"):\n    output = list()\n    for study_id, series_id, inst_id in df.groupby([\"study_id\", \"series_id\", \"instance_number\"]).first().index:\n        #print (study_id, series_id, inst_id)\n        \n        metadata = dict()\n        metadata['study_id'] = int(study_id)\n        metadata['series_id'] = int(series_id)\n        metadata['instance_number'] = int(inst_id)\n        \n        d = pydicom.dcmread(os.path.join(dirname, str(study_id), str(series_id), str(inst_id) + \".dcm\"))\n        \n        metadata[\"ImagePositionPatient\"] = d.ImagePositionPatient\n        metadata[\"ImageOrientationPatient\"] = d.ImageOrientationPatient\n        metadata[\"PixelSpacing\"] = d.PixelSpacing\n        metadata[\"height\"], metadata[\"width\"] = d.pixel_array.shape\n        output.append(metadata)\n        \n    return df.merge(pd.DataFrame(output), how=\"inner\", on=['study_id', 'series_id', 'instance_number'])\n\nstudy_cords = connect_meta(study_cords)\n\nstudy_cords.head(2)","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:49.526164Z","iopub.execute_input":"2024-09-25T15:24:49.526556Z","iopub.status.idle":"2024-09-25T15:24:49.737772Z","shell.execute_reply.started":"2024-09-25T15:24:49.526513Z","shell.execute_reply":"2024-09-25T15:24:49.736303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def D_points_to_planes(axi_points, axi_planes, sag_planes):\n    assert len(axi_points) == len(axi_planes)\n    output_shape = [len(axi_points), len(sag_planes)]\n    D = np.zeros(output_shape)\n    cords = np.zeros(output_shape + [7])\n    \n    for i, (axi_point, (from_id, axi_plane)) in enumerate(zip(axi_points, axi_planes.iterrows())):\n        for j, (to_id, sag_plane) in enumerate(sag_planes.iterrows()):\n            d, point = project_ax2sag(sag_plane, axi_plane, axi_point[0], axi_point[1])\n            if np.all(point):\n                D[i, j] = d\n                cords[i, j] = i, j, point[0], point[1], d, from_id, to_id\n            else:\n                D[i, j] = np.inf\n    \n    return D, cords\n\ndef create_row(i, j, x, y, d, from_id, to_id, condition, level, planes):\n    row = planes[[\"study_id\", \"series_id\", \"instance_number\"]].iloc[j].to_dict()\n    row[\"condition\"] = condition\n    row[\"level\"] = level\n    row[\"x\"] = x\n    row[\"y\"] = y\n    row[\"distance\"] = d\n    row[\"from_id\"] = from_id\n    row[\"to_id\"] = to_id\n    return row\n\ndef project_points_to_planes(points_with_planes, other_planes, threshold):\n    D, points = D_points_to_planes(points_with_planes[[\"x\", \"y\"]].values, points_with_planes, other_planes)\n    output = list()\n    \n    #print (D)\n    for i, j, x, y, d, from_id, to_id in points[D <= threshold]:\n        condition, level = points_with_planes[[\"condition\", \"level\"]].iloc[int(i)] \n        output.append(create_row(int(i), int(j), x, y, d, from_id, to_id, condition, level, other_planes))\n    \n    return pd.DataFrame(output)\n\n# Axial T2 , Sagittal T1 , Sagittal T2/STIR\ndef project_combinations_planes(study_with_meta, meta_df, combinations, threshold=100):\n    output = list()\n    for from_plane, to_plane in combinations:\n        #print (from_plane, to_plane)\n        study_from = study_with_meta[study_with_meta.series_description == from_plane]\n        study_to = meta_df[meta_df.series_description == to_plane]\n        new_projections = project_points_to_planes(study_from, study_to, threshold)\n        if not new_projections.empty:\n            new_projections[\"series_description\"] = to_plane\n            output.append(new_projections)\n    if output:\n        return pd.concat(output, axis=0)\n    else:\n        return None\n\ncombinations = [\n    (\"Axial T2\", \"Sagittal T1\"),\n    (\"Axial T2\", \"Sagittal T2/STIR\"),\n    (\"Sagittal T1\", \"Axial T2\"),\n    (\"Sagittal T2/STIR\", \"Axial T2\")\n]\n\nanswers = project_combinations_planes(study_cords, study_cords, combinations, threshold=5)\n\nanswers","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:49.739326Z","iopub.execute_input":"2024-09-25T15:24:49.739797Z","iopub.status.idle":"2024-09-25T15:24:50.03679Z","shell.execute_reply.started":"2024-09-25T15:24:49.739753Z","shell.execute_reply":"2024-09-25T15:24:50.035547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\" | \".join(answers.iloc[0][[\"series_description\", \"condition\", \"level\"]].values)","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:24:50.038563Z","iopub.execute_input":"2024-09-25T15:24:50.039834Z","iopub.status.idle":"2024-09-25T15:24:50.050498Z","shell.execute_reply.started":"2024-09-25T15:24:50.03977Z","shell.execute_reply":"2024-09-25T15:24:50.049174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_dicom_from_plane(study_id, series_id, inst_id, dicom_ending=\".dcm\"):\n    root_dir = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images\"\n    d = pydicom.dcmread(os.path.join(root_dir, str(study_id), str(series_id), str(inst_id) + dicom_ending))\n    return d.pixel_array, str(d.InstanceNumber)\n    \ndef visualize_correspondence(corresp_df, origin_df, origin_meta, draw_rows=5, iter_stop=20):\n    count = 0\n    for i, (_, row) in enumerate(corresp_df.iterrows()):\n        if count == iter_stop:\n            break\n        if count % draw_rows == 0:\n            fig, ax = plt.subplots(draw_rows, 2)\n            fig.set_size_inches((20, 50))\n        from_id = int(row[\"from_id\"])\n        to_id = int(row[\"to_id\"])\n        \n        orig_dicom, inst_id = load_dicom_from_plane(*origin_df.loc[from_id][[\"study_id\", \"series_id\", \"instance_number\"]].values)\n        \n        #print (type(inst_id))\n        ax[i % draw_rows, 0].imshow(orig_dicom, cmap=plt.cm.gray)\n        x, y = origin_df.loc[from_id][[\"x\", \"y\"]].values\n        ax[i % draw_rows, 0].plot(x, y, 'x', markersize=15, markeredgewidth=3, color=\"green\")\n        ax[i % draw_rows, 0].set_title(\" | \".join(origin_df.loc[from_id][[\"series_description\", \"condition\", \"level\"]].to_list() + [inst_id]))\n\n        refer_dicom, inst_id = load_dicom_from_plane(*origin_meta.loc[to_id][[\"study_id\", \"series_id\", \"instance_number\"]].values)\n        ax[i % draw_rows, 1].imshow(refer_dicom, cmap=plt.cm.gray)\n        ax[i % draw_rows, 1].plot(row[\"x\"], row[\"y\"], 'x', markersize=15, markeredgewidth=3, color=\"green\")\n        ax[i % draw_rows, 1].set_title(\" | \".join(row[[\"series_description\", \"condition\", \"level\"]].to_list() + [inst_id]))\n        \n        count += 1\n\nvisualize_correspondence(answers, study_cords, study_cords)","metadata":{"execution":{"iopub.status.busy":"2024-09-25T15:31:11.589326Z","iopub.execute_input":"2024-09-25T15:31:11.589781Z","iopub.status.idle":"2024-09-25T15:31:32.837975Z","shell.execute_reply.started":"2024-09-25T15:31:11.589729Z","shell.execute_reply":"2024-09-25T15:31:32.836479Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}