{"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":"I provide a 256x256 image dataset for the initial prototyping. Based on the size of the detected features, I'd expect that the appropriate tile size for this data should be 1024x1024. However, it would be an overshoot for the initial model development. Therefore, I reduced the images by 4 times to 256x256.\n\n* The corresponding dataset: https://www.kaggle.com/iafoss/hubmap-256x256\n* The dataset with 512x512 tiles (reduction of 2): https://www.kaggle.com/iafoss/hubmap-512x512\n* The dataset with 1024x1024 tiles (no resolution reduction): https://www.kaggle.com/iafoss/hubmap-1024x1024\n\n\n* Update (11/17): remove gray background tiles based on saturation check\n* Update (11/19): fix a bug with dimentions in cv2.resize, please use the latest version\n* Update (3/9): rerun with the new data (load images part by part because of insufficient RAM for one of the new images)","metadata":{}},{"cell_type":"code","source":"#derived from iafoss hubmap kernel\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom PIL import Image\nimport tifffile as tiff\nimport cv2\nimport os\nfrom tqdm.notebook import tqdm\nimport zipfile\nimport rasterio\nfrom rasterio.windows import Window\nfrom torch.utils.data import Dataset\nimport gc","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-08-09T09:35:46.465743Z","iopub.execute_input":"2022-08-09T09:35:46.466081Z","iopub.status.idle":"2022-08-09T09:35:46.472918Z","shell.execute_reply.started":"2022-08-09T09:35:46.466049Z","shell.execute_reply":"2022-08-09T09:35:46.471494Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sz = 256   #the size of tiles\nreduce = 4 #reduce the original images by 4 times \nMASKS = '../input/hubmap-kidney-segmentation/train.csv'\nDATA = '../input/hubmap-kidney-segmentation/train/'\nOUT_TRAIN = 'train.zip'\nOUT_MASKS = 'masks.zip'","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:35:50.355353Z","iopub.execute_input":"2022-08-09T09:35:50.356069Z","iopub.status.idle":"2022-08-09T09:35:50.360716Z","shell.execute_reply.started":"2022-08-09T09:35:50.356025Z","shell.execute_reply":"2022-08-09T09:35:50.359813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sz = 256   #the size of tiles\nreduce = 4 #reduce the original images by 4 times \nMASKS = '../input/mayo-clinic-strip-ai/train.csv'\nDATA = '../input/mayo-clinic-strip-ai/train/'\nOUT_TRAIN = 'train.zip'","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:22:54.535786Z","iopub.execute_input":"2022-08-09T11:22:54.536241Z","iopub.status.idle":"2022-08-09T11:22:54.542633Z","shell.execute_reply.started":"2022-08-09T11:22:54.536202Z","shell.execute_reply":"2022-08-09T11:22:54.540936Z"},"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).set_index('image_id')\ndf_masks.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:23:03.966607Z","iopub.execute_input":"2022-08-09T11:23:03.967178Z","iopub.status.idle":"2022-08-09T11:23:04.00935Z","shell.execute_reply.started":"2022-08-09T11:23:03.967142Z","shell.execute_reply":"2022-08-09T11:23:04.007717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import time\n#time.sleep(3660)","metadata":{"execution":{"iopub.status.busy":"2022-08-08T15:23:08.921999Z","iopub.execute_input":"2022-08-08T15:23:08.922437Z","iopub.status.idle":"2022-08-08T16:17:20.964283Z","shell.execute_reply.started":"2022-08-08T15:23:08.922403Z","shell.execute_reply":"2022-08-08T16:17:20.962678Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#(1024-1224%256)%256\n1224%256\n\n(256-1224%256)%256","metadata":{"execution":{"iopub.status.busy":"2022-08-08T16:23:18.747419Z","iopub.execute_input":"2022-08-08T16:23:18.747799Z","iopub.status.idle":"2022-08-08T16:23:18.753692Z","shell.execute_reply.started":"2022-08-08T16:23:18.747766Z","shell.execute_reply":"2022-08-08T16:23:18.752765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#a=np.array([[1,2],[3,4]])\n#np.pad(a,((1,0),(1,0))),a\na=np.arange(0,np.power(300,2))\na=a.reshape(300,300)\n\nnp.pad(a,((12,12),(12,12))).shape\nif (np.pad(a,((12,12),(12,12))).shape)!=(324,324 ):\n    print('x')","metadata":{"execution":{"iopub.status.busy":"2022-08-09T10:36:45.67608Z","iopub.execute_input":"2022-08-09T10:36:45.676807Z","iopub.status.idle":"2022-08-09T10:36:45.684278Z","shell.execute_reply.started":"2022-08-09T10:36:45.676766Z","shell.execute_reply":"2022-08-09T10:36:45.683425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import math","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:35:48.406261Z","iopub.execute_input":"2022-08-09T11:35:48.406988Z","iopub.status.idle":"2022-08-09T11:35:48.411307Z","shell.execute_reply.started":"2022-08-09T11:35:48.406949Z","shell.execute_reply":"2022-08-09T11:35:48.410386Z"},"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\nreduce=4\n\n#class MayoDataset(Dataset):\nclass HuBMAPDataset(Dataset):\n    def __init__(self, idx, sz=sz, reduce=reduce, encs=None,ext='tiff'):\n        self.data = rasterio.open(os.path.join(DATA,idx+'.'+ext),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        print('self.shape',self.shape,self.sz)\n        self.pad0 = (self.sz - self.shape[0]%self.sz)%self.sz # if shape is divisible then pad0 should be 0 so we take reminder\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        print('pad',self.pad0,self.pad1,self.n0max,self.n1max)\n        #self.shape (31278, 25794) 1024\n        #pad 466 830 31 26\n        #ds length 806\n        #xo -233 -415 31278 25794\n        if encs is not None :\n            \n            self.mask = enc2mask(encs,(self.shape[1],self.shape[0])) if encs is not None else None\n        else:\n            self.mask=None\n    def __len__(self):\n        return self.n0max*self.n1max\n    \n    def pad_tile(self,tmp):\n        if tmp.shape[0]<self.sz:\n            #pad=self.pad0//2+1 if self.pad0%2!=0 else self.pad0//2\n            pad=(int(self.pad0//2+1),int(self.pad0//2)) if self.pad0%2!=0  else (self.pad0//2,self.pad0//2)\n        \n            tmp=np.pad(tmp,( pad,(0,0),(0,0)))\n        if tmp.shape[1]<self.sz:\n            #print(int(math.ceil(self.pad1/2)) )\n            pad=(int(self.pad1//2+1),int(self.pad1//2)) if self.pad1%2!=0 else (self.pad1//2,self.pad1//2)\n            #print(pad)\n            tmp=np.pad(tmp,((0,0), pad,(0,0)))\n        return tmp\n    \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        #print('no,n1',n0,n1) #indexer moving from row to row\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        x0,y0 =  n0*self.sz,  n1*self.sz\n        #print('xo,y0',x0,y0,'shape',self.shape[0],self.shape[1])\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        #print('ps',p00,p01,p10,p11)\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            if np.moveaxis(self.data.read([1,2,3],\n                window=Window.from_slices((p00,p01),(p10,p11))), 0, -1).shape !=(self.sz,self.sz,3):\n                \n                tmp=np.moveaxis(self.data.read([1,2,3],\n                window=Window.from_slices((p00,p01),(p10,p11))), 0, -1)\n                \n                print('tmp',tmp.shape)\n               \n                img=self.pad_tile(tmp)\n                print('shape not eq',img.shape)\n            else:\n                \n                img = np.moveaxis(self.data.read([1,2,3],\n                    window=Window.from_slices((p00,p01),(p10,p11))), 0, -1)\n            #print((p00 ),(p01 ),(p10 ),(p11 ) \n            #       )\n            \n        else:\n            for i,layer in enumerate(self.layers):\n                if  layer.read(1,window=Window.from_slices((p00,p01),(p10,p11))).shape !=(self.sz,self.sz,3):\n                \n                    tmp=layer.read(1,window=Window.from_slices((p00,p01),(p10,p11)))\n                    img=self.pad_tile(tmp)\n                #print('tmp',tmp.shape)\n               \n                 \n                    print('shape not eq',img.shape)\n                else:\n                    \n                    img  =  layer.read(1,window=Window.from_slices((p00,p01),(p10,p11)))\n        if self.mask is not None: \n            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            if self.mask is not None: \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        if self.mask is not None: \n            \n            return img, mask, (-1 if (s>s_th).sum() <= p_th or img.sum() <= p_th else idx)\n        else:\n            return img, None, (-1 if (s>s_th).sum() <= p_th or img.sum() <= p_th else idx)\n            ","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:57:17.354699Z","iopub.execute_input":"2022-08-09T11:57:17.355151Z","iopub.status.idle":"2022-08-09T11:57:17.396495Z","shell.execute_reply.started":"2022-08-09T11:57:17.355114Z","shell.execute_reply":"2022-08-09T11:57:17.39526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with_mask=False\nx_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=None,ext='tif')\n        print('ds length',len(ds))\n         \n        for i in range(len(ds)):\n            im,m,idx = ds[i]\n            #break\n            if idx < 0: \n                #print('idx',idx)\n                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            \n            if with_mask:\n                m = cv2.imencode('.png',m)[1]\n                \n                mask_out.writestr(f'{index}_{idx:04d}.png', m)\n        #break\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":{"execution":{"iopub.status.busy":"2022-08-09T11:57:21.340163Z","iopub.execute_input":"2022-08-09T11:57:21.340565Z","iopub.status.idle":"2022-08-09T11:58:29.568463Z","shell.execute_reply.started":"2022-08-09T11:57:21.340529Z","shell.execute_reply":"2022-08-09T11:58:29.567141Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"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            if with_mask:\n                \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            if with_mask:\n                \n                plt.imshow(Image.fromarray(mask), alpha=0.2)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-09T11:59:35.190778Z","iopub.execute_input":"2022-08-09T11:59:35.191162Z","iopub.status.idle":"2022-08-09T11:59:36.23563Z","shell.execute_reply.started":"2022-08-09T11:59:35.191126Z","shell.execute_reply":"2022-08-09T11:59:36.233769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}