{"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":"# Creates 9 Sagittal midline JPG images for each case in the training dataset.\n\n- Windows were too dark in prior versions\n- WindowCenter and WindowLevel used from metadata rather than fixed values\n- Fastai DataFrame.from_dicoms used to create metadata csv","metadata":{}},{"cell_type":"code","source":"#Some of the data is compressed and needs GDCM\n! pip install pylibjpeg -q\n! pip install python-gdcm -q\n! pip install pylibjpeg-libjpeg -q","metadata":{"execution":{"iopub.status.busy":"2022-08-23T17:10:28.465923Z","iopub.execute_input":"2022-08-23T17:10:28.466651Z","iopub.status.idle":"2022-08-23T17:11:07.235755Z","shell.execute_reply.started":"2022-08-23T17:10:28.466589Z","shell.execute_reply":"2022-08-23T17:11:07.234272Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from fastai.vision.all import *\nfrom fastai.medical.imaging import *\nfrom fastcore.all import *\nimport pandas as pd\nimport pydicom\nimport numpy as np\n\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\n%matplotlib inline","metadata":{"execution":{"iopub.status.busy":"2022-08-23T17:11:07.23996Z","iopub.execute_input":"2022-08-23T17:11:07.240366Z","iopub.status.idle":"2022-08-23T17:11:11.135531Z","shell.execute_reply.started":"2022-08-23T17:11:07.24033Z","shell.execute_reply":"2022-08-23T17:11:11.134269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_DIR = Path(\"../input/rsna-2022-cervical-spine-fracture-detection/\")\ncases = (DATA_DIR/\"train_images\").glob('*')","metadata":{"execution":{"iopub.status.busy":"2022-08-23T17:11:11.13697Z","iopub.execute_input":"2022-08-23T17:11:11.13752Z","iopub.status.idle":"2022-08-23T17:11:11.142863Z","shell.execute_reply.started":"2022-08-23T17:11:11.137486Z","shell.execute_reply":"2022-08-23T17:11:11.141934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Get Metadata for first image in each study using fastai","metadata":{}},{"cell_type":"code","source":"first_images = []\nfor case in cases:\n    fns = list(case.glob('*'))\n    first_images.append(fns[0])\n\nlen(first_images)","metadata":{"execution":{"iopub.status.busy":"2022-08-23T17:49:36.456319Z","iopub.execute_input":"2022-08-23T17:49:36.456845Z","iopub.status.idle":"2022-08-23T17:51:11.001028Z","shell.execute_reply.started":"2022-08-23T17:49:36.456805Z","shell.execute_reply":"2022-08-23T17:51:11.000024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Fastai dataframe extension for dicom metadata\nmetadata = pd.DataFrame.from_dicoms(first_images)\nmetadata.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-23T17:51:11.002964Z","iopub.execute_input":"2022-08-23T17:51:11.003598Z","iopub.status.idle":"2022-08-23T17:51:50.285761Z","shell.execute_reply.started":"2022-08-23T17:51:11.00356Z","shell.execute_reply":"2022-08-23T17:51:50.284366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"metadata.to_csv('/kaggle/working/df.csv',index=False)","metadata":{"execution":{"iopub.status.busy":"2022-08-23T17:51:50.287346Z","iopub.execute_input":"2022-08-23T17:51:50.287853Z","iopub.status.idle":"2022-08-23T17:51:50.340746Z","shell.execute_reply.started":"2022-08-23T17:51:50.287808Z","shell.execute_reply":"2022-08-23T17:51:50.339063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Columns in metadata","metadata":{}},{"cell_type":"code","source":"metadata.columns","metadata":{"execution":{"iopub.status.busy":"2022-08-23T18:08:34.580347Z","iopub.execute_input":"2022-08-23T18:08:34.581372Z","iopub.status.idle":"2022-08-23T18:08:34.588543Z","shell.execute_reply.started":"2022-08-23T18:08:34.581332Z","shell.execute_reply":"2022-08-23T18:08:34.587325Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Distribution of Window Center and Width","metadata":{}},{"cell_type":"code","source":"wc = metadata.WindowCenter.values\nww = metadata.WindowWidth.values\n\nplt.hist(wc)\nplt.hist(ww)","metadata":{"execution":{"iopub.status.busy":"2022-08-23T18:23:37.073017Z","iopub.execute_input":"2022-08-23T18:23:37.073426Z","iopub.status.idle":"2022-08-23T18:23:37.312877Z","shell.execute_reply.started":"2022-08-23T18:23:37.073372Z","shell.execute_reply":"2022-08-23T18:23:37.311742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### The mode of the Window Center and Width can be used for the default values","metadata":{}},{"cell_type":"code","source":"metadata.WindowCenter.mode(), metadata.WindowWidth.mode()","metadata":{"execution":{"iopub.status.busy":"2022-08-23T18:31:03.872496Z","iopub.execute_input":"2022-08-23T18:31:03.873448Z","iopub.status.idle":"2022-08-23T18:31:03.88421Z","shell.execute_reply.started":"2022-08-23T18:31:03.873403Z","shell.execute_reply":"2022-08-23T18:31:03.883286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Adapted from Pydicom 'Load CT slices and plot axial, sagittal and coronal images'","metadata":{}},{"cell_type":"code","source":"#Adapted from Pydicom: 'Load CT slices and plot axial, sagittal and coronal images'\ndef create_3d(case):\n    files = []\n    fns = get_dicom_files(case)\n    for fn in fns:\n        files.append(pydicom.dcmread(fn))\n        \n    # skip files with no InstanceNumber (eg. Scout)\n    slices = []\n    skipcount = 0\n    for f in files:\n        if hasattr(f, 'InstanceNumber'):\n            slices.append(f)\n        else:\n            skipcount = skipcount + 1\n\n    if skipcount > 0:\n        print(\"skipped, no InstanceNumber: {}\".format(skipcount))\n\n    # ensure they are in the correct order\n    slices = sorted(slices, key=lambda s: s.ImagePositionPatient[2])\n    dcm = slices[0]\n\n    # pixel aspects, assuming all slices are the same\n    ps = dcm.PixelSpacing\n    ss = dcm.SliceThickness\n    #ax_aspect = ps[1]/ps[0]\n    sag_aspect = ps[1]/ss\n    #cor_aspect = ss/ps[0]\n\n    # create 3D array\n    img_shape = list(dcm.pixel_array.shape)\n    img_shape.append(len(slices))\n    img3d = np.zeros(img_shape)\n    \n    # fill 3D array with the images from the files\n    for i, s in enumerate(slices):\n        img2d = dcm_apply_windows(s)\n        img3d[:, :, i] = img2d\n        \n    return img3d, sag_aspect, fn.parent.name\n\ndef save_sag(img3d, sag_aspect, folder_name, interval=5):\n    img_shape = img3d.shape\n    for i in range(-4 * interval, 5 * interval, interval):\n        arr = img3d[:, img_shape[1]//2 + i, :] * sag_aspect\n        #flip and rotate so that C1 is at the top\n        arr = np.fliplr(arr)\n        arr = np.transpose(arr)\n        #array is a tensor from 0 - 1.0, multiply by 256 to get 8 bit image\n        arr *= 256\n        if not os.path.isdir(f\"/kaggle/working/sag/{folder_name}\"):\n             os.makedirs(f\"/kaggle/working/sag/{folder_name}\", exist_ok=True)\n        im = Image.fromarray(arr)\n        #to avoid error from F mode\n        if im.mode != 'RGB':\n            im = im.convert('RGB')\n        save_fn = f'/kaggle/working/sag/{folder_name}/{folder_name}_{i}.jpg'\n        im.save(save_fn)\n    ","metadata":{"execution":{"iopub.status.busy":"2022-08-23T18:33:04.766097Z","iopub.execute_input":"2022-08-23T18:33:04.766589Z","iopub.status.idle":"2022-08-23T18:33:04.783845Z","shell.execute_reply.started":"2022-08-23T18:33:04.766552Z","shell.execute_reply":"2022-08-23T18:33:04.782456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_window_from_dicom(dcm, default = (2000, 500)):\n    \"\"\"\n    Returns window width and window center values or first example if MultiValue\n    Strips comma from value if present (seen in a different dataset)\n    If no window width/level is provided or available, returns default.\n    \"\"\"\n    width, level = default\n\n    if \"WindowWidth\" in dcm:\n        width = dcm.WindowWidth\n        if isinstance(width, pydicom.multival.MultiValue):\n            width = float(width[0])\n        else:\n            width = float(str(width).replace(',', ''))\n\n    if \"WindowCenter\" in dcm:\n        level = dcm.WindowCenter\n        if isinstance(level, pydicom.multival.MultiValue):\n            level = float(level[0])\n        else:\n            level = float(str(level).replace(',', ''))\n            \n    return width, level\n\n\ndef dcm_apply_windows(dcm):\n    \"\"\"\n    Applies Intercept/Slope and Window Center/Width\n    \"\"\"\n    arr = dcm.pixel_array\n    #slope, intercept\n    slope = 1\n    intercept = 0\n    if \"RescaleIntercept\" in dcm and \"RescaleSlope\" in dcm:\n        intercept = int(dcm.RescaleIntercept)\n        slope = int(dcm.RescaleSlope)\n        \n    arr = arr * slope + intercept\n    \n    #window\n    width,level = get_window_from_dicom(dcm)\n    if width is not None and level is not None:\n        arr = np.clip(arr, level - width // 2, level + width // 2)\n        \n     #scale\n    arr = (arr - np.min(arr)) / np.max(arr)\n        \n    return arr","metadata":{"execution":{"iopub.status.busy":"2022-08-23T18:35:30.596288Z","iopub.execute_input":"2022-08-23T18:35:30.596803Z","iopub.status.idle":"2022-08-23T18:35:30.610193Z","shell.execute_reply.started":"2022-08-23T18:35:30.596767Z","shell.execute_reply":"2022-08-23T18:35:30.60881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def show_case(folder_name, interval=5):\n    plt.figure(figsize=(16, 12))\n    count = 1\n    for i in range(-4 * interval, 5 * interval, interval):\n        fn = f'/kaggle/working/sag/{folder_name}/{folder_name}_{i}.jpg'\n        img = mpimg.imread(fn)\n        plt.subplot(3,3,count)\n        count += 1\n        plt.imshow(img, cmap='bone')\n        plt.axis(\"off\")\n    plt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"case = DATA_DIR/'train_images/1.2.826.0.1.3680043.2487'\n\nimg3d, sag_aspect, folder_name = create_3d(case)\nsave_sag( img3d, sag_aspect, folder_name)\nshow_case('1.2.826.0.1.3680043.2487')","metadata":{"execution":{"iopub.status.busy":"2022-08-23T18:35:32.185842Z","iopub.execute_input":"2022-08-23T18:35:32.186536Z","iopub.status.idle":"2022-08-23T18:35:42.953983Z","shell.execute_reply.started":"2022-08-23T18:35:32.186496Z","shell.execute_reply":"2022-08-23T18:35:42.952898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def to_sag(case):\n    img3d, sag_aspect, folder_name = create_3d(case)\n    save_sag( img3d, sag_aspect, folder_name)\n\n    \ncases = (DATA_DIR/\"train_images\").glob('*')\nparallel(to_sag, list(cases))","metadata":{"execution":{"iopub.status.busy":"2022-08-12T17:57:27.335904Z","iopub.execute_input":"2022-08-12T17:57:27.336468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 16bit case with <300 slices and severe lordosis (curvature)\nshow_case('1.2.826.0.1.3680043.17625')","metadata":{"execution":{"iopub.status.busy":"2022-08-16T16:32:48.85669Z","iopub.execute_input":"2022-08-16T16:32:48.857105Z","iopub.status.idle":"2022-08-16T16:32:49.617389Z","shell.execute_reply.started":"2022-08-16T16:32:48.857072Z","shell.execute_reply":"2022-08-16T16:32:49.616165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}