{"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":"code","source":"%matplotlib inline\n# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport glob\nimport os\nimport matplotlib.pyplot as plt\nimport pydicom # for reading dicom files at train, test datasets\nimport tensorflow as tf","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-10-13T14:11:24.96866Z","iopub.execute_input":"2022-10-13T14:11:24.969253Z","iopub.status.idle":"2022-10-13T14:11:30.733085Z","shell.execute_reply.started":"2022-10-13T14:11:24.969148Z","shell.execute_reply":"2022-10-13T14:11:30.731621Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# classes\nclass Patient:\n    \"\"\"class for patient information\"\"\"\n    def __init__(self, name, imnumber):\n        self.name=name\n        self.imnumber=imnumber\n    \n    def makeArray(self, baseDir = 'train'):\n        \"\"\"extract the DICOM image and load it into a numpy array\n        \"\"\"\n        #create path to file\n        if baseDir == 'test':\n            base = tests+'/'+self.name\n        else:\n            base = trains+'/'+self.name\n        pass_dicom = self.imnumber+\".dcm\"\n\n        #extract file\n        filename = pydicom.data.data_manager.get_files(base, pass_dicom)[0]\n        ds = pydicom.dcmread(filename)\n        #extract the array representing the DICOM image\n        return ds.pixel_array\n        ","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:11:30.735874Z","iopub.execute_input":"2022-10-13T14:11:30.736805Z","iopub.status.idle":"2022-10-13T14:11:30.748876Z","shell.execute_reply.started":"2022-10-13T14:11:30.736752Z","shell.execute_reply":"2022-10-13T14:11:30.747263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# functions\ndef loadPatient(name, imgNum, baseDir = \"train\"):\n    \"\"\"create a patient object given a UID and image number\n    \"\"\"\n    # choose path based on whether we want to look in the training or testing data\n    if baseDir == \"test\":\n        path = tests\n    else:\n        path = trains\n    # extract the folder specific to the name to get image numbers\n    path = path + '/{}'.format(name)\n    img_numbers = [s.split('/')[-1][:-4] for s in \n                   glob.glob(path+'/*.dcm')]\n    # create a patient object if the imgNum is in the set of image numbers\n    if imgNum in img_numbers:\n        return Patient(name, imgNum)\n    \ndef idVertebrae(name):\n    '''takes an index and returns the affected vertebrae\n    '''\n    # locate record under the name\n    inst = training[training['StudyInstanceUID']==name]\n    # get list of vertebrae\n    vertebrae = list(training.columns[2:])\n    # map vertebrae names to 0-1 values\n    df = inst[vertebrae]\n    df = pd.DataFrame(np.where(df == 1, df.columns, np.nan), columns=df.columns)\n    # create list naming affected vertebrae\n    return df.loc[0, :].dropna().values.flatten().tolist()\n\ndef synthesizeArray(name, number, baseDir = 'train'):\n    '''synthesize a numpy array for a given UID\n    '''\n    # create patient object\n    image = loadPatient(name, str(number))\n    # build an array of the images\n    try:\n        return image.makeArray(baseDir)\n        \n    except:\n        return np.nan\n    \ndef pictureBreak(name, num, x,y,w,h, baseDir = 'train'):\n    \"\"\"produce an image of the fracture in the bone\n    \"\"\"\n    try:\n        Img = synthesizeArray(name,num, baseDir)\n        a,b,c,d = int(x), int(w+x), int(y), int(y+h)\n        return Img[a:b,c:d]\n    except:\n        return np.nan\n    \ndef retrieve_slices(name, k, baseDir = 'train'):\n    '''return the first k slice numbers for a given patient UID\n    '''\n    if baseDir == 'test':\n        h = glob.glob(tests+'/'+name+'/*.dcm')\n    else:\n        h = glob.glob(trains+'/'+name+'/*.dcm')\n    b = [int(s.split('/')[-1][:-4]) for s in h]\n    return b[:k]\n\ndef make_label(name):\n    '''label UID as healthy or fractured\n    '''\n    myDict = dict(zip(training['StudyInstanceUID'], training['patient_overall']))\n    if myDict[name] == 0:\n        return myDict[name], 'healthy :-)'\n    else:\n        return myDict[name], '-'.join(idVertebrae(name))\n    \ndef normalize_img(data, label):\n    \"\"\"normalize the images in the training data\n    \"\"\"\n    tensor = tf.cast(data, dtype=tf.float32)\n    tensor = tf.math.divide_no_nan(\n            tf.subtract(tensor, tf.reduce_min(tensor)), \n            tf.reduce_max(tensor)\n        )\n    return tensor, label \n\ndef prediction(img_array):\n    \"\"\"predict presence of break using the model\n    \"\"\"\n    #predict probability of vertebrae break using model\n    predictions = model.predict(img_array)\n    #convert prediction to probability\n    score = tf.nn.softmax(predictions[0])\n    #provide score and result\n    return np.max(score), np.argmax(score)\n\ndef make_coords():\n    \"\"\"generate coordinates\n    \"\"\"\n    a = np.random.randint(min(bounds['x']),max(bounds['x']))\n    b = np.random.randint(min(bounds['y']),max(bounds['y']))\n    c = np.random.randint(min(bounds['width']),max(bounds['width']))\n    d = np.random.randint(min(bounds['height']),max(bounds['height']))\n    return a,b,c,d\n\ndef testingImages(name, baseDir = 'train', k=100):\n    '''retrieve and format testing data images for prediction\n    '''\n    pictures = []\n    # pick out first 100 slices\n    d = retrieve_slices(name, k, baseDir)\n    for sl in d:\n        # generate coordinates of the break location\n        x,y,width,height = make_coords()\n        # image of break\n        j = pictureBreak(name, sl, x,y,width,height, baseDir)\n        # return resized if possible\n        try:\n            pictures.append(st.resize(j, (64,64,1)))\n        except:\n            continue\n    return pictures\n\ndef get_best(testimgs):\n    '''get the best prediction from the model\n    '''\n    preDict = {}\n    for img in testimgs:\n        try:\n            img = np.array([img])\n            a,b = prediction(img)\n            preDict[a] = b\n        except:\n            continue\n\n    #m = max(preDict.keys())\n    return preDict\n\ndef get_label(my_dict):\n    '''use the model to make predictions\n    '''\n    # run the model on images in named directory\n    #my_dict = get_best(testingImages(name, baseDir, k))\n    #invert the dictionary to provide label names and average probability\n    my_inverted_dict = {}\n    for k,v in my_dict.items():\n        my_inverted_dict.setdefault(translator[v], list()).append(k)\n    # smooth results\n    stats = {k:np.mean(my_inverted_dict[k]) for k in my_inverted_dict.keys()}\n    # only state that the bone is healthy if there are no predictions indicating fractures\n    if list(stats.keys()) == ['0, healthy :-)']:\n        f = list(stats.keys())[0]\n        return f, stats[f]\n    else:\n        del stats['0, healthy :-)']\n        f = max(stats, key=stats.get)\n        return f, stats[f]\n    \ndef predict_new(name):\n    '''use model to predict on new data outside the training set'''\n    C = retrieve_slices(name,10, baseDir=\"test\")\n    bob = []\n    for v in C:\n        \n        try:\n            a,b,c,d = make_coords()\n            x = Patient(name,str(v)).makeArray(baseDir=\"test\")\n            y= x[a:a+c, b:b+d]\n            bob.append(st.resize(y, (64,64,1)))\n        except:\n            continue\n    try:\n        return get_label(get_best(np.array(bob)))\n    except:\n        return 0\n\ndef predict_affected(name, affected_vertebrae):\n    '''predict if a case is broken and if the stated vertebrae is affected\n    '''\n    vert = ['C1','C2','C3','C4','C5','C6','C7']\n    try:\n        # extract probability and label\n        b, p = predict_new(name)\n        b1 = int(b.split(', ')[0])\n        b2 = b.split(', ')[1].split('-')\n        if affected_vertebrae in vert:\n            # check if the vertebrae is affected by the break\n            if affected_vertebrae in b2:\n                return p\n            else:\n                return 0\n        else:\n            # return the overall probability\n            if b2 == ['healthy :-)']:\n                return 0\n            else:\n                return p\n    except:\n        return 0.0","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:27:53.309685Z","iopub.execute_input":"2022-10-13T14:27:53.310191Z","iopub.status.idle":"2022-10-13T14:27:53.351349Z","shell.execute_reply.started":"2022-10-13T14:27:53.310148Z","shell.execute_reply":"2022-10-13T14:27:53.350113Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load csv files\nbounds = pd.read_csv('../input/rsna-2022-cervical-spine-fracture-detection/train_bounding_boxes.csv')\ntraining = pd.read_csv('../input/rsna-2022-cervical-spine-fracture-detection/train.csv')\ntesting = pd.read_csv('../input/rsna-2022-cervical-spine-fracture-detection/test.csv')\n\n#paths to train/test images\ntrains = '../input/rsna-2022-cervical-spine-fracture-detection/train_images'\ntests = '../input/rsna-2022-cervical-spine-fracture-detection/test_images'","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:11:30.793857Z","iopub.execute_input":"2022-10-13T14:11:30.794335Z","iopub.status.idle":"2022-10-13T14:11:30.853596Z","shell.execute_reply.started":"2022-10-13T14:11:30.794292Z","shell.execute_reply":"2022-10-13T14:11:30.851949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get patients whose data appears in the list of those with bounding boxes\nbbox_p = list(set(bounds['StudyInstanceUID']))\n# get patients whose data appears in segmentation\nseglist = glob.glob('../input/rsna-2022-cervical-spine-fracture-detection/segmentations/*')\nseg_p = [s.split('/')[-1][:-4] for s in seglist]","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:11:30.858223Z","iopub.execute_input":"2022-10-13T14:11:30.859165Z","iopub.status.idle":"2022-10-13T14:11:30.896998Z","shell.execute_reply.started":"2022-10-13T14:11:30.859106Z","shell.execute_reply":"2022-10-13T14:11:30.895545Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get names and slices of healthy vertebrae\nsegs = training[training['StudyInstanceUID'].isin(seg_p)]\nhv = segs[segs['patient_overall']==0]['StudyInstanceUID'].to_list()\n# map the names to the slices\nname2slice = {hv[i]:sorted(retrieve_slices(hv[i], 150)) for i in range(len(hv))}\n# create the dataframe for healthy patients\nnames = []\nslicenums = []\nfor k in name2slice:\n    for v in name2slice[k]:\n        names.append(k)\n        slicenums.append(v)\nhealthy = pd.DataFrame({'StudyInstanceUID':names, \n                        'x': bounds['x'].to_list()[:len(names)],\n                        'y': bounds['y'].to_list()[:len(names)], \n                        'width': bounds['width'].to_list()[:len(names)], \n                        'height': bounds['height'].to_list()[:len(names)], \n                        'slice_number':slicenums})","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:11:30.898679Z","iopub.execute_input":"2022-10-13T14:11:30.899083Z","iopub.status.idle":"2022-10-13T14:11:32.118948Z","shell.execute_reply.started":"2022-10-13T14:11:30.899047Z","shell.execute_reply":"2022-10-13T14:11:32.117225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#draw samples of the control and treatment groups\nA, B = healthy.sample(frac=0.45), bounds.sample(frac=0.45)\n#combine the resulting dataframes\ndf = pd.concat([A, B], axis=0).dropna()","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:11:32.120889Z","iopub.execute_input":"2022-10-13T14:11:32.121429Z","iopub.status.idle":"2022-10-13T14:11:32.142702Z","shell.execute_reply.started":"2022-10-13T14:11:32.121374Z","shell.execute_reply":"2022-10-13T14:11:32.141251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#create arrays of images and labels\na = df.apply(lambda x: pictureBreak(x['StudyInstanceUID'],x[\"slice_number\"], x[\"x\"], x[\"y\"], x[\"width\"], x[\"height\"]),axis=1)\nb = df['StudyInstanceUID'].apply(make_label)\n# assemble into data frame with NaN elements removed\nmyDict = {'img_array': a.tolist(), 'label_array': b.tolist()}\ndf = pd.DataFrame(myDict).dropna()","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:11:32.144679Z","iopub.execute_input":"2022-10-13T14:11:32.145194Z","iopub.status.idle":"2022-10-13T14:13:18.870158Z","shell.execute_reply.started":"2022-10-13T14:11:32.145142Z","shell.execute_reply":"2022-10-13T14:13:18.868943Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import skimage.transform as st\n#prepare feature data for load\nx_train = np.array(df['img_array'])\nx_train = np.array([st.resize(val,(64,64)) for val in x_train])\n#prepare target data for load\ny_train = df[\"label_array\"].tolist()\ny_train = [str(y[0])+\", \"+y[1] for y in y_train]","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:13:18.872357Z","iopub.execute_input":"2022-10-13T14:13:18.872867Z","iopub.status.idle":"2022-10-13T14:13:22.709666Z","shell.execute_reply.started":"2022-10-13T14:13:18.872819Z","shell.execute_reply":"2022-10-13T14:13:22.7083Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"unquies = list(set(y_train))\nnumDict={unquies[i]:i for i in range(len(unquies))}\ntranslator = {v:k for k,v in numDict.items()}\nnumz = [numDict[y] for y in y_train]\ny_train=numz","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:13:22.711319Z","iopub.execute_input":"2022-10-13T14:13:22.712217Z","iopub.status.idle":"2022-10-13T14:13:22.723875Z","shell.execute_reply.started":"2022-10-13T14:13:22.712179Z","shell.execute_reply":"2022-10-13T14:13:22.722238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# data augmentation layer\ndata_augmentation = tf.keras.Sequential([\n  tf.keras.layers.RandomFlip(\"horizontal_and_vertical\"),\n  tf.keras.layers.RandomRotation(0.3),\n  tf.keras.layers.RandomZoom(0.1)\n])\n\ndef acceptor(img, label):\n    '''function to format the data for augmentation\n    '''\n    image = tf.cast(tf.expand_dims(img, -1), tf.float32)\n    return image, label","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:13:22.725653Z","iopub.execute_input":"2022-10-13T14:13:22.726247Z","iopub.status.idle":"2022-10-13T14:13:23.968321Z","shell.execute_reply.started":"2022-10-13T14:13:22.726199Z","shell.execute_reply":"2022-10-13T14:13:23.967094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#create tensorflow data set (the fun part)\nmy_ds = tf.data.Dataset.from_tensor_slices((x_train, y_train))\n#shuffle the data\nmy_ds = my_ds.shuffle(len(df.index))\n#normalize the image so that everything is between 0 and 1\nmy_ds = my_ds.map(normalize_img, num_parallel_calls=tf.data.AUTOTUNE)\n# augment the data\nmy_ds = my_ds.map(acceptor, num_parallel_calls=tf.data.AUTOTUNE)\nmy_ds = my_ds.map(lambda x, y: (data_augmentation(x), y),  num_parallel_calls=4)","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:13:23.969833Z","iopub.execute_input":"2022-10-13T14:13:23.970382Z","iopub.status.idle":"2022-10-13T14:13:24.759788Z","shell.execute_reply.started":"2022-10-13T14:13:23.970347Z","shell.execute_reply":"2022-10-13T14:13:24.75875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#build a simple model\nmodel = tf.keras.Sequential([\n    tf.keras.layers.Conv2D(16, 3, padding='same', activation='relu'),\n    tf.keras.layers.MaxPooling2D(),\n    tf.keras.layers.Conv2D(32, 3, padding='same', activation='relu'),\n    tf.keras.layers.MaxPooling2D(),\n    tf.keras.layers.Conv2D(64, 3, padding='same', activation='relu'),\n    tf.keras.layers.MaxPooling2D(),\n    tf.keras.layers.Dropout(0.15),\n    tf.keras.layers.Flatten(),\n    tf.keras.layers.Dense(128, activation='relu'),\n    tf.keras.layers.Dense(max(numDict.values())+1)\n])\n#compile the model\nmodel.compile(optimizer='adam',\n              loss=tf.keras.losses.SparseCategoricalCrossentropy(from_logits=True),\n              metrics=['accuracy'])","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:13:24.761892Z","iopub.execute_input":"2022-10-13T14:13:24.76314Z","iopub.status.idle":"2022-10-13T14:13:24.811378Z","shell.execute_reply.started":"2022-10-13T14:13:24.763093Z","shell.execute_reply":"2022-10-13T14:13:24.810228Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#separate data set into validation and training sets\nvalid_ds_f = my_ds.take(int(0.2*4432))\ntrain_ds_f = my_ds.skip(int(0.2*4432))\n#batch datasets\nvalid_ds = valid_ds_f.batch(32)\ntrain_ds = train_ds_f.batch(32)","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:13:24.815283Z","iopub.execute_input":"2022-10-13T14:13:24.815641Z","iopub.status.idle":"2022-10-13T14:13:24.829557Z","shell.execute_reply.started":"2022-10-13T14:13:24.81561Z","shell.execute_reply":"2022-10-13T14:13:24.828399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#set the number of epochs\nepochs=100\n#fit the model\noutput = model.fit(train_ds,\n          validation_data=valid_ds,\n          epochs=epochs)","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:13:24.831201Z","iopub.execute_input":"2022-10-13T14:13:24.831545Z","iopub.status.idle":"2022-10-13T14:26:33.915985Z","shell.execute_reply.started":"2022-10-13T14:13:24.831515Z","shell.execute_reply":"2022-10-13T14:26:33.914275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# extract entries in the testing data\nM = glob.glob(tests+\"/*\")\nN = [M[i].split('/')[-1] for i in range(len(M))]","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:28:01.324255Z","iopub.execute_input":"2022-10-13T14:28:01.324747Z","iopub.status.idle":"2022-10-13T14:28:01.333488Z","shell.execute_reply.started":"2022-10-13T14:28:01.324711Z","shell.execute_reply":"2022-10-13T14:28:01.331443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Fix mismatch with test_images folder\ntesting = pd.DataFrame(columns = ['row_id','StudyInstanceUID','prediction_type'])\nfor i in N:\n    for j in ['C1','C2','C3','C4','C5','C6','C7','patient_overall']:\n        testing = testing.append({'row_id':i+'_'+j,'StudyInstanceUID':i,'prediction_type':j},ignore_index=True)\n\ntesting['fractured'] = testing.apply(lambda x: predict_affected(x['StudyInstanceUID'], x['prediction_type']),axis=1)\ntesting = testing.dropna()","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:28:02.998961Z","iopub.execute_input":"2022-10-13T14:28:02.999482Z","iopub.status.idle":"2022-10-13T14:28:15.868218Z","shell.execute_reply.started":"2022-10-13T14:28:02.999443Z","shell.execute_reply":"2022-10-13T14:28:15.866988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# save outout to file\nsub = testing[['row_id', 'fractured']]\nsub.to_csv(\"submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2022-10-13T14:28:15.870286Z","iopub.execute_input":"2022-10-13T14:28:15.870992Z","iopub.status.idle":"2022-10-13T14:28:15.879178Z","shell.execute_reply.started":"2022-10-13T14:28:15.870913Z","shell.execute_reply":"2022-10-13T14:28:15.877977Z"},"trusted":true},"execution_count":null,"outputs":[]}]}