{"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":"# Notebook 1 - Create sagittal slices for segmentations and original images","metadata":{}},{"cell_type":"markdown","source":"This is notebook one out of two. It creates 16 jpgs and masks for each of the segmented cases. The second notebook uses the sagittal images and masks to train a model and predict on nonsegmented cases.\n\n- The CSpine helper notebook is used to import libraries if the internet is off\n- The dicom metadata Intercept/Slope and Window Center/Width is used to convert the raw dicom pixel array\n- The images and masks are rotated to both be oriented sagitally","metadata":{}},{"cell_type":"code","source":"import sys\n\n!{sys.executable} -m pip install '../input/cspine-helper/fastai-2.7.9-py3-none-any.whl' --upgrade --no-deps --ignore-installed -q \n!{sys.executable} -m pip install '../input/cspine-helper/fastcore-1.5.21-py3-none-any.whl' -Uqq\n!{sys.executable} -m pip install '../input/cspine-helper/pylibjpeg-1.4.0-py3-none-any.whl' -q\n!{sys.executable} -m pip install '../input/cspine-helper/pylibjpeg_libjpeg-1.3.1-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl' -q\n!{sys.executable} -m pip install '../input/cspine-helper/python_gdcm-3.0.15-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl' -q\n#!{sys.executable} -m pip install '../input/cspine-helper/timm-0.6.7-py3-none-any.whl' -q","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:47:16.500996Z","iopub.execute_input":"2022-09-01T15:47:16.501763Z","iopub.status.idle":"2022-09-01T15:48:11.434628Z","shell.execute_reply.started":"2022-09-01T15:47:16.501635Z","shell.execute_reply":"2022-09-01T15:48:11.433203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ! pip install pylibjpeg -q\n# ! pip install python-gdcm -q\n# ! pip install pylibjpeg-libjpeg -q","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:11.438217Z","iopub.execute_input":"2022-09-01T15:48:11.438698Z","iopub.status.idle":"2022-09-01T15:48:11.443809Z","shell.execute_reply.started":"2022-09-01T15:48:11.43865Z","shell.execute_reply":"2022-09-01T15:48:11.442585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import fastai\nfastai.__version__","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:11.445599Z","iopub.execute_input":"2022-09-01T15:48:11.446032Z","iopub.status.idle":"2022-09-01T15:48:11.467035Z","shell.execute_reply.started":"2022-09-01T15:48:11.445987Z","shell.execute_reply":"2022-09-01T15:48:11.465567Z"},"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.image as mpimg\nimport cv2\nimport nibabel as nib","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:11.470203Z","iopub.execute_input":"2022-09-01T15:48:11.471369Z","iopub.status.idle":"2022-09-01T15:48:15.410709Z","shell.execute_reply.started":"2022-09-01T15:48:11.471323Z","shell.execute_reply":"2022-09-01T15:48:15.409555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train.csv\")\ntest_df = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/test.csv\")\ntrain_df.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:15.412671Z","iopub.execute_input":"2022-09-01T15:48:15.413593Z","iopub.status.idle":"2022-09-01T15:48:15.452122Z","shell.execute_reply.started":"2022-09-01T15:48:15.413548Z","shell.execute_reply":"2022-09-01T15:48:15.45099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Get a list of the study_ids for the cases with segmentation data","metadata":{}},{"cell_type":"code","source":"seg_path = Path('../input/rsna-2022-cervical-spine-fracture-detection/segmentations')\norig_path = Path('../input/rsna-2022-cervical-spine-fracture-detection/train_images')\n\nsegmentations = list(seg_path.iterdir())\nstudy_ids = [o.stem for o in segmentations]\n\nsegmentations[:3], study_ids[0]","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:15.454111Z","iopub.execute_input":"2022-09-01T15:48:15.454963Z","iopub.status.idle":"2022-09-01T15:48:15.494443Z","shell.execute_reply.started":"2022-09-01T15:48:15.454918Z","shell.execute_reply":"2022-09-01T15:48:15.493145Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Functions to create a 3D array from the dicom image directory\n\n- Reads Dicom and applies meta data values Intercept/Slope and Window Center/Width\n- Creates a 3D array ","metadata":{}},{"cell_type":"code","source":"## Original\n\n#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        \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, ps, ss\n\ndef 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-09-01T15:48:15.496546Z","iopub.execute_input":"2022-09-01T15:48:15.497343Z","iopub.status.idle":"2022-09-01T15:48:15.51463Z","shell.execute_reply.started":"2022-09-01T15:48:15.4973Z","shell.execute_reply":"2022-09-01T15:48:15.513413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# index = 0\n# seg = segmentations[index]\n# study_id = seg.stem\n# dirname = orig_path/study_id\n# study_id","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:15.516212Z","iopub.execute_input":"2022-09-01T15:48:15.517047Z","iopub.status.idle":"2022-09-01T15:48:15.541466Z","shell.execute_reply.started":"2022-09-01T15:48:15.516993Z","shell.execute_reply":"2022-09-01T15:48:15.539994Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## ","metadata":{}},{"cell_type":"markdown","source":"## Create matching sagittal views of the segmentation masks and original images at intervals near center of the volume","metadata":{}},{"cell_type":"code","source":"def get_sag(study_id, n = 16, spread = 160):\n    thickness = int(spread/n)\n    orig = []\n    segs = []\n    #orig\n    dirname = orig_path/study_id\n\n    img, _,_ = create_3d(dirname)\n    slice = int(img.shape[1]/2 - spread/2)\n    for i in range(n):\n        arr = img[:, slice + (i*thickness),:] \n        #flip and rotate so that C1 is at the top\n        arr = np.transpose(arr)\n        arr = np.flip(arr, 0)\n        orig.append(arr)\n    \n    #seg\n    seg = nib.load(seg_path/(study_id + '.nii')).get_fdata()\n\n    slice = int(img.shape[0]/2 - spread/2)\n    for i in range(n):\n        #use first dimension to get sagittal views\n        arr = seg[slice + (i*thickness),:,:]\n        #flip and rotate so that C1 is at the top\n        arr = np.flip(arr, 0)\n        arr = np.rot90(arr)\n        segs.append(arr)\n        \n    return orig, segs\n\ndef get_sag_for_predict(study_id, n = 16, spread = 160):\n    thickness = int(spread/n)\n    orig = []\n    dirname = orig_path/study_id\n\n    img,_,_ = create_3d(dirname)\n    slice = int(img.shape[1]/2 - spread/2)\n    for i in range(n):\n        arr = img[:, slice + (i*thickness),:] \n        #flip and rotate so that C1 is at the top\n        arr = np.transpose(arr)\n        arr = np.flip(arr, 0)\n        orig.append(arr)    \n    return orig","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:15.548453Z","iopub.execute_input":"2022-09-01T15:48:15.549275Z","iopub.status.idle":"2022-09-01T15:48:15.575112Z","shell.execute_reply.started":"2022-09-01T15:48:15.549225Z","shell.execute_reply":"2022-09-01T15:48:15.573206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Save the images","metadata":{}},{"cell_type":"code","source":"save_dir = Path('/kaggle/working/')\n\ndef save_sag(study_id):\n    orig, segs = get_sag(study_id)\n\n    orig = np.array(orig)  * 256\n    segs =  np.array(segs)\n\n    slices = segs.shape[0]\n    #print(slices, segs.shape)\n\n    for i in range(slices):\n        np.save(    f'{save_dir/study_id}_{i}.npy', segs[i,:,:])\n        cv2.imwrite(f'{save_dir/study_id}_{i}.jpg', orig[i,:,:])\n        \ndef save_sag_for_predict(study_id):\n    orig = get_sag_for_predict(study_id)\n    orig = np.array(orig)  * 256\n    slices = orig.shape[0]\n\n    for i in range(slices):\n        cv2.imwrite(f'{save_dir/study_id}_{i}.jpg', orig[i,:,:])","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:15.591829Z","iopub.execute_input":"2022-09-01T15:48:15.594162Z","iopub.status.idle":"2022-09-01T15:48:15.622544Z","shell.execute_reply.started":"2022-09-01T15:48:15.594101Z","shell.execute_reply":"2022-09-01T15:48:15.618929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def show_case(study_id, ext = '.jpg', aspect = 2):\n    fig, axs = plt.subplots(4,4, figsize=(16, 12))\n    axs = axs.flatten()\n    for i in range(len(axs)):\n        fn = f'/kaggle/working/{study_id}_{i}{ext}'\n        if ext == '.npy':\n            img = np.load(fn)\n            axs[i].imshow(img)\n        else:\n            img = mpimg.imread(fn)\n            axs[i].imshow(img, cmap='bone')\n     \n        axs[i].axis(\"off\")\n        axs[i].set_aspect(aspect)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:15.62686Z","iopub.execute_input":"2022-09-01T15:48:15.630173Z","iopub.status.idle":"2022-09-01T15:48:15.648497Z","shell.execute_reply.started":"2022-09-01T15:48:15.630107Z","shell.execute_reply":"2022-09-01T15:48:15.646095Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Look at a sample","metadata":{}},{"cell_type":"code","source":"study_id = '1.2.826.0.1.3680043.12833'\norig, segs = get_sag(study_id)\n\nsave_sag(study_id)\nshow_case(study_id, ext='.jpg')","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:15.650675Z","iopub.execute_input":"2022-09-01T15:48:15.651491Z","iopub.status.idle":"2022-09-01T15:48:28.732366Z","shell.execute_reply.started":"2022-09-01T15:48:15.651447Z","shell.execute_reply":"2022-09-01T15:48:28.731191Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"show_case(study_id, ext='.npy')","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:28.733509Z","iopub.execute_input":"2022-09-01T15:48:28.737083Z","iopub.status.idle":"2022-09-01T15:48:29.828858Z","shell.execute_reply.started":"2022-09-01T15:48:28.737039Z","shell.execute_reply":"2022-09-01T15:48:29.827758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#show_case(study_id)","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:29.830559Z","iopub.execute_input":"2022-09-01T15:48:29.831931Z","iopub.status.idle":"2022-09-01T15:48:29.837598Z","shell.execute_reply.started":"2022-09-01T15:48:29.831885Z","shell.execute_reply":"2022-09-01T15:48:29.835913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Convert all studies","metadata":{}},{"cell_type":"code","source":"from tqdm.auto import tqdm\nfrom joblib import Parallel, delayed\n\n_= Parallel(n_jobs=4)(delayed(save_sag)(study_id) for study_id in tqdm(study_ids))\n\n# for study_id in study_ids:\n#     save_sag(study_id)","metadata":{"execution":{"iopub.status.busy":"2022-09-01T15:48:29.839607Z","iopub.execute_input":"2022-09-01T15:48:29.84032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# list of images\nimage_list = []\nn = 16\nfor study_id in study_ids:\n    for i in range(n):\n        image_list.append(f'{save_dir/study_id}_{i}.jpg') \n        \nlen(image_list), image_list[0]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}