{"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":[{"sourceType":"competition","sourceId":71549,"databundleVersionId":8561470}],"dockerImageVersionId":30786,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip -q install natsort","metadata":{"execution":{"iopub.status.busy":"2024-11-27T07:47:29.432681Z","iopub.execute_input":"2024-11-27T07:47:29.43308Z","iopub.status.idle":"2024-11-27T07:47:42.678091Z","shell.execute_reply.started":"2024-11-27T07:47:29.433046Z","shell.execute_reply":"2024-11-27T07:47:42.676273Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# code flow\n# https://www.kaggle.com/code/hengck23/2d-to-3d-projection-for-dicom/notebook\n# https://www.kaggle.com/code/hengck23/ver-1-demo-workflow-2-stage-approach/notebook?scriptVersionId=191553260\n# https://www.kaggle.com/code/hengck23/ver-2-more-magic-single-stage-model\n# https://www.kaggle.com/code/hengck23/place-multi-series-level-prediction-for-axial-t2\n# https://www.kaggle.com/code/hengck23/post-lhw-v24-ensemble-add-heng\n# https://www.kaggle.com/code/hengck23/ver-0-nonmax-supression-for-multi-point-detector","metadata":{"execution":{"iopub.status.busy":"2024-11-18T16:06:25.773444Z","iopub.execute_input":"2024-11-18T16:06:25.773868Z","iopub.status.idle":"2024-11-18T16:06:25.779011Z","shell.execute_reply.started":"2024-11-18T16:06:25.773817Z","shell.execute_reply":"2024-11-18T16:06:25.777759Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Import libraries\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport warnings\nimport os\nfrom tqdm import tqdm\nimport glob\nimport pydicom\nimport cv2\nfrom natsort import natsorted","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-11-27T07:47:42.68102Z","iopub.execute_input":"2024-11-27T07:47:42.681487Z","iopub.status.idle":"2024-11-27T07:47:44.949538Z","shell.execute_reply.started":"2024-11-27T07:47:42.681444Z","shell.execute_reply":"2024-11-27T07:47:44.948162Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Read training data\ntrain = pd.read_csv(\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-11-18T11:09:14.406148Z","iopub.execute_input":"2024-11-18T11:09:14.40667Z","iopub.status.idle":"2024-11-18T11:09:14.444523Z","shell.execute_reply.started":"2024-11-18T11:09:14.406628Z","shell.execute_reply":"2024-11-18T11:09:14.443503Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Total cases: \", len(train))","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:19:40.357643Z","iopub.execute_input":"2024-10-27T13:19:40.358064Z","iopub.status.idle":"2024-10-27T13:19:40.363913Z","shell.execute_reply.started":"2024-10-27T13:19:40.358021Z","shell.execute_reply":"2024-10-27T13:19:40.362653Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Explore columns\ntrain.columns","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:19:40.637783Z","iopub.execute_input":"2024-10-27T13:19:40.638273Z","iopub.status.idle":"2024-10-27T13:19:40.648633Z","shell.execute_reply.started":"2024-10-27T13:19:40.638202Z","shell.execute_reply":"2024-10-27T13:19:40.647259Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# View Sample\ntrain.head()","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:19:41.023878Z","iopub.execute_input":"2024-10-27T13:19:41.024348Z","iopub.status.idle":"2024-10-27T13:19:41.066851Z","shell.execute_reply.started":"2024-10-27T13:19:41.024304Z","shell.execute_reply":"2024-10-27T13:19:41.065349Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Checking Distrbution for Foraminal, Subatricular and Canal","metadata":{}},{"cell_type":"code","source":"# Create distribution plot\nfigure, axis = plt.subplots(nrows = 1, ncols = 3, figsize = (20, 5))\nfor idx, d in enumerate([\"foraminal\", \"subarticular\", \"canal\"]):\n    diagnosis = list(filter(lambda x: x.find(d) > -1, train.columns))\n    dff = train[diagnosis]\n    with warnings.catch_warnings():\n        warnings.simplefilter(action = 'ignore', category = FutureWarning)\n        value_counts = dff.apply(pd.value_counts).fillna(0).T\n    value_counts.plot(kind = \"bar\", stacked = True, ax = axis[idx])\n    axis[idx].set_title(f'{d} distribution')","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:19:41.894355Z","iopub.execute_input":"2024-10-27T13:19:41.894789Z","iopub.status.idle":"2024-10-27T13:19:43.505988Z","shell.execute_reply.started":"2024-10-27T13:19:41.894749Z","shell.execute_reply":"2024-10-27T13:19:43.504572Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Loading Images","metadata":{}},{"cell_type":"code","source":"# List out all studies\npart_1 = os.listdir('/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images')","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:19:43.507925Z","iopub.execute_input":"2024-10-27T13:19:43.508363Z","iopub.status.idle":"2024-10-27T13:19:43.689777Z","shell.execute_reply.started":"2024-10-27T13:19:43.508319Z","shell.execute_reply":"2024-10-27T13:19:43.68799Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Load metadata\ndf_meta_f = pd.read_csv(\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_series_descriptions.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:19:43.69145Z","iopub.execute_input":"2024-10-27T13:19:43.691891Z","iopub.status.idle":"2024-10-27T13:19:43.709772Z","shell.execute_reply.started":"2024-10-27T13:19:43.691851Z","shell.execute_reply":"2024-10-27T13:19:43.708549Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Create metadata for each scan\np1 = [(x, f\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/{x}\") for x in part_1]\nmeta_obj = { p[0]: { 'folder_path': p[1], \n                    'SeriesInstanceUIDs': [] \n                   } \n            for p in p1 }","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:19:47.471488Z","iopub.execute_input":"2024-10-27T13:19:47.472085Z","iopub.status.idle":"2024-10-27T13:19:47.483056Z","shell.execute_reply.started":"2024-10-27T13:19:47.472033Z","shell.execute_reply":"2024-10-27T13:19:47.481449Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Get SeriesInstanceUIDs in meta_obj\nfor m in meta_obj:\n    meta_obj[m]['SeriesInstanceUIDs'] = os.listdir(meta_obj[m]['folder_path'])","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:19:48.271331Z","iopub.execute_input":"2024-10-27T13:19:48.272548Z","iopub.status.idle":"2024-10-27T13:19:56.205327Z","shell.execute_reply.started":"2024-10-27T13:19:48.272493Z","shell.execute_reply":"2024-10-27T13:19:56.203746Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Grab the corresponding series instance description from train_series_descriptions\nfor study in tqdm(meta_obj):\n    for series_instance in meta_obj[study]['SeriesInstanceUIDs']:\n        if 'SeriesDescriptions' not in meta_obj[study]:\n            meta_obj[study]['SeriesDescriptions'] = []\n        meta_obj[study]['SeriesDescriptions'].append(df_meta_f[(df_meta_f['study_id'] == int(study)) & (df_meta_f['series_id'] == int(series_instance))]['series_description'].iloc[0])","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:19:56.207442Z","iopub.execute_input":"2024-10-27T13:19:56.207896Z","iopub.status.idle":"2024-10-27T13:20:00.329231Z","shell.execute_reply.started":"2024-10-27T13:19:56.207847Z","shell.execute_reply":"2024-10-27T13:20:00.327761Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"meta_obj['4003253']","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:20:00.331893Z","iopub.execute_input":"2024-10-27T13:20:00.332487Z","iopub.status.idle":"2024-10-27T13:20:00.342078Z","shell.execute_reply.started":"2024-10-27T13:20:00.332426Z","shell.execute_reply":"2024-10-27T13:20:00.340114Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Pull up images for one patient","metadata":{}},{"cell_type":"code","source":"patient = train.iloc[1]","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:20:00.343978Z","iopub.execute_input":"2024-10-27T13:20:00.34447Z","iopub.status.idle":"2024-10-27T13:20:00.352907Z","shell.execute_reply.started":"2024-10-27T13:20:00.344425Z","shell.execute_reply":"2024-10-27T13:20:00.351619Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"ptobj = meta_obj[str(patient['study_id'])]","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:20:00.731167Z","iopub.execute_input":"2024-10-27T13:20:00.73165Z","iopub.status.idle":"2024-10-27T13:20:00.737947Z","shell.execute_reply.started":"2024-10-27T13:20:00.731605Z","shell.execute_reply":"2024-10-27T13:20:00.736212Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(ptobj)","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:20:01.342927Z","iopub.execute_input":"2024-10-27T13:20:01.343413Z","iopub.status.idle":"2024-10-27T13:20:01.351874Z","shell.execute_reply.started":"2024-10-27T13:20:01.343367Z","shell.execute_reply":"2024-10-27T13:20:01.348789Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Get data into the format\n\"\"\"\nim_list_dcm = {\n    '{SeriesInstanceUID}': {\n        'images': [\n            {'SOPInstanceUID': ...,\n             'dicom': PyDicom object\n            },\n            ...,\n        ],\n        'description': # SeriesDescription\n    },\n    ...\n}\n\"\"\"\nim_list_dcm = {}\nfor idx, i in enumerate(ptobj['SeriesInstanceUIDs']):\n    im_list_dcm[i] = {'images': [], 'description': ptobj['SeriesDescriptions'][idx]}\n    images = glob.glob(f\"{ptobj['folder_path']}/{ptobj['SeriesInstanceUIDs'][idx]}/*.dcm\")\n    for j in sorted(images, key=lambda x: int(x.split('/')[-1].replace('.dcm', ''))):\n        im_list_dcm[i]['images'].append({\n            'SOPInstanceUID': j.split('/')[-1].replace('.dcm', ''), \n            'dicom': pydicom.dcmread(j) })","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:20:02.912271Z","iopub.execute_input":"2024-10-27T13:20:02.912787Z","iopub.status.idle":"2024-10-27T13:20:03.907652Z","shell.execute_reply.started":"2024-10-27T13:20:02.912741Z","shell.execute_reply":"2024-10-27T13:20:03.906313Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Function to display images\ndef display_images(images, title, max_images_per_row=4):\n    # Calculate the number of rows needed\n    num_images = len(images)\n    num_rows = (num_images + max_images_per_row - 1) // max_images_per_row  # Ceiling division\n\n    # Create a subplot grid\n    fig, axes = plt.subplots(num_rows, max_images_per_row, figsize=(5, 1.5 * num_rows))\n    \n    # Flatten axes array for easier looping if there are multiple rows\n    if num_rows > 1:\n        axes = axes.flatten()\n    else:\n        axes = [axes]  # Make it iterable for consistency\n\n    # Plot each image\n    for idx, image in enumerate(images):\n        ax = axes[idx]\n        ax.imshow(image, cmap='gray')  # Assuming grayscale for simplicity, change cmap as needed\n        ax.axis('off')  # Hide axes\n\n    # Turn off unused subplots\n    for idx in range(num_images, len(axes)):\n        axes[idx].axis('off')\n    fig.suptitle(title, fontsize=16)\n\n    plt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:22:33.289113Z","iopub.execute_input":"2024-10-27T13:22:33.290508Z","iopub.status.idle":"2024-10-27T13:22:33.300516Z","shell.execute_reply.started":"2024-10-27T13:22:33.290454Z","shell.execute_reply":"2024-10-27T13:22:33.298903Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"for i in im_list_dcm:\n    display_images([x['dicom'].pixel_array for x in im_list_dcm[i]['images']], \n                   im_list_dcm[i]['description'])","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:22:33.673675Z","iopub.execute_input":"2024-10-27T13:22:33.674167Z","iopub.status.idle":"2024-10-27T13:22:40.761416Z","shell.execute_reply.started":"2024-10-27T13:22:33.674123Z","shell.execute_reply":"2024-10-27T13:22:40.759667Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Plot co-ordinates of Pathologies for patients","metadata":{}},{"cell_type":"code","source":"df_coor = pd.read_csv(\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_label_coordinates.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:23:19.641102Z","iopub.execute_input":"2024-10-27T13:23:19.641725Z","iopub.status.idle":"2024-10-27T13:23:19.771377Z","shell.execute_reply.started":"2024-10-27T13:23:19.641665Z","shell.execute_reply":"2024-10-27T13:23:19.769736Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df_coor.head()","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:23:26.508798Z","iopub.execute_input":"2024-10-27T13:23:26.509297Z","iopub.status.idle":"2024-10-27T13:23:26.528613Z","shell.execute_reply.started":"2024-10-27T13:23:26.509222Z","shell.execute_reply":"2024-10-27T13:23:26.526965Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def display_coor_on_img(c, i, title):\n    center_coordinates = (int(c['x']), int(c['y']))\n    radius = 10\n    color = (255, 0, 0)  # Red color in BGR\n    thickness = 2\n    IMG = i['dicom'].pixel_array\n    IMG_normalized = cv2.normalize(IMG, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)\n    \n    IMG_with_circle = cv2.circle(IMG_normalized.copy(), center_coordinates, radius, color, thickness)\n    \n    # Convert the image from BGR to RGB for correct color display in matplotlib\n    IMG_with_circle = cv2.cvtColor(IMG_with_circle, cv2.COLOR_BGR2RGB)\n    \n    # Display the image\n    plt.imshow(IMG_with_circle)\n    plt.axis('off')  # Turn off axis numbers and ticks\n    plt.title(title)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:26:59.091835Z","iopub.execute_input":"2024-10-27T13:26:59.092803Z","iopub.status.idle":"2024-10-27T13:26:59.102477Z","shell.execute_reply.started":"2024-10-27T13:26:59.092715Z","shell.execute_reply":"2024-10-27T13:26:59.101143Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"coor_entries = df_coor[df_coor['study_id'] == int(patient['study_id'])]","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:27:11.01249Z","iopub.execute_input":"2024-10-27T13:27:11.013054Z","iopub.status.idle":"2024-10-27T13:27:11.021238Z","shell.execute_reply.started":"2024-10-27T13:27:11.012988Z","shell.execute_reply":"2024-10-27T13:27:11.019875Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Only showing severe cases for this patient\")\nfor idc, c in coor_entries.iterrows():\n    for i in im_list_dcm[str(c['series_id'])]['images']:\n        if int(i['SOPInstanceUID']) == int(c['instance_number']):\n            try:\n                patient_severity = patient[\n                    f\"{c['condition'].lower().replace(' ', '_')}_{c['level'].lower().replace('/', '_')}\"\n                ]\n            except Exception as e:\n                patient_severity = \"unknown severity\"\n            title = f\"{i['SOPInstanceUID']} \\n{c['level']}, {c['condition']}: {patient_severity} \\n{c['x']}, {c['y']}\"\n            if patient_severity == 'Severe':\n                display_coor_on_img(c, i, title)","metadata":{"execution":{"iopub.status.busy":"2024-10-27T13:27:18.893299Z","iopub.execute_input":"2024-10-27T13:27:18.893745Z","iopub.status.idle":"2024-10-27T13:27:19.366448Z","shell.execute_reply.started":"2024-10-27T13:27:18.8937Z","shell.execute_reply":"2024-10-27T13:27:19.365216Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Data Preprocessing","metadata":{}},{"cell_type":"code","source":"kaggle_dir = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification\"","metadata":{"execution":{"iopub.status.busy":"2024-11-27T06:15:50.594922Z","iopub.execute_input":"2024-11-27T06:15:50.596123Z","iopub.status.idle":"2024-11-27T06:15:50.601545Z","shell.execute_reply.started":"2024-11-27T06:15:50.596068Z","shell.execute_reply":"2024-11-27T06:15:50.600294Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#helper\nclass dotdict(dict):\n    __setattr__ = dict.__setitem__\n    __delattr__ = dict.__delitem__\n\n    def __getattr__(self, name):\n        try:\n            return self[name]\n        except KeyError:\n            raise AttributeError(name)","metadata":{"execution":{"iopub.status.busy":"2024-11-27T07:48:08.989339Z","iopub.execute_input":"2024-11-27T07:48:08.989742Z","iopub.status.idle":"2024-11-27T07:48:08.996164Z","shell.execute_reply.started":"2024-11-27T07:48:08.989707Z","shell.execute_reply":"2024-11-27T07:48:08.994799Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Dicom reader\n# Normalize image to 8-bit grayscale format\ndef normalise_to_8bit(x, lower = 0.1, upper = 99.9):\n    # Get 0.1 and 99.9 percentile\n    lower, upper = np.percentile(x, (lower, upper))\n    # Clip values below 0.1 and above 99.9 percentile\n    x = np.clip(x, lower, upper)\n    # Shift x so min becomes 0\n    x = x - np.min(x)\n    # Scale values between 0 and 1\n    x = x / np.max(x)\n    return (x * 255).astype(np.uint8)","metadata":{"execution":{"iopub.status.busy":"2024-11-27T06:15:53.133158Z","iopub.execute_input":"2024-11-27T06:15:53.133579Z","iopub.status.idle":"2024-11-27T06:15:53.140781Z","shell.execute_reply.started":"2024-11-27T06:15:53.133545Z","shell.execute_reply":"2024-11-27T06:15:53.139287Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# discontinous volume will be splitted to continous chunks\ndef load_and_split_mri_from_dicom_dir(\n        study_id,\n        series_id,\n        series_description,\n    ):\n    dicom_dir= f'{kaggle_dir}/train_images/{study_id}/{series_id}'\n\n    #---\n    dicom_file = natsorted(glob.glob( f'{dicom_dir}/*.dcm'))\n    instance_number = [int(f.split('/')[-1].split('.')[0]) for f in dicom_file]\n    dicom = [pydicom.dcmread(f) for f in dicom_file]\n\n    dicom_df = []\n    for i,d in zip(instance_number,dicom):\n        dicom_df.append(\n            dotdict(\n                study_id=study_id,\n                series_id=series_id,\n                series_description=series_description,\n                instance_number=i,\n                ImagePositionPatient=tuple([float(v) for v in d.ImagePositionPatient]),  # not continous\n                ImageOrientationPatient=tuple([float(v) for v in d.ImageOrientationPatient]),\n                PixelSpacing=tuple([float(v) for v in d.PixelSpacing]),\n                SpacingBetweenSlices=float(d.SpacingBetweenSlices),\n                SliceThickness=float(d.SliceThickness),\n            )\n        )\n    dicom_df = pd.DataFrame(dicom_df)\n    # All images may or may not have the same orientation\n    dicom_df = [d for _,d in dicom_df.groupby('ImageOrientationPatient')]\n\n    #sort slice\n    mri=[]\n    for df in dicom_df:\n        position = np.array(df['ImagePositionPatient'].values.tolist())\n        orientation = np.array(df['ImageOrientationPatient'].values.tolist())\n        normal = np.cross(orientation[:,:3], orientation[:,3:])\n        projection = np.sum(normal*position,1) #np.dot(normal, position)\n        df.loc[:,'projection'] = projection\n        df = df.sort_values('projection')\n\n        #todo: assert all slices are continous\n        assert len(df.SliceThickness.unique())==1\n        assert len(df.ImageOrientationPatient.unique())==1\n        assert len(df.ImageOrientationPatient.unique())==1\n        assert len(df.SpacingBetweenSlices.unique())==1\n\n        volume =[\n            dicom[instance_number.index(i)].pixel_array for i in df.instance_number\n        ]\n        volume = np.stack(volume)\n        volume = normalise_to_8bit(volume)\n        mri.append(dotdict(\n            df=df,\n            volume=volume,\n        ))\n    return mri","metadata":{"execution":{"iopub.status.busy":"2024-11-27T06:15:54.791267Z","iopub.execute_input":"2024-11-27T06:15:54.791707Z","iopub.status.idle":"2024-11-27T06:15:54.804533Z","shell.execute_reply.started":"2024-11-27T06:15:54.791671Z","shell.execute_reply":"2024-11-27T06:15:54.803243Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from mpl_toolkits.mplot3d import Axes3D\n\nimage_3d = load_and_split_mri_from_dicom_dir(\"4003253\", \"702807833\", \"Sagittal T2/STIR\")[0].volume\n\n# Visualizing a slice from the 3D image\nplt.imshow(image_3d[image_3d.shape[0] // 2], cmap='gray')  # Display a middle slice\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-27T06:15:55.580666Z","iopub.execute_input":"2024-11-27T06:15:55.581144Z","iopub.status.idle":"2024-11-27T06:15:56.578209Z","shell.execute_reply.started":"2024-11-27T06:15:55.581104Z","shell.execute_reply":"2024-11-27T06:15:56.576914Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Convert 2D (x,y) to 3D (x,y,z) for point in train_label_coordinates\ndef add_XYZ_to_label_df(study_id_df):\n    # Add Columns Width and Height of image\n    for col in ['W','H']:\n        study_id_df[col] = 0\n    # Add World coordinates xx, yy, zz\n    for col in ['xx','yy','zz']:\n        study_id_df[col] = 0.0\n\n    for t,d in study_id_df.iterrows():\n        dicom_file = f'/{kaggle_dir}/train_images/{d.study_id}/{d.series_id}/{d.instance_number}.dcm'\n        dicom = pydicom.dcmread(dicom_file)\n        H, W = dicom.pixel_array.shape\n        sx, sy, sz = [float(v) for v in dicom.ImagePositionPatient]\n        o0, o1, o2, o3, o4, o5 = [float(v) for v in dicom.ImageOrientationPatient]\n        delx, dely = dicom.PixelSpacing\n        xx =  o0*delx*d.x + o3*dely*d.y + sx\n        yy =  o1*delx*d.x + o4*dely*d.y + sy\n        zz =  o2*delx*d.x + o5*dely*d.y + sz\n        study_id_df.loc[t,'W'] = W\n        study_id_df.loc[t,'H'] = H\n        study_id_df.loc[t,'xx'] = xx\n        study_id_df.loc[t,'yy'] = yy\n        study_id_df.loc[t,'zz'] = zz\n    return study_id_df","metadata":{"execution":{"iopub.status.busy":"2024-11-27T06:15:56.58032Z","iopub.execute_input":"2024-11-27T06:15:56.580713Z","iopub.status.idle":"2024-11-27T06:15:56.592889Z","shell.execute_reply.started":"2024-11-27T06:15:56.580676Z","shell.execute_reply":"2024-11-27T06:15:56.591481Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Read all MRI for one patient\ndef load_for_one(study_id_df):\n    group = study_id_df.groupby(['series_id', 'series_description']).agg('first').index\n    mri = []\n    for series_id, series_description in group:\n        mri += load_and_split_mri_from_dicom_dir(\n            study_id = study_id,\n            series_id = series_id,\n            series_description = series_description)\n    return mri","metadata":{"execution":{"iopub.status.busy":"2024-11-27T06:15:57.193945Z","iopub.execute_input":"2024-11-27T06:15:57.194345Z","iopub.status.idle":"2024-11-27T06:15:57.201241Z","shell.execute_reply.started":"2024-11-27T06:15:57.194311Z","shell.execute_reply":"2024-11-27T06:15:57.199625Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Backproject 3D to 2D\ndef backproject_XYZ(xx, yy, zz, mri):\n    r = mri\n    d0 = r.df.iloc[0]\n    \n    sx, sy, sz = [float(v) for v in d0.ImagePositionPatient]\n    o0, o1, o2, o3, o4, o5 = [float(v) for v in d0.ImageOrientationPatient]\n    delx, dely = d0.PixelSpacing\n    delz = d0.SpacingBetweenSlices\n    \n    ax = np.array([o0, o1, o2])\n    ay = np.array([o3, o4, o5])\n    az = np.cross(ax, ay)\n    \n    p = np.array([xx-sx, yy-sy, zz-sz])\n    x = np.dot(ax, p)/delx\n    y = np.dot(ay, p)/dely\n    z = np.dot(az, p)/delz\n    x = round(x)\n    y = round(y)\n    z = round(z)\n    \n    D, H, W = r.volume.shape\n    inside = (x >= 0) & (x < W) & (y >= 0) & (y < H) & (z >= 0) & (z < D)\n    \n    if not inside:\n        return False, 0, 0, 0, 0\n    \n    n = r.df.instance_number.values[z]\n    return True, x, y, z, n","metadata":{"execution":{"iopub.status.busy":"2024-11-27T06:15:58.809335Z","iopub.execute_input":"2024-11-27T06:15:58.809745Z","iopub.status.idle":"2024-11-27T06:15:58.819861Z","shell.execute_reply.started":"2024-11-27T06:15:58.809711Z","shell.execute_reply":"2024-11-27T06:15:58.818409Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"desc_df = pd.read_csv(\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_series_descriptions.csv\")\nlabel_df = pd.read_csv(\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_label_coordinates.csv\")\nlabel_df = pd.merge(label_df, desc_df, on = ['study_id', 'series_id'], how = \"left\")","metadata":{"execution":{"iopub.status.busy":"2024-11-27T06:15:59.939286Z","iopub.execute_input":"2024-11-27T06:15:59.939747Z","iopub.status.idle":"2024-11-27T06:16:00.094099Z","shell.execute_reply.started":"2024-11-27T06:15:59.939709Z","shell.execute_reply":"2024-11-27T06:16:00.092898Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"label_df.loc[label_df['study_id'] == 4646740]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T06:32:29.420926Z","iopub.execute_input":"2024-11-27T06:32:29.421431Z","iopub.status.idle":"2024-11-27T06:32:29.44222Z","shell.execute_reply.started":"2024-11-27T06:32:29.42134Z","shell.execute_reply":"2024-11-27T06:32:29.440982Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"study_id = 4003253\nstudy_id_df = label_df[label_df.study_id==study_id]\nstudy_id_df = add_XYZ_to_label_df(study_id_df)\nprint(study_id_df.iloc[0])\n\nmri = load_for_one(study_id_df)\nprint('len(mri): ', len(mri))","metadata":{"execution":{"iopub.status.busy":"2024-11-27T06:16:01.435915Z","iopub.execute_input":"2024-11-27T06:16:01.436339Z","iopub.status.idle":"2024-11-27T06:16:03.187965Z","shell.execute_reply.started":"2024-11-27T06:16:01.4363Z","shell.execute_reply":"2024-11-27T06:16:03.186559Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#verify code\n\n'''\nstudy_id\tseries_id\tseries_description\n0\t4003253\t702807833\tSagittal T2/STIR\n1\t4003253\t1054713880\tSagittal T1\n2\t4003253\t2448190387\tAxial T2\n'''\n\n# backproject within the same view as truth point:\n# - given x,y,z of view v from label.csv\n# - project to 3d xx,yy,zz\n# - for the same view v, backproject from xx,yy,zz to vx,vy,vz\n# - check (vx,vy,vz) must be same as (x,y,z)\n\nprint('START VERIFICATION !!!')\nfor t, d in study_id_df.iterrows():\n    print('=================================')\n    print(d)\n    print('')\n    \n    xx, yy, zz = d.xx, d.yy, d.zz\n    found = 0\n    for r in mri:\n        d0 = r.df.iloc[0]\n        if not (d0.study_id == d.study_id) & (d0.series_id == d.series_id) & (d0.series_description == d.series_description):\n            continue\n        inside, x, y, z, n = backproject_XYZ(xx, yy, zz, r)\n        found += 1\n        print(\"Truth: \", d.instance_number, d.x, d.y)\n        print(\"Predict: \", n, x, y, f'inside={inside}', f'array index z={z}')\n    print(f'found={found}')\n    print('')    ","metadata":{"execution":{"iopub.status.busy":"2024-11-18T16:06:50.122999Z","iopub.execute_input":"2024-11-18T16:06:50.123944Z","iopub.status.idle":"2024-11-18T16:06:50.184463Z","shell.execute_reply.started":"2024-11-18T16:06:50.123876Z","shell.execute_reply":"2024-11-18T16:06:50.183219Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(mri[0].volume.shape)","metadata":{"execution":{"iopub.status.busy":"2024-11-18T16:07:05.578423Z","iopub.execute_input":"2024-11-18T16:07:05.578844Z","iopub.status.idle":"2024-11-18T16:07:05.584983Z","shell.execute_reply.started":"2024-11-18T16:07:05.578803Z","shell.execute_reply":"2024-11-18T16:07:05.583852Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Example of correct cross view 3d to 2d projection\nd = study_id_df.iloc[0] \nprint(d)\nprint('')\n\n\nfor r in mri:\n    d0 = r.df.iloc[0]\n    print(d0.series_id)\n    print(d0.series_description)\n    inside,x,y,z,n = backproject_XYZ(d.xx,d.yy,d.zz,r) \n    print(\"Truth: \", d.instance_number, d.x, d.y)\n    print('predict:', n, x, y, f'inside={inside}', f'array index z={z}')\n \n    \n    if inside:\n        slice = r.volume[z].copy()\n        slice[y] = 255 #slice[y] // 2 + 127\n        slice[:,x] = 255 #slice[:,x] // 2 + 127\n    else:\n        D,H,W = r.volume.shape\n        slice = np.zeros((H,W),dtype=np.uint8)\n    plt.imshow(slice, cmap='gray')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-18T16:19:16.5584Z","iopub.execute_input":"2024-11-18T16:19:16.55885Z","iopub.status.idle":"2024-11-18T16:19:18.299423Z","shell.execute_reply.started":"2024-11-18T16:19:16.558807Z","shell.execute_reply":"2024-11-18T16:19:18.29814Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#example of wrong cross view 3d to 2d projection ?????\nstudy_id = 4096820034\nstudy_id_df = label_df[label_df.study_id==study_id]\nstudy_id_df = add_XYZ_to_label_df(study_id_df)\n#print(study_id_df.iloc[0])\n\nmri = load_for_one(study_id_df)\nprint('len(mri): ', len(mri))\nprint('')\n\n\nfor r in mri:\n    d0 = r.df.iloc[0]\n    print(d0.series_id)\n    print(d0.series_description)\n    if not (d0.study_id == d.study_id) & (d0.series_id == d.series_id) & (d0.series_description == d.series_description):\n            continue\n    inside,x,y,z,n = backproject_XYZ(d.xx,d.yy,d.zz,r) \n    print(\"Truth: \", d.instance_number, d.x, d.y)\n    print('predict:', n, x, y, f'inside={inside}', f'array index z={z}')\n \n    \n    if inside:\n        slice = r.volume[z].copy()\n        slice[y] = 255 #slice[y] // 2 + 127\n        slice[:,x] = 255 #slice[:,x] // 2 + 127\n    else:\n        D,H,W = r.volume.shape\n        slice = np.zeros((H,W),dtype=np.uint8)\n    plt.imshow(slice, cmap='gray')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-18T16:18:52.843178Z","iopub.execute_input":"2024-11-18T16:18:52.844204Z","iopub.status.idle":"2024-11-18T16:18:55.791568Z","shell.execute_reply.started":"2024-11-18T16:18:52.844155Z","shell.execute_reply":"2024-11-18T16:18:55.790497Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#here is the magic !!!\n#plot axial slices\n\ncolor_table = [\n    [1,0,0],\n    [0,1,0],\n    [0,0,1],\n    [1,1,0],\n    [0,1,1],\n]\n\nstudy_id = 40745534 #88465004 #40745534 #4003253 #11340341 #4096820034\nstudy_id_df = label_df[label_df.study_id == study_id]\nstudy_id_df = add_XYZ_to_label_df(study_id_df)\nstudy_id_df = study_id_df.sort_values(['series_id','instance_number'])\nmri = load_for_one(study_id_df)\nprint('len(mri): ', len(mri))\nprint('')\n\n# For this study_id there are 3 series_id but 7 mri series\n# Remember MRI series are grouped according to ImageOrientationPatient\n# There are 5 axial slices and 2 saggital slices\n\nfig = plt.figure()\nax = fig.add_subplot(projection='3d')\nax.set_aspect('equal')\n\ncolor_i=0\nfor j,r in enumerate(mri):\n    D, H, W = r.volume.shape\n    if not 'Axial' in r.df.iloc[0].series_description: continue\n\n    color = color_table[color_i]\n    color_i += 1\n\n    N = len(r.df)\n    for i in range(N):\n        d0 = r.df.iloc[i]\n\n        o0, o1, o2, o3, o4, o5 = d0.ImageOrientationPatient\n        ox = np.array([o0, o1, o2])\n        oy = np.array([o3, o4, o5])\n        sx,sy,sz = d0.ImagePositionPatient\n        s = np.array([sx,sy,sz])\n        delx, dely = d0.PixelSpacing\n\n        p0 = s\n        p1 = s + W*delx*ox\n        p2 = s + H*dely*oy\n        p3 = s + H*dely*oy + W*delx*ox\n\n        grid = np.stack([p0, p1, p2, p3]).reshape(2,2,3)\n        gx = grid[:,:,0]\n        gy = grid[:, :, 1]\n        gz = grid[:, :, 2]\n\n        if i==0:\n            ax.plot_surface(gx, gy, gz, alpha=0.7, color=color)\n            ax.scatter([sx], [sy], [sz], color='black')\n\n        else:\n            ax.plot_surface(gx, gy, gz, alpha=0.3, color=color)\n            ax.scatter([sx], [sy], [sz], alpha=0.1, color='black')\n\n    #---\n\n    #plt.show()\n\nxLabel = ax.set_xlabel('x')\nyLabel = ax.set_ylabel('y')\nzLabel = ax.set_zlabel('z')\nplt.show()\nzz=0","metadata":{"execution":{"iopub.status.busy":"2024-11-18T16:32:10.299373Z","iopub.execute_input":"2024-11-18T16:32:10.29989Z","iopub.status.idle":"2024-11-18T16:32:13.577352Z","shell.execute_reply.started":"2024-11-18T16:32:10.299843Z","shell.execute_reply":"2024-11-18T16:32:13.576245Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#demo with real kaggle data and visualisation\n\n\nDATA_KAGGLE_DIR = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification'\n\ndef np_dot(a,b):\n    return np.sum(a * b, 1)\n\ndef normalise_to_8bit(x, lower=0.1, upper=99.9): \n    lower, upper = np.percentile(x, (lower, upper))\n    x = np.clip(x, lower, upper)\n    x = x - np.min(x)\n    x = x / np.max(x)\n    return (x * 255).astype(np.uint8)\n\ndef read_series(study_id,series_id,series_description):\n    error_code = ''\n    \n    data_kaggle_dir = DATA_KAGGLE_DIR\n    dicom_dir = f'{data_kaggle_dir}/train_images/{study_id}/{series_id}'\n\n    # read dicom file\n    dicom_file = natsorted(glob.glob(f'{dicom_dir}/*.dcm'))\n    instance_number = [int(f.split('/')[-1].split('.')[0]) for f in dicom_file]\n    dicom = [pydicom.dcmread(f) for f in dicom_file]\n\n    # make dicom header df\n    dicom_df = []\n    for i, d in zip(instance_number, dicom):  # d__.dict__\n        dicom_df.append(\n            dotdict(\n                study_id=study_id,\n                series_id=series_id,\n                series_description=series_description,\n                instance_number=i,\n                # InstanceNumber = d.InstanceNumber,\n                ImagePositionPatient=[float(v) for v in d.ImagePositionPatient],\n                ImageOrientationPatient=[float(v) for v in d.ImageOrientationPatient],\n                PixelSpacing=[float(v) for v in d.PixelSpacing],\n                SpacingBetweenSlices=float(d.SpacingBetweenSlices),\n                SliceThickness=float(d.SliceThickness),\n                grouping=str([round(float(v), 3) for v in d.ImageOrientationPatient]),\n                H=d.pixel_array.shape[0],\n                W=d.pixel_array.shape[1],\n            )\n        )\n    dicom_df = pd.DataFrame(dicom_df)\n    # dicom_df.to_csv('dicom_df.csv',index=False)\n    # exit(0)\n\n    #----\n    if ((dicom_df.W.nunique()!=1) or (dicom_df.H.nunique()!=1)):\n        error_code = '[multi-shape]'\n    Wmax = dicom_df.W.max()\n    Hmax = dicom_df.H.max()\n\n    # sort slices\n    dicom_df = [d for _, d in dicom_df.groupby('grouping')]\n\n    data = []\n    sort_data_by_group = []\n    for df in dicom_df:\n        position = np.array(df['ImagePositionPatient'].values.tolist())\n        orientation = np.array(df['ImageOrientationPatient'].values.tolist())\n        normal = np.cross(orientation[:, :3], orientation[:, 3:])\n        projection = np_dot(normal, position)\n        df.loc[:, 'projection'] = projection\n        df = df.sort_values('projection')\n\n\n        # todo: assert all slices are continous ??\n        # use  (position[-1]-position[0])/N = SpacingBetweenSlices ??\n        assert len(df.SliceThickness.unique()) == 1\n        #assert len(df.SpacingBetweenSlices.unique()) == 1\n\n\n        volume = []\n        for i in df.instance_number:\n            v = dicom[instance_number.index(i)].pixel_array\n            if error_code.find('multi-shape')!=-1:\n                H,W = v.shape\n                v=np.pad(v,[(0,Hmax-H),(0,Wmax-W)],'reflect')\n            volume.append(v)\n\n        volume = np.stack(volume)\n        volume = normalise_to_8bit(volume)\n\n        data.append(dotdict(\n            df=df,\n            volume=volume,\n        ))\n\n        if 'sagittal' in series_description.lower():\n            sort_data_by_group.append(position[0, 0])  # x\n        if 'axial' in series_description.lower():\n            sort_data_by_group.append(position[0, 2])  # z\n\n    data = [r for _, r in sorted(zip(sort_data_by_group, data))]\n    for i, r in enumerate(data):\n        r.df.loc[:, 'group'] = i\n\n    df = pd.concat([r.df for r in data])\n    df.loc[:, 'z'] = np.arange(len(df))\n    volume = np.concatenate([r.volume for r in data])\n    return volume, df, error_code\n\ndef do_resize_and_center(\n    image, reference_size\n):\n   \n    H, W = image.shape[:2]\n    if (W==reference_size) & (H==reference_size):\n        return image, (1,0,0)\n\n    s = reference_size / max(H, W)\n    m = cv2.resize(image, dsize=None, fx=s, fy=s)\n    h, w = m.shape[:2]\n    padx0 = (reference_size-w)//2\n    padx1 = reference_size-w-padx0\n    pady0 = (reference_size-h)//2\n    pady1 = reference_size-h-pady0\n\n    m = np.pad(m, [[pady0, pady1], [padx0, padx1], [0, 0]], mode='constant', constant_values=0)\n    #p = point * s +[[padx0,pady0]]\n    scale_param = s,padx0,pady0\n    return m, scale_param\n\n\n#read kaggle data----------------------------------------------------------\nstudy_id = 267842058\nseries_id = 894248358\nseries_description = 'sagittal_t1'\n\n#ground truth\ntruth_grade=[\n    'Normal/Mild',\n    'Moderate',\n    'Severe',\n    'Moderate',\n    'Severe',\n    'Normal/Mild',\n    'Normal/Mild',\n    'Moderate',\n    'Severe',\n    'Severe',\n]\n\n\n\nvolume, dicom_df, _ = read_series(study_id,series_id,series_description)\nprint('volume', volume.shape)\n\n# image = np.ascontiguousarray(volume.transpose(1,2,0))\n# image, scale_param = do_resize_and_center(\n#     image, reference_size=320\n# )\n# image = np.ascontiguousarray(image.transpose(2,0,1))\n# print('image', image.shape)\n\n# batch = {\n#     'D': [len(image)], #only one volume\n#     'image': torch.from_numpy(image).byte(), \n# }\n\n# net = Net(pretrained=False, cfg=None)#\n# state_dict = torch.load(\n#     '/kaggle/input/rnsa2024-single-stage-model/00034142.pth',\n#     map_location=lambda storage, loc: storage, weights_only=True)['state_dict']\n# print(net.load_state_dict(state_dict, strict=False))  # True\n\n# #net = net.cuda()\n# net = net.eval()\n# net.output_type = ['infer']\n\n# with torch.no_grad():\n#     with torch.cuda.amp.autocast(enabled=True):\n#         output = net(batch)\n\n# zxy_mask = output['zxy_mask'].data.cpu().numpy()\n# xy = output['xy'][0].data.cpu().numpy()\n# z  = output['z'][0].data.cpu().numpy()\n# grade = output['grade'][0].data.cpu().numpy()\n# #print(grade) \n# print('predict ok!')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T07:48:43.47889Z","iopub.execute_input":"2024-11-27T07:48:43.479379Z","iopub.status.idle":"2024-11-27T07:48:43.70326Z","shell.execute_reply.started":"2024-11-27T07:48:43.479339Z","shell.execute_reply":"2024-11-27T07:48:43.701793Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}