{"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":"<div class='alert alert-info' style='text-align: center'><h1>DICOM Full Range Pixels as CNN Input</h1>\n    - yet another chest x-ray processing notebook -\n</div>\n\n**This notebook uses the full range of DICOM pixel values as input for a CNN and a Sequential model, rather than exporting 8 bit JPGs.**\n\n- It does not address MONOCHROME1 (inverted pixel intensities). Obviously this would have to be added.\n- It resizes the images to a managable size (256,256) by default. Which is lossy of course.\n- It ignores JPG compressed Transfer Syntaxes and only gets images that are in Explicit VR LE (to avoid pylibjpg install).\n- It only uses a few images and a few epochs to demonstrate the thought process. Clearly, it won't be accurate at all.\n- I didn't bother normalizing or otherwise apply any kind of processing to the data.\n- I'm sure the pros already do this, but I'm still trying to wrap my brain around it.\n\n**It seems that adding a VOI LUT, or bit plane slicing on the full range images would produce more usable results for the model.**","metadata":{}},{"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore')\nimport os\nimport cv2\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom skimage.transform import resize\nfrom sklearn.model_selection import train_test_split\nfrom skimage.color import gray2rgb\nfrom tensorflow.keras import Model\nfrom tensorflow.keras.optimizers import Adam, RMSprop\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import Dense, Dropout, BatchNormalization, Conv2D, MaxPooling2D, GlobalMaxPooling2D, Activation, Flatten\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator\nfrom tensorflow.keras.callbacks import ModelCheckpoint, EarlyStopping\nfrom tensorflow.keras.applications.inception_v3 import InceptionV3","metadata":{"execution":{"iopub.status.busy":"2021-07-06T04:17:47.161975Z","iopub.execute_input":"2021-07-06T04:17:47.162725Z","iopub.status.idle":"2021-07-06T04:17:55.185014Z","shell.execute_reply.started":"2021-07-06T04:17:47.162603Z","shell.execute_reply":"2021-07-06T04:17:55.184101Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Path to the SIIM Covid19 dataset\nbase_path = \"/kaggle/input/siim-covid19-detection/\"\n\n# Final image size for model input\nimg_size = 256\n\n# number of images in each of the two classes to get into the train set. Keep it real small for this demo.\nnum_images = 100\n\n# Number of epochs\nepochs = 10\n\n# number of batches\nbatch_size = 10","metadata":{"execution":{"iopub.status.busy":"2021-07-06T04:17:59.301281Z","iopub.execute_input":"2021-07-06T04:17:59.301667Z","iopub.status.idle":"2021-07-06T04:17:59.308214Z","shell.execute_reply.started":"2021-07-06T04:17:59.301632Z","shell.execute_reply":"2021-07-06T04:17:59.30711Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the data\nstudies_df = pd.read_csv(os.path.join(base_path,\"train_study_level.csv\"))\nimages_df = pd.read_csv(os.path.join(base_path,\"train_image_level.csv\"))","metadata":{"execution":{"iopub.status.busy":"2021-07-06T04:18:02.543181Z","iopub.execute_input":"2021-07-06T04:18:02.54358Z","iopub.status.idle":"2021-07-06T04:18:02.625244Z","shell.execute_reply.started":"2021-07-06T04:18:02.543548Z","shell.execute_reply":"2021-07-06T04:18:02.624486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"images_df.head(10)","metadata":{"execution":{"iopub.status.busy":"2021-07-06T04:24:08.857814Z","iopub.execute_input":"2021-07-06T04:24:08.858225Z","iopub.status.idle":"2021-07-06T04:24:08.872671Z","shell.execute_reply.started":"2021-07-06T04:24:08.858193Z","shell.execute_reply":"2021-07-06T04:24:08.87158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Loop through the dfs and remove the \"_study\" and \"_image\" from the ids so we can use them as keys to join on later\nstudies_df['id'] = studies_df['id'].map(lambda x: x.rstrip('_study'))\nimages_df['id'] = images_df['id'].map(lambda x: x.rstrip('_image'))","metadata":{"execution":{"iopub.status.busy":"2021-06-25T15:15:34.321862Z","iopub.execute_input":"2021-06-25T15:15:34.32236Z","iopub.status.idle":"2021-06-25T15:15:34.344956Z","shell.execute_reply.started":"2021-06-25T15:15:34.32231Z","shell.execute_reply":"2021-06-25T15:15:34.343876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Merge the study and image df's on StudyInstanceUID\ndata_df = pd.merge(images_df, studies_df, how='inner', left_on='StudyInstanceUID', right_on='id')\ndata_df.drop(['id_y'], axis=1, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2021-06-25T15:15:37.267421Z","iopub.execute_input":"2021-06-25T15:15:37.267852Z","iopub.status.idle":"2021-06-25T15:15:37.295996Z","shell.execute_reply.started":"2021-06-25T15:15:37.267819Z","shell.execute_reply":"2021-06-25T15:15:37.295014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_df.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-25T15:19:10.522183Z","iopub.execute_input":"2021-06-25T15:19:10.522701Z","iopub.status.idle":"2021-06-25T15:19:10.539278Z","shell.execute_reply.started":"2021-06-25T15:19:10.522655Z","shell.execute_reply":"2021-06-25T15:19:10.538138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"new_df = data_df[data_df['Atypical Appearance'] == 1]\n\nnew_df2 = new_df[new_df['Indeterminate Appearance'] == 1]\nnew_df2.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-25T15:22:36.764685Z","iopub.execute_input":"2021-06-25T15:22:36.76551Z","iopub.status.idle":"2021-06-25T15:22:36.779887Z","shell.execute_reply.started":"2021-06-25T15:22:36.765467Z","shell.execute_reply":"2021-06-25T15:22:36.779118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Split the studies dataframe into negative and positive dataframes\npositive_images_df = data_df[data_df['Negative for Pneumonia'] == 0]\nnegative_images_df = data_df[data_df['Negative for Pneumonia'] == 1]","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:42:29.95939Z","iopub.execute_input":"2021-06-01T00:42:29.959823Z","iopub.status.idle":"2021-06-01T00:42:29.970421Z","shell.execute_reply.started":"2021-06-01T00:42:29.959779Z","shell.execute_reply":"2021-06-01T00:42:29.969013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Make some helper functions","metadata":{}},{"cell_type":"code","source":"# This function finds the first image in a StudyInstanceUID directory and returns its path\ndef get_image_by_study_id(study_id):\n    study_path = base_path + \"train/\" + study_id + \"/\"\n    for subdir, dirs, files in os.walk(study_path):\n        for file in files:     \n            image = os.path.join(subdir, file)\n            if os.path.isfile(image):\n                return image\n    return \"none\"","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:42:29.972558Z","iopub.execute_input":"2021-06-01T00:42:29.973342Z","iopub.status.idle":"2021-06-01T00:42:29.9802Z","shell.execute_reply.started":"2021-06-01T00:42:29.973292Z","shell.execute_reply":"2021-06-01T00:42:29.978876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to resize pixels\ndef resize_image(pixels_in): \n    return resize(pixels_in, (img_size, img_size), anti_aliasing=True).astype(float)","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:42:29.982023Z","iopub.execute_input":"2021-06-01T00:42:29.982913Z","iopub.status.idle":"2021-06-01T00:42:29.992402Z","shell.execute_reply.started":"2021-06-01T00:42:29.982856Z","shell.execute_reply":"2021-06-01T00:42:29.990808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Visualize an image","metadata":{}},{"cell_type":"code","source":"# Take a look at a random image\nimg_file = get_image_by_study_id(\"00c74279c5b7\")\nimg = pydicom.dcmread(img_file)\npixels = img.pixel_array\n\nprint(\"Pixel range: \" + str(np.amin(pixels)) + \" - \" + str(np.amax(pixels)))\nplt.imshow(pixels,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:42:29.996321Z","iopub.execute_input":"2021-06-01T00:42:29.996843Z","iopub.status.idle":"2021-06-01T00:42:30.945106Z","shell.execute_reply.started":"2021-06-01T00:42:29.996796Z","shell.execute_reply":"2021-06-01T00:42:30.943935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Get the pixels","metadata":{}},{"cell_type":"code","source":"# Iterate through images in the positive set and extract pixels. This takes a while with a large set.\n# Get only LE Explicit DICOMs and ignore other transfer syntaxes so we don't have to deal with pylibjpg\n\nX_train_data = []\ny_train_data = []\n\n# Iterate through the rows of the 'positive' DF\ncount = 0\nfor index, row in positive_images_df.iterrows():\n    img_file = get_image_by_study_id(row['StudyInstanceUID'])\n    img = pydicom.dcmread(img_file)\n    \n    # Get only Explicit VR LE Transfer Syntax\n    if img.file_meta.TransferSyntaxUID == \"1.2.840.10008.1.2.1\":\n        pixels = resize_image(img.pixel_array)\n        X_train_data.append(pixels)\n        y_train_data.append(row['Negative for Pneumonia'])\n        if (count == num_images):\n            break\n        count += 1\n    \nprint(\"Done getting \" + str(count) + \" negative images\")","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:42:30.947193Z","iopub.execute_input":"2021-06-01T00:42:30.947788Z","iopub.status.idle":"2021-06-01T00:43:39.381424Z","shell.execute_reply.started":"2021-06-01T00:42:30.947693Z","shell.execute_reply":"2021-06-01T00:43:39.380109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extract pixels from images in the negative set\ncount = 0\nfor index, row in negative_images_df.iterrows():\n    img_file = get_image_by_study_id(row['StudyInstanceUID'])\n    img = pydicom.dcmread(img_file)\n    if img.file_meta.TransferSyntaxUID == \"1.2.840.10008.1.2.1\":\n        pixels = resize_image(img.pixel_array)\n        X_train_data.append(pixels)\n        y_train_data.append(row['Negative for Pneumonia'])\n        if (count == num_images):\n            break\n        count += 1 \nprint(\"Done getting \" + str(count) + \" positive images\")","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:43:39.388324Z","iopub.execute_input":"2021-06-01T00:43:39.392393Z","iopub.status.idle":"2021-06-01T00:44:38.94583Z","shell.execute_reply.started":"2021-06-01T00:43:39.392314Z","shell.execute_reply":"2021-06-01T00:44:38.944646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Split the train/test data\nX_train, X_test, y_train, y_test = train_test_split(X_train_data, y_train_data, test_size = 0.3, random_state = 82, shuffle=True)\nprint(\"X_train len: \" + str(len(X_train)))\nprint(\"y_train len: \" + str(len(y_train)))\nprint(\"X_test len: \" + str(len(X_test)))\nprint(\"y_test len: \" + str(len(y_test)))","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:44:38.947433Z","iopub.execute_input":"2021-06-01T00:44:38.948181Z","iopub.status.idle":"2021-06-01T00:44:38.961847Z","shell.execute_reply.started":"2021-06-01T00:44:38.948129Z","shell.execute_reply":"2021-06-01T00:44:38.960583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Convert the lists to arrays\nX_train = np.asarray(X_train)\nX_test = np.asarray(X_test)\ny_train = np.asarray(y_train)\ny_test = np.asarray(y_test)\n\n# Reshape the images to 3 channels\nX_train = gray2rgb(X_train)\nX_test = gray2rgb(X_test)","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:44:38.963735Z","iopub.execute_input":"2021-06-01T00:44:38.965341Z","iopub.status.idle":"2021-06-01T00:44:39.189112Z","shell.execute_reply.started":"2021-06-01T00:44:38.965283Z","shell.execute_reply":"2021-06-01T00:44:39.188087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Build the models","metadata":{}},{"cell_type":"code","source":"# Use InceptionV3 for transfer learning and add our data as the last layer.\npre_trained_model = InceptionV3(input_shape = (img_size, img_size, 3), include_top = False, weights = \"imagenet\")\n\n# Freeze the lower layers\nfor layer in pre_trained_model.layers:\n     layer.trainable = False\n        \nlast_layer = pre_trained_model.get_layer('mixed7')\nlast_output = last_layer.output\n\n# Add our layer\nlayer = Flatten()(last_output)\nlayer = Dense(1024, activation='relu')(layer)\nlayer = Dropout(0.2)(layer)                  \nlayer = Dense(1, activation='sigmoid')(layer)           \n\nmodel = Model(pre_trained_model.input, layer) \nmodel.compile(optimizer = RMSprop(lr=0.0001), loss = 'binary_crossentropy', metrics = ['acc'])\n\n# Fit the data\nhistory=model.fit(X_train,y_train,epochs=epochs,verbose=1,validation_data=(X_test,y_test))","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:44:39.19054Z","iopub.execute_input":"2021-06-01T00:44:39.190927Z","iopub.status.idle":"2021-06-01T00:44:52.792516Z","shell.execute_reply.started":"2021-06-01T00:44:39.190876Z","shell.execute_reply":"2021-06-01T00:44:52.791546Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot training history\nplt.plot(history.history['loss'], label='train')\nplt.plot(history.history['val_loss'], label='test')\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:44:52.794797Z","iopub.execute_input":"2021-06-01T00:44:52.795237Z","iopub.status.idle":"2021-06-01T00:44:52.968491Z","shell.execute_reply.started":"2021-06-01T00:44:52.795193Z","shell.execute_reply":"2021-06-01T00:44:52.967248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Build a simple Sequential model with a couple layers","metadata":{}},{"cell_type":"code","source":"model2 = Sequential()\nmodel2.add(Conv2D(32, (3, 3), activation=\"relu\", input_shape = X_train.shape[1:]))\nmodel2.add(BatchNormalization())\nmodel2.add(MaxPooling2D((2, 2)))\nmodel2.add(Dropout(0.2))\n\nmodel2.add(Conv2D(64, (3, 3), activation=\"relu\"))\nmodel2.add(BatchNormalization())\nmodel2.add(MaxPooling2D((2, 2)))\nmodel2.add(Dropout(0.3))\n\nmodel2.add(GlobalMaxPooling2D())\nmodel2.add(Dense(256, activation=\"relu\"))\nmodel2.add(Dropout(0.5))\nmodel2.add(Dense(1, activation=\"sigmoid\"))\n\nopt = Adam(learning_rate=0.001)\nmodel2.compile(loss='binary_crossentropy',optimizer=opt,metrics=['accuracy'])\nearly_stopping = EarlyStopping(monitor='val_loss', min_delta=0.001, patience=10)\n\nhistory = model2.fit(X_train,y_train,batch_size=batch_size,epochs=epochs,verbose=1,validation_data=(X_test, y_test),callbacks=[early_stopping])","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:44:52.970092Z","iopub.execute_input":"2021-06-01T00:44:52.970689Z","iopub.status.idle":"2021-06-01T00:45:02.862952Z","shell.execute_reply.started":"2021-06-01T00:44:52.970641Z","shell.execute_reply":"2021-06-01T00:45:02.861863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot training history\nplt.plot(history.history['loss'], label='train')\nplt.plot(history.history['val_loss'], label='test')\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-01T00:45:02.865386Z","iopub.execute_input":"2021-06-01T00:45:02.865866Z","iopub.status.idle":"2021-06-01T00:45:03.046499Z","shell.execute_reply.started":"2021-06-01T00:45:02.865816Z","shell.execute_reply":"2021-06-01T00:45:03.045375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Here are some other processing notebooks I made:\n- Lung Segmentation Without CNN -> https://www.kaggle.com/davidbroberts/lung-segmentation-without-cnn\n- Applying filters to x-rays -> https://www.kaggle.com/davidbroberts/applying-filters-to-chest-x-rays\n- Rib supression on Chest X-Rays -> https://www.kaggle.com/davidbroberts/rib-suppression-poc\n- Manual DICOM VOI LUT -> https://www.kaggle.com/davidbroberts/manual-dicom-voi-lut\n- Apply Unsharp Mask to Chest X-Rays -> https://www.kaggle.com/davidbroberts/unsharp-masking-chest-x-rays\n- Cropping Chest X-Rays -> https://www.kaggle.com/davidbroberts/cropping-chest-x-rays\n- Bounding Boxes on Cropped Images -> https://www.kaggle.com/davidbroberts/bounding-boxes-on-cropped-images\n- Visualizing Chest X-Ray bit planes -> https://www.kaggle.com/davidbroberts/visualizing-chest-x-ray-bitplanes\n- Standardizing Chest X-Ray Dataset Exports -> https://www.kaggle.com/davidbroberts/standardizing-cxr-datasets","metadata":{}}]}