{"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":"## 768×768 dataset with Stain Normalization by StainNet\n\n combination of [Converting to 256x256](https://www.kaggle.com/code/thedevastator/converting-to-256x256) and [Stain Normalization: StainGAN, StainNet](https://www.kaggle.com/code/nghihuynh/stain-normalization-staingan-stainnet)\n \n if you want a dataset for multi-segmentation, set multi_seg = True.","metadata":{}},{"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":2.622755,"end_time":"2022-07-04T13:11:37.002731","exception":false,"start_time":"2022-07-04T13:11:34.379976","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-08-18T11:55:09.013717Z","iopub.execute_input":"2022-08-18T11:55:09.014737Z","iopub.status.idle":"2022-08-18T11:55:11.445819Z","shell.execute_reply.started":"2022-08-18T11:55:09.014681Z","shell.execute_reply":"2022-08-18T11:55:11.444604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"OUT_TRAIN = 'train.zip'\nOUT_MASKS = 'masks.zip'\nsz = 768   # the size of tiles\nreduce = 1 # reduce the original images by 4 times \nmulti_seg = False # if True, creat a dataset for multi segmentation \nMASKS = '../input/hubmap-organ-segmentation/train.csv'\nDATA = '../input/hubmap-organ-segmentation/train_images'","metadata":{"papermill":{"duration":0.011739,"end_time":"2022-07-04T13:11:37.018315","exception":false,"start_time":"2022-07-04T13:11:37.006576","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-08-18T11:55:11.447684Z","iopub.execute_input":"2022-08-18T11:55:11.448255Z","iopub.status.idle":"2022-08-18T11:55:11.456097Z","shell.execute_reply.started":"2022-08-18T11:55:11.448216Z","shell.execute_reply":"2022-08-18T11:55:11.452833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# functions to convert encoding to mask and mask to encoding\ndef enc2mask(mask_rle, shape):\n    img = np.zeros(shape[0]*shape[1], dtype=np.uint8)\n    s = mask_rle.split()\n    starts, lengths = [np.asarray(x, dtype=int) for x in (s[0::2], s[1::2])]\n    starts -= 1\n    ends = starts + lengths\n    for lo, hi in zip(starts, ends):\n        img[lo:hi] = 1\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\nmask_map = dict(\n    kidney=1,\n    prostate=2,\n    largeintestine=3,\n    spleen=4,\n    lung=5)#add\n\ndf_masks = pd.read_csv(MASKS)[['id', 'organ', 'rle']].set_index('id')#add\ndf_masks.head()","metadata":{"papermill":{"duration":0.342618,"end_time":"2022-07-04T13:11:37.363696","exception":false,"start_time":"2022-07-04T13:11:37.021078","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-08-18T11:55:11.457350Z","iopub.execute_input":"2022-08-18T11:55:11.459484Z","iopub.status.idle":"2022-08-18T11:55:11.777496Z","shell.execute_reply.started":"2022-08-18T11:55:11.459447Z","shell.execute_reply":"2022-08-18T11:55:11.776444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s_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        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        #print(self.shape, self.reduce, self.sz, self.pad0, self.pad1, self.n0max, self.n1max)\n        \n    def __len__(self):\n        return self.n0max*self.n1max\n    \n    def __getitem__(self, idx):\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        \n    \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.02608,"end_time":"2022-07-04T13:11:37.39294","exception":false,"start_time":"2022-07-04T13:11:37.36686","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-08-18T11:55:11.779659Z","iopub.execute_input":"2022-08-18T11:55:11.780045Z","iopub.status.idle":"2022-08-18T11:55:11.798854Z","shell.execute_reply.started":"2022-08-18T11:55:11.780006Z","shell.execute_reply":"2022-08-18T11:55:11.797670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sys\nsys.path.append('../input/stainnet/')\nfrom models import StainNet, ResnetGenerator\nimport torch\n\ndef norm(image):\n    image = np.array(image).astype(np.float32)\n    image = image.transpose((2, 0, 1))\n    image = ((image / 255) - 0.5) / 0.5\n    image=image[np.newaxis, ...]\n    image=torch.from_numpy(image)\n    return image\n\ndef un_norm(image):\n    image = image.cpu().detach().numpy()[0]\n    image = ((image * 0.5 + 0.5) * 255).astype(np.uint8).transpose((1,2,0))\n    return image\n\nmodel_Net = StainNet().cuda()\nmodel_Net.load_state_dict(torch.load(\"../input/stainnet/checkpoints/aligned_histopathology_dataset/StainNet-Public_layer3_ch32.pth\"))\nmodel_Net.eval()\n\ndef stain_normalize(source, verbose=False):\n    with torch.no_grad():\n        img_net=model_Net(norm(source).cuda())\n        img_net=un_norm(img_net)\n        if verbose: plt.imshow(img_net); plt.show()\n        return img_net","metadata":{"execution":{"iopub.status.busy":"2022-08-18T11:58:28.698053Z","iopub.execute_input":"2022-08-18T11:58:28.698591Z","iopub.status.idle":"2022-08-18T11:58:28.724667Z","shell.execute_reply.started":"2022-08-18T11:58:28.698551Z","shell.execute_reply":"2022-08-18T11:58:28.722764Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x_tot,x2_tot = [],[]\nmask_ = [] #add\nwith zipfile.ZipFile(OUT_TRAIN, 'w') as img_out,\\\n zipfile.ZipFile(OUT_MASKS, 'w') as mask_out:\n    for index,(organ, encs) in tqdm(df_masks.iterrows(),total=len(df_masks)):#add\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            im = stain_normalize(im, verbose=False)\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            if multi_seg:\n                m = mask_map[organ]*m #add\n            for e in np.unique(m): #add\n                mask_.append(e) #add\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)\nprint('classes of mask', set(mask_)) #add","metadata":{"papermill":{"duration":187.425795,"end_time":"2022-07-04T13:14:44.822066","exception":false,"start_time":"2022-07-04T13:11:37.396271","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-08-18T12:05:51.775862Z","iopub.execute_input":"2022-08-18T12:05:51.776275Z","iopub.status.idle":"2022-08-18T12:05:56.172621Z","shell.execute_reply.started":"2022-08-18T12:05:51.776236Z","shell.execute_reply":"2022-08-18T12:05:56.171446Z"},"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.3)\nplt.show()","metadata":{"papermill":{"duration":1.729573,"end_time":"2022-07-04T13:14:46.555048","exception":false,"start_time":"2022-07-04T13:14:44.825475","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-08-18T11:41:01.134946Z","iopub.execute_input":"2022-08-18T11:41:01.135352Z","iopub.status.idle":"2022-08-18T11:41:05.559337Z","shell.execute_reply.started":"2022-08-18T11:41:01.135300Z","shell.execute_reply":"2022-08-18T11:41:05.558235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"papermill":{"duration":0.020686,"end_time":"2022-07-04T13:14:46.597035","exception":false,"start_time":"2022-07-04T13:14:46.576349","status":"completed"},"tags":[]},"execution_count":null,"outputs":[]}]}