{"cells":[{"metadata":{},"cell_type":"markdown","source":"### Pipeline for Stain Normilized Training (Keras)\nWith this kernel you can train a Deep Neural Network with stain normalized images. The basic idea is to convert train and test images to the similar color space. All source images are transformed on the basis of a randomly selected target image. During the transformation, however, I had problems with the some images leading to a singular matrix during the conversion. I found these images in the train and test set. Here are some papers describing Stain Normalization:\n\n- [Neural Stain Normalization and Unsupervised Classification of Cell Nuclei in Histopathological Breast Cancer Images](http://https://arxiv.org/abs/1811.03815)\n- [The importance of stain normalization in colorectal tissue classification with convolutional networks](http://https://arxiv.org/abs/1702.05931)\n- [Stain normalization of histopathology images using generative adversarial networks](http://https://ieeexplore.ieee.org/document/8363641)\n\nI took and changed the main functions from https://github.com/Peter554/StainTools\n"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"from sklearn.utils import shuffle\nimport pandas as pd\nimport os\n\n# Train \ndf = pd.read_csv(\"../input/histopathologic-cancer-detection/train_labels.csv\")\n\n# Test\ndf_test = pd.read_csv('../input/histopathologic-cancer-detection/sample_submission.csv')\n\n# Test cleaned\ndf_test_cleaned = pd.read_csv('../input/histopathologic-cancer-detection/sample_submission.csv')\n\ndf = shuffle(df,random_state=123)\n\nprint(len(df))\nprint(len(df_test))\nprint(len(df_test_cleaned))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"aa8a97c6c11e2b9d7390e8f0521884f1e6142edd","_kg_hide-input":false},"cell_type":"code","source":"# Remove error mages from the training set\n\ndf_train_error1 = pd.read_csv('../input/train-error-images1/train_error_images1.csv')\ndf_train_error2 = pd.read_csv('../input/train-error-images2/train_error_images2.csv')\ndf_train_error3 = pd.read_csv('../input/train-error-images3/train_error_images3.csv')\ndf_train_error3 = df_train_error3.drop(df_train_error3.columns[0], axis=1)\ndf_train_error1.columns = ['id','label']\ndf_train_error2.columns = ['id','label']\ndf_train_error3.columns = ['id','label']\ndf_train_error = pd.concat([df_train_error1,df_train_error1,df_train_error1])\ndf_train_error = df_train_error.drop_duplicates(subset='id', keep='first')\nfor i in range(len(df_train_error)):\n     df = df[df['id'] != df_train_error.iloc[i,0]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7ceaed5cf0ef1ece34685f48d2fac71b2d9073df","_kg_hide-input":false},"cell_type":"code","source":"# Remove error images from the test set\n\ndf_test_error = pd.read_csv('../input/test-error-images1/test_error_images.csv')\ndf_test_error = df_test_error.drop(df_test_error.columns[0], axis=1)\ndf_test_error = df_test_error.drop_duplicates(subset='id', keep='first')\n\nfor error in df_test_error['id']:\n     df_test_cleaned = df_test_cleaned[df_test_cleaned['id'] != error]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"16f521f21a4a5fa6a0c0d8351c9e16143dd75651"},"cell_type":"code","source":"# Reduce train set for demonstration only\n\ndf = df[0:5000]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ef9a6aee05dec5ec4ba97fd5e8d2724ff4f743ae","_kg_hide-input":false},"cell_type":"code","source":"# Split data set  to train and validation sets\n\nfrom sklearn.model_selection import train_test_split\n\n# Use stratify= df['label'] to get balance ratio 1/1 in train and validation sets\ndf_train, df_val = train_test_split(df, test_size=0.1, stratify= df['label'], random_state=123)\n\n# Check balancing\nprint(\"Train data: \" + str(len(df_train[df_train[\"label\"] == 1]) + len(df_train[df_train[\"label\"] == 0])))\nprint(\"True positive in train data: \" +  str(len(df_train[df_train[\"label\"] == 1])))\nprint(\"True negative in train data: \" +  str(len(df_train[df_train[\"label\"] == 0])))\nprint(\"Valid data: \" + str(len(df_val[df_val[\"label\"] == 1]) + len(df_val[df_val[\"label\"] == 0])))\nprint(\"True positive in validation data: \" +  str(len(df_val[df_val[\"label\"] == 1])))\nprint(\"True negative in validation data: \" +  str(len(df_val[df_val[\"label\"] == 0])))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"89f0b8d6439c876464b9e44062852a9eb747426f","_kg_hide-input":false},"cell_type":"code","source":"# Train List\ntrain_list = df_train['id'].tolist()\ntrain_list = ['../input/histopathologic-cancer-detection/train/'+ name + \".tif\" for name in train_list]\n\n# Validation List\nval_list = df_val['id'].tolist()\nval_list = ['../input/histopathologic-cancer-detection/train/'+ name + \".tif\" for name in val_list]\n\n# Test list\ntest_list = df_test['id'].tolist()\ntest_list = ['../input/histopathologic-cancer-detection/test/'+ name + \".tif\" for name in test_list]\n\n# Test cleaned\ntest_cleaned_list = df_test_cleaned['id'].tolist()\ntest_cleaned_list = ['../input/histopathologic-cancer-detection/test/'+ name + \".tif\" for name in test_cleaned_list]\n\n# Test error\ntest_error_list = df_test_error['id'].tolist()\ntest_error_list = ['../input/histopathologic-cancer-detection/test/'+ name + \".tif\" for name in test_error_list]\n\n# id dictionary\nid_label_map = {k:v for k,v in zip(df.id.values, df.label.values)}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"- Demonstration of images where I could not perform stain normalization because they lead to a singular matrix during the conversion."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Example of images where you will have trouble with stain normalization\n\nimport cv2 as cv\nimport matplotlib.pyplot as plt\n\nfig = plt.figure()\nfig, ax = plt.subplots(1,3, figsize=(20,20))\n\nax[0].imshow(cv.imread(test_error_list[10]))\nax[0].set_title(\"Image from the Test Set\",fontsize=14)\nax[1].imshow(cv.imread(test_error_list[15]))\nax[1].set_title(\"Image from the Test Set\",fontsize=14)\nax[2].imshow(cv.imread(test_error_list[16]))\nax[2].set_title(\"Image from the Test Set\",fontsize=14)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5d5639a7f9ebeb8fa044f5743bb6237d73f7fe09"},"cell_type":"code","source":"# Functions for generators\n\nfrom imgaug import augmenters as iaa\nimport imgaug as ia\nimport numpy as np \nimport cv2\n\ndef get_id_from_file_path(file_path):\n    return file_path.split(os.path.sep)[-1].replace('.tif', '')\n\ndef chunker(seq, size):\n    return (seq[pos:pos + size] for pos in range(0, len(seq), size))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"96052304f0140065cbf996c54ea2b0a76cb98e5b"},"cell_type":"code","source":"# Augmentation\n\nfrom imgaug import augmenters as iaa\nimport imgaug as ia\nimport numpy as np \nimport cv2\n\ndef augmentation():\n    \n    sometimes = lambda aug: iaa.Sometimes(0.5, aug)\n    \n    seq = iaa.Sequential([\n        # Horizontal flips 50% images  \n        iaa.Fliplr(0.5),   \n        \n        # Hrizontal flips 50% images\n        iaa.Flipud(0.5),\n        \n        # Crop some of the images by 0-10% of their height/width\n        #sometimes(iaa.Crop(percent=(0, 0.1))),\n        \n        #iaa.Sometimes(0.5, iaa.GaussianBlur(sigma=(0, 0.5))), # gaussian blur\n        \n        #iaa.Sometimes(0.3, iaa.AdditiveGaussianNoise(scale=(0, 0.05*255))) # gaussian noise\n        \n        # Translate images by -20 to +20% on x- and y-axis independently:\n        #iaa.Affine(translate_percent={\"x\": -0.20}, mode=ia.ALL, cval=(0, 255)),\n        \n        # Rotate images by -45 to 45 degrees:\n        #iaa.Affine(rotate=(-45, 45)),\n        \n         ], random_order=True) # apply augmenters in random order\n    \n    return seq","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5053f430de595e2d8d07d95a92242b557ccb1b8d"},"cell_type":"code","source":"!pip install spams","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e32fb867eea587c9ac2aa0ba1883223d76f3325c"},"cell_type":"code","source":"# STAIN NORMALIZATION FUNCTIONS\n\nimport spams\n\nclass TissueMaskException(Exception):\n    pass\n\n######################################################################################################\n\ndef is_uint8_image(I):\n    if not is_image(I):\n        return False\n    if I.dtype != np.uint8:\n        return False\n    return True\n######################################################################################################\n\ndef is_image(I):\n    if not isinstance(I, np.ndarray):\n        return False\n    if not I.ndim == 3:\n        return False\n    return True\n######################################################################################################\n\ndef get_tissue_mask(I, luminosity_threshold=0.8):\n    I_LAB = cv.cvtColor(I, cv.COLOR_RGB2LAB)\n    L = I_LAB[:, :, 0] / 255.0  # Convert to range [0,1].\n    mask = L < luminosity_threshold\n\n    # Check it's not empty\n    if mask.sum() == 0:\n        raise TissueMaskException(\"Empty tissue mask computed\")\n\n    return mask\n\n######################################################################################################\n\ndef convert_RGB_to_OD(I):\n    mask = (I == 0)\n    I[mask] = 1\n    \n\n    #return np.maximum(-1 * np.log(I / 255), 1e-6)\n    return np.maximum(-1 * np.log(I / 255), np.zeros(I.shape) + 0.1)\n\n######################################################################################################\n\ndef convert_OD_to_RGB(OD):\n    \n    assert OD.min() >= 0, \"Negative optical density.\"\n    \n    OD = np.maximum(OD, 1e-6)\n    \n    return (255 * np.exp(-1 * OD)).astype(np.uint8)\n\n######################################################################################################\n\ndef normalize_matrix_rows(A):\n    return A / np.linalg.norm(A, axis=1)[:, None]\n\n######################################################################################################\n\n\ndef get_concentrations(I, stain_matrix, regularizer=0.01):\n    OD = convert_RGB_to_OD(I).reshape((-1, 3))\n    return spams.lasso(X=OD.T, D=stain_matrix.T, mode=2, lambda1=regularizer, pos=True).toarray().T\n\n######################################################################################################\n\ndef get_stain_matrix(I, luminosity_threshold=0.8, angular_percentile=99):\n\n    # Convert to OD and ignore background\n    tissue_mask = get_tissue_mask(I, luminosity_threshold=luminosity_threshold).reshape((-1,))\n    OD = convert_RGB_to_OD(I).reshape((-1, 3))\n    \n    OD = OD[tissue_mask]\n\n    # Eigenvectors of cov in OD space (orthogonal as cov symmetric)\n    _, V = np.linalg.eigh(np.cov(OD, rowvar=False))\n\n    # The two principle eigenvectors\n    V = V[:, [2, 1]]\n\n    # Make sure vectors are pointing the right way\n    if V[0, 0] < 0: V[:, 0] *= -1\n    if V[0, 1] < 0: V[:, 1] *= -1\n\n    # Project on this basis.\n    That = np.dot(OD, V)\n\n    # Angular coordinates with repect to the prinicple, orthogonal eigenvectors\n    phi = np.arctan2(That[:, 1], That[:, 0])\n\n    # Min and max angles\n    minPhi = np.percentile(phi, 100 - angular_percentile)\n    maxPhi = np.percentile(phi, angular_percentile)\n\n    # the two principle colors\n    v1 = np.dot(V, np.array([np.cos(minPhi), np.sin(minPhi)]))\n    v2 = np.dot(V, np.array([np.cos(maxPhi), np.sin(maxPhi)]))\n\n    # Order of H and E.\n    # H first row.\n    if v1[0] > v2[0]:\n        HE = np.array([v1, v2])\n    else:\n        HE = np.array([v2, v1])\n\n    return normalize_matrix_rows(HE)\n\n######################################################################################################\n\ndef mapping(target,source):\n    \n    stain_matrix_target = get_stain_matrix(target)\n    target_concentrations = get_concentrations(target,stain_matrix_target)\n    maxC_target = np.percentile(target_concentrations, 99, axis=0).reshape((1, 2))\n    stain_matrix_target_RGB = convert_OD_to_RGB(stain_matrix_target) \n    \n    stain_matrix_source = get_stain_matrix(source)\n    source_concentrations = get_concentrations(source, stain_matrix_source)\n    maxC_source = np.percentile(source_concentrations, 99, axis=0).reshape((1, 2))\n    source_concentrations *= (maxC_target / maxC_source)\n    tmp = 255 * np.exp(-1 * np.dot(source_concentrations, stain_matrix_target))\n    return tmp.reshape(source.shape).astype(np.uint8)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3859fffec087e71acda3cc992290245152759829"},"cell_type":"code","source":"# Show example stain transformation\nimport cv2 as cv\nimport matplotlib.pyplot as plt\n\ntarget = cv.imread(test_list[10])\nsource = cv.imread(train_list[555])\n\n# Convert from cv2 standard of BGR to our convention of RGB.\ntarget = cv.cvtColor(target, cv.COLOR_BGR2RGB)\nsource = cv.cvtColor(source, cv.COLOR_BGR2RGB)\n\n# Perform stain normalization\ntransformed = mapping(target,source)\n\nfig = plt.figure()\nfig, ax = plt.subplots(1,3, figsize=(20,20))\n\nax[0].imshow(source)\nax[0].set_title(\"Source Image\",fontsize=14)\nax[1].imshow(target)\nax[1].set_title(\"Target Image\",fontsize=14)\nax[2].imshow(transformed)\nax[2].set_title(\"Transformed Image\",fontsize=14)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"94bfd61b2dfcd0f4616e47ede87ab9ddb1a2f83e"},"cell_type":"code","source":"# Import Pretrained Models\nfrom keras.applications.densenet import DenseNet169, preprocess_input\nfrom keras.layers import Dense, Input, Dropout, MaxPooling2D, Concatenate, GlobalAveragePooling2D, GlobalMaxPooling2D, Flatten, Concatenate\nfrom keras.models import Model\nfrom keras.preprocessing.image import ImageDataGenerator\nfrom keras.layers.normalization import BatchNormalization\nfrom keras.layers import Activation\nfrom keras.models import Sequential\nfrom keras import applications\nfrom keras.layers import Dense, Input, Dropout, MaxPooling2D, Concatenate, GlobalAveragePooling2D, GlobalMaxPooling2D, Flatten, Concatenate\n\n# Load pretrained model\ninputTensor = Input((96,96,3))\nmodel_DenseNet169 = DenseNet169(include_top=False, weights='imagenet')\n\n# Concatenate Pretrained Models\n######\n#models = [model_NASNet,model_DenseNet201]\nmodels = [model_DenseNet169]\n######\n\noutputTensors = [m(inputTensor) for m in models]\nif len(models) > 1:\n    output = Concatenate()(outputTensors) \nelse:\n    output = outputTensors[0]\n    \n# Classifier \nout1 = GlobalMaxPooling2D()(output)\nout2 = GlobalAveragePooling2D()(output)\nout = Concatenate(axis=-1)([out1, out2])\nout = Dropout(0.8)(out)\nout = Dense(1, activation=\"sigmoid\")(out)\nmodel = Model(inputTensor,out)\nmodel.summary()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5b8619e486a52a3d7589e7ee1eef1d87ace6c856"},"cell_type":"code","source":"# Read and convert images to rgb\nimport cv2 as cv\nimport os\n\ndef read_image(path):\n    im = cv.imread(path)\n    # Convert from cv2 standard of BGR to our convention of RGB.\n    im = cv.cvtColor(im, cv.COLOR_BGR2RGB)\n    return im","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Stain Normalization on the fly\n\nfrom random import randint\nimport random\n\nimport cv2 as cv\n\n\ndef stain_normalization(sources,target_list = test_list):\n    \n    X = []\n    # Chose any random target iamge for stain normalization\n    target = cv.imread(train_list[random.randint(0,len(train_list))])\n    target = cv.cvtColor(target, cv.COLOR_BGR2RGB)\n    \n    for source in sources:     \n        \n        # Perform stain normalization\n        transformed = mapping(target,source)\n        \n        X.append(transformed)\n        \n    return np.asarray(X)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f5eb4e06ec6fdcd61241ae2a4ca4576620e7759c","scrolled":false},"cell_type":"code","source":"# Train Generator\n\ndef train_gen(list_files, id_label_map, batch_size, augment = False, stain = False):\n    \n    while True:\n        \n        shuffle(list_files)\n        \n        for batch in chunker(list_files, batch_size):\n            \n            X = [read_image(x) for x in batch]\n            Y = [id_label_map[get_id_from_file_path(x)] for x in batch]      \n            \n            if stain:\n                X = stain_normalization(X)\n            \n            if augment:\n                X = augmentation().augment_images(X)   \n                \n            X = [preprocess_input(x) for x in X]   \n            yield np.array(X), np.array(Y)\n   ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"112181d488af001b0f1f485a021c8afb0e48da20","scrolled":true},"cell_type":"code","source":"# Train the model\nfrom keras.optimizers import Adam\nfrom keras.callbacks import ModelCheckpoint\n\nbatch_size = 64\nstain = True\n\n\nmodel.compile(loss='binary_crossentropy', optimizer=Adam(0.0001),metrics=['accuracy'])\n\ncallbacks =  ModelCheckpoint('model.h5', monitor='val_loss', save_best_only=True, mode='min',verbose=1)\n\nhistory = model.fit_generator(train_gen(train_list, id_label_map, batch_size, augment = False, stain = stain),\n                              validation_data=train_gen(val_list, id_label_map, batch_size, augment = False, stain = stain),\n                              epochs = 10,\n                              steps_per_epoch = len(train_list) // batch_size + 1,\n                              validation_steps = len(val_list) // batch_size + 1,\n                              callbacks=[callbacks],\n                              verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1f4037334ea852917034daf037bc976e17fed9ac","scrolled":true},"cell_type":"code","source":"# Plot validation and accuracies over epochs\nval_loss = history.history['val_loss']\nloss = history.history['loss']\n\nplt.plot(range(len(val_loss)),val_loss,'c',label='Validation loss')\nplt.plot(range(len(loss)),loss,'m',label='Train loss')\n\nplt.title('Training and validation losses')\nplt.legend()\nplt.xlabel('epochs')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"141fb90b20395685f6c426e8e5c5991b0d09b749"},"cell_type":"code","source":"# Predict cleaned stain normalized test data\nmodel.load_weights('model.h5')\n\npreds = []\nids = []\n\ni = len(test_cleaned_list)\nfor batch in chunker(test_cleaned_list[:10], batch_size):\n    X = [read_image(x) for x in batch]\n    X = stain_normalization(X)\n    X = [preprocess_input(x) for x in X]\n    ids_batch = [get_id_from_file_path(x) for x in batch]\n    X = np.array(X)\n    preds_batch = model.predict(X).ravel().tolist()\n    preds += preds_batch\n    ids += ids_batch\n    print(i)\n    i = i - 32\n    \ndf_subm_cleaned = pd.DataFrame({'id':ids, 'label':preds})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dc867647c6779ab525b2300c47e80c6f8975809a"},"cell_type":"code","source":"# Predict error test data without stain normalization\n\npreds = []\nids = []\n\nfor batch in chunker(test_error_list, batch_size):\n    X = [read_image(x) for x in batch]\n    X = [preprocess_input(x) for x in X]\n    ids_batch = [get_id_from_file_path(x) for x in batch]\n    X = np.array(X)\n    preds_batch = model.predict(X).ravel().tolist()\n    preds += preds_batch\n    ids += ids_batch\n\ndf_subm_error = pd.DataFrame({'id':ids, 'label':preds})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_subm = pd.concat([df_subm_cleaned,df_subm_error])\ndf_subm.to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_subm.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}