--- a +++ b/uNet_FullImage_Segmentation_FullDirectory.py @@ -0,0 +1,454 @@ +# %% importing packages + +import numpy as np +import tensorflow as tf +from tensorflow import keras +from tensorflow.keras import layers +from tensorflow.keras import mixed_precision +from tensorflow.python.ops.numpy_ops import np_config +np_config.enable_numpy_behavior() +from skimage import measure +from skimage import morphology +from scipy import ndimage +import cv2 as cv +import os +import matplotlib.pyplot as plt +import tqdm +from natsort import natsorted +plt.rcParams['figure.figsize'] = [50, 150] + + +# %% Citations +############################################################# +############################################################# + + +# Defining Functions +############################################################# +############################################################# + +class EncoderBlock(layers.Layer): + '''This function returns an encoder block with two convolutional layers and + an option for returning both a max-pooled output with a stride and pool + size of (2,2) and the output of the second convolution for skip + connections implemented later in the network during the decoding + section. All padding is set to "same" for cleanliness. + + When initializing it receives the number of filters to be used in both + of the convolutional layers as well as the kernel size and stride for + those same layers. It also receives the trainable variable for use with + the batch normalization layers.''' + + def __init__(self, + filters, + kernel_size=(3,3), + strides=(1,1), + trainable=True, + name='encoder_block', + **kwargs): + + super(EncoderBlock,self).__init__(trainable, name, **kwargs) + # When initializing this object receives a trainable parameter for + # freezing the convolutional layers. + + # including the image normalization within the network for easier image + # processing during inference + self.image_normalization = layers.Normalization() + + # below creates the first of two convolutional layers + self.conv1 = layers.Conv2D(filters=filters, + kernel_size=kernel_size, + strides=strides, + padding='same', + name='encoder_conv1', + trainable=trainable) + + # second of two convolutional layers + self.conv2 = layers.Conv2D(filters=filters, + kernel_size=kernel_size, + strides=strides, + padding='same', + name='encoder_conv2', + trainable=trainable) + + # creates the max-pooling layer for downsampling the image. + self.enc_pool = layers.MaxPool2D(pool_size=(2,2), + strides=(2,2), + padding='same', + name='enc_pool') + + # ReLU layer for activations. + self.ReLU = layers.ReLU() + + # both batch normalization layers for use with their corresponding + # convolutional layers. + self.batch_norm1 = tf.keras.layers.BatchNormalization() + self.batch_norm2 = tf.keras.layers.BatchNormalization() + + def call(self,input,training=True,include_pool=True): + + # first conv of the encoder block + x = self.image_normalization(input) + x = self.conv1(x) + x = self.batch_norm1(x,training=training) + x = self.ReLU(x) + + # second conv of the encoder block + x = self.conv2(x) + x = self.batch_norm2(x,training=training) + x = self.ReLU(x) + + # calculate and include the max pooling layer if include_pool is true. + # This output is used for the skip connections later in the network. + if include_pool: + pooled_x = self.enc_pool(x) + return(x,pooled_x) + + else: + return(x) + + +############################################################# + +class DecoderBlock(layers.Layer): + '''This function returns a decoder block that when called receives both an + input and a "skip connection". The input is passed to the + "up convolution" or transpose conv layer to double the dimensions before + being concatenated with its associated skip connection from the encoder + section of the network. All padding is set to "same" for cleanliness. + The decoder block also has an option for including an additional + "segmentation" layer, which is a (1,1) convolution with 4 filters, which + produces the logits for the one-hot encoded ground truth. + + When initializing it receives the number of filters to be used in the + up convolutional layer as well as the other two forward convolutions. + The received kernel_size and stride is used for the forward convolutions, + with the up convolution kernel and stride set to be (2,2).''' + def __init__(self, + filters, + trainable=True, + kernel_size=(3,3), + strides=(1,1), + name='DecoderBlock', + **kwargs): + + super(DecoderBlock,self).__init__(trainable, name, **kwargs) + + # creating the up convolution layer + self.up_conv = layers.Conv2DTranspose(filters=filters, + kernel_size=(2,2), + strides=(2,2), + padding='same', + name='decoder_upconv', + trainable=trainable) + + # the first of two forward convolutional layers + self.conv1 = layers.Conv2D(filters=filters, + kernel_size=kernel_size, + strides=strides, + padding='same', + name ='decoder_conv1', + trainable=trainable) + + # second convolutional layer + self.conv2 = layers.Conv2D(filters=filters, + kernel_size=kernel_size, + strides=strides, + padding='same', + name ='decoder_conv2', + trainable=trainable) + + # this creates the output prediction logits layer. + self.seg_out = layers.Conv2D(filters=7, + kernel_size=(1,1), + name='conv_feature_map') + + # ReLU for activation of all above layers + self.ReLU = layers.ReLU() + + # the individual batch normalization layers for their respective + # convolutional layers. + self.batch_norm1 = tf.keras.layers.BatchNormalization() + self.batch_norm2 = tf.keras.layers.BatchNormalization() + + + def call(self,input,skip_conn,training=True,segmentation=False): + + up = self.up_conv(input) # perform image up convolution + # concatenate the input and the skip_conn along the features axis + concatenated = layers.concatenate([up,skip_conn],axis=-1) + + # first convolution + x = self.conv1(concatenated) + x = self.batch_norm1(x,training=training) + x = self.ReLU(x) + + # second convolution + x = self.conv2(x) + x = self.batch_norm2(x,training=training) + x = self.ReLU(x) + + # if segmentation is True, then run the segmentation (1,1) convolution + # and use the Softmax to produce a probability distribution. + if segmentation: + seg = self.seg_out(x) + # deliberately set as "float32" to ensure proper calculation if + # switching to mixed precision for efficiency + prob = layers.Softmax(dtype='float32')(seg) + return(prob) + + else: + return(x) + + +############################################################# + +class uNet(keras.Model): + '''This is a sub-classed model that uses the encoder and decoder blocks + defined above to create a custom unet. The differences from the original + paper include a variable filter scalar (filter_multiplier), batch + normalization between each convolutional layer and the associated ReLU + activation, as well as feature normalization implemented in the first + layer of the network.''' + def __init__(self,filter_multiplier=2,**kwargs): + super(uNet,self).__init__() + + # Defining encoder blocks + self.encoder_block1 = EncoderBlock(filters=2*filter_multiplier, + name='Enc1') + self.encoder_block2 = EncoderBlock(filters=4*filter_multiplier, + name='Enc2') + self.encoder_block3 = EncoderBlock(filters=8*filter_multiplier, + name='Enc3') + self.encoder_block4 = EncoderBlock(filters=16*filter_multiplier, + name='Enc4') + self.encoder_block5 = EncoderBlock(filters=32*filter_multiplier, + name='Enc5') + + # Defining decoder blocks. The names are in reverse order to make it + # (hopefully) easier to understand which skip connections are associated + # with which decoder layers. + self.decoder_block4 = DecoderBlock(filters=16*filter_multiplier, + name='Dec4') + self.decoder_block3 = DecoderBlock(filters=8*filter_multiplier, + name='Dec3') + self.decoder_block2 = DecoderBlock(filters=4*filter_multiplier, + name='Dec2') + self.decoder_block1 = DecoderBlock(filters=2*filter_multiplier, + name='Dec1') + + + def call(self,inputs,training): + + # encoder + enc1,enc1_pool = self.encoder_block1(input=inputs,training=training) + enc2,enc2_pool = self.encoder_block2(input=enc1_pool,training=training) + enc3,enc3_pool = self.encoder_block3(input=enc2_pool,training=training) + enc4,enc4_pool = self.encoder_block4(input=enc3_pool,training=training) + enc5 = self.encoder_block5(input=enc4_pool, + include_pool=False, + training=training) + + # decoder + dec4 = self.decoder_block4(input=enc5,skip_conn=enc4,training=training) + dec3 = self.decoder_block3(input=dec4,skip_conn=enc3,training=training) + dec2 = self.decoder_block2(input=dec3,skip_conn=enc2,training=training) + seg_logits_out = self.decoder_block1(input=dec2, + skip_conn=enc1, + segmentation=True, + training=training) + + return(seg_logits_out) + +############################################################# + +def get_image_blocks(image,tile_distance=512,tile_size=1024): + '''Receives an image as well as a minimum distance between tiles. + Returns the name of the image processed, the image dimensions, and a list + of tile centers evenly distributed across the tissue surface.''' + image_dimensions = image.shape + + safe_mask = np.zeros([image_dimensions[0],image_dimensions[1]]) + safe_mask[int(tile_size/2):image_dimensions[0]-int(tile_size/2), + int(tile_size/2):image_dimensions[1]-int(tile_size/2)] = 1 + + grid_0 = np.arange(0,image_dimensions[0],tile_distance) + grid_1 = np.arange(0,image_dimensions[1],tile_distance) + + + + center_indexes = [] + + for grid0 in grid_0: + for grid1 in grid_1: + if safe_mask[grid0,grid1]: + center_indexes.append([grid0,grid1]) + + return([image_dimensions,center_indexes]) + +############################################################# + +def get_reduced_tile_indexes(tile_center,returned_size=1024): + start_0 = int(tile_center[0] - returned_size/2) + end_0 = int(tile_center[0] + returned_size/2) + + start_1 = int(tile_center[1] - returned_size/2) + end_1 = int(tile_center[1] + returned_size/2) + + return([start_0,end_0],[start_1,end_1]) + +############################################################# + +def segment_tiles(unet,center_indexes,image,scaling_factor=1,tile_size=1024): + + m,n,z = image.shape + segmentation = np.zeros((m,n)) + + for idx in tqdm.tqdm(range(len(center_indexes))): + center = center_indexes[idx] + dim0, dim1 = get_reduced_tile_indexes(center,tile_size) + sub_sectioned_tile = image[dim0[0]:dim0[1],dim1[0]:dim1[1]] + + full_tile_dim0,full_tile_dim1,z = sub_sectioned_tile.shape + + color_tile = sub_sectioned_tile[:,:,0:3] + + if scaling_factor > 1: + height = color_tile.shape[0] + width = color_tile.shape[1] + + height2 = int(height/scaling_factor) + width2 = int(width/scaling_factor) + + color_tile = cv.resize(color_tile,[height2,width2],cv.INTER_AREA) + + color_tile = color_tile[None,:,:,:] + + prediction = unet.predict(color_tile,verbose=0) + + prediction_tile = np.squeeze(np.asarray(tf.argmax(prediction,axis=-1)).astype(np.float32).copy()) + + if scaling_factor > 1: + prediction_tile = cv.resize(prediction_tile,[full_tile_dim0,full_tile_dim1],cv.INTER_NEAREST) + + + dim0, dim1 = get_reduced_tile_indexes(center,returned_size=512) + + # fix this hard coding of the tile indexes for the prediction + segmentation[dim0[0]:dim0[1],dim1[0]:dim1[1]] = prediction_tile[256:768,256:768] + + return(segmentation) + +############################################################# + +def segment_directory(JPG_directory, + unet,tile_size=2048, + tile_distance=512, + scaling_factor=2, + HeartID='0', + ): + os.chdir(JPG_directory) + + out_directory = './../uNet_Segmentations/' + + # create the directory for saving if it doesn't already exist + if not os.path.isdir(out_directory): + os.mkdir(out_directory) + + os.chdir(out_directory) + + file_names = tf.io.gfile.glob(JPG_directory + HeartID + '*.jpg') + + for idx,file in enumerate(file_names): + print(f'segmenting file {idx} of {len(file_names)}') + + file_id = file.split('/')[-1].split('.')[0] + + image = cv.imread(file,cv.IMREAD_UNCHANGED) + image = cv.copyMakeBorder(image,4000,4000,4000,4000,cv.BORDER_REPLICATE) + + dimensions,center_indexes = get_image_blocks(image, + tile_distance=tile_distance, + tile_size=tile_size + ) + try: + + segmentation = segment_tiles(unet, + center_indexes, + image, + scaling_factor=scaling_factor, + tile_size=tile_size) + + except Exception as e: + print(file) + + cv.imwrite( + file_id + + f'_uNetSegmentation.png', + segmentation + ) + + return() + + +############################################################# + +def double_check_produced_dataset(new_directory,image_idx=0): + '''this function samples a random image from a given directory, crops off + the ground truth from the 4th layer, and displays the color image to + verify they work.''' + os.chdir(new_directory) + file_names = tf.io.gfile.glob('./*.png') + file_names = natsorted(file_names) + # pick a random image index number + if image_idx == 0: + image_idx = int(np.random.random()*len(file_names)) + else: + pass + + print(image_idx) + # reading specific file from the random index + segmentation = cv.imread(file_names[image_idx],cv.IMREAD_UNCHANGED) + # changing the color for the tile from BGR to RGB + print(file_names[image_idx]) + # plotting the images next to each other + plt.imshow(segmentation,vmin=0, vmax=6) + print(np.unique(segmentation)) + plt.show() + + +############################################################# +############################################################# +# %% +tile_size = 1024 +unet_directory = '/home/briancottle/Research/Semantic_Segmentation/dataset_shards_4/' +os.chdir(unet_directory) + +sample_data = np.zeros((1,1024,1024,3)).astype(np.int8) +unet = uNet(filter_multiplier=12) +out = unet(sample_data) +unet.summary() + +unet.load_weights('./unet_seg_weights.49-0.52-0.94-0.92.h5') + +# %% + + + +JPG_directory = '/var/confocaldata/HumanNodal/HeartData/14/01/JPG/' + +segment_directory(JPG_directory, + unet, + tile_size=tile_size, + tile_distance=512, + scaling_factor=1, + HeartID='14', + ) + + +# %% + +# double_check_produced_dataset('/var/confocaldata/HumanNodal/HeartData/08/02/uNet_Segmentations', +# image_idx=0) + +# %%