{"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":"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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"orig = 1024\nsz = 512 #128 #256 #the size of tiles\nreduce = orig//sz  #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*(sz//256)**2 #threshold for the minimum number of pixels\n\ntop_n = 5 # 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":"code","source":"def 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('id')\ndf_masks.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class HuBMAPDataset(Dataset):\n    def __init__(self, idx, sz=sz, reduce=reduce, encs=None):\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('tmp.tiff', tiff.imread(os.path.join(DATA,idx+'.tiff')), photometric='rgb')\n            self.data = rasterio.open('tmp.tiff',num_threads='all_cpus')\n            gc.collect()\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 + 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        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\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: 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":{"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":"code","source":"os.makedirs('train')\n\nx_tot,x2_tot = [],[]\nshift = 0 \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)\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            \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    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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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\nmini_size = 64\ndef _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\ndef load_dataset(filenames):\n    dataset = tf.data.TFRecordDataset(filenames)\n    dataset = dataset.map(lambda ex: _parse_image_function(ex))\n    return dataset\n\nN = 8\ndef get_dataset(FILENAME):\n    dataset = load_dataset(FILENAME)\n    dataset = dataset.batch(N*N)\n    return dataset","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nfrom skimage.segmentation import mark_boundaries\nfor imgs, masks in get_dataset(train_images[0]).take(1):\n    pass\n\nplt.figure(figsize = (N,N))\ngs1 = gridspec.GridSpec(N,N)\n\nfor 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\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"os.makedirs('train2')\n\nx_tot,x2_tot = [],[]\nshift = orig\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)\n\n    filename = 'train2/'+ 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            \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,'train2/'+ index + '-'+str(cnt) +'.tfrec')\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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train2_images = glob.glob('train2/*.tfrec')\nctrain2i = count_data_items(train2_images)\nprint(f'Num train2 images: {ctrain2i}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for imgs, masks in get_dataset(train2_images[0]).take(1):\n    pass\nplt.figure(figsize = (N,N))\nfor 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\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"p = pathlib.Path('../input/hubmap-kidney-segmentation')\nidentity = rasterio.Affine(1, 0, 0, 0, 1, 0)\nos.makedirs('test', exist_ok = True)\n\nfor 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":"code","source":"test_images = glob.glob('test/*.tfrec')\nctesti = count_data_items(test_images)\nprint(f'Num test images: {ctesti}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DIM = sz\nmini_size = 64\ndef _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\ndef load_dataset(filenames):\n    dataset = tf.data.TFRecordDataset(filenames)\n    dataset = dataset.map(lambda ex: _parse_image_function(ex))\n    return dataset\n\nN = 8\ndef get_dataset(FILENAME):\n    dataset = load_dataset(FILENAME)\n    dataset = dataset.batch(N*N)\n    return dataset","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\nplt.figure(figsize = (N,N))\ngs1 = gridspec.GridSpec(N,N)\n\nfor 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\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}