{"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":"## Refferences:\n* @iafoss https://www.kaggle.com/iafoss/256x256-images (many thanks - huge part of a code presented below is COPIED from this notebook, kindly please upvote original notebook and dataset)\n* @cdeotte https://www.kaggle.com/cdeotte/how-to-create-tfrecords\n* https://www.tensorflow.org/tutorials/load_data/tfrecord\n* @leighplt https://www.kaggle.com/leighplt/pytorch-fcn-resnet50 (another tiling idea with make_grid() and rasterio - useful for inference)\n\n## Table of contents:\n* Introduction\n* Init\n    * Packages\n    * Parameters\n    * Functions\n* Train\n* Train2\n* Inference\n\n# Introduction\nA lot of users asks for this notebook so here it is.\nThis notebooks presents how to create TFrecords for next steps (training and inference) in HuBMAP - Hacking the Kidney competition. So far I've created datatasets with different tilings listed (for train and train2) below:\n* [1024-512](https://www.kaggle.com/wrrosa/hubmap-tfrecords-1024-512)\n* [1024-256](https://www.kaggle.com/wrrosa/hubmap-tfrecords-1024-256)\n* [1024-128](https://www.kaggle.com/wrrosa/hubmap-tfrecords-1024-128)\n* [1024-64](https://www.kaggle.com/wrrosa/hubmap-tfrecords-1024-64)\n* [1536-512](https://www.kaggle.com/wrrosa/hubmap-tfrecords-1536-512)\n* [768-384](https://www.kaggle.com/wrrosa/hubmap-tfrecords-768-384)\n\n\n## Tensorflow HuBMAP competition starter kit:\n* this notebook (how to create training and inference tfrecords)\n* https://www.kaggle.com/wrrosa/hubmap-tf-with-tpu-efficientunet-512x512-train (training pipeline)\n* https://www.kaggle.com/wrrosa/hubmap-tf-with-tpu-efficientunet-512x512-subm (inference with submission)\n\n\n# Init\n## Packages","metadata":{}},{"cell_type":"code","source":"import 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 tensorflow as tf\nimport gc\nimport rasterio\nfrom rasterio.windows import Window\nfrom torch.utils.data import Dataset","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from kaggle_datasets import KaggleDatasets\n# GCS_DS_PATH_MAIN = KaggleDatasets().get_gcs_path('hubmap-kidney-segmentation')\n# GCS_DS_PATH_MODEL = KaggleDatasets().get_gcs_path('hubmap-fold-model')\n# print(GCS_DS_PATH_MAIN, GCS_DS_PATH_MODEL)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Parameters","metadata":{}},{"cell_type":"code","source":"patch_size = 1024\npad_offset = 64\nreduce = 2  #reduce the original images by 'reduce' times \nMASKS = '../input/hubmap-kidney-segmentation/train.csv'\nDATA = '../input/hubmap-kidney-segmentation/train/'\ns_th = 40  #saturation blancking threshold\np_th = 1000*((patch_size//reduce)//512)**2 #threshold for the minimum number of pixels\nSHIFT = 257\ntop_n = 25 # only first 5 tiff files for train, train2 and test will be processed due to output 20gb limit","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Functions","metadata":{}},{"cell_type":"code","source":"#functions to convert encoding to mask and mask to encoding\nimport pandas as pd\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])\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('id')\ndf_masks.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(df_masks)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!mkdir /kaggle/tmp","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Thank you @iafoss\nimport os\nclass HuBMAPDataset(Dataset):\n    def __init__(self, idx, patch_size=patch_size, reduce=reduce, encs=None, shift=0, pad_offset=pad_offset):\n        self.data = rasterio.open(os.path.join(DATA,idx+'.tiff'),num_threads='all_cpus')\n        # some images have issues with format \n        # and must be saved correctly before reading with rasterio\n        if self.data.count == 1:\n            tiff.imwrite('/kaggle/tmp/tmp.tiff', tiff.imread(os.path.join(DATA,idx+'.tiff')), photometric='rgb')\n            self.data = rasterio.open('/kaggle/tmp/tmp.tiff',num_threads='all_cpus')\n            gc.collect()\n        self.shape = self.data.shape\n        self.reduce = reduce\n        self.sz = patch_size\n        self.pad_offset = pad_offset  \n        self.pad0 = (self.sz - self.shape[0]%self.sz)%self.sz + shift\n        self.pad1 = (self.sz - self.shape[1]%self.sz)%self.sz + shift\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        \n\n        img = np.zeros((self.sz +self.pad_offset*2,self.sz + self.pad_offset*2,3),np.uint8)\n        mask = np.zeros((self.sz +self.pad_offset*2,self.sz + self.pad_offset*2),np.uint8)\n        \n        x0,y0 = -self.pad0//2 + n0*self.sz - self.pad_offset, -self.pad1//2 + n1*self.sz - self.pad_offset\n        # make sure that the region to read is within the image\n        p00,p01 = max(0,x0 ), min(\n            x0 + self.sz + self.pad_offset*2 ,self.shape[0])\n        p10,p11 = max(0,y0 ), min(\n            y0 + self.sz + self.pad_offset*2, self.shape[1])\n        # mapping the loade region to the tile\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 self.mask is not None: \n#             x0,y0 = -self.pad0//2 + n0*self.sz, -self.pad1//2 + n1*self.sz\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            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 + self.pad_offset*2)//reduce,(self.sz + self.pad_offset*2)//reduce),\n                             interpolation = cv2.INTER_AREA)\n            mask = cv2.resize(mask,((self.sz + self.pad_offset*2)//reduce,(self.sz + self.pad_offset*2)//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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The following function can be used to convert a value to a type compatible\n# with tf.train.Example.\n\ndef _bytes_feature(value):\n  \"\"\"Returns a bytes_list from a string / byte.\"\"\"\n  if isinstance(value, type(tf.constant(0))):\n    value = value.numpy() # BytesList won't unpack a string from an EagerTensor.\n  return tf.train.Feature(bytes_list=tf.train.BytesList(value=[value]))\n\ndef serialize_example(image, mask):\n  \"\"\"\n  Creates a tf.train.Example message ready to be written to a file.\n  \"\"\"\n  # Create a dictionary mapping the feature name to the tf.train.Example-compatible\n  # data type.\n  feature = {\n      'image': _bytes_feature(image),\n      'mask': _bytes_feature(mask),\n  }\n\n  # Create a Features message using tf.train.Example.\n\n  example_proto = tf.train.Example(features=tf.train.Features(feature=feature))\n  return example_proto.SerializeToString()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train","metadata":{}},{"cell_type":"code","source":"os.makedirs('train', exist_ok=True)\n\nx_tot,x2_tot = [],[]\n\nfor index, encs in tqdm(df_masks.head(top_n).iterrows()):\n    print(index)\n    #read image and generate the mask\n    ds = HuBMAPDataset(index,encs=encs, shift=SHIFT)\n\n    filename = 'train/'+ index + '.tfrec'\n    cnt = 0\n    with tf.io.TFRecordWriter(filename) as writer:\n        \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            #write data   \n            im = cv2.cvtColor(im, cv2.COLOR_RGB2BGR)\n            example = serialize_example(im.tobytes(),m.tobytes())\n            writer.write(example)\n            cnt +=1\n            \n    os.rename(filename,'train/'+ index + '-'+str(cnt) +'.tfrec')\n    print('filename done')\n    gc.collect()        \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)\nos.remove('/kaggle/tmp/tmp.tiff')\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check","metadata":{}},{"cell_type":"code","source":"import re\nimport glob\ndef count_data_items(filenames):\n    n = [int(re.compile(r\"-([0-9]*)\\.\").search(filename).group(1)) for filename in filenames]\n    return np.sum(n)\ntrain_images = glob.glob('train/*.tfrec')\nctraini = count_data_items(train_images)\nprint(f'Num train images: {ctraini}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# DIM = sz\n# mini_size = 64\n# def _parse_image_function(example_proto):\n#     image_feature_description = {\n#         'image': tf.io.FixedLenFeature([], tf.string),\n#         'mask': tf.io.FixedLenFeature([], tf.string)\n#     }\n#     single_example = tf.io.parse_single_example(example_proto, image_feature_description)\n#     image = tf.reshape( tf.io.decode_raw(single_example['image'],out_type=np.dtype('uint8')), (DIM,DIM, 3))\n#     mask =  tf.reshape(tf.io.decode_raw(single_example['mask'],out_type='bool'),(DIM,DIM,1))\n    \n#     image = tf.image.resize(image,(mini_size,mini_size))/255.0\n#     mask = tf.image.resize(tf.cast(mask,'uint8'),(mini_size,mini_size))\n#     return image, mask\n\n\n# def load_dataset(filenames):\n#     dataset = tf.data.TFRecordDataset(filenames)\n#     dataset = dataset.map(lambda ex: _parse_image_function(ex))\n#     return dataset\n\n# N = 8\n# def get_dataset(FILENAME):\n#     dataset = load_dataset(FILENAME)\n#     dataset = dataset.batch(N*N)\n#     return dataset\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import matplotlib.pyplot as plt\n# import matplotlib.gridspec as gridspec\n# from skimage.segmentation import mark_boundaries\n# for imgs, masks in get_dataset(train_images[0]).take(1):\n#     pass\n\n# plt.figure(figsize = (N,N))\n# gs1 = gridspec.GridSpec(N,N)\n\n# for i in range(N*N):\n#    # i = i + 1 # grid spec indexes from 0\n#     ax1 = plt.subplot(gs1[i])\n#     plt.axis('on')\n#     ax1.set_xticklabels([])\n#     ax1.set_yticklabels([])\n#     ax1.set_aspect('equal')\n#     ax1.imshow(mark_boundaries(imgs[i], masks[i].numpy().squeeze().astype('bool')))\n\n# plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train 2\nSame idea as train, but shifted by padding.","metadata":{}},{"cell_type":"markdown","source":"## Check","metadata":{}},{"cell_type":"code","source":"# for imgs, masks in get_dataset(train2_images[0]).take(1):\n#     pass\n# plt.figure(figsize = (N,N))\n# for i in range(N*N):\n#    # i = i + 1 # grid spec indexes from 0\n#     ax1 = plt.subplot(gs1[i])\n#     plt.axis('on')\n#     ax1.set_xticklabels([])\n#     ax1.set_yticklabels([])\n#     ax1.set_aspect('equal')\n#     ax1.imshow(mark_boundaries(imgs[i], masks[i].numpy().squeeze().astype('bool')))\n\n# plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Inference\nThe approach presented above is superfast and elegant, but it does not contain the coordinates of the image (x1, y1) so (without modification) it is useless for inference. Now @leighplt https://www.kaggle.com/leighplt/pytorch-fcn-resnet50 presented how to use rasterio - but it doesn't support batching. So idea is to create tfrecords using rasterio and use them in inference - should be faster.","metadata":{}},{"cell_type":"code","source":"WINDOW = orig #1024\nMIN_OVERLAP = 300\nNEW_SIZE = sz #512\n\nimport numpy as np\nimport pandas as pd\nimport os\nimport glob\nimport gc\n\nimport rasterio\nfrom rasterio.windows import Window\n\nimport pathlib\nfrom tqdm.notebook import tqdm\nimport cv2\n\nimport tensorflow as tf\n\ndef make_grid(shape, window=256, min_overlap=32):\n    \"\"\"\n        Return Array of size (N,4), where N - number of tiles,\n        2nd axis represente slices: x1,x2,y1,y2 \n    \"\"\"\n    x, y = shape\n    nx = x // (window - min_overlap) + 1\n    x1 = np.linspace(0, x, num=nx, endpoint=False, dtype=np.int64)\n    x1[-1] = x - window\n    x2 = (x1 + window).clip(0, x)\n    ny = y // (window - min_overlap) + 1\n    y1 = np.linspace(0, y, num=ny, endpoint=False, dtype=np.int64)\n    y1[-1] = y - window\n    y2 = (y1 + window).clip(0, y)\n    slices = np.zeros((nx,ny, 4), dtype=np.int64)\n    \n    for i in range(nx):\n        for j in range(ny):\n            slices[i,j] = x1[i], x2[i], y1[j], y2[j]    \n    return slices.reshape(nx*ny,4)\n\ndef _bytes_feature(value):\n  \"\"\"Returns a bytes_list from a string / byte.\"\"\"\n  if isinstance(value, type(tf.constant(0))):\n    value = value.numpy() # BytesList won't unpack a string from an EagerTensor.\n  return tf.train.Feature(bytes_list=tf.train.BytesList(value=[value]))\n\ndef _int64_feature(value):\n  \"\"\"Returns an int64_list from a bool / enum / int / uint.\"\"\"\n  return tf.train.Feature(int64_list=tf.train.Int64List(value=[value]))\n\ndef serialize_example(image, x1, y1):\n  feature = {\n      'image': _bytes_feature(image),\n      'x1': _int64_feature(x1),\n      'y1': _int64_feature(y1)\n  }\n  example_proto = tf.train.Example(features=tf.train.Features(feature=feature))\n  return example_proto.SerializeToString()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# p = pathlib.Path('../input/hubmap-kidney-segmentation')\n# identity = rasterio.Affine(1, 0, 0, 0, 1, 0)\n# os.makedirs('test', exist_ok = True)\n\n# for i, filename in tqdm(enumerate(p.glob('test/*.tiff')), \n#                         total = len(list(p.glob('test/*.tiff')))):\n    \n#     print(f'{i+1} Creating tfrecords for image: {filename.stem}')\n#     dataset = rasterio.open(filename.as_posix(), transform = identity)\n#     slices = make_grid(dataset.shape, window=WINDOW, min_overlap=MIN_OVERLAP)\n    \n#     if dataset.count != 3:\n#         layers = [rasterio.open(subd) for subd in dataset.subdatasets]\n    \n#     print(slices.shape[0])\n#     cnt = 0\n#     part = 0 \n#     fname = f'test/{filename.stem}-part{part}.tfrec'\n#     writer = tf.io.TFRecordWriter(fname) \n#     for (x1,x2,y1,y2) in slices:\n#         if cnt>999:\n#             writer.close()\n#             os.rename(fname, f'test/{filename.stem}-part{part}-{cnt}.tfrec')\n#             part += 1\n#             fname = f'test/{filename.stem}-part{part}.tfrec'\n#             writer = tf.io.TFRecordWriter(fname)\n#             cnt = 0\n        \n#         if dataset.count == 3:\n#             image = dataset.read([1,2,3],\n#                         window=Window.from_slices((x1,x2),(y1,y2)))\n#             image = np.moveaxis(image, 0, -1)\n#         else:\n#             image = np.zeros((WINDOW, WINDOW, 3), dtype=np.uint8)\n#             for fl in range(3):\n#                 image[:,:,fl] = layers[fl].read(window=Window.from_slices((x1,x2),(y1,y2)))\n                \n#         image = cv2.resize(image, (NEW_SIZE, NEW_SIZE),interpolation = cv2.INTER_AREA)\n#         image = cv2.cvtColor(image, cv2.COLOR_RGB2BGR)\n#         example = serialize_example(image.tobytes(),x1,y1)\n#         writer.write(example)\n#         cnt+=1\n#     writer.close()\n#     del writer\n#     os.rename(fname, f'test/{filename.stem}-part{part}-{cnt}.tfrec')\n#     gc.collect();","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check","metadata":{}},{"cell_type":"code","source":"# test_images = glob.glob('test/*.tfrec')\n# ctesti = count_data_items(test_images)\n# print(f'Num test images: {ctesti}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# DIM = sz\n# mini_size = 64\n# def _parse_image_function(example_proto):\n#     image_feature_description = {\n#         'image': tf.io.FixedLenFeature([], tf.string),\n#         'x1': tf.io.FixedLenFeature([], tf.int64),\n#         'y1': tf.io.FixedLenFeature([], tf.int64)\n#     }\n#     single_example = tf.io.parse_single_example(example_proto, image_feature_description)\n#     image = tf.reshape( tf.io.decode_raw(single_example['image'],out_type=np.dtype('uint8')), (DIM,DIM, 3))\n#     x1 = single_example['x1']\n#     y1 = single_example['y1']\n#     image = tf.image.resize(image,(mini_size,mini_size))/255.0\n#     return image, x1, y1\n\n\n# def load_dataset(filenames):\n#     dataset = tf.data.TFRecordDataset(filenames)\n#     dataset = dataset.map(lambda ex: _parse_image_function(ex))\n#     return dataset\n\n# N = 8\n# def get_dataset(FILENAME):\n#     dataset = load_dataset(FILENAME)\n#     dataset = dataset.batch(N*N)\n#     return dataset\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# for imgs, x1, y1 in get_dataset(test_images[1]).take(2):\n#     pass\n\n# plt.figure(figsize = (N,N))\n# gs1 = gridspec.GridSpec(N,N)\n\n# for i in range(N*N):\n#    # i = i + 1 # grid spec indexes from 0\n#     ax1 = plt.subplot(gs1[i])\n#     plt.axis('on')\n#     ax1.set_xticklabels([])\n#     ax1.set_yticklabels([])\n#     ax1.set_aspect('equal')\n#     ax1.set_title(f'{x1[i]}; {y1[i]}', fontsize=6)\n#     ax1.imshow(imgs[i])\n\n# plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}