{"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":"
[WSI-HuBMAP] Dataset creation 256x256
\n\n\n\n
Created by: Nghi Huynh
","metadata":{}},{"cell_type":"markdown","source":"#
1. Imports
","metadata":{}},{"cell_type":"code","source":"import gc\nimport os\nimport cv2\nimport zipfile\nimport rasterio\nimport numpy as np\nimport pandas as pd\nimport tifffile as tiff\nfrom PIL import Image\nimport matplotlib.pyplot as plt\nfrom rasterio.windows import Window\nfrom torch.utils.data import Dataset\nfrom tqdm.notebook import tqdm\nimport math\nimport re","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-09-20T06:07:18.875907Z","iopub.execute_input":"2022-09-20T06:07:18.876549Z","iopub.status.idle":"2022-09-20T06:07:18.886291Z","shell.execute_reply.started":"2022-09-20T06:07:18.876497Z","shell.execute_reply":"2022-09-20T06:07:18.884846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#
2. Configs
","metadata":{}},{"cell_type":"code","source":"config = {\n 'resize': (768,768),\n 'resolution': (1024,1024),\n 'DATA': '../input/hubmap-organ-segmentation/train_images',\n 'MASKS' : '../input/hubmap-organ-segmentation/train.csv',\n 'Window' : (250,2670),\n 'bs': 64,\n 'nfolds': 4,\n 'fold': 0,\n 'NUM_WORKERS': 4,\n 'OUT_TRAIN' : 'train.zip',\n 'OUT_MASKS' : 'masks.zip'\n}","metadata":{"execution":{"iopub.status.busy":"2022-09-20T05:51:37.574578Z","iopub.execute_input":"2022-09-20T05:51:37.575451Z","iopub.status.idle":"2022-09-20T05:51:37.581630Z","shell.execute_reply.started":"2022-09-20T05:51:37.575407Z","shell.execute_reply":"2022-09-20T05:51:37.580745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#
3. Helper functions
","metadata":{}},{"cell_type":"code","source":"# functions to convert encoding to mask and mask to encoding\ndef enc2mask(encs, shape):\n '''\n Args:\n encs: list of rle masks\n shape: mask shape\n '''\n img = np.zeros(shape[0]*shape[1], dtype=np.uint8)\n for m,enc in enumerate(encs):\n if isinstance(enc,np.float32) 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","metadata":{"execution":{"iopub.status.busy":"2022-09-20T05:51:39.953461Z","iopub.execute_input":"2022-09-20T05:51:39.953902Z","iopub.status.idle":"2022-09-20T05:51:39.968421Z","shell.execute_reply.started":"2022-09-20T05:51:39.953863Z","shell.execute_reply":"2022-09-20T05:51:39.966716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def visualize_batch(df, organ, slide):\n count = 0\n fig = plt.figure(figsize=(30,30))\n for index, row in df.iterrows():\n if count == 25: break\n \n if row['organ'] == organ:\n img = rasterio.open(os.path.join(config['DATA'],str(index)+'.tiff'), num_threads='all_cpus')\n mask = enc2mask([row['rle']],(img.shape[1],img.shape[0]))\n \n #read window slide\n img = img.read([1,2,3],window=Window.from_slices(slide, slide)) \n mask = mask[slide[0]:slide[1], slide[0]:slide[1]] \n \n plt.subplot(5,5, count+1)\n plt.imshow(np.transpose(img, (1,2,0)))\n plt.imshow(mask, cmap='seismic', alpha=0.4)\n plt.title(f'{index}', size=40)\n plt.axis('off')\n \n count += 1\n del img, mask\n \n plt.suptitle(f'{organ}', size=50, weight=\"bold\", y = 1.0) \n plt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-20T05:58:52.443966Z","iopub.execute_input":"2022-09-20T05:58:52.444436Z","iopub.status.idle":"2022-09-20T05:58:52.456766Z","shell.execute_reply.started":"2022-09-20T05:58:52.444391Z","shell.execute_reply":"2022-09-20T05:58:52.455733Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_masks = pd.read_csv(config['MASKS'])[['id', 'rle','organ']].set_index('id')\ndf_masks.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-20T05:51:45.211509Z","iopub.execute_input":"2022-09-20T05:51:45.212013Z","iopub.status.idle":"2022-09-20T05:51:45.722799Z","shell.execute_reply.started":"2022-09-20T05:51:45.211974Z","shell.execute_reply":"2022-09-20T05:51:45.721543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#
4. Create HPA dataset class
\n\nBased on the information provided by the host about tissue area within HPA image around 2500 x 2500 pixels, my approach is to create whole image dataset contains only tissue area. \n\n* [256x256](https://www.kaggle.com/datasets/nghihuynh/hubmap-2022-wsi-256x256)\n* [512x512](https://www.kaggle.com/datasets/nghihuynh/hubmap-2022-wsi-512x512)\n* [768x768](https://www.kaggle.com/datasets/nghihuynh/hubmap-2022-wsi-768x768)\n\n`Window.from_slices((row_start, row_stop), (col_start, col_stop))`","metadata":{}},{"cell_type":"code","source":"class HPADataset(Dataset):\n def __init__(self, idx, resize=config['resize'], slide=config['Window'], encs=None):\n self.data = rasterio.open(os.path.join(config['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 \n self.shape = self.data.shape\n self.slide = slide\n self.resize = resize\n self.idx = idx\n \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 1\n \n def __getitem__(self, idx):\n # read img (RGB), mask (grayscale) from window slide\n # img, mask: uint8\n img = self.data.read([1,2,3],window=Window.from_slices(self.slide,self.slide)) \n mask = self.mask[self.slide[0]:self.slide[1], self.slide[0]:self.slide[1]] \n \n # resize\n img = cv2.resize(np.transpose(img,(1,2,0)),(self.resize[0],self.resize[1]),\n interpolation = cv2.INTER_AREA)\n mask = cv2.resize(mask,(self.resize[0], self.resize[1]),\n interpolation = cv2.INTER_NEAREST)\n \n return img, mask, self.idx","metadata":{"execution":{"iopub.status.busy":"2022-09-20T05:51:48.503796Z","iopub.execute_input":"2022-09-20T05:51:48.504822Z","iopub.status.idle":"2022-09-20T05:51:48.518715Z","shell.execute_reply.started":"2022-09-20T05:51:48.504752Z","shell.execute_reply":"2022-09-20T05:51:48.517496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#
5. Visualize batches
\n\nVisualize batch of images based on organ with cropped tissue area.","metadata":{}},{"cell_type":"code","source":"visualize_batch(df_masks, 'lung', config['Window'])","metadata":{"execution":{"iopub.status.busy":"2022-09-20T05:59:36.688884Z","iopub.execute_input":"2022-09-20T05:59:36.689560Z","iopub.status.idle":"2022-09-20T06:00:12.571482Z","shell.execute_reply.started":"2022-09-20T05:59:36.689493Z","shell.execute_reply":"2022-09-20T06:00:12.569889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#
6. Dataset creation
","metadata":{}},{"cell_type":"code","source":"x_tot,x2_tot = [],[]\nwith zipfile.ZipFile(config['OUT_TRAIN'], 'w') as img_out,\\\n zipfile.ZipFile(config['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 = HPADataset(index,encs=encs)\n if index == 31800: continue\n for i in range(len(ds)):\n img, m, idx = ds[i]\n\n x_tot.append((img/255.0).reshape(-1,3).mean(0))\n x2_tot.append(((img/255.0)**2).reshape(-1,3).mean(0))\n\n #write data \n img = cv2.imencode('.png',cv2.cvtColor(img, cv2.COLOR_RGB2BGR))[1]\n img_out.writestr(f'{index}.png', img)\n m = cv2.imencode('.png',m)[1]\n mask_out.writestr(f'{index}.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":{"execution":{"iopub.status.busy":"2022-09-20T06:02:37.886130Z","iopub.execute_input":"2022-09-20T06:02:37.886633Z","iopub.status.idle":"2022-09-20T06:04:56.861289Z","shell.execute_reply.started":"2022-09-20T06:02:37.886598Z","shell.execute_reply":"2022-09-20T06:04:56.859397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#
7. Sanity check
","metadata":{}},{"cell_type":"code","source":"cols, rows = 4,4\nidx0 = 20\nfig=plt.figure(figsize=(cols*4, rows*4))\nwith zipfile.ZipFile(config['OUT_TRAIN'], 'r') as img_arch, \\\n zipfile.ZipFile(config['OUT_MASKS'], 'r') as msk_arch:\n fnames = sorted(img_arch.namelist())[8:]\n #print(fnames)\n for i in range(rows):\n for j in range(cols):\n idx = i+j*cols\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 fig.add_subplot(rows, cols, idx+1)\n plt.axis('off')\n plt.imshow(Image.fromarray(img))\n plt.imshow(Image.fromarray(mask), alpha=0.2)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-20T06:07:23.410712Z","iopub.execute_input":"2022-09-20T06:07:23.412251Z","iopub.status.idle":"2022-09-20T06:07:27.219366Z","shell.execute_reply.started":"2022-09-20T06:07:23.412197Z","shell.execute_reply":"2022-09-20T06:07:27.218406Z"},"trusted":true},"execution_count":null,"outputs":[]}]}