{"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":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":0.022738,"end_time":"2022-11-25T01:43:50.180306","exception":false,"start_time":"2022-11-25T01:43:50.157568","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-01-31T10:57:50.83499Z","iopub.execute_input":"2023-01-31T10:57:50.835344Z","iopub.status.idle":"2023-01-31T10:57:50.84087Z","shell.execute_reply.started":"2023-01-31T10:57:50.835313Z","shell.execute_reply":"2023-01-31T10:57:50.839913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from PIL import Image\nimport tensorflow as tf\nimport tensorflow.keras.preprocessing as image\nimport tensorflow_io as tfio\nfrom sklearn.model_selection import train_test_split\nimport cv2\nimport tifffile as tifi\nimport gc\nimport os\nimport openslide\nfrom openslide import OpenSlide\nimport math\nfrom keras.layers.merge import concatenate\nfrom keras.layers import Input\nimport tensorflow as tf\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator\nfrom tensorflow.keras.layers import Conv2D, MaxPooling2D,GlobalAveragePooling2D, Flatten, Dense, Dropout, BatchNormalization\nfrom tensorflow.keras import layers\nfrom tensorflow.keras.applications.densenet import DenseNet169\nfrom openslide import open_slide\nfrom matplotlib import pyplot as plt\nimport pandas as pd\ninp_size=512","metadata":{"papermill":{"duration":7.712145,"end_time":"2022-11-25T01:43:57.896485","exception":false,"start_time":"2022-11-25T01:43:50.18434","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-01-31T10:57:50.842926Z","iopub.execute_input":"2023-01-31T10:57:50.843821Z","iopub.status.idle":"2023-01-31T10:57:50.852879Z","shell.execute_reply.started":"2023-01-31T10:57:50.843785Z","shell.execute_reply":"2023-01-31T10:57:50.851935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_train_file(x):\n    return \"../input/mayo-clinic-strip-ai/train/\" + x + \".tif\"\n\ndef make_test_file(x):\n    return x + \".tif\"","metadata":{"papermill":{"duration":0.012515,"end_time":"2022-11-25T01:43:57.912488","exception":false,"start_time":"2022-11-25T01:43:57.899973","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-01-31T10:57:50.854815Z","iopub.execute_input":"2023-01-31T10:57:50.855819Z","iopub.status.idle":"2023-01-31T10:57:50.861693Z","shell.execute_reply.started":"2023-01-31T10:57:50.855784Z","shell.execute_reply":"2023-01-31T10:57:50.860907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv('../input/mayo-clinic-strip-ai/train.csv')\ntrain.head()\ntrain_data = pd.DataFrame({'image_id': train.image_id.apply(make_train_file), 'label': train.label})","metadata":{"papermill":{"duration":0.03399,"end_time":"2022-11-25T01:43:57.950017","exception":false,"start_time":"2022-11-25T01:43:57.916027","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-01-31T10:57:50.864962Z","iopub.execute_input":"2023-01-31T10:57:50.865338Z","iopub.status.idle":"2023-01-31T10:57:50.877845Z","shell.execute_reply.started":"2023-01-31T10:57:50.865303Z","shell.execute_reply":"2023-01-31T10:57:50.876861Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_data.head()","metadata":{"papermill":{"duration":0.023583,"end_time":"2022-11-25T01:43:57.977105","exception":false,"start_time":"2022-11-25T01:43:57.953522","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-01-31T10:57:50.879307Z","iopub.execute_input":"2023-01-31T10:57:50.879676Z","iopub.status.idle":"2023-01-31T10:57:50.889032Z","shell.execute_reply.started":"2023-01-31T10:57:50.879641Z","shell.execute_reply":"2023-01-31T10:57:50.887958Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_data.label.value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-01-31T10:57:50.891044Z","iopub.execute_input":"2023-01-31T10:57:50.89144Z","iopub.status.idle":"2023-01-31T10:57:50.901265Z","shell.execute_reply.started":"2023-01-31T10:57:50.891368Z","shell.execute_reply":"2023-01-31T10:57:50.899933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# path = '../input/mayo-clinic-strip-ai/train/'","metadata":{"papermill":{"duration":0.01179,"end_time":"2022-11-25T01:43:57.992571","exception":false,"start_time":"2022-11-25T01:43:57.980781","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-01-31T10:57:50.90289Z","iopub.execute_input":"2023-01-31T10:57:50.903227Z","iopub.status.idle":"2023-01-31T10:57:50.907679Z","shell.execute_reply.started":"2023-01-31T10:57:50.903194Z","shell.execute_reply":"2023-01-31T10:57:50.906631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"id=0\nsample = train_data.image_id[id]\nlabel = train_data.label[id]\nprint(label)\nslide = open_slide(sample)\nprint(slide.properties)\nprint(slide.dimensions)\nprint(slide.level_dimensions)\n#in the slide object, we have only 1 level bcoz it is the way it was stored originally\n#so that is why we can only view images of original resolution\nslide=slide.read_region((1300,1900),0,(5000,5000))\nslide=slide.convert('RGB')\nslide=np.array(slide)\nplt.imshow(slide)","metadata":{"execution":{"iopub.status.busy":"2023-01-31T10:57:50.909346Z","iopub.execute_input":"2023-01-31T10:57:50.910157Z","iopub.status.idle":"2023-01-31T10:57:55.120561Z","shell.execute_reply.started":"2023-01-31T10:57:50.910121Z","shell.execute_reply":"2023-01-31T10:57:55.11952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from openslide.deepzoom import DeepZoomGenerator\nslide = open_slide(sample)\ntiles=DeepZoomGenerator(slide,tile_size=inp_size,overlap=0,limit_bounds=False)","metadata":{"execution":{"iopub.status.busy":"2023-01-31T10:57:55.121931Z","iopub.execute_input":"2023-01-31T10:57:55.122374Z","iopub.status.idle":"2023-01-31T10:57:55.134823Z","shell.execute_reply.started":"2023-01-31T10:57:55.122328Z","shell.execute_reply":"2023-01-31T10:57:55.134002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"number of levels : \" , tiles.level_count)\nprint(\"dimensions of each level : \" , tiles.level_dimensions)\nprint('tot num of tiles: ' , tiles.tile_count)\nprint(\"get only last level  :--------- \")\nprint(\"grid size : \" , tiles.level_tiles[tiles.level_count-1])\nprint(\"each tile dim : \" , tiles.get_tile_dimensions(tiles.level_count-1,(1,4)))","metadata":{"execution":{"iopub.status.busy":"2023-01-31T10:57:55.138408Z","iopub.execute_input":"2023-01-31T10:57:55.138677Z","iopub.status.idle":"2023-01-31T10:57:55.145243Z","shell.execute_reply.started":"2023-01-31T10:57:55.138653Z","shell.execute_reply":"2023-01-31T10:57:55.144304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def norm_HnE(img, Io=250, alpha=1, beta=0.15):\n\n\n    ######## Step 1: Convert RGB to OD ###################\n    ## reference H&E OD matrix.\n    #Can be updated if you know the best values for your image. \n    #Otherwise use the following default values. \n    #Read the above referenced papers on this topic. \n    HERef = np.array([[0.5626, 0.2159],\n                      [0.7201, 0.8012],\n                      [0.4062, 0.5581]])\n    ### reference maximum stain concentrations for H&E\n    maxCRef = np.array([1.9705, 1.0308])\n    \n    \n    # extract the height, width and num of channels of image\n    h, w, c = img.shape\n    \n    # reshape image to multiple rows and 3 columns.\n    #Num of rows depends on the image size (wxh)\n    img = img.reshape((-1,3))\n    \n    # calculate optical density\n    # OD = −log10(I)  \n    #OD = -np.log10(img+0.004)  #Use this when reading images with skimage\n    #Adding 0.004 just to avoid log of zero. \n    \n    OD = -np.log10((img.astype(np.float)+1)/Io) #Use this for opencv imread\n    #Add 1 in case any pixels in the image have a value of 0 (log 0 is indeterminate)\n    \n    \n    ############ Step 2: Remove data with OD intensity less than β ############\n    # remove transparent pixels (clear region with no tissue)\n    ODhat = OD[~np.any(OD < beta, axis=1)] #Returns an array where OD values are above beta\n    #Check by printing ODhat.min()\n    \n    ############# Step 3: Calculate SVD on the OD tuples ######################\n    #Estimate covariance matrix of ODhat (transposed)\n    # and then compute eigen values & eigenvectors.\n    eigvals, eigvecs = np.linalg.eigh(np.cov(ODhat.T))\n    \n    \n    ######## Step 4: Create plane from the SVD directions with two largest values ######\n    #project on the plane spanned by the eigenvectors corresponding to the two \n    # largest eigenvalues    \n    That = ODhat.dot(eigvecs[:,1:3]) #Dot product\n    \n    ############### Step 5: Project data onto the plane, and normalize to unit length ###########\n    ############## Step 6: Calculate angle of each point wrt the first SVD direction ########\n    #find the min and max vectors and project back to OD space\n    phi = np.arctan2(That[:,1],That[:,0])\n    \n    minPhi = np.percentile(phi, alpha)\n    maxPhi = np.percentile(phi, 100-alpha)\n    \n    vMin = eigvecs[:,1:3].dot(np.array([(np.cos(minPhi), np.sin(minPhi))]).T)\n    vMax = eigvecs[:,1:3].dot(np.array([(np.cos(maxPhi), np.sin(maxPhi))]).T)\n    \n    \n    # a heuristic to make the vector corresponding to hematoxylin first and the \n    # one corresponding to eosin second\n    if vMin[0] > vMax[0]:    \n        HE = np.array((vMin[:,0], vMax[:,0])).T\n        \n    else:\n        HE = np.array((vMax[:,0], vMin[:,0])).T\n    \n    \n    # rows correspond to channels (RGB), columns to OD values\n    Y = np.reshape(OD, (-1, 3)).T\n    \n    # determine concentrations of the individual stains\n    C = np.linalg.lstsq(HE,Y, rcond=None)[0]\n    \n    # normalize stain concentrations\n    maxC = np.array([np.percentile(C[0,:], 99), np.percentile(C[1,:],99)])\n    tmp = np.divide(maxC,maxCRef)\n    C2 = np.divide(C,tmp[:, np.newaxis])\n    \n    ###### Step 8: Convert extreme values back to OD space\n    # recreate the normalized image using reference mixing matrix \n    \n    Inorm = np.multiply(Io, np.exp(-HERef.dot(C2)))\n    Inorm[Inorm>255] = 254\n    Inorm = np.reshape(Inorm.T, (h, w, 3)).astype(np.uint8)  \n    \n    # Separating H and E components\n    \n    H = np.multiply(Io, np.exp(np.expand_dims(-HERef[:,0], axis=1).dot(np.expand_dims(C2[0,:], axis=0))))\n    H[H>255] = 254\n    H = np.reshape(H.T, (h, w, 3)).astype(np.uint8)\n    \n    E = np.multiply(Io, np.exp(np.expand_dims(-HERef[:,1], axis=1).dot(np.expand_dims(C2[1,:], axis=0))))\n    E[E>255] = 254\n    E = np.reshape(E.T, (h, w, 3)).astype(np.uint8)\n    \n    return (Inorm, H, E)","metadata":{"execution":{"iopub.status.busy":"2023-01-31T10:57:55.14664Z","iopub.execute_input":"2023-01-31T10:57:55.147143Z","iopub.status.idle":"2023-01-31T10:57:55.166393Z","shell.execute_reply.started":"2023-01-31T10:57:55.147108Z","shell.execute_reply":"2023-01-31T10:57:55.16552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#get and print a single tile\n\ntile=tiles.get_tile(tiles.level_count-1,(16,0))\ntile=tile.convert(\"RGB\")\ntile=np.array(tile)\nnorm_img,H,E= norm_HnE(tile)\nplt.figure(figsize=(12,12))\nplt.subplot(221)\nplt.title(\"original image\")\nplt.imshow(tile)\n\nplt.subplot(222)\nplt.title(\"norm image\")\nplt.imshow(norm_img)\n\nplt.subplot(223)\nplt.title(\"H image\")\nplt.imshow(H)\n\nplt.subplot(224)\nplt.title(\"E image\")\nplt.imshow(E)\nprint(tile.mean())\nprint(tile.std())\nprint(norm_img.mean())\nprint(norm_img.std())","metadata":{"execution":{"iopub.status.busy":"2023-01-31T10:57:55.167556Z","iopub.execute_input":"2023-01-31T10:57:55.167911Z","iopub.status.idle":"2023-01-31T10:57:55.948402Z","shell.execute_reply.started":"2023-01-31T10:57:55.167878Z","shell.execute_reply":"2023-01-31T10:57:55.94731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import shutil\n# shutil.rmtree(\"/kaggle/working/\")\n# # os.remove(\"/kaggle/working/download.zip\")","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:12:29.189132Z","iopub.execute_input":"2023-01-31T11:12:29.189511Z","iopub.status.idle":"2023-01-31T11:12:29.193658Z","shell.execute_reply.started":"2023-01-31T11:12:29.189478Z","shell.execute_reply":"2023-01-31T11:12:29.192614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#save all tiles of a image in to directory\n\n\nbase_dir = '/kaggle/working/DATASET/'\nCE_dir = 'ce/'\nLAA_dir = 'laa/'\n\noriginal_path = base_dir+'original/'\nnormalized_path= base_dir+'norm/'\n\n\n\nif not os.path.exists(normalized_path+LAA_dir):\n    os.makedirs(normalized_path+LAA_dir)\n    \nif not os.path.exists(normalized_path+CE_dir):\n    os.makedirs(normalized_path+CE_dir)\n    \n\nif not os.path.exists(original_path+LAA_dir):\n    os.makedirs(original_path+LAA_dir)\n    \nif not os.path.exists(original_path+CE_dir):\n    os.makedirs(original_path+CE_dir)","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:12:34.718671Z","iopub.execute_input":"2023-01-31T11:12:34.719022Z","iopub.status.idle":"2023-01-31T11:12:34.726474Z","shell.execute_reply.started":"2023-01-31T11:12:34.71899Z","shell.execute_reply":"2023-01-31T11:12:34.725542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import traceback\nce_imgs=0\nlaa_imgs=0\n\nfor x in range(int(train_data.shape[0])):\n    img_path = train_data.image_id[x]\n    label = train_data.label[x]\n    print(x,label)\n    label=label.lower()\n    slide = open_slide(img_path)\n    tiles=DeepZoomGenerator(slide,tile_size=inp_size,overlap=0,limit_bounds=False)\n    cols,rows = tiles.level_tiles[tiles.level_count-1]\n    \n    count=0\n    if label=='ce':\n        thresh=6\n    else:\n        thresh=22\n    \n    for row in range(0,rows,5):\n        for col in range(0,cols,5):\n            file_name = str(x)+'_'+ str(col)+'_'+str(row)+'.png'\n            tile=tiles.get_tile(tiles.level_count-1,(col,row))\n            tile=tile.convert(\"RGB\")\n            tile=np.array(tile)\n            try:\n                if tile.mean()<180 and tile.std()>50:\n                    norm_img,H,E= norm_HnE(tile)\n                    plt.imsave(original_path+label+'/'+file_name,tile)\n                    plt.imsave(normalized_path+label+'/'+file_name,norm_img)\n                    if label=='ce':ce_imgs+=1\n                    else:laa_imgs+=1\n                    count+=1\n                    if count>thresh:break\n                    \n            except:\n                print('exception')\n                traceback.print_exc()\n                plt.imshow(tile)\n                plt.imshow(norm_img)\n                pass\n        if count>thresh:break\n    print(x,ce_imgs,laa_imgs)","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:12:42.983917Z","iopub.execute_input":"2023-01-31T11:12:42.984273Z","iopub.status.idle":"2023-01-31T11:14:32.103552Z","shell.execute_reply.started":"2023-01-31T11:12:42.984241Z","shell.execute_reply":"2023-01-31T11:14:32.102294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import shutil\n# shutil.make_archive(\"mean_180-std_50-20_gap_dataset\", 'zip', '/kaggle/working/')\n\n# !tar -czf dataset.tar.gz /kaggle/working","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:03:18.641907Z","iopub.status.idle":"2023-01-31T11:03:18.642646Z","shell.execute_reply.started":"2023-01-31T11:03:18.642383Z","shell.execute_reply":"2023-01-31T11:03:18.642408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from IPython.display import FileLink\nFileLink(r'/kaggle/working/dataset.tar.gz')","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:03:18.64401Z","iopub.status.idle":"2023-01-31T11:03:18.644735Z","shell.execute_reply.started":"2023-01-31T11:03:18.644475Z","shell.execute_reply":"2023-01-31T11:03:18.6445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:15:01.100394Z","iopub.execute_input":"2023-01-31T11:15:01.100754Z","iopub.status.idle":"2023-01-31T11:15:01.110602Z","shell.execute_reply.started":"2023-01-31T11:15:01.100722Z","shell.execute_reply":"2023-01-31T11:15:01.109447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}