{"cells":[{"metadata":{},"cell_type":"markdown","source":"IMPORTS"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport numpy as np \nimport pandas as pd \nimport os\nfrom pydicom import dcmread # for dcm files\nimport pydicom\nimport cv2\nimport csv\n\nfrom sklearn.metrics import confusion_matrix, accuracy_score\nfrom plotly.offline import iplot, init_notebook_mode\nfrom keras.losses import categorical_crossentropy\nfrom keras.optimizers import Adadelta\nimport plotly.graph_objs as go\nfrom matplotlib.pyplot import cm\nfrom keras.models import Model\nimport numpy as np\nimport keras\nimport h5py\n\n\nimport tensorflow as tf\nfrom skimage import measure, morphology\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# VISUALIZE SINGLE DCM FILE"},{"metadata":{},"cell_type":"markdown","source":"def visualize_single_dcm_file(path)"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"\ndef visualize_single_dcm_file(path):\n    \n    dataset = dcmread(path)\n    print(\"SOP name.....:\", dataset.SOPClassUID)\n    print(\"Series name.....:\", dataset.SeriesInstanceUID)\n\n    if 'PixelData' in dataset:\n        rows = int(dataset.Rows)\n        cols = int(dataset.Columns)\n        print(\"Image size.......: {rows:d} x {cols:d}, {size:d} bytes\".format(\n            rows=rows, cols=cols, size=len(dataset.PixelData)))\n        if 'PixelSpacing' in dataset:\n            print(\"Pixel spacing....:\", dataset.PixelSpacing)\n\n\n    print(\"Image position...:\", dataset.ImagePositionPatient )\n    cut_constant = 50\n    pixel_arrays_shape = (dataset.pixel_array).shape\n    crop = dataset.pixel_array[:,cut_constant:pixel_arrays_shape[1]-1-cut_constant]\n    crop = crop[cut_constant:pixel_arrays_shape[0]-1-cut_constant,:]\n    image = np.asarray(dataset.pixel_array)\n\n    print(image.shape)\n    #print(crop.shape)\n    inter_area_resize_256 = cv2.resize(dataset.pixel_array , (256, 256),  interpolation = cv2.INTER_AREA) \n    plt.imshow(image, cmap=plt.cm.bone)\n    plt.show()\n    plt.imshow(crop, cmap=plt.cm.bone)\n    plt.show()\n    plt.imshow(dataset.pixel_array, cmap=plt.cm.bone)\n    plt.show()\n    plt.imshow(inter_area_resize_256, cmap=plt.cm.bone)\n    plt.show()\n    inter_area_resize_64 = cv2.resize(crop , (64, 64),  interpolation = cv2.INTER_AREA)\n    plt.imshow(inter_area_resize_64, cmap=plt.cm.bone)\n    plt.show()\n    #return dataset.pixel_array\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Testing visualize_single_dcm_file(path) function"},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_single_dcm_file('../input/rsna-str-pulmonary-embolism-detection/train/0003b3d648eb/d2b2960c2bbf/00ac73cfc372.dcm')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# PATHS OF SERIES IN THE TRAINING SET"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_path = '../input/rsna-str-pulmonary-embolism-detection/train'\nstudies_in_train = os.listdir(train_path)\nprint(len(studies_in_train))\n\n\nstudies_in_train.sort()\n#print(train_path+'/'+studies_in_train[0]+'/'+ os.listdir(train_path+'/'+studies_in_train[0])[0])\nseries_paths = []\nfor x in range(len(studies_in_train)):\n    temp_series_path = train_path+'/'+studies_in_train[x]+'/'+ os.listdir(train_path+'/'+studies_in_train[x])[0]\n    series_paths.append(temp_series_path)\n\nprint(len(series_paths))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# PREPARE Y LABELS"},{"metadata":{"trusted":true},"cell_type":"code","source":"with open ('../input/rsna-str-pulmonary-embolism-detection/train.csv') as csvfile:  \n    readCSVFeatures = csv.reader(csvfile, delimiter=',')\n    y = list(csv.reader(csvfile))\n    \ny = np.asarray(y)\nstudy_instances = y[:,0]\nunique_study_instances, indexes_study_instances = np.unique(study_instances, return_index=True)\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nprint(unique_study_instances[0])\nunique_study_instances =np.asarray(unique_study_instances)\n\ntrain_y = y[indexes_study_instances]\nprint(len(train_y))\ntrain_y_cleaned = []\nfor i in range(0,len(studies_in_train)):\n    ind = np.where(unique_study_instances==studies_in_train[i])\n    train_y_cleaned.append(train_y[ind][0])\nprint(len(train_y_cleaned))    \ntrain_y_cleaned= np.asarray(train_y_cleaned)\ny_train = train_y_cleaned[:,3:len(train_y_cleaned[0])]\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preprocessing of a single serie\n "},{"metadata":{},"cell_type":"markdown","source":"* Load Scans in given folder path\n* Get pixels of slices \n* Plot 3d \nReference : https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load the scans in given folder path\ndef load_scan(path):\n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices\n\n\n\ndef get_pixels_hu(slices):\n    image = np.stack([s.pixel_array for s in slices])\n    image_64 = np.ones((len(slices),64,64))\n    # Convert to int16 (from sometimes int16), \n    # should be possible as values should always be low enough (<32k)\n    image = image.astype(np.int16)\n\n    # Set outside-of-scan pixels to 0\n    # The intercept is usually -1024, so air is approximately 0\n    image[image == -2000] = 0\n    \n    # Convert to Hounsfield units (HU)\n    for slice_number in range(len(slices)):\n        \n        intercept = slices[slice_number].RescaleIntercept\n        slope = slices[slice_number].RescaleSlope\n        \n        if slope != 1:\n            image[slice_number] = slope * image[slice_number].astype(np.float64)\n            image[slice_number] = image[slice_number].astype(np.int16)\n            \n        image[slice_number] += np.int16(intercept)\n        cut_constant = 50\n        pixel_arrays_shape = (image[slice_number]).shape\n        temp = image[slice_number][:,cut_constant:pixel_arrays_shape[1]-1-cut_constant]\n        temp = temp[cut_constant:pixel_arrays_shape[0]-1-cut_constant,:]\n        image_64[slice_number] = cv2.resize(temp , (64, 64),  interpolation = cv2.INTER_AREA)\n        #image[slice_number] = inter_area_resize_64\n        \n        \n    \n\n    #return np.array(image, dtype=np.int16),np.array(image_64, dtype=np.int16)\n    return np.array(image_64, dtype=np.int16)\n\n\n\n\ndef plot_3d(image, threshold=-300):\n    \n    # Position the scan upright, \n    # so the head of the patient would be at the top facing the camera\n    p = image.transpose(2,1,0)\n    \n    verts, faces, norm, val = measure.marching_cubes_lewiner(p, threshold)\n\n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n\n    # Fancy indexing: `verts[faces]` to generate a collection of triangles\n    mesh = Poly3DCollection(verts[faces], alpha=0.70)\n    face_color = [0.45, 0.45, 0.75]\n    mesh.set_facecolor(face_color)\n    ax.add_collection3d(mesh)\n\n    ax.set_xlim(0, p.shape[0])\n    ax.set_ylim(0, p.shape[1])\n    ax.set_zlim(0, p.shape[2])\n\n    plt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# TRY \n* def load_scan\n* def get_pixels_hu\n* plot_3d(first_patient_pixels, 400)    "},{"metadata":{"trusted":true},"cell_type":"code","source":"slices = load_scan('../input/rsna-str-pulmonary-embolism-detection/train/005a0dbcb4b7/4ceaee66edc8')\nfirst_patient_pixels= get_pixels_hu(slices)\nplt.hist(first_patient_pixels.flatten(), bins=80, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()\n\n# Show some slice in the middle\nplt.imshow(first_patient_pixels[80], cmap=plt.cm.gray)\nplt.show()\n\n#plot_3d(first_patient_pixels, 400)    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# RESAMPLING "},{"metadata":{"trusted":true},"cell_type":"code","source":"def resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    #return np.array(image, dtype=np.int16), np.array([slices[0].SliceThickness, slices[0].PixelSpacing[0], slices[0].PixelSpacing[1]], dtype=np.float32)\n    return np.array(image, dtype=np.int16), np.array([1, 1, 1], dtype=np.float32)\n\n#pix_resampled,spacing = resample(first_patient_pixels,slices, [1,1,1])\n#plot_3d(pix_resampled, 400)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# PLOT 3D image of a given series path"},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_scans_3d(path): \n    number_of_slices_to_be_used = 300\n    slices = load_scan(path)\n    first_patient_pixels = get_pixels_hu(slices)\n    plt.hist(first_patient_pixels.flatten(), bins=80, color='c')\n    plt.xlabel(\"Hounsfield Units (HU)\")\n    plt.ylabel(\"Frequency\")\n    plt.show()\n    pix_resampled,spacing = resample(first_patient_pixels,slices, [1,1,1])\n    # Show some slice in the middle\n    if pix_resampled.shape[0] > number_of_slices_to_be_used:\n        pix_resampled_slice_editted = pix_resampled[len(pix_resampled)-number_of_slices_to_be_used:len(pix_resampled)]\n    elif pix_resampled.shape[0] < number_of_slices_to_be_used:       \n\n        dummy_slice = -1024*np.ones((number_of_slices_to_be_used-pix_resampled.shape[0],512, 512), dtype=int)\n        pix_resampled_slice_editted = np.concatenate((dummy_slice, pix_resampled), axis=0)\n            \n    \n    \n    plt.imshow(first_patient_pixels[10], cmap=plt.cm.gray)\n    plt.show()\n    plot_3d(first_patient_pixels, 400) \n    plot_3d(pix_resampled, 400) \n    plot_3d(pix_resampled_slice_editted, 400) \n    \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_scans_3d(series_paths[6222])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Normalization**\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"\nMIN_BOUND = -1000.0\nMAX_BOUND = 400.0\n    \ndef normalize(image):\n    image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)\n    image[image>1] = 1.\n    image[image<0] = 0.\n    return image\n\nPIXEL_MEAN = 0.25\ndef zero_center(image):\n    image = image - PIXEL_MEAN\n    return image\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# PREPROCESS CT SCANS OF A GIVEN PATH OF A SERIES"},{"metadata":{"trusted":true},"cell_type":"code","source":"def preprocess(path):\n    number_of_slices_to_be_used = 50\n    slices = load_scan(path)\n    pixels_of_serie = get_pixels_hu(slices)\n    pixels_of_serie_resampled,spacing = resample(pixels_of_serie,slices, [1,1,1])\n    if pixels_of_serie_resampled.shape[0] > number_of_slices_to_be_used:\n        pixels_of_serie_resampled = pixels_of_serie_resampled[len(pixels_of_serie_resampled)-number_of_slices_to_be_used:len(pixels_of_serie_resampled)]\n    elif pixels_of_serie_resampled.shape[0] < number_of_slices_to_be_used:       \n\n        dummy_slice = -1024*np.ones((number_of_slices_to_be_used-pixels_of_serie_resampled.shape[0],64, 64), dtype=int)\n        pixels_of_serie_resampled = np.concatenate((dummy_slice, pixels_of_serie_resampled), axis=0)\n                \n    preprocessed_image = normalize(pixels_of_serie_resampled)\n    preprocessed_image = zero_center(preprocessed_image)\n    preprocessed_image = preprocessed_image.reshape( number_of_slices_to_be_used, 64, 64, 1)\n    #preprocessed_image = preprocessed_image.reshape( 1,number_of_slices_to_be_used, 64, 64, 1)\n\n    return preprocessed_image\n    \n    \n    \n    \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# TEST \ndef preprocess(path)"},{"metadata":{"trusted":true},"cell_type":"code","source":"number_of_slices_to_be_used = 50\nx_train = np.ones((4,number_of_slices_to_be_used, 64, 64, 1),  dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(x_train[0].shape)\ndene = preprocess(series_paths[0])\nprint(dene.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\nfor i in range(0,4):\n    try :         \n        x_train[i]= preprocess(series_paths[i])\n    except:\n        print(i)\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(y_train.shape)\nprint(x_train.shape)\n\nprint(x_train.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#delete_indexes = [5,21,25,48,99,107,115,153,185,193]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#x_train  = np.delete(x_train , delete_indexes)\ny_train = y_train[0:4]\n#y_train = np.delete(y_train , delete_indexes)\n\nprint(y_train.shape)\nprint(x_train.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# The inputs are 28x28x28 volumes with a single channel, and the\n# batch size is 4\nnumber_of_slices_to_be_used = 50\n#input_shape =( 240, 64, 64, 1)\ninput_layer = tf.keras.layers.Input(( number_of_slices_to_be_used, 64, 64, 1))\n\nconv_layer1 = tf.keras.layers.Conv3D(filters=8, kernel_size=(3, 3, 3), activation='relu')(input_layer)\nconv_layer2 = tf.keras.layers.Conv3D(filters=16, kernel_size=(3, 3, 3), activation='relu')(conv_layer1)\n\npooling_layer1 = tf.keras.layers.MaxPool3D(pool_size=(2, 2, 2))(conv_layer2)\n\nconv_layer3 = tf.keras.layers.Conv3D(filters=32, kernel_size=(3, 3, 3), activation='relu')(pooling_layer1)\nconv_layer4 = tf.keras.layers.Conv3D(filters=64, kernel_size=(3, 3, 3), activation='relu')(conv_layer3)\npooling_layer2 = tf.keras.layers.MaxPool3D(pool_size=(2, 2, 2))(conv_layer4)\n\n## perform batch normalization on the convolution outputs before feeding it to MLP architecture\npooling_layer2 = tf.keras.layers.BatchNormalization()(pooling_layer2)\nflatten_layer = tf.keras.layers.Flatten()(pooling_layer2)\n\nprint(flatten_layer.shape)\n\ndense_layer1 = tf.keras.layers.Dense(units=2048, activation='relu')(flatten_layer)\ndense_layer1 = tf.keras.layers.Dropout(0.4)(dense_layer1)\ndense_layer2 = tf.keras.layers.Dense(units=512, activation='relu')(dense_layer1)\ndense_layer2 = tf.keras.layers.Dropout(0.4)(dense_layer2)\noutput_layer = tf.keras.layers.Dense(units=14, activation='softmax')(dense_layer2)\nprint(output_layer.shape)\n## define the model with input layer and output layer\nmodel = Model(inputs=input_layer, outputs=output_layer)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ny_train = y_train.astype(np.float)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.compile(loss=categorical_crossentropy, optimizer=Adadelta(lr=0.1), metrics=['acc'])\nmodel.fit(x=x_train, y=y_train, batch_size=2, epochs=50, validation_split=0.2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# CREATE 3D CONVOLUTIONAL NEURAL NETWORK"}],"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":4,"nbformat_minor":4}