{"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":"# Use Fastai to label images containing CSpine\n\nThere are extra slices above and below the cspine which don't contribute to predictions. It would be helpful to only focus on the pertinent images.","metadata":{}},{"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-15T23:46:15.742034Z","iopub.execute_input":"2022-09-15T23:46:15.742512Z","iopub.status.idle":"2022-09-15T23:46:55.406242Z","shell.execute_reply.started":"2022-09-15T23:46:15.742416Z","shell.execute_reply":"2022-09-15T23:46:55.404489Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Libraries","metadata":{}},{"cell_type":"code","source":"from fastai.vision.all import *\nfrom fastai.medical.imaging import *\nfrom fastcore.all import *\n\nimport pandas as pd\nimport pydicom\nimport numpy as np\n\nimport matplotlib.image as mpimg\nimport cv2\nimport nibabel as nib\n\nfrom tqdm.auto import tqdm","metadata":{"execution":{"iopub.status.busy":"2022-09-15T23:46:55.409604Z","iopub.execute_input":"2022-09-15T23:46:55.410111Z","iopub.status.idle":"2022-09-15T23:46:59.5776Z","shell.execute_reply.started":"2022-09-15T23:46:55.410058Z","shell.execute_reply":"2022-09-15T23:46:59.576704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Gather .csv files","metadata":{}},{"cell_type":"code","source":"root = Path('../input/rsna-2022-cervical-spine-fracture-detection')\ntrain_folder = root/'train_images'\ntest_folder = root/'test_images'\n\nseg_folder = root/'segmentations'\nsegmentations = list(seg_folder.iterdir())\nseg_ids = [o.stem for o in segmentations]\n\ntrain_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\")\nbbdf = pd.read_csv('../input/rsna-2022-cervical-spine-fracture-detection/train_bounding_boxes.csv')\n\ntrain_df.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-15T23:46:59.578755Z","iopub.execute_input":"2022-09-15T23:46:59.580572Z","iopub.status.idle":"2022-09-15T23:46:59.65813Z","shell.execute_reply.started":"2022-09-15T23:46:59.580532Z","shell.execute_reply":"2022-09-15T23:46:59.656762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load segmentation files and get top and bottom slices with CSpine labels\n\n### CSpine labels are 1-7, 0 is background","metadata":{}},{"cell_type":"code","source":"def get_index_boundaries(study_id):\n    seg = nib.load(seg_folder/(study_id + '.nii')).get_fdata()\n    c = {1., 2., 3., 4., 5., 6., 7.}\n    last = 0\n    first = 100000\n    count = seg.shape[2]\n    for i in range(count):\n        arr = seg[:,:,i]\n        rowvals = set(np.unique(arr))\n        if len(c.intersection(rowvals)) > 0:\n            if i < first:\n                first = i\n            if i > last:\n                last = i\n    return first, last, count        \n            \ndef get_slice_boundaries_from_sag(study_id):\n    seg = nib.load(seg_folder/(study_id + '.nii')).get_fdata()\n    c = {1., 2., 3., 4., 5., 6., 7.}\n    n=16\n    spread=160\n    thickness = int(spread/n)\n    slice = int(seg.shape[0]/2 - spread/2)\n\n    max = 0\n    min = 100000\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        vals = set(np.unique(arr))\n        if (c.intersection(vals) == c):\n            #ignore first and last slice because of issues with '1.2.826.0.1.3680043.1363'\n            for i in range(1, arr.shape[0] - 1):\n                row = arr[i]\n                rowvals = set(np.unique(row))\n                if len(c.intersection(rowvals)) > 0:\n                    if i > max:\n                        max = i\n                    if i < min:\n                        min = i\n    return (min, max, seg.shape)","metadata":{"execution":{"iopub.status.busy":"2022-09-15T17:31:12.305305Z","iopub.execute_input":"2022-09-15T17:31:12.305685Z","iopub.status.idle":"2022-09-15T17:31:12.319063Z","shell.execute_reply.started":"2022-09-15T17:31:12.305651Z","shell.execute_reply":"2022-09-15T17:31:12.318158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Create a dataframe of the max and min slices containing any cspine label\n\n### Slices above the cspine are labeled 'head' and below are labeled 'tspine'","metadata":{}},{"cell_type":"code","source":"def get_sorted_fns(study_path):\n    files = []\n    fns = get_dicom_files(study_path)   \n    fn_position_dict = {}\n    for fn in fns:\n        ds = pydicom.dcmread(fn, stop_before_pixels = True)\n        fn_position_dict[fn.stem] = ds.ImagePositionPatient[2]\n        files.append(ds)\n    \n    #sort reverse to get head first, the z-axis is increasing towards the cranial (as opposed to caudal) end of the patient for the neck\n    fn_position_dict = {key: val for key, val in sorted(fn_position_dict.items(), key = lambda ele: ele[1],reverse=True)}\n    sorted_fns = list(fn_position_dict.keys())\n     #returns list of image names sorted by patient position   \n    return sorted_fns\n\ndef process_seg(seg_ids):\n    dict_for_df = []\n\n    for study_id in tqdm(seg_ids):\n        min, max, shape = get_slice_boundaries_from_sag(study_id)\n        fns = get_sorted_fns(train_folder/study_id)\n        fns = [train_folder/study_id/(o + '.dcm') for o in fns]\n        \n        #ignore first and last slice because of issues with '1.2.826.0.1.3680043.1363'\n        for i in range (1,shape[2]-1):\n            label = 'cspine'\n            if i < min:\n                label = 'head'\n            elif i > max:\n                label = 'tspine'\n            fn = Path(fns[i])\n            row = {'study_id':study_id, 'slice':fn.stem,'label':label }\n            dict_for_df.append(row)\n\n    df = pd.DataFrame(dict_for_df)\n    return df\n\n\nif Path('../input/cspine-data/spine_labels.csv').exists():\n    df = pd.read_csv('../input/cspine-data/spine_labels.csv')\nelse:\n    df = process_seg(seg_ids)\n    df.to_csv('spine_labels.csv', index=False)\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-15T18:18:40.642709Z","iopub.execute_input":"2022-09-15T18:18:40.643256Z","iopub.status.idle":"2022-09-15T18:18:40.694167Z","shell.execute_reply.started":"2022-09-15T18:18:40.64321Z","shell.execute_reply":"2022-09-15T18:18:40.692575Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.label.value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-09-15T14:57:05.249775Z","iopub.execute_input":"2022-09-15T14:57:05.25082Z","iopub.status.idle":"2022-09-15T14:57:05.269184Z","shell.execute_reply.started":"2022-09-15T14:57:05.250776Z","shell.execute_reply":"2022-09-15T14:57:05.268132Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_fn(study_id, slice):\n    return train_folder/study_id/f'{slice}.dcm'","metadata":{"execution":{"iopub.status.busy":"2022-09-16T00:11:16.98443Z","iopub.execute_input":"2022-09-16T00:11:16.985015Z","iopub.status.idle":"2022-09-16T00:11:16.990959Z","shell.execute_reply.started":"2022-09-16T00:11:16.984962Z","shell.execute_reply":"2022-09-16T00:11:16.989821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Modified Fastai PILDicom class to get around mode errors","metadata":{}},{"cell_type":"code","source":"#PILDicom2 adapted from https://forums.fast.ai/t/fastai2-problems-with-medical-images-dicom/76138/4 by @deep_derping\nclass PILDicom2(PILBase):\n    \"same as PILDicom but changed pixel type  as np.int16 cannot be handled by PIL\"\n    _open_args,_tensor_cls,_show_args = {},TensorDicom,TensorDicom._show_args\n    @classmethod\n    def create(cls, fn:(Path,str,bytes), mode=None)->None:\n        \"Open a `DICOM file` from path `fn` or bytes `fn` and load it as a `PIL Image`\"\n        if isinstance(fn,bytes): im = Image.fromarray(dcmread_to8bit(pydicom.filebase.DicomBytesIO(fn)))\n        if isinstance(fn,(Path,str)): im = Image.fromarray(dcmread_to8bit(fn))\n        im.load()\n        im = im._new(im.im)\n        return cls(im.convert(mode) if mode else im)\n    \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# apply slope/intercept and W/L\ndef dcmread_to8bit(fn):\n    dcm = fn.dcmread()\n    arr=dcm.pixel_array\n    \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 = intercept + arr * slope\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    arr = (arr * 255).astype(\"uint8\")\n    \n    if \"PhotometricInterpretation\" in dcm and dcm.PhotometricInterpretation == \"MONOCHROME1\":\n        arr = 255 - arr\n        \n    return arr","metadata":{"execution":{"iopub.status.busy":"2022-09-15T14:57:15.610037Z","iopub.execute_input":"2022-09-15T14:57:15.610469Z","iopub.status.idle":"2022-09-15T14:57:15.623501Z","shell.execute_reply.started":"2022-09-15T14:57:15.610418Z","shell.execute_reply":"2022-09-15T14:57:15.622124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fastai datablock","metadata":{}},{"cell_type":"code","source":"datablock = DataBlock(blocks=(ImageBlock(cls=PILDicom2), CategoryBlock),\n                   get_x=lambda x:train_folder/x[0]/f'{x[1]}.dcm',\n                   get_y=lambda x:x[-1],\n                   item_tfms=[Resize(224)],\n                   batch_tfms=[*aug_transforms(size=224),Normalize.from_stats(*imagenet_stats)])\n\ndls = datablock.dataloaders(df.values, num_workers=4)\ndls.vocab","metadata":{"execution":{"iopub.status.busy":"2022-09-14T13:24:03.708759Z","iopub.execute_input":"2022-09-14T13:24:03.710656Z","iopub.status.idle":"2022-09-14T13:24:14.236973Z","shell.execute_reply.started":"2022-09-14T13:24:03.710618Z","shell.execute_reply":"2022-09-14T13:24:14.235652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fastai vision_learner","metadata":{}},{"cell_type":"code","source":"learn = vision_learner(dls, resnet34, metrics=accuracy).to_fp16()\nlearn.fine_tune(10)","metadata":{"execution":{"iopub.status.busy":"2022-09-14T13:24:14.242022Z","iopub.execute_input":"2022-09-14T13:24:14.242615Z","iopub.status.idle":"2022-09-14T14:36:22.60846Z","shell.execute_reply.started":"2022-09-14T13:24:14.242568Z","shell.execute_reply":"2022-09-14T14:36:22.607409Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Save model","metadata":{}},{"cell_type":"code","source":"learn.save('spine_label.h5')","metadata":{"execution":{"iopub.status.busy":"2022-09-14T14:36:22.611829Z","iopub.execute_input":"2022-09-14T14:36:22.612284Z","iopub.status.idle":"2022-09-14T14:36:23.081763Z","shell.execute_reply.started":"2022-09-14T14:36:22.612247Z","shell.execute_reply":"2022-09-14T14:36:23.080671Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# !mkdir  'models'\n# !cp  '../input/cspine-data/spine_label.h5.pth' 'models'","metadata":{"execution":{"iopub.status.busy":"2022-09-14T14:36:23.083401Z","iopub.execute_input":"2022-09-14T14:36:23.083777Z","iopub.status.idle":"2022-09-14T14:36:23.088099Z","shell.execute_reply.started":"2022-09-14T14:36:23.083737Z","shell.execute_reply":"2022-09-14T14:36:23.087112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# learn.load('spine_label.h5')","metadata":{"execution":{"iopub.status.busy":"2022-09-14T14:36:23.08971Z","iopub.execute_input":"2022-09-14T14:36:23.090088Z","iopub.status.idle":"2022-09-14T14:36:23.099157Z","shell.execute_reply.started":"2022-09-14T14:36:23.090047Z","shell.execute_reply":"2022-09-14T14:36:23.098045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#learn.show_results(max_n=16)","metadata":{"execution":{"iopub.status.busy":"2022-09-14T14:36:23.100617Z","iopub.execute_input":"2022-09-14T14:36:23.101144Z","iopub.status.idle":"2022-09-14T14:36:23.109703Z","shell.execute_reply.started":"2022-09-14T14:36:23.101108Z","shell.execute_reply":"2022-09-14T14:36:23.108737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# interp = Interpretation.from_learner(learn)\n# interp.plot_top_losses(2)","metadata":{"execution":{"iopub.status.busy":"2022-09-14T14:36:23.111352Z","iopub.execute_input":"2022-09-14T14:36:23.111974Z","iopub.status.idle":"2022-09-14T14:36:23.119693Z","shell.execute_reply.started":"2022-09-14T14:36:23.111936Z","shell.execute_reply":"2022-09-14T14:36:23.118819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Create dataframe of the remaining training images for prediction","metadata":{}},{"cell_type":"code","source":"#get all image names\ncsv = Path('../input/cspine-data/dcm_filenames.csv')\nif csv.exists():\n    fndf = pd.read_csv('../input/cspine-data/dcm_filenames.csv')\nelse:\n    fn_glob = train_folder.glob('**/*.dcm')\n    fns = [fn for fn in fn_glob]\n    fndf = pd.DataFrame({'fn':fns})\n    fndf.to_csv('dcm_filenames.csv', index=False)\n\nfndf['study_id'] = fndf.fn.apply(lambda x: Path(x).parent.name)\nfndf['slice'] = fndf.fn.apply(lambda x: Path(x).stem)\nfndf = fndf.drop('fn', axis=1)\nfndf.slice = fndf.slice.astype(int)\nfndf = fndf.sort_values(by = ['study_id', 'slice'])\nfndf.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-16T00:22:42.726917Z","iopub.execute_input":"2022-09-16T00:22:42.72751Z","iopub.status.idle":"2022-09-16T00:22:57.985056Z","shell.execute_reply.started":"2022-09-16T00:22:42.727469Z","shell.execute_reply":"2022-09-16T00:22:57.983945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## List of images in study at intervals of 50/count. Take first and last 12.\n\nThis assumes that the central 50% of the volume includes CSpine and make predictions on the front and back 25%.","metadata":{}},{"cell_type":"code","source":"def get_indices(count):\n    step = max(int(count/50), 1)\n    arr = np.arange(1,count,step)\n    retarr = arr[:12].copy()\n    retarr = np.append(retarr, arr[-12:])\n    return retarr\n\ndef get_indices_for_pred(group):\n    count = group.shape[0]\n    indices = get_indices(count)\n    ret = [o in indices for o in range(count)]\n    return ret\n\nfndf['should_predict'] = fndf.groupby('study_id').transform(get_indices_for_pred)\nfndf.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-15T17:31:29.696751Z","iopub.execute_input":"2022-09-15T17:31:29.697294Z","iopub.status.idle":"2022-09-15T17:31:35.337985Z","shell.execute_reply.started":"2022-09-15T17:31:29.697246Z","shell.execute_reply":"2022-09-15T17:31:35.336871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Only predict on images not in segmentation group","metadata":{}},{"cell_type":"code","source":"train_ids = train_df.StudyInstanceUID.unique()\nstudy_ids = list(set(train_ids) - set(seg_ids))\n\ntdf = fndf[(fndf.should_predict == True) & (fndf.study_id.isin(study_ids))]\ntdf = tdf[['study_id','slice']]\ntdf.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-15T15:24:44.769739Z","iopub.execute_input":"2022-09-15T15:24:44.770141Z","iopub.status.idle":"2022-09-15T15:24:44.824204Z","shell.execute_reply.started":"2022-09-15T15:24:44.770111Z","shell.execute_reply":"2022-09-15T15:24:44.822926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tdf.to_csv('nonseg_image_paths.csv', index = False)","metadata":{"execution":{"iopub.status.busy":"2022-09-15T15:26:48.262827Z","iopub.execute_input":"2022-09-15T15:26:48.263305Z","iopub.status.idle":"2022-09-15T15:26:48.338785Z","shell.execute_reply.started":"2022-09-15T15:26:48.263271Z","shell.execute_reply":"2022-09-15T15:26:48.337733Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Predict","metadata":{}},{"cell_type":"code","source":"test_dl = learn.dls.test_dl(tdf.values, num_workers=1)\npreds,_ = learn.get_preds( dl = test_dl)","metadata":{"execution":{"iopub.status.busy":"2022-09-14T14:39:04.514789Z","iopub.execute_input":"2022-09-14T14:39:04.515287Z","iopub.status.idle":"2022-09-14T14:59:13.681681Z","shell.execute_reply.started":"2022-09-14T14:39:04.515257Z","shell.execute_reply":"2022-09-14T14:59:13.679178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Convert argmax of predictions to labels","metadata":{}},{"cell_type":"code","source":"vocab = learn.dls.vocab\npred_labels = np.array(np.argmax(preds, axis = 1))\nlabels = [vocab[o] for o in pred_labels] \n\npred_labels[:3], labels[:3]","metadata":{"execution":{"iopub.status.busy":"2022-09-14T14:59:13.685589Z","iopub.execute_input":"2022-09-14T14:59:13.686671Z","iopub.status.idle":"2022-09-14T14:59:13.973225Z","shell.execute_reply.started":"2022-09-14T14:59:13.686623Z","shell.execute_reply":"2022-09-14T14:59:13.972293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Add labels to dataframe","metadata":{}},{"cell_type":"code","source":"tdf['label'] = labels\n\nlabel_dict = {'head':0, 'cspine':1, 'tspine':0}\ntdf['label_int'] = tdf.label.apply(lambda x: label_dict[x])\ntdf.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-14T14:59:13.984182Z","iopub.execute_input":"2022-09-14T14:59:13.985181Z","iopub.status.idle":"2022-09-14T14:59:14.039129Z","shell.execute_reply.started":"2022-09-14T14:59:13.985144Z","shell.execute_reply":"2022-09-14T14:59:14.038048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#fn = train_folder/f'1.2.826.0.1.3680043.1479/1.dcm'\n# fn = train_folder/f'1.2.826.0.1.3680043.16206/57.dcm'\n# im = PILDicom2.create(fn)\n# plt.imshow(im)","metadata":{"execution":{"iopub.status.busy":"2022-09-14T14:59:14.040411Z","iopub.execute_input":"2022-09-14T14:59:14.040846Z","iopub.status.idle":"2022-09-14T14:59:14.336925Z","shell.execute_reply.started":"2022-09-14T14:59:14.04081Z","shell.execute_reply":"2022-09-14T14:59:14.33595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tdf.to_csv('preds.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2022-09-14T14:59:14.338655Z","iopub.execute_input":"2022-09-14T14:59:14.339021Z","iopub.status.idle":"2022-09-14T14:59:14.435494Z","shell.execute_reply.started":"2022-09-14T14:59:14.33897Z","shell.execute_reply":"2022-09-14T14:59:14.434474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Smooth out predictions\n\nThe predictions go through 2 rounds of simple template smoothing. If a CSpine is between 2 non-Cspine images, that image gets demoted and if a non-CSpine image is between two positives, it gets promoted.\n\nThe large central 50% gap between slices is expanded on either side if there are CSpine images. This prevents false positives at first and last slice to be ignored, since those would not be smoothed out.","metadata":{}},{"cell_type":"code","source":"window_size = 3\n  \n# Convert array of integers to pandas series\n\ndef smooth_step(numbers_series):\n    # Get the window of series\n    # of observations of specified window size\n    windows = numbers_series.rolling(window_size)\n\n    for w in windows:\n        if len(w == 3):\n            i = w.keys().start + 1\n            if np.array_equal(w.values, [1, 0, 1]):\n                #print('0 to 1 ', i)\n                numbers_series[i] = 1\n            if np.array_equal(w.values, [0, 1, 0]):\n                #print('1 to 0 ', i)\n                numbers_series[i] = 0\n    return numbers_series\n\ndef smooth(arr):\n    numbers_series = pd.Series(arr)\n    numbers_series = smooth_step(numbers_series)\n    numbers_series = smooth_step(numbers_series)\n    return numbers_series\n\ndef get_range(df_group):\n    study_assess = pd.Series(dtype='object')\n\n    #above\n    above = df_group.label_int.values[:12].copy()\n    above = smooth(above)\n    not_cspine = np.where(above != 1)\n    i = not_cspine[0]\n    if len(i) == 0:\n        i = 0\n    else:\n        i = i.max()\n    \n    slices = df_group.slice.values\n    top_range = slices[i] + 1\n    \n    #below\n    below = df_group.label_int.values[12:].copy()\n    below = smooth(below)\n    not_cspine = np.where(below != 1)\n    i = not_cspine[0]\n    if len(i) == 0:\n        i = len(below)-1\n    else:\n        i = i.min()\n\n    if (len(above) + i) >= len(slices):\n        print(study_id, len(above) + i, len(slices), len(above), i, len(below))\n        bottom_range = slices[-1]\n    else:\n        bottom_range = slices[len(above) + i] - 1\n\n    study_assess['top_range'] = top_range\n    study_assess['bottom_range'] = bottom_range\n    return study_assess\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-15T23:48:11.087498Z","iopub.execute_input":"2022-09-15T23:48:11.08794Z","iopub.status.idle":"2022-09-15T23:48:11.103751Z","shell.execute_reply.started":"2022-09-15T23:48:11.087902Z","shell.execute_reply":"2022-09-15T23:48:11.102339Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred_range_df = tdf.groupby('study_id').apply(get_range).reset_index()\npred_range_df.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-15T23:48:11.364194Z","iopub.execute_input":"2022-09-15T23:48:11.364612Z","iopub.status.idle":"2022-09-15T23:48:33.44057Z","shell.execute_reply.started":"2022-09-15T23:48:11.364574Z","shell.execute_reply":"2022-09-15T23:48:33.439078Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Range of images in segmented set","metadata":{}},{"cell_type":"code","source":"# segs\ndf = pd.read_csv('../input/cspine-data/spine_labels.csv')","metadata":{"execution":{"iopub.status.busy":"2022-09-16T00:24:34.693784Z","iopub.execute_input":"2022-09-16T00:24:34.695007Z","iopub.status.idle":"2022-09-16T00:24:34.73676Z","shell.execute_reply.started":"2022-09-16T00:24:34.694952Z","shell.execute_reply":"2022-09-16T00:24:34.735514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"seg_range_df = df[df.label == 'cspine'].groupby('study_id')['slice'].agg(top_range = 'min',bottom_range = 'max').reset_index()\nseg_range_df.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-16T00:24:35.239489Z","iopub.execute_input":"2022-09-16T00:24:35.240747Z","iopub.status.idle":"2022-09-16T00:24:35.271306Z","shell.execute_reply.started":"2022-09-16T00:24:35.240687Z","shell.execute_reply":"2022-09-16T00:24:35.270073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Concat both predicted and segmented and save to .csv","metadata":{}},{"cell_type":"code","source":"final_df = pd.concat([pred_range_df, seg_range_df])\nfinal_df.to_csv('study_ranges_with_cspine.csv', index = False)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T00:24:36.987091Z","iopub.execute_input":"2022-09-16T00:24:36.987981Z","iopub.status.idle":"2022-09-16T00:24:37.005347Z","shell.execute_reply.started":"2022-09-16T00:24:36.987927Z","shell.execute_reply":"2022-09-16T00:24:37.004302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"final_df.shape","metadata":{"execution":{"iopub.status.busy":"2022-09-15T20:36:56.320335Z","iopub.execute_input":"2022-09-15T20:36:56.321328Z","iopub.status.idle":"2022-09-15T20:36:56.330104Z","shell.execute_reply.started":"2022-09-15T20:36:56.321272Z","shell.execute_reply":"2022-09-15T20:36:56.328904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Show a few examples with predicted ranges","metadata":{}},{"cell_type":"markdown","source":"### List of filenames for 3D","metadata":{}},{"cell_type":"code","source":"def filename_list_from_range(study_id, first, last, image_folder):\n    return [image_folder/study_id/f'{o}.dcm' for o in range(first, last + 1)]","metadata":{"execution":{"iopub.status.busy":"2022-09-15T23:48:33.500072Z","iopub.execute_input":"2022-09-15T23:48:33.500882Z","iopub.status.idle":"2022-09-15T23:48:33.507813Z","shell.execute_reply.started":"2022-09-15T23:48:33.500818Z","shell.execute_reply":"2022-09-15T23:48:33.50652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_dicom(dcm, voi_lut = True, fix_monochrome = True):\n    '''ref: https://www.kaggle.com/code/raddar/convert-dicom-to-np-array-the-correct-way\n    '''\n    \n    # VOI LUT (if available by DICOM device) is used to transform raw DICOM data to \"human-friendly\" view\n    if voi_lut:\n        data = pydicom.pixel_data_handlers.util.apply_voi_lut(dcm.pixel_array, dcm)\n    else:\n        data = dcm.pixel_array\n               \n    # depending on this value, X-ray may look inverted - fix that:\n    if fix_monochrome and dcm.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n        \n    data = data - np.min(data)\n    data = data / np.max(data)\n    data = (data * 255).astype(np.uint8)\n    return data","metadata":{"execution":{"iopub.status.busy":"2022-09-15T23:50:36.940263Z","iopub.execute_input":"2022-09-15T23:50:36.940686Z","iopub.status.idle":"2022-09-15T23:50:36.949356Z","shell.execute_reply.started":"2022-09-15T23:50:36.940653Z","shell.execute_reply":"2022-09-15T23:50:36.947641Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Adapted from Pydicom: 'Load CT slices and plot axial, sagittal and coronal images'\ndef create_3d_from_filenames(filenames):\n    files = [pydicom.dcmread(fn) for fn in filenames]\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 = read_dicom(s, voi_lut = False)\n        img3d[:, :, i] = img2d\n        \n    return img3d, ps, ss","metadata":{"execution":{"iopub.status.busy":"2022-09-15T23:54:27.202521Z","iopub.execute_input":"2022-09-15T23:54:27.20441Z","iopub.status.idle":"2022-09-15T23:54:27.217159Z","shell.execute_reply.started":"2022-09-15T23:54:27.204359Z","shell.execute_reply":"2022-09-15T23:54:27.216128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#df = pred_range_df\ndef get_center_sag(img3d):\n    i = int(img3d.shape[1]/2)\n    \n    arr = img3d[:,i,:] \n    #flip and rotate so that C1 is at the top\n    arr = np.transpose(arr)\n    arr = np.flip(arr, 0)\n   \n    return arr\n\ndef create_center_image(index, range_df):\n    row = range_df.iloc[index]\n    study_id = row.study_id\n    first = row.top_range\n    last = row.bottom_range\n    filenames = filename_list_from_range(study_id, first, last, train_folder)\n\n    img3d, _, _ = create_3d_from_filenames(filenames)\n    arr = get_center_sag(img3d)\n    im = Image.fromarray(arr)\n    return im, first, last, study_id","metadata":{"execution":{"iopub.status.busy":"2022-09-16T00:05:31.332908Z","iopub.execute_input":"2022-09-16T00:05:31.333779Z","iopub.status.idle":"2022-09-16T00:05:31.34272Z","shell.execute_reply.started":"2022-09-16T00:05:31.333734Z","shell.execute_reply":"2022-09-16T00:05:31.341125Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def show_cases(indices, aspect = 2):\n    fig, axs = plt.subplots(4,4, figsize=(16, 12))\n    axs = axs.flatten()\n    for i in range(len(axs)):\n        img, first, last, study_id = create_center_image(indices[i], pred_range_df)\n        axs[i].set_title(f'{study_id}\\n {first} - {last}')\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-16T00:06:29.899687Z","iopub.execute_input":"2022-09-16T00:06:29.900986Z","iopub.status.idle":"2022-09-16T00:06:29.909383Z","shell.execute_reply.started":"2022-09-16T00:06:29.90094Z","shell.execute_reply":"2022-09-16T00:06:29.90834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Sample images\n\nThe images are getting the entire cspine although it may be prudent to give a small cushion above and below.","metadata":{}},{"cell_type":"code","source":"indices = [0,100, 20,30, 300, 500,550, 600, 650, 700, 750, 800, 850, 900, 950, 1000]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"show_cases(indices)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T00:06:32.393556Z","iopub.execute_input":"2022-09-16T00:06:32.394002Z","iopub.status.idle":"2022-09-16T00:07:50.124224Z","shell.execute_reply.started":"2022-09-16T00:06:32.393966Z","shell.execute_reply":"2022-09-16T00:07:50.122992Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Image Count: Before & After","metadata":{}},{"cell_type":"code","source":"num_images_before = fndf.shape[0]\n\ndff = final_df.copy()\ndff['total'] = dff.bottom_range - dff.top_range + 1\nnum_images_after = dff.total.sum()\n\nnum_images_before, num_images_after, num_images_before - num_images_after","metadata":{"execution":{"iopub.status.busy":"2022-09-16T00:28:42.56575Z","iopub.execute_input":"2022-09-16T00:28:42.566405Z","iopub.status.idle":"2022-09-16T00:28:42.582427Z","shell.execute_reply.started":"2022-09-16T00:28:42.566353Z","shell.execute_reply":"2022-09-16T00:28:42.580918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1/3 of the images are removed","metadata":{}},{"cell_type":"code","source":"num_images_after/num_images_before","metadata":{"execution":{"iopub.status.busy":"2022-09-16T00:29:00.201213Z","iopub.execute_input":"2022-09-16T00:29:00.201635Z","iopub.status.idle":"2022-09-16T00:29:00.210174Z","shell.execute_reply.started":"2022-09-16T00:29:00.201592Z","shell.execute_reply":"2022-09-16T00:29:00.208478Z"},"trusted":true},"execution_count":null,"outputs":[]}]}