{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## I want to focus on:\n- getting the segmentations\n- resampling them to a common spacing (let's say 1mm pixdim)\n- analysing the positioning of each class (C*) w.r.t. the Z coordinate.\n\nThat's a try to look if this is an approximate signal for which C* vertebra classification\n\nI'm using some tools from my previous work (https://www.kaggle.com/code/kretes/segmentations-fracture-zoom-in)","metadata":{"execution":{"iopub.status.busy":"2022-07-29T07:11:03.457277Z"}}},{"cell_type":"code","source":"# source: https://www.kaggle.com/code/ipythonx/cervical-spine-fracture-detection-quick-eda\n\n!pip install -q ../input/for-pydicom/pylibjpeg-1.4.0-py3-none-any.whl\n!pip install -q ../input/for-pydicom/python_gdcm-3.0.14-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl\n!pip install -q ../input/for-pydicom/pylibjpeg_libjpeg-1.3.1-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:52:41.438155Z","iopub.execute_input":"2022-08-09T09:52:41.438723Z","iopub.status.idle":"2022-08-09T09:54:19.942749Z","shell.execute_reply.started":"2022-08-09T09:52:41.438608Z","shell.execute_reply":"2022-08-09T09:54:19.9415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd \nimport matplotlib.pyplot as plt \n\nfrom path import Path\nimport os \nimport glob\n\nimport os ","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:54:19.94567Z","iopub.execute_input":"2022-08-09T09:54:19.946073Z","iopub.status.idle":"2022-08-09T09:54:19.965529Z","shell.execute_reply.started":"2022-08-09T09:54:19.946035Z","shell.execute_reply":"2022-08-09T09:54:19.964333Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:54:30.678223Z","iopub.execute_input":"2022-08-09T09:54:30.678618Z","iopub.status.idle":"2022-08-09T09:54:30.696418Z","shell.execute_reply.started":"2022-08-09T09:54:30.678586Z","shell.execute_reply":"2022-08-09T09:54:30.695026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Segmentations","metadata":{}},{"cell_type":"code","source":"segmentations_paths = list(Path(\"../input/rsna-2022-cervical-spine-fracture-detection/segmentations\").glob(\"*.nii\"))","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:54:34.209268Z","iopub.execute_input":"2022-08-09T09:54:34.209723Z","iopub.status.idle":"2022-08-09T09:54:34.229486Z","shell.execute_reply.started":"2022-08-09T09:54:34.209683Z","shell.execute_reply":"2022-08-09T09:54:34.228307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patients_with_seg = [p.stem for p in segmentations_paths]\nlen(patients_with_seg)","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:54:36.069277Z","iopub.execute_input":"2022-08-09T09:54:36.069735Z","iopub.status.idle":"2022-08-09T09:54:36.081335Z","shell.execute_reply.started":"2022-08-09T09:54:36.069697Z","shell.execute_reply":"2022-08-09T09:54:36.079925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Used notebook https://www.kaggle.com/code/muki2003/display-dicom-and-nifti-format-s","metadata":{}},{"cell_type":"code","source":"import cv2\nimport PIL\nimport pydicom as dicom\nimport nibabel as nib\nimport matplotlib.patches as patches","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:54:19.969385Z","iopub.execute_input":"2022-08-09T09:54:19.970237Z","iopub.status.idle":"2022-08-09T09:54:20.636619Z","shell.execute_reply.started":"2022-08-09T09:54:19.970181Z","shell.execute_reply":"2022-08-09T09:54:20.635479Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from typing import Optional, Tuple\n\ndef load_raw_data(patient_id):\n    \"\"\"\n    returns a tuple of dicom_data for 1st slice, nifti_data(optional), no_dcm_files\n    \"\"\"\n    patient_dir = Path(f\"../input/rsna-2022-cervical-spine-fracture-detection/train_images/{patient_id}\")\n    files = list(patient_dir.glob(\"*.dcm\"))\n    no_files = len(files)\n    image_path = files[0]\n    ds = dicom.dcmread(image_path)\n    \n    seg = None    \n    if patient_id in patients_with_seg:\n\n        patient_seg_path = Path(f\"../input/rsna-2022-cervical-spine-fracture-detection/segmentations/{patient_id}.nii\")\n        seg = nib.load(patient_seg_path)\n        \n    return ds, seg, no_files\n    \n\ndef load_case_data(patient_id):\n    \"\"\"\n    returns a tuple of dicom_data, image_data, vertebrae segmentation(optional)\n    \"\"\"\n    ds, seg_nib, no_files = load_raw_data(patient_id)\n    patient_dir = Path(f\"../input/rsna-2022-cervical-spine-fracture-detection/train_images/{patient_id}\")\n    \n    image = np.zeros((*ds.pixel_array.shape, no_files))\n    \n    for i in range(no_files):\n        image_path = patient_dir / f\"{i+1}.dcm\"\n        \n        if image_path.exists():\n    \n            image[:,:,i] = dicom.dcmread(image_path).pixel_array\n    \n    seg = None    \n    if seg_nib is not None:\n        seg = seg_nib.get_fdata()\n        # https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection/data\n        seg = np.flip(seg, axis=-1) # flip in Z\n        seg = np.rot90(seg) # rotate in XY \n        \n        assert seg.shape == image.shape\n        \n    return ds, image, seg\n    ","metadata":{"execution":{"iopub.status.busy":"2022-08-09T10:33:10.018255Z","iopub.execute_input":"2022-08-09T10:33:10.018677Z","iopub.status.idle":"2022-08-09T10:33:10.030116Z","shell.execute_reply.started":"2022-08-09T10:33:10.018645Z","shell.execute_reply":"2022-08-09T10:33:10.028792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def display_patient_img(img, n=42, offset=0):\n    \n    plt.figure(figsize=(20, n // 3))\n\n    # image\n    img_per_row = 7\n    rows = (n // img_per_row) + 1\n    for i in range(1, n+1):\n        ax = plt.subplot(rows, img_per_row, i)        \n        plt.axis('off')\n        slice_index = i+offset-1\n        plt.imshow(img[:,:,slice_index])\n        \n        \n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:35:37.403579Z","iopub.execute_input":"2022-08-09T11:35:37.404081Z","iopub.status.idle":"2022-08-09T11:35:37.412086Z","shell.execute_reply.started":"2022-08-09T11:35:37.404041Z","shell.execute_reply":"2022-08-09T11:35:37.411035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tqdm import tqdm\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2022-08-09T10:11:42.573869Z","iopub.execute_input":"2022-08-09T10:11:42.574391Z","iopub.status.idle":"2022-08-09T10:11:43.472612Z","shell.execute_reply.started":"2022-08-09T10:11:42.574354Z","shell.execute_reply":"2022-08-09T10:11:43.471502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This takes < 5 minutes to gather\n\ndata = []\n\nfor patient_id in tqdm(train_df[\"StudyInstanceUID\"]):\n    ds, seg_nib, no_slices = load_raw_data(patient_id)\n    data.append({\n        \"patient_id\": patient_id,\n        \"spacing_0\": ds.PixelSpacing[0], \n        \"spacing_1\": ds.PixelSpacing[1], \n        \"slice_thickness\": float(ds.SliceThickness),\n        \"orientation\": ds.ImageOrientationPatient,\n        \"no_slices\": no_slices\n    })\ndf = pd.DataFrame(data)\ndf.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:19:52.125822Z","iopub.execute_input":"2022-08-09T11:19:52.126269Z","iopub.status.idle":"2022-08-09T11:20:04.539848Z","shell.execute_reply.started":"2022-08-09T11:19:52.126234Z","shell.execute_reply":"2022-08-09T11:20:04.538817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# since I can't yet reorient the images I simply skip the ones that are non-standard:\ndf_all = df.copy()\ndf = df[df[\"orientation\"].apply(lambda x: x == [1,0,0,0,1,0])]\ndf.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:29:45.265586Z","iopub.execute_input":"2022-08-09T11:29:45.266022Z","iopub.status.idle":"2022-08-09T11:29:45.295735Z","shell.execute_reply.started":"2022-08-09T11:29:45.265988Z","shell.execute_reply":"2022-08-09T11:29:45.294396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df[(df[\"spacing_0\"] != df[\"spacing_1\"])]","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:29:56.249172Z","iopub.execute_input":"2022-08-09T11:29:56.250003Z","iopub.status.idle":"2022-08-09T11:29:56.262399Z","shell.execute_reply.started":"2022-08-09T11:29:56.249933Z","shell.execute_reply":"2022-08-09T11:29:56.261278Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Among those - spacing in X is equal to spacing in Y","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(12,12))\nsns.histplot(df[\"spacing_0\"], )\nprint(df[\"spacing_0\"].median())\nplt.figure(figsize=(12,12))\nsns.histplot(df[\"slice_thickness\"])\nprint(df[\"slice_thickness\"].median())","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:30:16.716899Z","iopub.execute_input":"2022-08-09T11:30:16.717399Z","iopub.status.idle":"2022-08-09T11:30:17.311601Z","shell.execute_reply.started":"2022-08-09T11:30:16.71736Z","shell.execute_reply":"2022-08-09T11:30:17.310054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# spacing in-slice is between 0.2 and 0.7 with median 0.312 mm\n\n# spacing between slices is in range 0.5 - 1.0 mm with median 0.625","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(16,16))\nsns.scatterplot(df[\"slice_thickness\"], df[\"spacing_0\"], size= df[\"no_slices\"])","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:30:21.72099Z","iopub.execute_input":"2022-08-09T11:30:21.722134Z","iopub.status.idle":"2022-08-09T11:30:22.191194Z","shell.execute_reply.started":"2022-08-09T11:30:21.72208Z","shell.execute_reply":"2022-08-09T11:30:22.189685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see above that the larger the tickness between slices - the less slices there are. \n\nIt seems that the scans are covering roughly the same physical extent - let's check that","metadata":{}},{"cell_type":"code","source":"df['physical_extent_z'] = df['slice_thickness'] * df[\"no_slices\"]","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:30:29.637058Z","iopub.execute_input":"2022-08-09T11:30:29.637483Z","iopub.status.idle":"2022-08-09T11:30:29.646409Z","shell.execute_reply.started":"2022-08-09T11:30:29.637449Z","shell.execute_reply":"2022-08-09T11:30:29.645021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(14,14))\nsns.histplot(df['physical_extent_z'])\ndf['physical_extent_z'].median()","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:30:30.02502Z","iopub.execute_input":"2022-08-09T11:30:30.025722Z","iopub.status.idle":"2022-08-09T11:30:30.347241Z","shell.execute_reply.started":"2022-08-09T11:30:30.025673Z","shell.execute_reply":"2022-08-09T11:30:30.345786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It seems that the physical extent is quite a tight distribution although with a large range","metadata":{"execution":{"iopub.status.busy":"2022-07-31T18:06:24.958425Z","iopub.execute_input":"2022-07-31T18:06:24.959072Z","iopub.status.idle":"2022-07-31T18:06:24.964342Z","shell.execute_reply.started":"2022-07-31T18:06:24.959033Z","shell.execute_reply":"2022-07-31T18:06:24.9635Z"}}},{"cell_type":"code","source":"df[\"physical_extent_z\"].min(), df[\"physical_extent_z\"].max()","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:30:37.574162Z","iopub.execute_input":"2022-08-09T11:30:37.574642Z","iopub.status.idle":"2022-08-09T11:30:37.58418Z","shell.execute_reply.started":"2022-08-09T11:30:37.574604Z","shell.execute_reply":"2022-08-09T11:30:37.583079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.sort_values(\"physical_extent_z\")[:3]","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:30:41.277992Z","iopub.execute_input":"2022-08-09T11:30:41.278407Z","iopub.status.idle":"2022-08-09T11:30:41.300058Z","shell.execute_reply.started":"2022-08-09T11:30:41.278375Z","shell.execute_reply":"2022-08-09T11:30:41.298865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1.2.826.0.1.3680043.20574","metadata":{}},{"cell_type":"code","source":"for patient_id in df.sort_values(\"physical_extent_z\")[:2][\"patient_id\"].values:\n    _, img, __ = load_case_data(patient_id)\n    print(patient_id, img.shape[-1])\n    display_patient_img(img, n = img.shape[-1])","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:35:44.414277Z","iopub.execute_input":"2022-08-09T11:35:44.41471Z","iopub.status.idle":"2022-08-09T11:36:09.421126Z","shell.execute_reply.started":"2022-08-09T11:35:44.414674Z","shell.execute_reply":"2022-08-09T11:36:09.41925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#  \t1.2.826.0.1.3680043.20574\nLooks like a thoracic scan without any cervical part","metadata":{}},{"cell_type":"code","source":"df.merge(train_df, left_on=\"patient_id\", right_on=\"StudyInstanceUID\").sort_values(\"physical_extent_z\")[:3]","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:33:58.462415Z","iopub.execute_input":"2022-08-09T11:33:58.462842Z","iopub.status.idle":"2022-08-09T11:33:58.496351Z","shell.execute_reply.started":"2022-08-09T11:33:58.462809Z","shell.execute_reply":"2022-08-09T11:33:58.495226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- 10 scans have non-standard orientation and they need to be re-oriented into same coordinate system before using\n- Looking at the non-standard orientations I can see that my way of reading data isn't preserving the right coordinate system","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:25:47.220698Z","iopub.execute_input":"2022-08-09T11:25:47.221195Z","iopub.status.idle":"2022-08-09T11:25:47.229201Z","shell.execute_reply.started":"2022-08-09T11:25:47.221159Z","shell.execute_reply":"2022-08-09T11:25:47.227802Z"}}},{"cell_type":"code","source":"df_all[df_all[\"orientation\"].apply(lambda x: x != [1,0,0,0,1,0])]","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:32:22.531327Z","iopub.execute_input":"2022-08-09T11:32:22.531748Z","iopub.status.idle":"2022-08-09T11:32:22.565933Z","shell.execute_reply.started":"2022-08-09T11:32:22.531714Z","shell.execute_reply":"2022-08-09T11:32:22.564694Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# TODO: understand how to resample into the same orientation","metadata":{},"execution_count":null,"outputs":[]}]}