{"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":"# 1. Description\nThere are two neural nets are used:  \n1) small scale, that detects glomeruli centers  \n2) big scale, that detects glomerulus edge  \n  \nDetails:  \n1) Small scale net is fully convolutional net without decoder, arhitecture:  \nc pcpc pcpc cccc  \nwhere c is convolution layer  \np - 2x2 max-polling layer  \nspace devide net in 4 block  \nin last 3 blocks pcpc pcpc cccc, each block has skip connection  \nThe skip connection is 2x2 mean-pooling and 1x1 convolution layer  \nLast skip connection for the last cccc block is only 1x1 convolution layer.  \nOut put has 1/16 resolution relative to input because there are 4 poolings in the net.  \nNet's effective receptive field has size that bigger than maximal size of glomeruli.  \nInput is original RGB image with resolution decreased with 0.25 factor.  \nOutput is 1-channel image where only centeres of mass of glomeruli.  \nEach center has gaussian shape with small sigma.  \nHere is a picture of the real predicition:  ","metadata":{}},{"cell_type":"code","source":"# code to draw pictures:\nimport matplotlib.pyplot as plt\nimport cv2\ndef imshow_from_file(file):\n    image_bgr = cv2.imread(file)\n    image = cv2.cvtColor(image_bgr, cv2.COLOR_BGR2RGB)\n    plt.figure(figsize=(25, 25))\n    plt.imshow(image)\n    plt.axis('off')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"imshow_from_file('../input/pics-j/small_scale.png')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"2) big scale net is convolutional net with arhitecture:  \nc pcpc pc pccuccc  \nu is x2 upsapling layer\nOut put has 1/8 resolution relative to input  \n1st net predicts center of glomerulus. Than it makes crop around this positions.  \nAster that crop is unrolled:  \n","metadata":{}},{"cell_type":"code","source":"imshow_from_file('../input/pics-j/big_scale.png')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"then unrolled image is go to the neural net and it predicts radius of the edge at each angle.  \nunrolled image is like cartesian to polar coordinates transformation. It done with cv2.remap function.  \nThe net's kernel sizes in convolution layers is organized in a way that the net can be considered as fully convolutional in angle direction and net with dence layer in radius direction.  \nThis makes translation invariation in angle direction.  \nThe input images is tiled 3 times in angle direction to emulate periodic condition.\nOtherwise it will have edge side-effects that can degrade prediction quality.  \nOutput shape is 1 x 1 x n_angles x 1 (batch_size x height x width x n_channels)  \nSo it predicts radii as regression.  \n  \nAdvatnages of this 2-nets method:  \n1) fast preidiction (about 4 minutes for 1 tiff image with GPU)  \n2) perspective method that can make high dice score (I started 2 month ago, not enough time for research, so only about 0.9 dice metrics was reached)  \n3) 1st net can be used to count number of glomeruli directly  \n4) pixels value in prediction of the 1st net can be used as probability of the detection\n","metadata":{}},{"cell_type":"markdown","source":"# 2. CODE","metadata":{}},{"cell_type":"code","source":"#!/usr/bin/env python3\n# -*- coding: utf-8 -*-\n\"\"\"\n\n@author: vedenev\n\"\"\"\n\n\n\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport cv2\nimport csv\nimport os\nimport gc\nimport rasterio\nimport tensorflow as tf\nimport glob\nimport json\nfrom rasterio.windows import Window\nfrom skimage.segmentation import watershed\n\n\n\nIMAGES_DIR = '../input/hubmap-kidney-segmentation/test'\nWEIGHTS_SMALL = '../input/weights-j/t198_best_model.ckpt'\nWEIGHTS_BIG = '../input/weights-j/t246_best_model.ckpt'\n\n\nOUTPUT_FILE = 'submission.csv'\n\n\nTHRESHOLD_SMALL = 0.6\n\nRESIZE_FACTOR_SMALL = 0.25\nROUNDING_ORIG_SMALL = 64\nOVERLAP_SMALL = 448\nBASE_SIZE_CROP_SMALL = 4096\nNET_SCALE_SMALL = 16\nN_INPUT_CHANNELS_SMALL = 3\nN_CLASSES_SMALL = 1\nN_FEATUREMAPS_SMALL = [N_INPUT_CHANNELS_SMALL, 4, 4, 8, 16, 24, 32, 32, 48, N_CLASSES_SMALL]\nKERNELS_X_SIZE_SMALL = 3\nKERNELS_Y_SIZE_SMALL = 3\nKERNELS_INCREASED_X_SIZE_SMALL = 7\nKERNELS_INCREASED_Y_SIZE_SMALL = 7\n\nBASE_SIZE_BIG = 900\nBASE_SIZE_HALF_BIG = BASE_SIZE_BIG // 2\nNET_SCALE_BIG = 8\n\nCENTER = (BASE_SIZE_BIG - 1) / 2\n\nN_PHI = 384 \nN_PHI_TILED = 3 * N_PHI\nN_PHI_OUT = N_PHI // NET_SCALE_BIG # number of angles in 2nd big scale net at output\nN_PHI_OUT_TILLED = 3 * N_PHI_OUT\nN_R = 128\nN_R_OUT = N_R // NET_SCALE_BIG\nR_MAX = BASE_SIZE_BIG / 2\nr_tmp_index = np.arange(N_R).astype(np.float32).reshape(N_R, 1)\nr_mean_max = (N_R - 1 + 0.5) / NET_SCALE_BIG - 0.5\n\nNET_SCALE_BIG = 8\nN_INPUT_CHANNELS_BIG = 3\nN_CLASSES_BIG = 1\nN_FEATUREMAPS_BIG = [N_INPUT_CHANNELS_BIG, 4, 4, 8, 16, 16, 32, 32, 32, N_CLASSES_BIG]\nKERNELS_X_SIZE_BIG = 3\nKERNELS_Y_SIZE_BIG = 3\nKERNELS_INCREASED_X_SIZE_BIG = 7\nKERNELS_INCREASED_Y_SIZE_BIG = 7\nKERNELS_INCREASED_X_SIZE_BIG_HALF = (KERNELS_INCREASED_X_SIZE_BIG - 1) / 2\n\ndef int_round(inp):\n    return int(np.round(inp))\n\nRESIZE_FACTOR_SMALL_SCALE = RESIZE_FACTOR_SMALL\n\n\ndef bool_to_u8(inp):\n    return 255 * (inp.astype(np.uint8))\n\ndef int_round(inp):\n    return int(np.round(inp))\n\n\ndef gauss_kernel(half_width):\n    x = np.arange(-half_width, half_width + 1, dtype=np.float32)\n    pseudo_sigma = 0.9 * half_width\n    y = np.exp(-(x / pseudo_sigma)**2)\n    y = y / np.sum(y)\n    return y\n\n#def moving_average(x, kernel):\n#    return np.convolve(x, kernel, mode='same')\n\ndef extend_contour(contour, n_extended):\n    contour_extended = np.zeros((contour.shape[0] + 2 * n_extended, 2), np.float32)\n    contour_extended[n_extended: -n_extended, :] = contour\n    contour_extended[0: n_extended, :] = contour[-n_extended:, :]\n    contour_extended[-n_extended:, :] = contour[0: n_extended, :]\n    return contour_extended\n\ndef smooth_contour(contour, kernel):\n    kernel_half_width = (kernel.size - 1) // 2\n    contour_extended = extend_contour(contour, kernel_half_width)\n    contour_smoothed = np.zeros((contour.shape[0], 2), np.float32)\n    contour_smoothed[:, 0] = np.convolve(contour_extended[:, 0], kernel, mode='valid')\n    contour_smoothed[:, 1] = np.convolve(contour_extended[:, 1], kernel, mode='valid')\n    return contour_smoothed\n    \n\ndef natural_parametrization(contour, points_per_pixel):\n    x = contour[:, 0]\n    y = contour[:, 1]\n    dx = np.diff(x)\n    dy = np.diff(y)\n    dd = np.sqrt(dx**2 + dy**2)\n    dd = np.concatenate(([0], dd))\n    d = np.cumsum(dd)\n    \n    x_extended = np.concatenate((x, [x[0]]))\n    y_extended = np.concatenate((y, [y[0]]))\n    dd_extention = np.sqrt((x[0] - x[-1])**2 + (y[0] - y[-1])**2)\n    d_extended = np.concatenate((d, [d[-1] + dd_extention]))\n    \n    #new_size = int_round((new_size_relative * x.size))\n    new_size = int_round(points_per_pixel * d_extended[-1])\n    d_uniform = np.linspace(0, d_extended[-1], new_size + 1)[0:-1]\n    x_iterpolated = np.interp(d_uniform, d_extended, x_extended)\n    y_iterpolated = np.interp(d_uniform, d_extended, y_extended)\n    contour_interpolated = np.zeros((x_iterpolated.size, 2), np.float32)\n    contour_interpolated[:, 0] = x_iterpolated\n    contour_interpolated[:, 1] = y_iterpolated\n    return contour_interpolated\n\n\ndef delete_files(dir_):\n    files = glob.glob(dir_ + '/*')\n    for file in files:\n        os.remove(file)\n\ndef prepare_dir(dir_):\n    if os.path.exists(dir_):\n        delete_files(dir_)\n    else:\n        os.makedirs(dir_)\n\ndef prepare_crop_indexes(x1, x2, size, size_crop):\n    if x1 < 0:\n        x1_src = 0\n        x2_src = x2\n        x1_dst = size_crop - x2\n        x2_dst = size_crop\n    elif x2 > size:\n        x1_src = x1\n        x2_src = size\n        x1_dst = 0\n        x2_dst = size_crop - (x2 - size)\n    else:\n        x1_src = x1\n        x2_src = x2\n        x1_dst = 0\n        x2_dst = size_crop\n    return x1_src, x2_src, x1_dst, x2_dst\n\ndef round_with_mult_higher(inp, mult):\n    return ((np.ceil(inp / mult)) * mult).astype(np.int64)\n\ndef crop_with_overlap_1d(size, overlap, base_size, rounding):\n    n_crops_non_integer = (size - overlap) / (base_size - overlap)\n    n_crops = int_round(n_crops_non_integer)\n    delimeters = np.linspace(0, size, n_crops + 1)\n    starts_0 = np.zeros(n_crops, np.float32)\n    starts_0[0] = 0.0\n    ends_0 = np.zeros(n_crops, np.float32)\n    ends_0[-1] = size\n    \n    \n    starts_0[1: ] = delimeters[1: -1]  - overlap / 2\n    ends_0[0: -1] = delimeters[1: -1]  + overlap / 2\n    \n    centers = (starts_0 + ends_0) / 2\n    widths_0 = ends_0 - starts_0\n    widths = round_with_mult_higher(widths_0, rounding)\n    starts = np.zeros(n_crops, np.int64)\n    starts[0] = 0\n    starts[1: -1] = np.round(centers[1: -1] - widths[1: -1] / 2).astype(np.int64)\n    base_size_rounded = round_with_mult_higher(base_size, rounding)\n    starts[-1] = size - base_size_rounded\n    ends = np.zeros(n_crops, np.int64)\n    ends[-1] = size\n    ends[0] = base_size_rounded\n    ends[1: -1] = starts[1: -1] + widths[1: -1]\n    \n    separators = np.zeros(n_crops + 1, np.float32)\n    separators[0] = 0.0\n    separators[-1] = size - 1\n    \n    separators[1:-1] = (ends[0: -1] - 1 + starts[1:]) / 2\n    \n    \n    \n    \n    return starts, ends, separators\n\n\n#t = crop_with_overlap_1d(50000, 448, 2048, 64)  \n##plt.plot(t[0], t[0] * 0, 'rx')\n##plt.plot(t[1], t[1] * 0, 'g+')\n#for i in range(t[0].size):\n#    plt.plot([t[0][i], t[1][i]], [i, i], 'k+-')\n\ndef weight_variable(name, shape):\n    \n    #return tf.Variable(tf.constant_initializer(0.01)(shape=shape), name=name)\n    return tf.Variable(tf.initializers.GlorotUniform()(shape=shape), name=name)\n\ndef weight_variable_for_skip_connection(name, shape):\n    return tf.Variable(tf.constant_initializer(1.0 / shape[2])(shape=shape), name=name)\n\ndef bias_variable(name, shape):\n    \n    return tf.Variable(tf.zeros_initializer()(shape=shape), name=name)\n\n    \ndef conv2d(x, W):\n    return tf.nn.conv2d(x, W, strides=[1,1,1,1], padding='SAME')\n\ndef max_pool_2x2(x):\n    return tf.nn.max_pool(x, ksize=[1,2,2,1], strides=[1,2,2,1], padding='SAME')\n\ndef avarage_pool_2x2(x):\n    return tf.nn.avg_pool(x, ksize=[1,2,2,1], strides=[1,2,2,1], padding='SAME')\n\ndef avarage_pool_4x4(x):\n    return tf.nn.avg_pool(x, ksize=[1,4,4,1], strides=[1,4,4,1], padding='SAME')\n\nclass NetModuleSmallScale(tf.Module):\n\n        def __init__(self):\n            self.W_conv1 = weight_variable(\"W_conv1\", [KERNELS_Y_SIZE_SMALL, KERNELS_X_SIZE_SMALL, N_FEATUREMAPS_SMALL[0], N_FEATUREMAPS_SMALL[1]])\n            self.b_conv1 = bias_variable(\"b_conv1\", [N_FEATUREMAPS_SMALL[1]])\n            \n            block_name = \"block1\"\n            n_featuremaps = N_FEATUREMAPS_SMALL[1: 3 + 1]\n            self.block1_W_conv1 = weight_variable(block_name + \"_\" + \"W_conv1\", [KERNELS_Y_SIZE_SMALL, KERNELS_X_SIZE_SMALL, n_featuremaps[0], n_featuremaps[1]])\n            self.block1_b_conv1 = bias_variable(block_name + \"_\" + \"b_conv1\", [n_featuremaps[1]])\n            self.block1_W_conv2 = weight_variable(block_name + \"_\" + \"W_conv2\", [KERNELS_Y_SIZE_SMALL, KERNELS_X_SIZE_SMALL, n_featuremaps[1], n_featuremaps[2]])\n            self.block1_b_conv2 = bias_variable(block_name + \"_\" + \"b_conv2\", [n_featuremaps[2]])\n            self.block1_W_conv3 = weight_variable_for_skip_connection(block_name + \"_\" + \"W_conv3\", [1, 1, n_featuremaps[0], n_featuremaps[2]])\n            \n            block_name = \"block2\"\n            n_featuremaps = N_FEATUREMAPS_SMALL[3: 5 + 1]\n            self.block2_W_conv1 = weight_variable(block_name + \"_\" + \"W_conv1\", [KERNELS_Y_SIZE_SMALL, KERNELS_X_SIZE_SMALL, n_featuremaps[0], n_featuremaps[1]])\n            self.block2_b_conv1 = bias_variable(block_name + \"_\" + \"b_conv1\", [n_featuremaps[1]])\n            self.block2_W_conv2 = weight_variable(block_name + \"_\" + \"W_conv2\", [KERNELS_Y_SIZE_SMALL, KERNELS_X_SIZE_SMALL, n_featuremaps[1], n_featuremaps[2]])\n            self.block2_b_conv2 = bias_variable(block_name + \"_\" + \"b_conv2\", [n_featuremaps[2]])\n            self.block2_W_conv3 = weight_variable_for_skip_connection(block_name + \"_\" + \"W_conv3\", [1, 1, n_featuremaps[0], n_featuremaps[2]])\n            \n            block_name = \"block3\"\n            n_featuremaps = N_FEATUREMAPS_SMALL[5: 10 + 1] # must be N_FEATUREMAPS[5: 9 + 1] but they have same result\n            self.block3_W_conv1 = weight_variable(block_name + \"_\" + \"W_conv1\", [KERNELS_INCREASED_Y_SIZE_SMALL, KERNELS_INCREASED_X_SIZE_SMALL, n_featuremaps[0], n_featuremaps[1]])\n            self.block3_b_conv1 = bias_variable(block_name + \"_\" + \"b_conv1\", [n_featuremaps[1]])\n            \n            self.block3_W_conv2 = weight_variable(block_name + \"_\" + \"W_conv2\", [KERNELS_INCREASED_Y_SIZE_SMALL, KERNELS_INCREASED_X_SIZE_SMALL, n_featuremaps[1], n_featuremaps[2]])\n            self.block3_b_conv2 = bias_variable(block_name + \"_\" + \"b_conv2\", [n_featuremaps[2]])\n            \n            self.block3_W_conv4 = weight_variable(block_name + \"_\" + \"W_conv4\", [KERNELS_INCREASED_Y_SIZE_SMALL, KERNELS_INCREASED_X_SIZE_SMALL, n_featuremaps[2], n_featuremaps[3]])\n            self.block3_b_conv4 = bias_variable(block_name + \"_\" + \"b_conv4\", [n_featuremaps[3]])\n            \n            self.block3_W_conv5 = weight_variable(block_name + \"_\" + \"W_conv5\", [KERNELS_INCREASED_Y_SIZE_SMALL, KERNELS_INCREASED_X_SIZE_SMALL, n_featuremaps[3], n_featuremaps[4]])\n            self.block3_b_conv5 = bias_variable(block_name + \"_\" + \"b_conv5\", [n_featuremaps[4]])\n            \n            self.block3_W_conv3 = weight_variable_for_skip_connection(block_name + \"_\" + \"W_conv3\", [1, 1, n_featuremaps[0], n_featuremaps[4]])\n            \n            \n            \n            self.var_list = [self.W_conv1, self.b_conv1,\n                        self.block1_W_conv1, self.block1_b_conv1, self.block1_W_conv2, self.block1_b_conv2, self.block1_W_conv3, \n                        self.block2_W_conv1, self.block2_b_conv1, self.block2_W_conv2, self.block2_b_conv2, self.block2_W_conv3,\n                        self.block3_W_conv1, self.block3_b_conv1, self.block3_W_conv2, self.block3_b_conv2, self.block3_W_conv3, self.block3_W_conv4, self.block3_b_conv4, self.block3_W_conv5, self.block3_b_conv5]\n                        \n                        \n                    \n        @tf.function(input_signature=[tf.TensorSpec(shape=[1, None, None, N_INPUT_CHANNELS_SMALL], dtype=tf.float32)])\n        def predict(self, x):\n            h_conv1 = tf.nn.relu(conv2d(x, self.W_conv1) + self.b_conv1)\n        \n        \n            featuremaps_input = h_conv1\n            featuremaps_max_pool_1 = max_pool_2x2(featuremaps_input)\n            block_h_conv1 = tf.nn.relu(conv2d(featuremaps_max_pool_1, self.block1_W_conv1) + self.block1_b_conv1)\n            featuremaps_max_pool_2 = max_pool_2x2(block_h_conv1)\n            block_h_conv2 = tf.nn.relu(conv2d(featuremaps_max_pool_2, self.block1_W_conv2) + self.block1_b_conv2)\n            # skip connection:\n            featuremaps_avarage_pool = avarage_pool_4x4(featuremaps_input)\n            block_h_conv3 = conv2d(featuremaps_avarage_pool, self.block1_W_conv3)\n            # sum main way and skip connection:\n            featuremaps_output = block_h_conv2 + block_h_conv3\n           \n            \n            featuremaps_input = featuremaps_output\n            featuremaps_max_pool_1 = max_pool_2x2(featuremaps_input)\n            block_h_conv1 = tf.nn.relu(conv2d(featuremaps_max_pool_1, self.block2_W_conv1) + self.block2_b_conv1)\n            featuremaps_max_pool_2 = max_pool_2x2(block_h_conv1)\n            block_h_conv2 = tf.nn.relu(conv2d(featuremaps_max_pool_2, self.block2_W_conv2) + self.block2_b_conv2)\n            # skip connection:\n            featuremaps_avarage_pool = avarage_pool_4x4(featuremaps_input)\n            block_h_conv3 = conv2d(featuremaps_avarage_pool, self.block2_W_conv3)\n            # sum main way and skip connection:\n            featuremaps_output = block_h_conv2 + block_h_conv3\n            \n            featuremaps_input = featuremaps_output\n            #featuremaps_max_pool_1 = max_pool_2x2(featuremaps_input)\n            featuremaps_max_pool_1 = featuremaps_input\n            block_h_conv1 = tf.nn.relu(conv2d(featuremaps_max_pool_1, self.block3_W_conv1) + self.block3_b_conv1)\n            block_h_conv2 = tf.nn.relu(conv2d(block_h_conv1, self.block3_W_conv2) + self.block3_b_conv2)\n            block_h_conv4 = tf.nn.relu(conv2d(block_h_conv2, self.block3_W_conv4) + self.block3_b_conv4)\n            block_h_conv5 = conv2d(block_h_conv4, self.block3_W_conv5) + self.block3_b_conv5\n            # skip connection:\n            #featuremaps_avarage_pool = avarage_pool_2x2(featuremaps_input)\n            featuremaps_avarage_pool = featuremaps_input\n            block_h_conv3 = conv2d(featuremaps_max_pool_1, self.block3_W_conv3)\n            # sum main way and skip connection:\n            output = tf.math.add(block_h_conv5, block_h_conv3, name='output')\n        \n            return output\n        \n        def load_weights(self, saved_wights_path):\n            ckpt = tf.train.Checkpoint(var_list=self.var_list)\n            ckpt.read(saved_wights_path)\n\nclass NetModuleBigScale(tf.Module):\n    def __init__(self):\n            self.W_conv1 = weight_variable(\"W_conv1\", [KERNELS_Y_SIZE_BIG, KERNELS_X_SIZE_BIG, N_FEATUREMAPS_BIG[0], N_FEATUREMAPS_BIG[1]])\n            self.b_conv1 = bias_variable(\"b_conv1\", [N_FEATUREMAPS_BIG[1]])\n            \n            block_name = \"block1\"\n            n_featuremaps = N_FEATUREMAPS_BIG[1: 3 + 1]\n            self.block1_W_conv1 = weight_variable(block_name + \"_\" + \"W_conv1\", [KERNELS_Y_SIZE_BIG, KERNELS_X_SIZE_BIG, n_featuremaps[0], n_featuremaps[1]])\n            self.block1_b_conv1 = bias_variable(block_name + \"_\" + \"b_conv1\", [n_featuremaps[1]])\n            self.block1_W_conv2 = weight_variable(block_name + \"_\" + \"W_conv2\", [KERNELS_Y_SIZE_BIG, KERNELS_X_SIZE_BIG, n_featuremaps[1], n_featuremaps[2]])\n            self.block1_b_conv2 = bias_variable(block_name + \"_\" + \"b_conv2\", [n_featuremaps[2]])\n            self.block1_W_conv3 = weight_variable_for_skip_connection(block_name + \"_\" + \"W_conv3\", [1, 1, n_featuremaps[0], n_featuremaps[2]])\n            \n            block_name = \"block2\"\n            n_featuremaps = N_FEATUREMAPS_BIG[3: 5 + 1]\n            self.block2_W_conv1 = weight_variable(block_name + \"_\" + \"W_conv1\", [KERNELS_Y_SIZE_BIG, KERNELS_X_SIZE_BIG, n_featuremaps[0], n_featuremaps[1]])\n            self.block2_b_conv1 = bias_variable(block_name + \"_\" + \"b_conv1\", [n_featuremaps[1]])\n            self.block2_W_conv2 = weight_variable(block_name + \"_\" + \"W_conv2\", [KERNELS_Y_SIZE_BIG, KERNELS_X_SIZE_BIG, n_featuremaps[1], n_featuremaps[2]])\n            self.block2_b_conv2 = bias_variable(block_name + \"_\" + \"b_conv2\", [n_featuremaps[2]])\n            self.block2_W_conv3 = weight_variable_for_skip_connection(block_name + \"_\" + \"W_conv3\", [1, 1, n_featuremaps[0], n_featuremaps[2]])\n            \n            block_name = \"block3\"\n            n_featuremaps = N_FEATUREMAPS_BIG[5: 10 + 1] # must be N_FEATUREMAPS[5: 9 + 1] but they have same result\n            self.block3_W_conv1 = weight_variable(block_name + \"_\" + \"W_conv1\", [KERNELS_INCREASED_Y_SIZE_BIG, KERNELS_INCREASED_X_SIZE_BIG, n_featuremaps[0], n_featuremaps[1]])\n            self.block3_b_conv1 = bias_variable(block_name + \"_\" + \"b_conv1\", [n_featuremaps[1]])\n            \n            self.block3_W_conv2 = weight_variable(block_name + \"_\" + \"W_conv2\", [N_R_OUT, KERNELS_INCREASED_X_SIZE_BIG, n_featuremaps[1], n_featuremaps[2]])\n            self.block3_b_conv2 = bias_variable(block_name + \"_\" + \"b_conv2\", [n_featuremaps[2]])\n            \n            self.block3_W_conv4 = weight_variable(block_name + \"_\" + \"W_conv4\", [1, KERNELS_INCREASED_X_SIZE_BIG, n_featuremaps[2], n_featuremaps[3]])\n            self.block3_b_conv4 = bias_variable(block_name + \"_\" + \"b_conv4\", [n_featuremaps[3]])\n            \n            self.block3_W_conv5 = weight_variable(block_name + \"_\" + \"W_conv5\", [1, KERNELS_INCREASED_X_SIZE_BIG, n_featuremaps[3], n_featuremaps[4]])\n            self.block3_b_conv5 = bias_variable(block_name + \"_\" + \"b_conv5\", [n_featuremaps[4]])\n            \n            self.block3_W_conv3 = weight_variable_for_skip_connection(block_name + \"_\" + \"W_conv3\", [N_R_OUT, 1, n_featuremaps[0], n_featuremaps[4]])\n            \n            \n            self.upsample_tile_parameter_tf_const = tf.constant([1, 1, 2, 1, 2, 1],  tf.int32)\n            self.upsample_reshape_parameter_ft_const = tf.constant([1, N_R_OUT, N_PHI_OUT_TILLED, N_FEATUREMAPS_BIG[-4]],  tf.int32)\n            \n            \n            self.var_list = [self.W_conv1, self.b_conv1,\n                        self.block1_W_conv1, self.block1_b_conv1, self.block1_W_conv2, self.block1_b_conv2, self.block1_W_conv3, \n                        self.block2_W_conv1, self.block2_b_conv1, self.block2_W_conv2, self.block2_b_conv2, self.block2_W_conv3,\n                        self.block3_W_conv1, self.block3_b_conv1, self.block3_W_conv2, self.block3_b_conv2, self.block3_W_conv3, self.block3_W_conv4, self.block3_b_conv4, self.block3_W_conv5, self.block3_b_conv5]\n                        \n                        \n    \n\n    @tf.function(input_signature=[tf.TensorSpec(shape=[1, N_R, N_PHI_TILED, N_INPUT_CHANNELS_BIG], dtype=tf.float32)])\n    def predict(self, x):\n        h_conv1 = tf.nn.relu(conv2d(x, self.W_conv1) + self.b_conv1)\n        \n        \n        featuremaps_input = h_conv1\n        featuremaps_max_pool_1 = max_pool_2x2(featuremaps_input)\n        block_h_conv1 = tf.nn.relu(conv2d(featuremaps_max_pool_1, self.block1_W_conv1) + self.block1_b_conv1)\n        featuremaps_max_pool_2 = max_pool_2x2(block_h_conv1)\n        block_h_conv2 = tf.nn.relu(conv2d(featuremaps_max_pool_2, self.block1_W_conv2) + self.block1_b_conv2)\n        # skip connection:\n        featuremaps_avarage_pool = avarage_pool_4x4(featuremaps_input)\n        block_h_conv3 = conv2d(featuremaps_avarage_pool, self.block1_W_conv3)\n        # sum main way and skip connection:\n        featuremaps_output = block_h_conv2 + block_h_conv3\n       \n        \n        featuremaps_input = featuremaps_output\n        featuremaps_max_pool_1 = max_pool_2x2(featuremaps_input)\n        block_h_conv1 = tf.nn.relu(conv2d(featuremaps_max_pool_1, self.block2_W_conv1) + self.block2_b_conv1)\n        \n        # skip connection end ealier now:\n        featuremaps_avarage_pool = avarage_pool_2x2(featuremaps_input)\n        block_h_conv3 = conv2d(featuremaps_avarage_pool, self.block2_W_conv3)\n        # sum main way and skip connection:\n        featuremaps_output = block_h_conv1 + block_h_conv3\n        \n        input_for_long_skip_connection = featuremaps_output\n        \n        featuremaps_max_pool_2 = max_pool_2x2(block_h_conv1)\n        block_h_conv2 = tf.nn.relu(conv2d(featuremaps_max_pool_2, self.block2_W_conv2) + self.block2_b_conv2)\n        \n        ## skip connection:\n        #featuremaps_avarage_pool = avarage_pool_2x2(featuremaps_input)\n        #block_h_conv3 = conv2d(featuremaps_avarage_pool, block2_W_conv3)\n        ## sum main way and skip connection:\n        #featuremaps_output = block_h_conv2 + block_h_conv3\n        featuremaps_output = block_h_conv2\n        \n        featuremaps_input = featuremaps_output\n        #featuremaps_max_pool_1 = max_pool_2x2(featuremaps_input)\n        featuremaps_max_pool_1 = featuremaps_input\n        block_h_conv1 = tf.nn.relu(conv2d(featuremaps_max_pool_1, self.block3_W_conv1) + self.block3_b_conv1)\n        \n        # upsample 2x by repeat values 4 times:\n        #upsample_tile_parameter_tf_const = tf.constant([1, 1, 2, 1, 2, 1],  tf.int32)\n        #upsample_reshape_parameter_ft_const = tf.constant([BATCH_SIZE, OUTPUT_SIZE, OUTPUT_SIZE, N_FEATUREMAPS[-4]],  tf.int32)\n        upsample_tmp_1 = tf.expand_dims(block_h_conv1, 2)\n        upsample_tmp_2 = tf.expand_dims(upsample_tmp_1, 4)\n        upsample_tmp_3 = tf.tile(upsample_tmp_2, self.upsample_tile_parameter_tf_const)\n        block_h_upsample2 = tf.reshape(upsample_tmp_3, self.upsample_reshape_parameter_ft_const)\n        \n        strides = [1,1,1,1]\n        padding = [[0, 0], [0, 0], [KERNELS_INCREASED_X_SIZE_BIG_HALF, KERNELS_INCREASED_X_SIZE_BIG_HALF], [0, 0]]\n        block_h_conv2 = tf.nn.relu(tf.nn.conv2d(block_h_upsample2, self.block3_W_conv2, strides=strides, padding=padding) + self.block3_b_conv2)\n        block_h_conv4 = tf.nn.relu(tf.nn.conv2d(block_h_conv2, self.block3_W_conv4, strides=strides, padding=padding) + self.block3_b_conv4)\n        block_h_conv5 = tf.nn.conv2d(block_h_conv4, self.block3_W_conv5, strides=strides, padding=padding) + self.block3_b_conv5\n        # skip connection:\n        #featuremaps_avarage_pool = avarage_pool_2x2(featuremaps_input)\n        padding = [[0, 0], [0, 0], [0, 0], [0, 0]]\n        block_h_conv3 = tf.nn.conv2d(input_for_long_skip_connection, self.block3_W_conv3, strides=strides, padding=padding)\n        # sum main way and skip connection:\n        output = tf.math.add(block_h_conv5, block_h_conv3, name='output')\n        #output =block_h_conv1\n    \n        return output\n    \n    def load_weights(self, saved_wights_path):\n        ckpt = tf.train.Checkpoint(var_list=self.var_list)\n        ckpt.read(saved_wights_path)\n    \n\nnet_module_small_scale = NetModuleSmallScale()\nnet_module_small_scale.load_weights(WEIGHTS_SMALL)\n\nnet_module_big_scale = NetModuleBigScale()\nnet_module_big_scale.load_weights(WEIGHTS_BIG)\n\n\n\nphi = np.linspace(0, 2*np.pi, N_PHI + 1)[0: -1]\nphi = phi.astype(np.float32)\nphi_out = phi.reshape(N_PHI // NET_SCALE_BIG, NET_SCALE_BIG)\nphi_out = np.mean(phi_out, axis=1)\nr = np.linspace(0, R_MAX, N_R + 1)[0: -1]\nr = r.astype(np.float32)\nr_indexes = np.arange(r.size).astype(np.float32)\nPHI, R = np.meshgrid(phi, r)\nmap_x = CENTER + R * np.cos(PHI) # maps for unroll with cv2.remap\nmap_y = CENTER + R * np.sin(PHI) # maps for unroll with cv2.remap\n        \n\n\n\n\nfiles = glob.glob(IMAGES_DIR + '/*.tiff')\n\noutput_file_fid = open(OUTPUT_FILE, 'w')\noutput_file_fid.write('id,predicted\\n')\n\nfor file_index in range(len(files)):\n\n    \n    file = files[file_index]\n    file_base = os.path.basename(file)\n    file_base_no_ext = os.path.splitext(file_base)[0]\n    \n    \n    image_id = file_base_no_ext\n    \n    \n\n    print(file_index, len(files) - 1)\n    \n\n    \n    \n    # 1st net inference:\n    \n    \n    fid = rasterio.open(IMAGES_DIR + '/' + image_id + '.tiff', 'r', num_threads='all_cpus')\n    \n    width = fid.shape[1]\n    height = fid.shape[0]\n    \n    if fid.count == 1:\n        layers = []\n        for subdataset_index in range(len(fid.subdatasets)):\n            subdataset = fid.subdatasets[subdataset_index]\n            layer = rasterio.open(subdataset)\n            layers.append(layer)\n    \n    starts_y, ends_y, separators_y = crop_with_overlap_1d(height, OVERLAP_SMALL, BASE_SIZE_CROP_SMALL, ROUNDING_ORIG_SMALL)\n    starts_x, ends_x, separators_x = crop_with_overlap_1d(width, OVERLAP_SMALL, BASE_SIZE_CROP_SMALL, ROUNDING_ORIG_SMALL)\n    \n    x_spot_all = []\n    y_spot_all = []\n    for y_index in range(starts_y.size):\n        y1_small = starts_y[y_index]\n        y2_small = ends_y[y_index]\n        separator_y1 = separators_y[y_index]\n        separator_y2 = separators_y[y_index + 1]\n        for x_index in range(starts_x.size):\n            x1_small = starts_x[x_index]\n            x2_small = ends_x[x_index]\n            separator_x1 = separators_x[x_index]\n            separator_x2 = separators_x[x_index + 1]\n            window = Window.from_slices((y1_small, y2_small),(x1_small, x2_small))\n            height_tmp = y2_small - y1_small\n            width_tmp = x2_small - x1_small\n            if fid.count == 1:\n                image_crop_small = np.zeros((height_tmp, width_tmp, 3), np.uint8)\n                for layer_index in range(len(layers)):\n                    layer = layers[layer_index]\n                    image_crop_small[:, :, layer_index] = layer.read(1, window=window)\n            else:\n                image_crop_small = np.moveaxis(fid.read([1, 2, 3], window=window), 0, -1)\n            \n            image_crop_small_resized = cv2.resize(image_crop_small, None, fx=RESIZE_FACTOR_SMALL, fy=RESIZE_FACTOR_SMALL, interpolation=cv2.INTER_CUBIC)\n            \n            \n                        \n            image_crop_small_resized_4d = np.expand_dims(image_crop_small_resized, 0).astype(np.float32) / 255.0\n            mask_crop_small = net_module_small_scale.predict(image_crop_small_resized_4d)\n            mask_crop_small = mask_crop_small.numpy()[0, :, :, 0]\n            \n            \n            heatmap_thresholded  = 255 * ((mask_crop_small >= THRESHOLD_SMALL).astype(np.uint8))\n            n_labels, labels, stats, centroids_tmp_0 = cv2.connectedComponentsWithStats(heatmap_thresholded, connectivity=4)\n            \n                                    \n            # (n_label - 1) - is number of detected glomeruli\n            if n_labels > 1:\n                \n                for max_ind in range(1, n_labels): # max_ind=0 is for background - not needed\n                \n                    y_bigest_spot, x_bigest_spot = np.where(labels == max_ind)\n                    weights = mask_crop_small[y_bigest_spot, x_bigest_spot]\n                    weights_sum = np.sum(weights)\n                    x_spot_out = np.sum(x_bigest_spot * weights) / weights_sum # center is calculated as weighted mean\n                    y_spot_out = np.sum(y_bigest_spot * weights) / weights_sum # center is calculated as weighted mean\n                    x_spot = (x_spot_out + 0.5) * NET_SCALE_SMALL / RESIZE_FACTOR_SMALL - 0.5\n                    y_spot = (y_spot_out + 0.5) * NET_SCALE_SMALL / RESIZE_FACTOR_SMALL - 0.5\n                    x_spot += x1_small\n                    y_spot += y1_small\n                    condition = True\n                    condition = condition and separator_x1 <= x_spot\n                    if x_index == (starts_x.size - 1):\n                        condition = condition and x_spot <= separator_x2\n                    else:\n                        condition = condition and x_spot < separator_x2\n                    condition = condition and separator_y1 <= y_spot\n                    if y_index == (starts_y.size - 1):\n                        condition = condition and y_spot <= separator_y2\n                    else:\n                        condition = condition and y_spot < separator_y2\n                    if condition:\n                        x_spot_all.append(x_spot)\n                        y_spot_all.append(y_spot)\n            \n            \n            \n\n            \n            \n            del image_crop_small_resized_4d\n            del image_crop_small\n            del image_crop_small_resized\n            del mask_crop_small\n            gc.collect()\n    \n    \n    if fid.count == 1:    \n        for layer in layers:\n            layer.close()\n        del layers\n        \n    fid.close()\n    gc.collect()\n    \n    #plt.plot(x_spot_all, y_spot_all, 'k.')\n    #plt.axis('equal')\n    #plt.gca().invert_yaxis()\n    \n    \n    \n    # 2nd net inference:\n    \n    fid = rasterio.open(file, 'r', num_threads='all_cpus')\n    if fid.count == 1:\n        layers = []\n        for subdataset_index in range(len(fid.subdatasets)):\n            subdataset = fid.subdatasets[subdataset_index]\n            layer = rasterio.open(subdataset)\n            layers.append(layer)\n            \n    mask_prediction = np.zeros((height, width), np.bool)\n    for spot_index in range(len(x_spot_all)):\n        x_spot = int_round(x_spot_all[spot_index])\n        y_spot = int_round(y_spot_all[spot_index])\n        x1_big = x_spot - BASE_SIZE_HALF_BIG\n        x2_big = x1_big + BASE_SIZE_BIG\n        y1_big = y_spot - BASE_SIZE_HALF_BIG\n        y2_big = y1_big + BASE_SIZE_BIG\n        \n        \n        x1_big_src, x2_big_src, x1_big_dst, x2_big_dst = prepare_crop_indexes(x1_big, x2_big, width, BASE_SIZE_BIG)\n        y1_big_src, y2_big_src, y1_big_dst, y2_big_dst = prepare_crop_indexes(y1_big, y2_big, height, BASE_SIZE_BIG)\n        \n        \n        image_crop_big = np.zeros((BASE_SIZE_BIG, BASE_SIZE_BIG, 3), np.uint8)\n        window = Window.from_slices((y1_big_src, y2_big_src),(x1_big_src, x2_big_src))\n        if fid.count == 1:\n            for layer_index in range(len(layers)):\n                layer = layers[layer_index]\n                image_crop_big[y1_big_dst: y2_big_dst, x1_big_dst: x2_big_dst, layer_index] = layer.read(1, window=window)\n        else:\n            image_crop_big[y1_big_dst: y2_big_dst, x1_big_dst: x2_big_dst, :] = np.moveaxis(fid.read([1, 2, 3], window=window), 0, -1)\n        \n        #plt.subplot(1, 2, 1)\n        #plt.imshow(image_crop_big)\n        \n        unrolled = cv2.remap(image_crop_big, map_x, map_y, cv2.INTER_CUBIC)\n        unrolled = np.tile(unrolled, (1, 3, 1))\n        \n        del image_crop_big\n        #gc.collect()\n        \n        unrolled_4d = np.expand_dims(unrolled, 0).astype(np.float32) / 255.0\n        del unrolled\n        #gc.collect()\n        big_predicted_tensor = net_module_big_scale.predict(unrolled_4d)\n        del unrolled_4d\n        #gc.collect()\n        r_mean = big_predicted_tensor.numpy()[0, 0, N_PHI_OUT: 2 * N_PHI_OUT, 0]\n        r_mean = r_mean * r_mean_max\n        r_mean = (r_mean + 0.5) * NET_SCALE_BIG - 0.5\n        #r_mean_x = (np.arange(r_mean.size) + 0.5) * NET_SCALE_BIG - 0.5\n        #phi_out\n        r_mean = np.clip(r_mean, 0.0, N_R - 1)\n        \n        r_mean = np.interp(r_mean, r_indexes, r)\n        \n        x_contour = CENTER + r_mean * np.cos(phi_out)\n        y_contour = CENTER + r_mean * np.sin(phi_out)\n        \n        contour_cv2 = np.zeros((x_contour.size, 1, 2), np.int32)\n        contour_cv2[:, 0, 0] = np.round(x_contour).astype(np.int32)\n        contour_cv2[:, 0, 1] = np.round(y_contour).astype(np.int32)\n        \n        mask_center_orig_u8 = np.zeros((BASE_SIZE_BIG, BASE_SIZE_BIG), np.uint8)\n        mask_center_orig_u8 = cv2.drawContours(mask_center_orig_u8, [contour_cv2], -1, (255, 255, 255), -1)\n        mask_center_orig = mask_center_orig_u8 >= 128\n        \n        #plt.subplot(1, 2, 2)\n        #plt.imshow(mask_center_orig)\n        #import sys\n        #sys.exit()\n\n        mask_prediction[y1_big_src: y2_big_src, x1_big_src: x2_big_src] = mask_prediction[y1_big_src: y2_big_src, x1_big_src: x2_big_src] | mask_center_orig[y1_big_dst: y2_big_dst, x1_big_dst: x2_big_dst]\n        \n        del mask_center_orig_u8\n        del mask_center_orig\n        del r_mean\n        gc.collect()\n        \n        \n        \n\n        \n    if fid.count == 1:    \n        for layer in layers:\n            layer.close()\n        del layers\n        \n    fid.close()\n    \n    \n    del fid\n    \n    gc.collect()\n    \n    \n    \n    \n    # convert predicted mask to RLE:\n    \n    y, x = np.nonzero(mask_prediction)\n    del mask_prediction\n    gc.collect()\n    indexes = np.sort(x * height + y + 1)\n    del x\n    gc.collect()\n    del y\n    gc.collect()\n    \n    \n    diff = np.diff(indexes)\n    gaps = np.nonzero(diff > 1)[0]\n    del diff\n    gc.collect()\n    \n    index_min = np.min(indexes)\n    index_max = np.max(indexes)\n    \n    starts = []\n    lengths = []\n\n    \n    start_tmp = index_min\n    for intermidieate_index in range(gaps.size):\n        starts.append(start_tmp)\n        lengths.append(indexes[gaps[intermidieate_index]] - start_tmp + 1)\n        start_tmp = indexes[gaps[intermidieate_index] + 1]\n    del indexes\n    gc.collect()\n    del gaps\n    gc.collect()\n    starts.append(start_tmp)\n    lengths.append(index_max - start_tmp + 1)\n    \n\n    line = image_id + ','\n    for segment_index in range(len(starts)):\n        start = starts[segment_index]\n        length = lengths[segment_index]\n        line = line + str(start) + ' ' + str(length)\n        if segment_index != (len(starts) - 1):\n            line = line + ' '\n    del starts\n    gc.collect()\n    del lengths\n    gc.collect()\n    output_file_fid.write(line + '\\n')\n    del line\n    gc.collect()\n    \n    \noutput_file_fid.close()    \n","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}