{"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 cv2\nimport matplotlib\nimport numpy as np\nimport os\nimport pandas as pd\nimport pickle\nimport pydicom as dicom\nimport tensorflow as tf\nimport time\n\nfrom keras import Sequential\nfrom keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau\nfrom keras.models import load_model\nfrom keras.optimizers import SGD\nfrom matplotlib import pyplot as plt\nfrom PIL import Image\nfrom scipy import ndimage\nfrom skimage import measure, morphology\nfrom sklearn.cluster import KMeans\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.utils import shuffle\nfrom statsmodels.formula.api import quantreg\nfrom tensorflow import keras\nfrom tensorflow.keras import layers","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# PART 1: Quantile Regression model\n# Model obtained from:\n# https://www.kaggle.com/titericz/tabular-simple-eda-linear-model \nos.chdir('/kaggle/input/osic-pulmonary-fibrosis-progression')\n\ntrain = pd.read_csv('train.csv')\ntest  = pd.read_csv('test.csv')\n\n# Marker to mark rows as from train or test csv\ntrain['traintest'] = 0\ntest ['traintest'] = 1\n\n# Generate format for submission.csv\nsubmission = pd.read_csv('sample_submission.csv')\nsubmission['Weeks'] = submission['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsubmission['Patient'] = submission['Patient_Week'].apply(lambda x: x.split('_')[0]) \n\ntrain = pd.concat((train, test))\ntrain.sort_values(['Patient','Weeks'], inplace=True)\n\n# Encode string value into categorical value\ntrain['Sex'] = pd.factorize(train['Sex'])[0]\ntrain['SmokingStatus'] = pd.factorize(train['SmokingStatus'])[0]\n\n# Standardize categorical values\ntrain['Percent']       = (train['Percent'] - train['Percent'].mean()) / train['Percent'].std()\ntrain['Age']           = (train['Age'] - train['Age'].mean()) / train['Age'].std()\ntrain['Sex']           = (train['Sex'] - train['Sex'].mean()) / train['Sex'].std()\ntrain['SmokingStatus'] = (train['SmokingStatus'] - train['SmokingStatus'].mean()) / train['SmokingStatus'].std()\n\n# Train model\nmodelL = quantreg('FVC ~ Weeks+Percent+Age+Sex+SmokingStatus', train).fit(q=0.15)\nmodel  = quantreg('FVC ~ Weeks+Percent+Age+Sex+SmokingStatus', train).fit(q=0.50)\nmodelH = quantreg('FVC ~ Weeks+Percent+Age+Sex+SmokingStatus', train).fit(q=0.85)\n\n# This part is used to evaluate model on train set\ntrain['ypredL'] = modelL.predict(train).values\ntrain['ypred']  = model.predict(train).values\ntrain['ypredH'] = modelH.predict(train).values\ntrain['ypredstd'] = 0.5 * np.abs(train['ypredH'] - train['ypred']) + \\\n                    0.5 * np.abs(train['ypred'] - train['ypredL'])\n\n# Prepare test dataframe for predicting FVC values\ndt = train.loc[train.traintest==1 ,['Patient','Percent','Age','Sex','SmokingStatus']]\ntest = pd.merge(submission, dt, on='Patient', how='left')\ntest.sort_values(['Patient','Weeks'], inplace=True)\n\n# Predict FVC values for test dataframe\ntest['ypredL'] = modelL.predict(test).values\ntest['FVC']    = model.predict(test).values\ntest['ypredH'] = modelH.predict(test).values\ntest['Confidence'] = np.abs(test['ypredH'] - test['ypredL']) / 2\n\n# Predict FVC from each Patient at Week = 0\n# Use FVC at Week = 0 to train CNN\ncnn_train = train.groupby('Patient').first()\ncnn_train['Week_0'] = 0\ncnn_train.rename(columns={'Weeks': 'Week_i', 'Week_0': 'Weeks'}, inplace=True)\n\ncnn_train['FVC_0'] = model.predict(cnn_train).values","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Preprocessing funtions\ndef processImage(image, slope, intercept):\n    # Crop image if not square\n    def crop(image):\n        if image.shape[0] != image.shape[1]:\n            height = image.shape[0]\n            width = image.shape[1]\n            \n            height_crop = (height - 512) // 2\n            width_crop = (width - 512) // 2\n            \n            image = image[height_crop:-height_crop, width_crop:-width_crop]\n        return image\n    \n    # Transform image values into Hounsfield Units with formula given in\n    # https://www.kaggle.com/avirdee/understanding-dicoms\n    def rescale(image):\n        rescaled = image * slope + intercept\n        return rescaled\n\n    # Normalize image values to between 0 and 1\n    def normalize(image):\n        min_value = -1000\n        max_value = 400\n\n        image[image < min_value] = min_value\n        image[image > max_value] = max_value\n\n        image = (image - min_value) / (max_value - min_value)\n        image = image.astype(\"float32\")\n\n        return image\n    \n    # Lung masking/segmentation\n    # https://www.kaggle.com/andradaolteanu/pulmonary-fibrosis-competition-eda-dicom-prep\n    def mask(image):\n        rows = image.shape[0]\n        cols = image.shape[1]\n\n        # Using K-Means to segment image\n        flat_image = np.reshape(image, [np.prod(image.shape), 1])\n        kmeans = KMeans(n_clusters=2).fit(flat_image)\n        centers = sorted(kmeans.cluster_centers_.flatten())\n        threshold = np.mean(centers)\n        segmented = np.where(image < threshold, 1.0, 0.0)\n\n        final_masks = []\n\n        # Access each segmented region individually\n        labels = measure.label(segmented) \n        regions = measure.regionprops(labels)\n        for region in regions:\n            B = region.bbox\n            # Hard-coded to ignore \"borders\" of the image\n            if (B[2] - B[0] < rows / 10 * 9 and \n                B[3] - B[1] < cols / 10 * 9 and \n                B[0] > rows / 10 and \n                B[2] < cols / 10 * 9):\n\n                final_masks.append(region.label)\n\n        # Final mask to segment out the lungs\n        mask = np.ndarray([rows, cols], dtype=np.int8)\n        mask[:] = 0\n        for N in final_masks:\n            mask = mask + np.where(labels == N, 1, 0)\n\n        return mask\n    \n    def resize(image):\n        resize_shape = (128, 128)\n        curr_shape = image.shape\n        \n        resized = ndimage.zoom(image, resize_shape[0] / curr_shape[0])\n        \n        return resized\n    \n    cropped = crop(image)\n    rescaled = rescale(cropped)\n    normalized = normalize(rescaled)\n    masked = mask(normalized)\n    resized = resize(masked)\n    \n    return resized\n\ndef processVolume(volume):\n    # Resize the 3D volume\n    def resize(volume):\n        resize_shape = (128, 128, 64)\n        \n        curr_shape = volume.shape\n        \n        height_ratio = resize_shape[0] / curr_shape[0]\n        width_ratio = resize_shape[1] / curr_shape[1]\n        depth_ratio = resize_shape[2] / curr_shape[2]\n        \n        resized = ndimage.zoom(volume, (height_ratio, width_ratio, depth_ratio))\n        \n        return resized\n    \n    volume = np.swapaxes(np.swapaxes(volume, 0, 1), 1, 2)\n    resized = resize(volume)\n    resized = resized.astype('int8')\n    \n    return resized","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extract 3D CT-scan from each patient\nos.chdir('/kaggle/input/osic-pulmonary-fibrosis-progression')\n\nCT_scans = {}\n\nfor patient in os.listdir('train'):\n    try:\n        patient_path = os.path.join('train', patient)\n\n        files = os.listdir(patient_path)\n        files = sorted(files, key=lambda x: int(x[:-4]))\n        \n        ratio = int(len(files) // 64) if len(files) >= 64 else 1\n        \n        volume = []\n        \n        # Extract CT-scan slices from patient's .dicom files\n        for index, file in enumerate(files):\n            if index % ratio == 0:\n                file_path = os.path.join(patient_path, file)\n\n                dcm = dicom.dcmread(file_path)\n\n                image = dcm.pixel_array\n                processed = processImage(image, dcm.RescaleSlope, dcm.RescaleIntercept)\n\n                volume.append(processed)\n        \n        volume = np.array(volume)        \n        volume = processVolume(volume)\n        \n        CT_scans[patient] = volume\n\n    except:\n        continue   ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Train / Test sets\nx_train = []\ny_train = []\n\n# Augment CT-scans\nfor item in CT_scans.items():\n    patient = item[0]\n    volume = item[1]\n    \n    volume[volume < 0] = 0\n    volume[volume > 1] = 1\n    \n    label = cnn_train[cnn_train.index == patient]['FVC_0'][0]\n    \n    def augmentVolume(volume, label):\n        exp_volume = np.expand_dims(volume, axis=3)\n\n        x_train.append(exp_volume)\n        y_train.append(label)      \n\n        # Rotation x5\n        for rotation in [30, 60, 90, 120, 150]:\n            augmented = ndimage.rotate(volume, rotation, reshape=False, mode='nearest')\n\n            augmented[augmented < 0] = 0\n            augmented[augmented > 1] = 1\n\n            exp_augmented = np.expand_dims(augmented, axis=3)\n\n            x_train.append(exp_augmented)\n            y_train.append(label)      \n\n    augmentVolume(volume, label)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prepare train / test sets\nx_train = np.array(x_train)\ny_train = np.array(y_train)\ny_train = StandardScaler().fit_transform(y_train.reshape(len(y_train), 1))[:, 0]\n\nx_train, x_test, y_train, y_test = train_test_split(x_train, y_train, test_size=0.25)\nx_train, y_train = shuffle(x_train, y_train)\nx_test, y_test = shuffle(x_test, y_test)\n\ntrain_loader = tf.data.Dataset.from_tensor_slices((x_train, y_train))\nvalidation_loader = tf.data.Dataset.from_tensor_slices((x_test, y_test))\n\nbatch_size = 8\n\ntrain_dataset = (\n    train_loader.shuffle(len(x_train))\n    .batch(batch_size)\n    .prefetch(8)\n)\n\nvalidation_dataset = (\n    validation_loader.shuffle(len(x_test))\n    .batch(batch_size)\n    .prefetch(8)\n)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 3D-CNN model\ndef build_model():\n\n    model = Sequential([\n        keras.Input((128, 128, 64, 1)),\n        layers.Conv3D(filters=64, kernel_size=3, activation=\"relu\"),\n        layers.MaxPool3D(pool_size=2),\n        layers.BatchNormalization(),\n\n        layers.Conv3D(filters=64, kernel_size=3, activation=\"relu\"),\n        layers.MaxPool3D(pool_size=2),\n        layers.BatchNormalization(),\n\n        layers.Conv3D(filters=128, kernel_size=3, activation=\"relu\"),\n        layers.MaxPool3D(pool_size=2),\n        layers.BatchNormalization(),\n\n        layers.Conv3D(filters=256, kernel_size=3, activation=\"relu\"),\n        layers.MaxPool3D(pool_size=2),\n        layers.BatchNormalization(),\n\n        layers.GlobalAveragePooling3D(),\n        layers.Dense(units=256, activation=\"relu\"),\n        layers.Dropout(0.3),\n\n        layers.Dense(units=3, activation=\"relu\", name='intermediate'),\n    \n        layers.Dense(units=1, activation=\"linear\")\n    ])\n\n    initial_learning_rate = 0.0001\n    lr_schedule = keras.optimizers.schedules.ExponentialDecay(\n        initial_learning_rate, decay_steps=100000, decay_rate=0.98, staircase=True\n    )\n    model.compile(loss='mean_squared_error', \n                  optimizer=keras.optimizers.Adam(learning_rate=lr_schedule),\n                  metrics=['mse'])\n    \n    return model\n\nmodel = build_model()\nmodel.summary()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Train model\ncallbacks = [EarlyStopping(monitor='val_loss', patience=20, verbose=1),\n             ModelCheckpoint('/kaggle/working/test.hdf5', \n                             monitor='val_loss', \n                             save_best_only=True, \n                             mode='min', \n                             verbose=1)]\n\nmodel.fit(train_dataset,\n          validation_data=validation_dataset,\n          epochs=300,\n          callbacks=callbacks)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = load_model('/kaggle/input/ctscans/test.hdf5')\n\n# Obtain feature extractor from trained model\nextractor = keras.Model(inputs=model.input,\n                        outputs=model.get_layer('intermediate').output)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extract features from each patient's CT-scan\nfeatures = {}\nf_1 = {}\nf_2 = {}\nf_3 = {}\n\nfor item in CT_scans.items():\n    patient = item[0]\n    volume = item[1]\n    feature = extractor.predict(np.array([volume]))\n    \n    features[patient] = feature[0]\n    f_1[patient] = feature[0][0]\n    f_2[patient] = feature[0][1]\n    f_3[patient] = feature[0][2]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# PART 2: Retrain Quantile Regression model using features\nos.chdir('/kaggle/input/osic-pulmonary-fibrosis-progression')\n\ntrain = pd.read_csv('train.csv')\ntest  = pd.read_csv('test.csv')\n\n# Marker to mark rows as from train or test csv\ntrain['traintest'] = 0\ntest ['traintest'] = 1\n\n# Generate format for submission.csv\nsubmission = pd.read_csv('sample_submission.csv')\nsubmission['Weeks'] = submission['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsubmission['Patient'] = submission['Patient_Week'].apply(lambda x: x.split('_')[0]) \n\ntrain = pd.concat((train, test))\ntrain.sort_values(['Patient','Weeks'], inplace=True)\n\n# Add feature columns\ntrain['f1'] = train.Patient.map(f_1)\ntrain['f2'] = train.Patient.map(f_2)\ntrain['f3'] = train.Patient.map(f_3)\n\n# Drop Patient ID00011637202177653955184\ntrain = train.dropna()\n\n# Encode string value into categorical value\ntrain['Sex'] = pd.factorize(train['Sex'])[0]\ntrain['SmokingStatus'] = pd.factorize(train['SmokingStatus'])[0]\n\n# Standardize categorical values\ntrain['Percent']       = (train['Percent'] - train['Percent'].mean()) / train['Percent'].std()\ntrain['Age']           = (train['Age'] - train['Age'].mean()) / train['Age'].std()\ntrain['Sex']           = (train['Sex'] - train['Sex'].mean()) / train['Sex'].std()\ntrain['SmokingStatus'] = (train['SmokingStatus'] - train['SmokingStatus'].mean()) / train['SmokingStatus'].std()\n\ntrain['f1']            = (train['f1'] - train['f1'].mean()) / train['f1'].std()\ntrain['f2']            = (train['f2'] - train['f2'].mean()) / train['f2'].std()\ntrain['f3']            = (train['f3'] - train['f3'].mean()) / train['f3'].std()\n\n# Train model\nmodelL = quantreg('FVC ~ Weeks+Percent+Age+Sex+SmokingStatus+f1+f2+f3', train).fit(q=0.15)\nmodel  = quantreg('FVC ~ Weeks+Percent+Age+Sex+SmokingStatus+f1+f2+f3', train).fit(q=0.50)\nmodelH = quantreg('FVC ~ Weeks+Percent+Age+Sex+SmokingStatus+f1+f2+f3', train).fit(q=0.85)\n\n# This part is used to evaluate model on train set\n# Can ignore for now\ntrain['ypredL'] = modelL.predict(train).values\ntrain['ypred']  = model.predict(train).values\ntrain['ypredH'] = modelH.predict(train).values\ntrain['ypredstd'] = 0.5 * np.abs(train['ypredH'] - train['ypred']) + \\\n                    0.5 * np.abs(train['ypred'] - train['ypredL'])\n\n# Prepare test dataframe for predicting FVC values\ndt = train.loc[train.traintest==1, ['Patient','Percent','Age','Sex','SmokingStatus','f1','f2','f3']]\ntest = pd.merge(submission, dt, on='Patient', how='left')\ntest.sort_values(['Patient','Weeks'], inplace=True)\n\n# Predict FVC values for test dataframe\ntest['ypredL'] = modelL.predict(test).values\ntest['FVC']    = model.predict(test).values\ntest['ypredH'] = modelH.predict(test).values\ntest['Confidence'] = np.abs(test['ypredH'] - test['ypredL']) / 2\n\n# Save submission file\nsubmission = test[['Patient_Week','FVC','Confidence']]\nsubmission.to_csv('/kaggle/working/submission.csv', index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Additional function used to make CT-scan gif\ndef make_gif():\n    os.chdir('/kaggle/input/osic-pulmonary-fibrosis-progression')\n    \n    # Make gif of CT scan from patient 0\n    patient = os.listdir('train')[0]\n    patient_path = os.path.join('train', patient)\n\n    files = os.listdir(patient_path)\n    files = sorted(files, key=lambda x: int(x[:-4]))\n\n    for index, file in enumerate(files):\n        # Select just 1/10 of the total slices\n        if index % 10 == 0:\n            file_path = os.path.join(patient_path, file)\n\n            dcm = dicom.dcmread(file_path)\n\n            image = dcm.pixel_array\n            processed = processImage(image, dcm.RescaleSlope, dcm.RescaleIntercept)\n            masked = mask(processed)\n\n            # Save slices as .png files first\n            gif_original_path = os.path.join('/kaggle/working/original', '{}.png'.format(index))\n            gif_masked_path = os.path.join('/kaggle/working/masked', '{}.png'.format(index))\n\n            matplotlib.image.imsave(gif_original_path, processed)\n            matplotlib.image.imsave(gif_masked_path, masked)\n\n    os.chdir('/kaggle/working')\n            \n    # Create gif for original scan\n    original_frames = []\n\n    originals = os.listdir('original')\n    originals = sorted(originals, key=lambda x: int(x[:-4]))\n\n    for original in originals:\n        image_path = os.path.join('original', original)\n        frame = Image.open(image_path)\n        original_frames.append(frame)\n\n    original_frames[0].save('original.gif', format='GIF',\n                            append_images=original_frames[1:],\n                            save_all=True,\n                            duration=150, \n                            loop=10)\n\n    # Create gif for masked scan\n    masked_frames = []\n\n    masks = os.listdir('masked')\n    masks = sorted(masks, key=lambda x: int(x[:-4]))\n\n    for mask in masks:\n        image_path = os.path.join('masked', mask)\n        frame = Image.open(image_path)\n        masked_frames.append(frame)\n\n    masked_frames[0].save('masked.gif', format='GIF',\n                           append_images=masked_frames[1:],\n                           save_all=True,\n                           duration=150, \n                           loop=10)\n    \nshow_gif(filename=\"masked.gif\", format='png', width=200, height=200)","metadata":{},"execution_count":null,"outputs":[]}]}