{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceType":"competition","sourceId":13451,"datasetId":654585,"databundleVersionId":1188070}],"dockerImageVersionId":30648,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pydicom\nimport os\nimport matplotlib.pyplot as plt\nimport collections\nimport cv2\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom sklearn.model_selection import ShuffleSplit","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:30:56.775898Z","iopub.execute_input":"2024-02-25T18:30:56.776657Z","iopub.status.idle":"2024-02-25T18:30:56.781962Z","shell.execute_reply.started":"2024-02-25T18:30:56.776622Z","shell.execute_reply":"2024-02-25T18:30:56.780861Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Input paths\ninput_path = \"../input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection/\"\ntest_images_dir = os.path.join(input_path, 'stage_2_test/')\ntrain_images_dir = os.path.join(input_path, 'stage_2_train/')","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:30:56.896621Z","iopub.execute_input":"2024-02-25T18:30:56.896993Z","iopub.status.idle":"2024-02-25T18:30:56.901646Z","shell.execute_reply.started":"2024-02-25T18:30:56.896965Z","shell.execute_reply":"2024-02-25T18:30:56.900716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def correct_dcm(dcm):\n    \"\"\"\n    Correct DICOM pixel values.\n    \"\"\"\n    x = dcm.pixel_array + 1000\n    px_mode = 4096\n    x[x >= px_mode] -= px_mode\n    dcm.PixelData = x.tobytes()\n    dcm.RescaleIntercept = -1000\n\ndef window_image(dcm, window_center, window_width):\n    \"\"\"\n    Apply specified window level and width to the DICOM image.\n    \"\"\"\n    if (dcm.BitsStored == 12) and (dcm.PixelRepresentation == 0) and (int(dcm.RescaleIntercept) > -100):\n        correct_dcm(dcm)\n    \n    img = dcm.pixel_array * dcm.RescaleSlope + dcm.RescaleIntercept\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    img = np.clip(img, img_min, img_max)\n\n    return img\n\ndef bsb_window(dcm):\n    \"\"\"\n    Apply Brain, Subdural, and Soft tissue windowing to the DICOM image.\n    \"\"\"\n    brain_img = window_image(dcm, 40, 80)\n    subdural_img = window_image(dcm, 80, 200)\n    soft_img = window_image(dcm, 40, 380)\n    \n    brain_img = (brain_img - 0) / 80\n    subdural_img = (subdural_img - (-20)) / 200\n    soft_img = (soft_img - (-150)) / 380\n    bsb_img = np.array([brain_img, subdural_img, soft_img]).transpose(1, 2, 0)\n\n    return bsb_img\n\n# Load a sample DICOM image for visualization\ndicom = pydicom.dcmread(os.path.join(train_images_dir, 'ID_5c8b5d701.dcm'))\n\n# Display the DICOM image after applying Brain, Subdural, and Soft tissue windowing\nplt.imshow(bsb_window(dicom), cmap=plt.cm.bone)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:30:57.025541Z","iopub.execute_input":"2024-02-25T18:30:57.026352Z","iopub.status.idle":"2024-02-25T18:30:57.266495Z","shell.execute_reply.started":"2024-02-25T18:30:57.026319Z","shell.execute_reply":"2024-02-25T18:30:57.265589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def window_with_correction(dcm, window_center, window_width):\n    \"\"\"\n    Apply windowing to DICOM image with correction for specific conditions.\n    \"\"\"\n    if (dcm.BitsStored == 12) and (dcm.PixelRepresentation == 0) and (int(dcm.RescaleIntercept) > -100):\n        correct_dcm(dcm)\n    img = dcm.pixel_array * dcm.RescaleSlope + dcm.RescaleIntercept\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    img = np.clip(img, img_min, img_max)\n    return img\n\ndef window_without_correction(dcm, window_center, window_width):\n    \"\"\"\n    Apply windowing to DICOM image without correction.\n    \"\"\"\n    img = dcm.pixel_array * dcm.RescaleSlope + dcm.RescaleIntercept\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    img = np.clip(img, img_min, img_max)\n    return img\n\ndef window_testing(img, window):\n    \"\"\"\n    Apply Brain, Subdural, and Soft tissue windowing to the DICOM image.\n    \"\"\"\n    brain_img = window(img, 40, 80)\n    subdural_img = window(img, 80, 200)\n    soft_img = window(img, 40, 380)\n    \n    brain_img = (brain_img - 0) / 80\n    subdural_img = (subdural_img - (-20)) / 200\n    soft_img = (soft_img - (-150)) / 380\n    bsb_img = np.array([brain_img, subdural_img, soft_img]).transpose(1, 2, 0)\n\n    return bsb_img\n\n# Load an example DICOM image for visualization\ndicom = pydicom.dcmread(os.path.join(train_images_dir, \"ID_036db39b7.dcm\"))\n\n# Plot original and corrected images side by side\nfig, ax = plt.subplots(1, 2)\nax[0].imshow(window_testing(dicom, window_without_correction), cmap=plt.cm.bone)\nax[0].set_title(\"Original\")\nax[1].imshow(window_testing(dicom, window_with_correction), cmap=plt.cm.bone)\nax[1].set_title(\"Corrected\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:30:57.268028Z","iopub.execute_input":"2024-02-25T18:30:57.268282Z","iopub.status.idle":"2024-02-25T18:30:58.050644Z","shell.execute_reply.started":"2024-02-25T18:30:57.268259Z","shell.execute_reply":"2024-02-25T18:30:58.049731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def _read(path, desired_size):\n    \"\"\"\n    Read DICOM file, apply windowing, and resize the image.\n    This function will be used in DataGenerator.\n    \"\"\"\n    dcm = pydicom.dcmread(path)\n    \n    try:\n        img = bsb_window(dcm)\n    except Exception as e:\n        print(f\"An error occurred: {e}\")\n        img = np.zeros(desired_size)\n    \n    img = cv2.resize(img, desired_size[:2], interpolation=cv2.INTER_LINEAR)\n    \n    return img\n\n# Another sanity check \nplt.imshow(\n    _read(train_images_dir + 'ID_5c8b5d701' + '.dcm', (128, 128)), cmap=plt.cm.bone\n)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:31:32.70902Z","iopub.execute_input":"2024-02-25T18:31:32.709692Z","iopub.status.idle":"2024-02-25T18:31:32.926207Z","shell.execute_reply.started":"2024-02-25T18:31:32.709665Z","shell.execute_reply":"2024-02-25T18:31:32.925238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class DataGenerator(keras.utils.Sequence):\n    def __init__(self, list_IDs, labels=None, batch_size=1, img_size=(512, 512, 1), \n                 img_dir=train_images_dir, *args, **kwargs):\n        self.list_IDs = list_IDs\n        self.labels = labels\n        self.batch_size = batch_size\n        self.img_size = img_size\n        self.img_dir = img_dir\n        self.on_epoch_end()\n\n    def __len__(self):\n        return int(np.ceil(len(self.indices) / self.batch_size))\n\n    def __getitem__(self, index):\n        indices = self.indices[index * self.batch_size:(index + 1) * self.batch_size]\n        list_IDs_temp = [self.list_IDs[k] for k in indices]\n        \n        if self.labels is not None:\n            X, Y = self.__data_generation(list_IDs_temp)\n            return X, Y\n        else:\n            X = self.__data_generation(list_IDs_temp)\n            return X\n\n    def on_epoch_end(self):\n        if self.labels is not None:\n            keep_prob = self.labels.iloc[:, 0].map({0: 0.35, 1: 0.5})\n            keep = (keep_prob > np.random.rand(len(keep_prob)))\n            self.indices = np.arange(len(self.list_IDs))[keep]\n            np.random.shuffle(self.indices)\n        else:\n            self.indices = np.arange(len(self.list_IDs))\n\n    def __data_generation(self, list_IDs_temp):\n        X = np.empty((self.batch_size, *self.img_size))\n        \n        if self.labels is not None:\n            Y = np.empty((self.batch_size, 6), dtype=np.float32)\n        \n            for i, ID in enumerate(list_IDs_temp):\n                X[i,] = _read(os.path.join(self.img_dir, ID + \".dcm\"), self.img_size)\n                Y[i,] = self.labels.loc[ID].values\n        \n            return X, Y\n        \n        else:\n            for i, ID in enumerate(list_IDs_temp):\n                X[i,] = _read(os.path.join(self.img_dir, ID + \".dcm\"), self.img_size)\n            \n            return X","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:31:35.821848Z","iopub.execute_input":"2024-02-25T18:31:35.822583Z","iopub.status.idle":"2024-02-25T18:31:35.838071Z","shell.execute_reply.started":"2024-02-25T18:31:35.82254Z","shell.execute_reply":"2024-02-25T18:31:35.837053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from keras import backend as K\n\ndef weighted_log_loss(y_true, y_pred):\n    class_weights = np.array([2., 1., 1., 1., 1., 1.])\n    eps = K.epsilon()\n    y_pred = K.clip(y_pred, eps, 1.0 - eps)\n    out = -(y_true * K.log(y_pred) * class_weights + (1.0 - y_true) * K.log(1.0 - y_pred) * class_weights)\n    return K.mean(out, axis=-1)\n\ndef _normalized_weighted_average(arr, weights=None):\n    if weights is not None:\n        scl = K.sum(weights)\n        weights = K.expand_dims(weights, axis=1)\n        return K.sum(K.dot(arr, weights), axis=1) / scl\n    return K.mean(arr, axis=1)\n\ndef weighted_loss(y_true, y_pred):\n    class_weights = K.variable([2., 1., 1., 1., 1., 1.])\n    eps = K.epsilon()\n    y_pred = K.clip(y_pred, eps, 1.0 - eps)\n    loss = -(y_true * K.log(y_pred) + (1.0 - y_true) * K.log(1.0 - y_pred))\n    loss_samples = _normalized_weighted_average(loss, class_weights)\n    return K.mean(loss_samples)\n\ndef weighted_log_loss_metric(trues, preds):\n    class_weights = [2., 1., 1., 1., 1., 1.]\n    epsilon = 1e-7\n    preds = np.clip(preds, epsilon, 1 - epsilon)\n    loss = trues * np.log(preds) + (1 - trues) * np.log(1 - preds)\n    loss_samples = np.average(loss, axis=1, weights=class_weights)\n    return -loss_samples.mean()","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:31:37.777283Z","iopub.execute_input":"2024-02-25T18:31:37.777675Z","iopub.status.idle":"2024-02-25T18:31:37.788938Z","shell.execute_reply.started":"2024-02-25T18:31:37.777647Z","shell.execute_reply":"2024-02-25T18:31:37.788021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class PredictionCheckpoint(keras.callbacks.Callback):\n    def __init__(self, test_df, valid_df, test_images_dir=test_images_dir, valid_images_dir=train_images_dir, \n                 batch_size=32, input_size=(224, 224, 3)):\n        self.test_df = test_df\n        self.valid_df = valid_df\n        self.test_images_dir = test_images_dir\n        self.valid_images_dir = valid_images_dir\n        self.batch_size = batch_size\n        self.input_size = input_size\n\n    def on_train_begin(self, logs={}):\n        self.test_predictions = []\n        self.valid_predictions = []\n\n    def on_epoch_end(self, batch, logs={}):\n        self.test_predictions.append(\n            self.model.predict_generator(\n                DataGenerator(self.test_df.index, None, self.batch_size, self.input_size, self.test_images_dir), \n                verbose=2)[:len(self.test_df)]\n        )\n        # Uncomment this part if you want to compute validation predictions and losses\n        # self.valid_predictions.append(\n        #     self.model.predict_generator(\n        #         DataGenerator(self.valid_df.index, None, self.batch_size, self.input_size, self.valid_images_dir), \n        #         verbose=2)[:len(self.valid_df)]\n        # )\n        # print(\"validation loss: %.4f\" % weighted_log_loss_metric(self.valid_df.values, np.average(self.valid_predictions, axis=0, weights=[2**i for i in range(len(self.valid_predictions))])))\n\nclass MyDeepModel:\n    def __init__(self, engine, input_dims, batch_size=5, num_epochs=4, learning_rate=1e-3, \n                 decay_rate=1.0, decay_steps=1, weights=\"imagenet\", verbose=1):\n        self.engine = engine\n        self.input_dims = input_dims\n        self.batch_size = batch_size\n        self.num_epochs = num_epochs\n        self.learning_rate = learning_rate\n        self.decay_rate = decay_rate\n        self.decay_steps = decay_steps\n        self.weights = weights\n        self.verbose = verbose\n        self._build()\n\n    def _build(self):\n        engine = self.engine(include_top=False, weights=self.weights, input_shape=self.input_dims,\n                             backend=keras.backend, layers=keras.layers, models=keras.models, utils=keras.utils)\n        x = keras.layers.GlobalAveragePooling2D(name='avg_pool')(engine.output)\n        out = keras.layers.Dense(6, activation=\"sigmoid\", name='dense_output')(x)\n        self.model = keras.models.Model(inputs=engine.input, outputs=out)\n        self.model.compile(loss=\"binary_crossentropy\", optimizer=keras.optimizers.Adam(), metrics=[weighted_loss])\n\n    def fit_and_predict(self, train_df, valid_df, test_df):\n        pred_history = PredictionCheckpoint(test_df, valid_df, input_size=self.input_dims)\n        scheduler = keras.callbacks.LearningRateScheduler(lambda epoch: self.learning_rate * pow(self.decay_rate, floor(epoch / self.decay_steps)))\n        self.model.fit_generator(\n            DataGenerator(train_df.index, train_df, self.batch_size, self.input_dims, train_images_dir),\n            epochs=self.num_epochs,\n            verbose=self.verbose,\n            use_multiprocessing=True,\n            workers=4,\n            callbacks=[pred_history, scheduler]\n        )\n        return pred_history\n\n    def save(self, path):\n        self.model.save_weights(path)\n\n    def load(self, path):\n        self.model.load_weights(path)","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:31:42.042987Z","iopub.execute_input":"2024-02-25T18:31:42.043345Z","iopub.status.idle":"2024-02-25T18:31:42.060287Z","shell.execute_reply.started":"2024-02-25T18:31:42.043319Z","shell.execute_reply":"2024-02-25T18:31:42.05921Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_testset(filename=input_path + \"stage_2_sample_submission.csv\"):\n    df = pd.read_csv(filename)\n    df[\"Image\"] = df[\"ID\"].str.slice(stop=12)\n    df[\"Diagnosis\"] = df[\"ID\"].str.slice(start=13)\n    \n    df = df.loc[:, [\"Label\", \"Diagnosis\", \"Image\"]]\n    df = df.set_index(['Image', 'Diagnosis']).unstack(level=-1)\n    \n    return df\n\ndef read_trainset(filename=input_path + \"stage_2_train.csv\"):\n    df = pd.read_csv(filename)\n    df[\"Image\"] = df[\"ID\"].str.slice(stop=12)\n    df[\"Diagnosis\"] = df[\"ID\"].str.slice(start=13)\n    \n    duplicates_to_remove = [\n        56346, 56347, 56348, 56349,\n        56350, 56351, 1171830, 1171831,\n        1171832, 1171833, 1171834, 1171835,\n        3705312, 3705313, 3705314, 3705315,\n        3705316, 3705317, 3842478, 3842479,\n        3842480, 3842481, 3842482, 3842483\n    ]\n    \n    df = df.drop(index=duplicates_to_remove)\n    df = df.reset_index(drop=True)\n    \n    df = df.loc[:, [\"Label\", \"Diagnosis\", \"Image\"]]\n    df = df.set_index(['Image', 'Diagnosis']).unstack(level=-1)\n    \n    return df\n\ntest_df = read_testset()\ndf = read_trainset()","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:31:58.927806Z","iopub.execute_input":"2024-02-25T18:31:58.928193Z","iopub.status.idle":"2024-02-25T18:32:13.558607Z","shell.execute_reply.started":"2024-02-25T18:31:58.928166Z","shell.execute_reply":"2024-02-25T18:32:13.557588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head(3)","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:32:18.400199Z","iopub.execute_input":"2024-02-25T18:32:18.400559Z","iopub.status.idle":"2024-02-25T18:32:18.413925Z","shell.execute_reply.started":"2024-02-25T18:32:18.400533Z","shell.execute_reply":"2024-02-25T18:32:18.413012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_df.head(3)","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:32:20.71661Z","iopub.execute_input":"2024-02-25T18:32:20.716949Z","iopub.status.idle":"2024-02-25T18:32:20.732303Z","shell.execute_reply.started":"2024-02-25T18:32:20.716925Z","shell.execute_reply":"2024-02-25T18:32:20.731304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns\n\n# Plotting count distribution of labels\nplt.figure(figsize=(10, 6))\nsns.countplot(df.Label)\nplt.title('Distribution of Labels')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:32:23.282672Z","iopub.execute_input":"2024-02-25T18:32:23.283388Z","iopub.status.idle":"2024-02-25T18:32:23.964221Z","shell.execute_reply.started":"2024-02-25T18:32:23.283355Z","shell.execute_reply":"2024-02-25T18:32:23.963257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 1. Histogram of Class Distribution\nclass_distribution = df['Label'].sum(axis=0)\nplt.figure(figsize=(10, 6))\nplt.bar(class_distribution.index, class_distribution.values)\nplt.title('Class Distribution')\nplt.xlabel('Class')\nplt.ylabel('Count')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:32:26.558251Z","iopub.execute_input":"2024-02-25T18:32:26.559009Z","iopub.status.idle":"2024-02-25T18:32:26.754284Z","shell.execute_reply.started":"2024-02-25T18:32:26.558976Z","shell.execute_reply":"2024-02-25T18:32:26.753297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def window_image(dcm, window_center, window_width, intercept, slope):\n    \"\"\"\n    Apply specified window level and width to the DICOM image.\n    \"\"\"\n    img = dcm.pixel_array * slope + intercept\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    img = np.clip(img, img_min, img_max)\n    \n    return img\n\ndef bsb_window(dcm):\n    \"\"\"\n    Apply Brain, Subdural, and Soft tissue windowing to the DICOM image.\n    \"\"\"\n    brain_img = window_image(dcm, 40, 80, dcm.RescaleIntercept, dcm.RescaleSlope)\n    subdural_img = window_image(dcm, 80, 200, dcm.RescaleIntercept, dcm.RescaleSlope)\n    soft_img = window_image(dcm, 40, 380, dcm.RescaleIntercept, dcm.RescaleSlope)\n    \n    brain_img = (brain_img - 0) / 80\n    subdural_img = (subdural_img - (-20)) / 200\n    soft_img = (soft_img - (-150)) / 380\n    bsb_img = np.array([brain_img, subdural_img, soft_img]).transpose(1, 2, 0)\n    \n    return bsb_img\n\n# Now, let's plot the sample DICOM images with their windowed versions\nsample_ids = df.sample(n=3).index\nplt.figure(figsize=(15, 5))\nfor i, img_id in enumerate(sample_ids):\n    plt.subplot(3, 2, 2*i+1)\n    dcm = pydicom.dcmread(os.path.join(train_images_dir, f'{img_id}.dcm'))\n    plt.imshow(dcm.pixel_array, cmap=plt.cm.bone)\n    plt.title('Original')\n    plt.axis('off')\n    \n    plt.subplot(3, 2, 2*i+2)\n    plt.imshow(bsb_window(dcm), cmap=plt.cm.bone)\n    plt.title('Windowed')\n    plt.axis('off')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-02-25T18:32:29.53875Z","iopub.execute_input":"2024-02-25T18:32:29.539622Z","iopub.status.idle":"2024-02-25T18:32:30.149763Z","shell.execute_reply.started":"2024-02-25T18:32:29.539593Z","shell.execute_reply":"2024-02-25T18:32:30.148888Z"},"trusted":true},"execution_count":null,"outputs":[]}]}