{"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":"import os\nimport requests\nimport glob\nimport cv2 as cv\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport plotly.graph_objects as go\nimport seaborn as sns\nimport random\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport nibabel as nib\nimport pydicom\nimport math\n\nfrom path import Path\nfrom tqdm import tqdm\nimport nibabel as nib\nimport pydicom as dicom\nfrom pydicom import dcmread\nfrom skimage import measure\nfrom skimage.segmentation import clear_border\nfrom skimage.measure import label,regionprops\nfrom scipy import ndimage as ndi\nfrom scipy.ndimage import measurements,center_of_mass,binary_dilation,zoom\nimport pickle\n\nfrom tensorflow import keras\nimport tensorflow as tf\nimport tensorflow_hub as hub\nfrom keras import layers\nfrom tensorflow.keras.utils import to_categorical\nfrom sklearn.model_selection import StratifiedKFold\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nfrom keras_preprocessing.image import load_img, img_to_array\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.layers import Input, Conv3D, MaxPooling3D, concatenate, Conv3DTranspose, BatchNormalization, Dropout, Lambda,Dense\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.metrics import MeanIoU\nkernel_initializer =  'he_uniform'","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train.csv\")\ntest = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/test.csv\")\ntrain_b = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train_bounding_boxes.csv\")\nsubmission = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/sample_submission.csv\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.sort_values('StudyInstanceUID',inplace=True,axis = 0)\ntrain_b.sort_values('StudyInstanceUID',inplace=True,axis = 0)\ntrain.reset_index(inplace=True,drop= True)\ntrain_b.reset_index(inplace=True,drop= True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_dir =\"../input/rsna-2022-cervical-spine-fracture-detection/train_images\"\ntest_dir =\"../input/rsna-2022-cervical-spine-fracture-detection/test_images\"\nseg_dir =\"../input/rsna-2022-cervical-spine-fracture-detection/segmentations\"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patients = os.listdir(train_dir)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# General EDA  ","metadata":{}},{"cell_type":"code","source":"images=[]\nfor p in patients[:1]:\n    dir = os.path.join(train_dir,p)\n    FileName = os.listdir(os.path.join(train_dir,p))\n    for index in range(150,168):\n        sec = str(index)+'.dcm'\n        image_path= os.path.join(dir,sec)\n        print(image_path)\n        images.append(image_path)\n        # print(len(FileName),images.pixel_array.shape)\n        # plt.pcolormesh(image.pixel_array)\n        # plt.colorbar()\n        # plt.show()\n    fig, axes = plt.subplots(nrows=3, ncols=6, figsize=(24,12))\n    fig.suptitle(f'ID: {p}', weight=\"bold\", size=20)\n\n    j=0\n    for i in range(150,168,1):\n        slice_no = i\n        # Plot the image\n        x = (j) // 6\n        y = (j) % 6\n        img = dicom.dcmread(images[j])\n        j+=1\n        axes[x, y].imshow(img.pixel_array, cmap=\"bone\")\n        axes[x, y].set_title(f\"Slice: {slice_no}\", fontsize=14, weight='bold')\n        axes[x, y].axis('off')\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom matplotlib import animation, rc; rc('animation', html='jshtml')\nimport re\nimport os\nimport cv2\nimport gc\nfrom glob import glob\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_animation(patient_id, drc='train_images', save_dir=\"D:\\Projects\\Kaggle\\RSNA\", save=True, fps=10):\n    # Get paths\n    base_path = \"D:\\Projects\\Kaggle\\RSNA\"\n    dcm_paths = glob(f\"{base_path}/{drc}/{patient_id}/*\")\n    def atoi(text):\n        return int(text) if text.isdigit() else text\n    def natural_keys(text):\n        return [atoi(c) for c in re.split(r'(\\d+)', text)]\n    dcm_paths.sort(key=natural_keys)\n    \n    # Get images\n    files = [pydicom.dcmread(path) for path in dcm_paths]\n    images = [apply_voi_lut(file.pixel_array, file) for file in files]\n    \n    # Optimise memory\n    images_small = [np.clip(cv2.resize(img,dsize=[256,256]),-32768,32767).astype('int16') for img in images]\n    \n    # Stack images\n    animation_arr = np.stack(images_small, axis=0)\n\n    del images, images_small\n    gc.collect()\n    \n    # Initialise plot\n    fig = plt.figure(figsize=(3,3))  # if size is too big then gif gets truncated\n    im = plt.imshow(animation_arr[0], cmap='bone')\n    plt.axis('off')\n    plt.title(f\"{patient_id}\", fontweight=\"bold\")\n    \n    # Load next frame\n    def animate_func(i):\n        im.set_array(animation_arr[i])\n        return [im]\n    plt.close()\n    \n    # Animation function\n    anim = animation.FuncAnimation(fig, animate_func, frames = animation_arr.shape[0], interval = 1000//fps)\n    \n    # Save\n    if save:\n        os.makedirs(save_dir, exist_ok=True)\n        anim.save(os.path.join(save_dir, f\"patient_{patient_id}.gif\"), fps=10, writer='imagemagick')\n        \n    return anim","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"create_animation('1.2.826.0.1.3680043.18659', fps=30)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"absolute_path = 'D:\\Projects\\Kaggle\\RSNA\\segmentations'\nbase_path = os.path.join(absolute_path)\nfull_directory = glob(os.path.join(absolute_path, '*'))\nprint(full_directory)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# NOT SO GOOD GIF CODE","metadata":{}},{"cell_type":"code","source":"def create_gif(input_image, title='.gif', filename='test.gif'):\n    # see example from matplotlib documentation\n    import imageio\n    import matplotlib.animation as animate\n    images = []\n    input_image_data = input_image.get_fdata()\n    fig = plt.figure()\n    for i in range(len(input_image_data)):\n        im = plt.imshow(input_image_data[i], animated=True)\n        images.append([im])\n    \n    ani = animate.ArtistAnimation(fig, images, interval=25,\\\n        blit=True, repeat_delay=500)\n    plt.title(title, fontsize=20)\n    plt.axis('off')\n    ani.save(filename)\n    plt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for val in full_directory[1:2]:\n    create_gif(nib.load(val),p)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pyglet\nanimation = pyglet.image.load_animation('test.gif')\nbin = pyglet.image.atlas.TextureBin()\nanimation.add_to_texture_bin(bin)\nsprite = pyglet.sprite.Sprite(img=animation)\n\nwindow = pyglet.window.Window()\n\n@window.event\ndef on_draw():\n    window.clear()\n    sprite.draw()\n\npyglet.app.run()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(20,6))\nplt.subplot(1,2,1)\nax1 = sns.countplot(data=train, x='patient_overall')\nfor container in ax1.containers:\n    ax1.bar_label(container)\nplt.title('Fractures by patient')\nplt.ylim([0,1300])\n\n# Unpivot train_df for plotting\ntrain_melt = pd.melt(train, id_vars = ['StudyInstanceUID', 'patient_overall'],\n             value_vars = ['C1','C2','C3','C4','C5','C6','C7'],\n             var_name=\"Vertebrae\",\n             value_name=\"Fractured\")\n\nplt.subplot(1,2,2)\nax2 = sns.countplot(data=train_melt, x='Vertebrae', hue='Fractured')\nfor container in ax2.containers:\n    ax2.bar_label(container)\nplt.title('Fractures by vertebrae')\nplt.ylim([0,2800])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12,5))\nax = sns.countplot(x = train[['C1','C2','C3','C4','C5','C6','C7']].sum(axis=1))\nfor container in ax.containers:\n    ax.bar_label(container)\nplt.title('Number of fractures by patient')\nplt.xlabel('Number of fractures')\nplt.ylim([0,1300])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(6,5))\nsns.heatmap(train[['C1','C2','C3','C4','C5','C6','C7']].corr(), cmap='bwr', vmin=-1, vmax=1)\nplt.title('Correlations')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"reminder Need to check the depth for each patient","metadata":{}},{"cell_type":"code","source":"train_b","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Testing the data cell ","metadata":{}},{"cell_type":"code","source":"path = os.path.join(train_dir,'1.2.826.0.1.3680043.10001')\n# print(path)\nfor s in os.listdir(path):\n    ds =dicom.read_file(path + '/' + s)\n    data = apply_voi_lut(ds.pixel_array, ds)\n\n    data = data - np.min(data)\n    data = data / np.max(data)\n    data = (data * 255).astype(np.uint8)\n# print(data)\n        \nslices = [dicom.dcmread(path + '/' + s) for s in os.listdir(path)]\nslices.sort( key = lambda x: int(x.ImagePositionPatient[2]))\n\n# print(slices)\nnew_slices = []\nslices = [cv.resize(np.array(each_slice.pixel_array),(128,128)) for each_slice in slices]\ndata=[]\nfor slice in slices:\n    data1 = slice - np.min(slice)\n    data1 = data1 / np.max(data1)\n    data1 = (data1 * 255).astype(np.uint8)\n    data.append(data1)\nprint(np.array(data).shape)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Chunks creator \nfor every patient there is diffirent size depth so I created chuck of size 16 to make the depth of equal size","metadata":{}},{"cell_type":"code","source":"def chunks(l, n):\n    # Credit: Ned Batchelder\n    # Link: http://stackoverflow.com/questions/312443/how-do-you-split-a-list-into-evenly-sized-chunks\n    \"\"\"Yield successive n-sized chunks from l.\"\"\"\n    for i in range(0, len(l), n):\n        yield l[i:i + n]\n\n\ndef mean(a):\n    return sum(a) / len(a)\n\ndef train_process_data(patient,labels_df,img_px_size=128, hm_slices=20, visualize=False):\n    cur_label = []\n    trainlabel = []\n    cur_label.append(int(labels_df.loc[labels_df['StudyInstanceUID']==patient,\"C1\"].any()))\n    cur_label.append(int(labels_df.loc[labels_df['StudyInstanceUID']==patient,\"C2\"].any()))\n    cur_label.append(int(labels_df.loc[labels_df['StudyInstanceUID']==patient,\"C3\"].any()))\n    cur_label.append(int(labels_df.loc[labels_df['StudyInstanceUID']==patient,\"C4\"].any()))\n    cur_label.append(int(labels_df.loc[labels_df['StudyInstanceUID']==patient,\"C5\"].any()))\n    cur_label.append(int(labels_df.loc[labels_df['StudyInstanceUID']==patient,\"C6\"].any()))\n    cur_label.append(int(labels_df.loc[labels_df['StudyInstanceUID']==patient,\"C7\"].any()))    \n    trainlabel += [cur_label]\n    # print(label)\n\n    path = os.path.join(train_dir,patient)\n    # print(path)\n    slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: int(x.ImagePositionPatient[2]))\n\n    new_slices = []\n    slices = [cv.resize(np.array(each_slice.pixel_array),(img_px_size,img_px_size)) for each_slice in slices]\n    # data=[]\n\n    # for slice in slices:\n    #     data1 = slice - np.min(slice)\n    #     data1 = data1 / np.max(data1)\n    #     data1 = (data1).astype(np.uint8)\n    #     data.append(data1)\n    # print(data)\n\n    chunk_sizes = math.ceil(len(slices) / hm_slices)\n    for slice_chunk in chunks(slices, chunk_sizes):\n        slice_chunk = list(map(mean, zip(*slice_chunk)))\n        new_slices.append(slice_chunk)\n\n    if len(new_slices) == hm_slices-1:\n            new_slices.append(new_slices[-1])\n\n    if len(new_slices) == hm_slices-2:\n            new_slices.append(new_slices[-1])\n            new_slices.append(new_slices[-1])\n    if len(new_slices) == hm_slices-3:\n            new_slices.append(new_slices[-1])\n            new_slices.append(new_slices[-1])\n            new_slices.append(new_slices[-1])\n\n    if len(new_slices) == hm_slices+2:\n            new_val = list(map(mean, zip(*[new_slices[hm_slices-1],new_slices[hm_slices],])))\n            del new_slices[hm_slices]\n            new_slices[hm_slices-1] = new_val\n\n    if len(new_slices) == hm_slices+1:\n            new_val = list(map(mean, zip(*[new_slices[hm_slices-1],new_slices[hm_slices],])))\n            del new_slices[hm_slices]\n            new_slices[hm_slices-1] = new_val\n        \n\n    return np.array(new_slices),np.array(trainlabel)\n\n#                                               stage 1 for real.\n\n\nmuch_data =[]\nfor num,patient in enumerate(patients):\n    if num % 100 == 0:\n        print(num)\n    try:\n        img_data,label = train_process_data(patient,train,img_px_size=128, hm_slices=16)\n        # print(img_data.shape,label)\n        much_data.append([img_data,label])\n    except KeyError as e:\n        print('This is unlabeled data!')\n        pass\n\nnp.save('muchdata-{}-{}-{}-label7.npy'.format(128,128,16), much_data)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def test_process_data(patient,labels_df,img_px_size=128, hm_slices=20, visualize=False):\n\n    \n    # print(label)\n\n    path = os.path.join(test_dir,patient)\n    # print(path)\n    slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: int(x.ImagePositionPatient[2]))\n\n    new_slices = []\n    slices = [cv.resize(np.array(each_slice.pixel_array),(img_px_size,img_px_size)) for each_slice in slices]\n#     data=[]\n\n#     for slice in slices:\n#         data1 = slice - np.min(slice)\n#         data1 = data1 / np.max(data1)\n#         data1 = (data1).astype(np.uint8)\n#         data.append(data1)\n    # print(data)\n\n    chunk_sizes = math.ceil(len(slices) / hm_slices)\n    for slice_chunk in chunks(slices, chunk_sizes):\n        slice_chunk = list(map(mean, zip(*slice_chunk)))\n        new_slices.append(slice_chunk)\n\n    if len(new_slices) == hm_slices-1:\n            new_slices.append(new_slices[-1])\n\n    if len(new_slices) == hm_slices-2:\n            new_slices.append(new_slices[-1])\n            new_slices.append(new_slices[-1])\n    if len(new_slices) == hm_slices-3:\n            new_slices.append(new_slices[-1])\n            new_slices.append(new_slices[-1])\n            new_slices.append(new_slices[-1])\n\n    if len(new_slices) == hm_slices+2:\n            new_val = list(map(mean, zip(*[new_slices[hm_slices-1],new_slices[hm_slices],])))\n            del new_slices[hm_slices]\n            new_slices[hm_slices-1] = new_val\n\n    if len(new_slices) == hm_slices+1:\n            new_val = list(map(mean, zip(*[new_slices[hm_slices-1],new_slices[hm_slices],])))\n            del new_slices[hm_slices]\n            new_slices[hm_slices-1] = new_val\n        \n\n    return np.array(new_slices)\n\n#                                               stage 1 for real.\n\n\ntest_data =[]\npatients = os.listdir(test_dir)\nfor num,patient in enumerate(patients):\n    \n    print(num)\n    try:\n        img_data = test_process_data(patient,train,img_px_size=128, hm_slices=16)\n        # print(img_data.shape,label)\n        test_data.append([img_data,label])\n    except KeyError as e:\n        # print('This is unlabeled data!')\n        pass\n\nnp.save('test-{}-{}-{}-label7.npy'.format(128,128,16), test_data)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"much_data = np.load('../input/chuck-file/muchdata-128-128-16-label7.npy',allow_pickle=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Normalization","metadata":{}},{"cell_type":"code","source":"len(much_data)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# `for i in range(len(much_data)):\n#     for j in range(len(much_data[i][0])):\n#         much_data[i][0][j] = much_data[i][0][j]/255`","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(len(much_data)):\n    if much_data[i][0].shape[0]!= 16:\n        print(much_data[i][0].shape)\n        print(i)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_data = much_data[:-300]\nval_data = much_data[-300:]\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_data[0][0].shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(len(train_data)):\n    if train_data[i][0].shape[0]!= 16:\n        print(train_data[i][0].shape)\n        print(i)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_data[1500][1][0].shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"DataGenerator for the batch size of 32","metadata":{}},{"cell_type":"code","source":"def DataGenerator(data,batch_size=32):\n    while True:\n        x, y= [],[] \n        for i in range(len(data)):\n            x.append(data[i][0])\n            y.append(data[i][1][0])\n            if len(y) == batch_size:                    \n                        yield np.array(x),np.array(y)\n                        x, y= [],[] ","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# model 1","metadata":{}},{"cell_type":"code","source":"def get_model():\n    \"\"\"Build a 3D convolutional neural network model.\"\"\"\n    inputs = keras.Input((None, None, None, 1))\n\n    x = layers.Conv3D(filters=64, kernel_size=3, activation=\"relu\",kernel_initializer=kernel_initializer, padding='same')(inputs)\n    x = layers.MaxPooling3D(pool_size=2)(x)\n    x = layers.BatchNormalization()(x)\n\n    x = layers.Conv3D(filters=64, kernel_size=3, activation=\"relu\",kernel_initializer=kernel_initializer, padding='same')(x)\n    x = layers.MaxPooling3D(pool_size=2)(x)\n    x = layers.BatchNormalization()(x)\n\n    x = layers.Conv3D(filters=128, kernel_size=3, activation=\"relu\",kernel_initializer=kernel_initializer, padding='same')(x)\n    x = layers.MaxPooling3D(pool_size=2)(x)\n    x = layers.BatchNormalization()(x)\n\n    x = layers.Conv3D(filters=256, kernel_size=3, activation=\"relu\",kernel_initializer=kernel_initializer, padding='same')(x)\n    x = layers.MaxPooling3D(pool_size=2)(x)\n    x = layers.BatchNormalization()(x)\n\n    x = layers.GlobalAveragePooling3D()(x)\n    x = layers.Dense(units=512, activation=\"relu\")(x)\n    x = layers.Dropout(0.3)(x)\n\n    outputs = layers.Dense(7, activation=\"sigmoid\")(x)\n    \n    model = keras.Model(inputs, outputs, name=\"3dcnn\")\n    return model","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = get_model()\nmodel.compile(optimizer = \"Adam\", loss=\"binary_crossentropy\", metrics=[\"accuracy\"])\nmodel.summary()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"checking if the model can overfit","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nhistory =model.fit(DataGenerator(train_data,32),validation_data=DataGenerator(val_data,32),epochs =2,callbacks = [keras.callbacks.EarlyStopping(monitor = 'loss', patience = 2, restore_best_weights = True)]\n,verbose = 1,\n validation_steps = len(val_data)//32 ,                        \nsteps_per_epoch = len(train_data)//32)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def TestGenerator(data,batch_size=32):\n    while True:\n        x= [] \n        for i in range(len(data)):\n            x.append(data[i][0])\n\n            if len(x) == batch_size:                    \n                        yield np.array(x)\n                        x = []","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_data= np.load('../input/chunk-test/test-128-128-16-label7.npy',allow_pickle=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(next(TestGenerator(test_data,3)))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(test_data)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prediction_type_mapping = test['prediction_type'].map({'C1': 0, 'C2': 1, 'C3': 2, 'C4': 3, 'C5': 4, 'C6': 5, 'C7': 6}).values","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preds = model.predict(TestGenerator(test_data,1), steps = len(test_data))\n                      \nnew_preds = []\nfor pred_idx in range(len(preds)):\n    new_preds.append(preds[pred_idx][prediction_type_mapping[pred_idx]])\n        # submission['fractured'] += preds[:, prediction_type_mapping] / 5\nsubmission['fractured'] += np.array(new_preds) / 5","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv('submission.csv', index = 0)","metadata":{},"execution_count":null,"outputs":[]}]}