{"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":"<center>\n    <h1>Converting the images to 256 x 256</h1>\n<center>","metadata":{"execution":{"iopub.status.busy":"2022-07-04T13:07:53.895023Z","iopub.execute_input":"2022-07-04T13:07:53.895609Z","iopub.status.idle":"2022-07-04T13:07:53.926544Z","shell.execute_reply.started":"2022-07-04T13:07:53.895491Z","shell.execute_reply":"2022-07-04T13:07:53.925036Z"}}},{"cell_type":"markdown","source":"<center>\n<img src=\"https://hubmapconsortium.org/wp-content/uploads/2019/01/HuBMAP-Retina-Logo-Color.png\">\n</center>","metadata":{}},{"cell_type":"markdown","source":"## Converting the images to 256 x 256\n\nThis notebook converts the images to 256x256 image dataset for the initial fast prototyping. Based on the size of the detected features, you can chose which dataset size you want to train with. \n\nI've made serveral options so you can experiment with them on your own:\n\n**Precomputed Datasets**\n\n#### [Dataset (512 x 512)](https://www.kaggle.com/datasets/thedevastator/hubmap-2022-512x512/)\n\n#### [Dataset (256 x 256)](https://www.kaggle.com/datasets/thedevastator/hubmap-2022-256x256/)\n\n#### [Dataset (128 x 128)](https://www.kaggle.com/datasets/thedevastator/hubmap-2022-128x128/settings)\n\n_____\n\n### Description \n\nThis notebook is part of a series providing a fast.ai starter pytorch pipeline based on a U-shape network (UneXt50) that was used on multiple competitions in the past and includes several tricks from the previous segmentation competitions.\nIt is [dividing the images into tiles](https://www.kaggle.com/code/thedevastator/converting-to-256x256), selection of tiles with tissue, evaluation of the predictions of multiple models with TTA, combining the tile masks back into image level masks, and conversion into RLE. The [inference](https://www.kaggle.com/code/thedevastator/inference-fastai-baseline) is performed based on models trained in the [fast.ai training notebook](https://www.kaggle.com/code/thedevastator/training-fastai-baseline).\n\n**Training & Dataset Creation**\n\n- #### Training Notebook [here](https://www.kaggle.com/code/thedevastator/training-fastai-baseline). \n- #### Inference Notebook [here](https://www.kaggle.com/code/thedevastator/inference-fastai-baseline). \n\n\n\n> **Credit**\n> #### Based on the excellent [notebook](https://www.kaggle.com/iafoss/256x256-images) by [iafoss](https://www.kaggle.com/iafoss) \n> All credit to belongs to the original author!","metadata":{"papermill":{"duration":0.007474,"end_time":"2021-03-12T06:31:41.563501","exception":false,"start_time":"2021-03-12T06:31:41.556027","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import gc\nimport os\nimport cv2\nimport zipfile\nimport rasterio\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image\nimport tifffile as tiff\nfrom tqdm.notebook import tqdm\nimport matplotlib.pyplot as plt\nfrom rasterio.windows import Window\nfrom torch.utils.data import Dataset","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":1.650158,"end_time":"2021-03-12T06:31:43.222111","exception":false,"start_time":"2021-03-12T06:31:41.571953","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-06-22T20:03:29.729153Z","iopub.execute_input":"2022-06-22T20:03:29.729802Z","iopub.status.idle":"2022-06-22T20:03:33.517128Z","shell.execute_reply.started":"2022-06-22T20:03:29.72968Z","shell.execute_reply":"2022-06-22T20:03:33.516013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"OUT_TRAIN = 'train.zip'\nOUT_MASKS = 'masks.zip'\nsz = 256   # the size of tiles\nreduce = 4 # reduce the original images by 4 times \nMASKS = '../input/hubmap-organ-segmentation/train.csv'\nDATA = '../input/hubmap-organ-segmentation/train_images'","metadata":{"papermill":{"duration":0.021573,"end_time":"2021-03-12T06:31:43.251307","exception":false,"start_time":"2021-03-12T06:31:43.229734","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-06-22T21:13:34.359567Z","iopub.execute_input":"2022-06-22T21:13:34.360348Z","iopub.status.idle":"2022-06-22T21:13:34.391588Z","shell.execute_reply.started":"2022-06-22T21:13:34.360241Z","shell.execute_reply":"2022-06-22T21:13:34.390298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# functions to convert encoding to mask and mask to encoding\ndef enc2mask(encs, shape):\n    img = np.zeros(shape[0]*shape[1], dtype=np.uint8)\n    for m,enc in enumerate(encs):\n        if isinstance(enc,np.float) and np.isnan(enc): continue\n        s = enc.split()\n        for i in range(len(s)//2):\n            start = int(s[2*i]) - 1\n            length = int(s[2*i+1])\n            img[start:start+length] = 1 + m\n    return img.reshape(shape).T\n\ndef mask2enc(mask, n=1):\n    pixels = mask.T.flatten()\n    encs = []\n    for i in range(1,n+1):\n        p = (pixels == i).astype(np.int8)\n        if p.sum() == 0: encs.append(np.nan)\n        else:\n            p = np.concatenate([[0], p, [0]])\n            runs = np.where(p[1:] != p[:-1])[0] + 1\n            runs[1::2] -= runs[::2]\n            encs.append(' '.join(str(x) for x in runs))\n    return encs\n\ndf_masks = pd.read_csv(MASKS)[['id', 'rle']].set_index('id')\ndf_masks.head()","metadata":{"papermill":{"duration":0.472979,"end_time":"2021-03-12T06:31:43.731052","exception":false,"start_time":"2021-03-12T06:31:43.258073","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-06-22T20:06:42.706435Z","iopub.execute_input":"2022-06-22T20:06:42.706789Z","iopub.status.idle":"2022-06-22T20:06:42.860288Z","shell.execute_reply.started":"2022-06-22T20:06:42.70672Z","shell.execute_reply":"2022-06-22T20:06:42.859244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# one of the new images cannot be loaded into 16GB RAM\n# use rasterio to load image part by part\n# using a dataset similar to my submission kernel\n\ns_th = 40  # saturation blancking threshold\np_th = 1000*(sz // 256) ** 2 # threshold for the minimum number of pixels\n\n\nclass HuBMAPDataset(Dataset):\n    def __init__(self, idx, sz=sz, reduce=reduce, encs=None):\n        self.data = rasterio.open(os.path.join(DATA,str(idx)+'.tiff'),num_threads='all_cpus')\n        # some images have issues with their format \n        # and must be saved correctly before reading with rasterio\n        if self.data.count != 3:\n            subdatasets = self.data.subdatasets\n            self.layers = []\n            if len(subdatasets) > 0:\n                for i, subdataset in enumerate(subdatasets, 0):\n                    self.layers.append(rasterio.open(subdataset))\n        self.shape = self.data.shape\n        self.reduce = reduce\n        self.sz = reduce*sz\n        self.pad0 = (self.sz - self.shape[0]%self.sz)%self.sz\n        self.pad1 = (self.sz - self.shape[1]%self.sz)%self.sz\n        self.n0max = (self.shape[0] + self.pad0)//self.sz\n        self.n1max = (self.shape[1] + self.pad1)//self.sz\n        self.mask = enc2mask(encs,(self.shape[1],self.shape[0])) if encs is not None else None\n        \n    def __len__(self):\n        return self.n0max*self.n1max\n    \n    def __getitem__(self, idx):\n        # the code below may be a little bit difficult to understand,\n        # but the thing it does is mapping the original image to\n        # tiles created with adding padding (like in the previous version of the kernel)\n        # then the tiles are loaded with rasterio\n        # n0,n1 - are the x and y index of the tile (idx = n0*self.n1max + n1)\n        n0,n1 = idx//self.n1max, idx%self.n1max\n        # x0,y0 - are the coordinates of the lower left corner of the tile in the image\n        # negative numbers correspond to padding (which must not be loaded)\n        x0,y0 = -self.pad0//2 + n0*self.sz, -self.pad1//2 + n1*self.sz\n\n        # make sure that the region to read is within the image\n        p00,p01 = max(0,x0), min(x0+self.sz,self.shape[0])\n        p10,p11 = max(0,y0), min(y0+self.sz,self.shape[1])\n        img = np.zeros((self.sz,self.sz,3),np.uint8)\n        mask = np.zeros((self.sz,self.sz),np.uint8)\n        # mapping the loade region to the tile\n        if self.data.count == 3:\n            img[(p00-x0):(p01-x0),(p10-y0):(p11-y0)] = np.moveaxis(self.data.read([1,2,3],\n                window=Window.from_slices((p00,p01),(p10,p11))), 0, -1)\n        else:\n            for i,layer in enumerate(self.layers):\n                img[(p00-x0):(p01-x0),(p10-y0):(p11-y0),i] =\\\n                  layer.read(1,window=Window.from_slices((p00,p01),(p10,p11)))\n        if self.mask is not None: mask[(p00-x0):(p01-x0),(p10-y0):(p11-y0)] = self.mask[p00:p01,p10:p11]\n        \n        if self.reduce != 1:\n            img = cv2.resize(img,(self.sz//reduce,self.sz//reduce),\n                             interpolation = cv2.INTER_AREA)\n            mask = cv2.resize(mask,(self.sz//reduce,self.sz//reduce),\n                             interpolation = cv2.INTER_NEAREST)\n        #check for empty imges\n        hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)\n        h,s,v = cv2.split(hsv)\n        #return -1 for empty images\n        return img, mask, (-1 if (s>s_th).sum() <= p_th or img.sum() <= p_th else idx)","metadata":{"papermill":{"duration":0.042498,"end_time":"2021-03-12T06:31:43.781178","exception":false,"start_time":"2021-03-12T06:31:43.73868","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-06-22T20:06:44.481788Z","iopub.execute_input":"2022-06-22T20:06:44.48215Z","iopub.status.idle":"2022-06-22T20:06:44.508378Z","shell.execute_reply.started":"2022-06-22T20:06:44.482116Z","shell.execute_reply":"2022-06-22T20:06:44.507004Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x_tot,x2_tot = [],[]\nwith zipfile.ZipFile(OUT_TRAIN, 'w') as img_out,\\\n zipfile.ZipFile(OUT_MASKS, 'w') as mask_out:\n    for index, encs in tqdm(df_masks.iterrows(),total=len(df_masks)):\n        #image+mask dataset\n        ds = HuBMAPDataset(index,encs=encs)\n        for i in range(len(ds)):\n            im,m,idx = ds[i]\n            if idx < 0: continue\n                \n            x_tot.append((im/255.0).reshape(-1,3).mean(0))\n            x2_tot.append(((im/255.0)**2).reshape(-1,3).mean(0))\n            \n            #write data   \n            im = cv2.imencode('.png',cv2.cvtColor(im, cv2.COLOR_RGB2BGR))[1]\n            img_out.writestr(f'{index}_{idx:04d}.png', im)\n            m = cv2.imencode('.png',m)[1]\n            mask_out.writestr(f'{index}_{idx:04d}.png', m)\n        \n#image stats\nimg_avr =  np.array(x_tot).mean(0)\nimg_std =  np.sqrt(np.array(x2_tot).mean(0) - img_avr**2)\nprint('mean:',img_avr, ', std:', img_std)","metadata":{"papermill":{"duration":827.081467,"end_time":"2021-03-12T06:45:30.870295","exception":false,"start_time":"2021-03-12T06:31:43.788828","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-06-22T20:06:52.988971Z","iopub.execute_input":"2022-06-22T20:06:52.989393Z","iopub.status.idle":"2022-06-22T20:11:00.366494Z","shell.execute_reply.started":"2022-06-22T20:06:52.989363Z","shell.execute_reply":"2022-06-22T20:11:00.3646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"columns, rows = 4,4\nidx0 = 20\nfig=plt.figure(figsize=(columns*4, rows*4))\nwith zipfile.ZipFile(OUT_TRAIN, 'r') as img_arch, \\\n     zipfile.ZipFile(OUT_MASKS, 'r') as msk_arch:\n    fnames = sorted(img_arch.namelist())[8:]\n    for i in range(rows):\n        for j in range(columns):\n            idx = i+j*columns\n            img = cv2.imdecode(np.frombuffer(img_arch.read(fnames[idx0+idx]), \n                                             np.uint8), cv2.IMREAD_COLOR)\n            img = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)\n            mask = cv2.imdecode(np.frombuffer(msk_arch.read(fnames[idx0+idx]), \n                                              np.uint8), cv2.IMREAD_GRAYSCALE)\n    \n            fig.add_subplot(rows, columns, idx+1)\n            plt.axis('off')\n            plt.imshow(Image.fromarray(img))\n            plt.imshow(Image.fromarray(mask), alpha=0.2)\nplt.show()","metadata":{"papermill":{"duration":1.929113,"end_time":"2021-03-12T06:45:32.808849","exception":false,"start_time":"2021-03-12T06:45:30.879736","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-06-22T20:11:00.369207Z","iopub.execute_input":"2022-06-22T20:11:00.369646Z","iopub.status.idle":"2022-06-22T20:11:02.570218Z","shell.execute_reply.started":"2022-06-22T20:11:00.369615Z","shell.execute_reply":"2022-06-22T20:11:02.56844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"papermill":{"duration":0.049358,"end_time":"2021-03-12T06:45:32.907031","exception":false,"start_time":"2021-03-12T06:45:32.857673","status":"completed"},"tags":[]},"execution_count":null,"outputs":[]}]}