{"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":"# CIS 3115 - Learning Notebook for Unit 8\n\nThis is a notebook for my undergrad CIS 3115 Machine Learning students. This provides an introduction to encoder-decoder models and the U-net model. As we will see in this unit, these are generative models that: \n\n1. Take an image as input\n1. Encodes that image into some internal represenation\n1. Then decodes the internal representation into an output image\n\n## This is a work in progress... \n- Version 1.8 -- Debugging -- Still generating null predictions\n- Version 1.7 -- Brought back call backs, spefically learning rate reduction, but that did not help\n- Version 1.6 -- Only select random windows with text in them--set sum threshold on output window between 0 and 10,000\n- Version 1.5 -- Some clean up \n- Version 1.4 -- Prediction basically working but network converges on predicting zero's for every pixel\n - TODO -- Add threshold to prediction, remove padding, and generate RLE\n - TODO -- Look into why network learns to prediction all zeros, possibly train more on interesting windows where sum > threshold\n- Version 1.3 -- Starting predictions\n- Version 1.2 -- Fixed issue with 16 bit pixels on tiff input channels\n - The Tiff input images seem to be 16 bit so the max is 65025 instead of 255. Need to verify this on Kaggle\n - This fix did not improve performance, still maxed after 1 epoch at \"loss: 0.1629 - dice_coef: 0.8371\"\n- Version 1.1 -- Issue with input channels not getting normalized from 0-255 to 0-1, but IR channel is normalized\n - TODO -- Fix issue with normalizing input channels\n- Version 1.1 -- Added IR image as a input channel temporarily -- little improvement\n - Results with IR \"loss: 0.1480 - dice_coef: 0.8520\" \n - Results withour IR \"loss: 0.1299 - dice_coef: 0.8701\"\n- Version 1.0 -- U-net model with Dice loss, but not training well\n - TODO -- try adding IR layer to see if that helps convere--temp fix since IR layer not available in testing -- Does not seem to help\n - TODO -- try more channels -- 20 channels fit in RAM fine\n - TODO -- check if Dice loss is accurate -- Updated and now I think it is good\n- Version 0.9 -- Cleaned up version\n- Version 0.7 -- Custom Dice loss function, training on 10 random channels of Fragment 1\n- Version 0.6 -- Training on 10 channels of Fragment 1, wrong loss function\n- Version 0.3 -- Trained on multiple channels of Fragment 1, wrong loss function\n- Version 0.2 -- Basic U-net model trains on a single channel of Fragment 1, wrong loss function","metadata":{}},{"cell_type":"markdown","source":"# Part 0 - Setup","metadata":{}},{"cell_type":"code","source":"# First, we'll import pandas and numpy, two data processing libraries\nimport pandas as pd\nimport numpy as np\n\n# We'll also import seaborn and matplot, twp Python graphing libraries\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n# Import the needed sklearn libraries\nfrom sklearn.preprocessing import MinMaxScaler, StandardScaler\nfrom sklearn import datasets\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.neighbors import KNeighborsClassifier\nfrom sklearn.svm import SVC\nfrom sklearn.decomposition import PCA\nfrom sklearn.preprocessing import LabelEncoder\n\n# The Keras library provides support for neural networks and deep learning\n# Use the updated Keras library from Tensorflow -- provides support for neural networks and deep learning\nimport tensorflow as tf\nimport tensorflow.keras as keras\nfrom keras import backend as K # TODO -- check if these should be tensorflow.keras\nfrom keras.utils import Sequence\nfrom tensorflow.keras import layers\nfrom tensorflow.keras.models import Sequential, Model\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator\nfrom tensorflow.keras.layers import Input, Dense, Dropout, Activation, Lambda, Flatten, LSTM\nfrom tensorflow.keras.layers import Conv2D, Convolution2D, Conv2DTranspose\nfrom tensorflow.keras.layers import MaxPooling2D, AveragePooling2D, GlobalAveragePooling2D, concatenate\nfrom tensorflow.keras.optimizers import Adam, RMSprop\nfrom tensorflow.keras.utils import to_categorical\n#from keras.utils import np_utils\n\n# Other libraries we need\nfrom pathlib import Path\nfrom tqdm.notebook import tqdm, trange\nfrom PIL import Image\nimport random\nimport os\n\nprint (\"All libraries imported\")","metadata":{"execution":{"iopub.status.busy":"2023-03-19T00:43:34.476722Z","iopub.execute_input":"2023-03-19T00:43:34.477257Z","iopub.status.idle":"2023-03-19T00:43:46.998362Z","shell.execute_reply.started":"2023-03-19T00:43:34.477209Z","shell.execute_reply":"2023-03-19T00:43:46.996826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"base_path = Path(\"/kaggle/input/vesuvius-challenge-ink-detection/\")\nTRAIN_PATH = base_path / \"train\"\nall_fragments = sorted([f.name for f in TRAIN_PATH.iterdir()])\nprint(\"All fragments:\", all_fragments)\n# Due to limited memory on Kaggle, we can only load 1 full fragment\ntrain_fragments = [TRAIN_PATH / fragment_name for fragment_name in [\"1\"]]\ntrain_fragments\n\nimage_id = \"1\"\nfrag1 = TRAIN_PATH / image_id / \"surface_volume\" \nfrag1_layer1 = TRAIN_PATH / image_id/ \"surface_volume\" / \"00.tif\"\noutput_image_file = TRAIN_PATH / image_id /\"inklabels.png\"\nir_image_file = TRAIN_PATH / image_id /\"ir.png\"\n\nprint (\"frag1_layer1 = \",frag1_layer1 )\n#TRAIN_PATH = Path(\"/kaggle/input/vesuvius-challenge-ink-detection/train\")\n","metadata":{"execution":{"iopub.status.busy":"2023-03-19T00:43:47.000672Z","iopub.execute_input":"2023-03-19T00:43:47.001975Z","iopub.status.idle":"2023-03-19T00:43:47.014340Z","shell.execute_reply.started":"2023-03-19T00:43:47.001927Z","shell.execute_reply":"2023-03-19T00:43:47.012753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Part 1 - Visualization","metadata":{}},{"cell_type":"markdown","source":"This code is take from [ENRIC DOMINGO's Visualization of files notebook](https://www.kaggle.com/code/edomingo/visualization-of-files/notebook)\n","metadata":{}},{"cell_type":"code","source":"# Plotting the first layer, the Ink Mask and IR images\ndef displayImage(image_id):\n plt.figure(figsize=(15,6))\n print(f\"Displaying image {image_id}\")\n\n img = Image.open(TRAIN_PATH / image_id / \"surface_volume\" / \"00.tif\")\n print(f\"Surface scanned image size: ({img.size[0]}, {img.size[1]})\")\n plt.subplot(1, 3, 1)\n plt.title(\"First layer\")\n plt.imshow(img)\n\n img_labels = Image.open(TRAIN_PATH / image_id /\"inklabels.png\")\n print(f\"Ink Labels Mask image size: ({img.size[0]}, {img.size[1]})\")\n plt.subplot(1, 3, 2)\n plt.imshow(img_labels)\n plt.title(\"Ink Labels Mask\")\n\n img_ir = Image.open(TRAIN_PATH / image_id / \"ir.png\")\n print(f\"IR image size: ({img.size[0]}, {img.size[1]})\")\n plt.subplot(1, 3, 3)\n plt.imshow(img_ir)\n plt.title(\"IR Image\")\n \ndisplayImage(\"1\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting the second fragment -- first layer, the Ink Mask and IR images\ndisplayImage(\"2\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting the third fragment -- first layer, the Ink Mask and IR images\ndisplayImage(\"3\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting every 5 images (layers) of the first fagment:\n\nfig, axs = plt.subplots(nrows=4, ncols=3, figsize=(12,16))\n\nfor i, ax in tqdm(enumerate(axs.flat)):\n img = Image.open(TRAIN_PATH / \"1\" / \"surface_volume\" / f\"{(i*5):02d}.tif\")\n ax.imshow(img)\n ax.set_title(f\"{(i*5):02d}.tif\")\n \nfig.tight_layout()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Part 2 - Preprocessing Image","metadata":{}},{"cell_type":"markdown","source":"From FRANCOIS LEMARCHAND's High-res samples into multi-input CNN (Keras) \nhttps://www.kaggle.com/code/frlemarchand/high-res-samples-into-multi-input-cnn-keras\n\n## Custom input generator\n\nWe create a custom generator that will randomly take samples from each full-scale image. It ensures there is actually some kind of content in the returned sample and give the associated groundtruth.\n\n\n\nImage from https://www.kaggle.com/code/jpposma/vesuvius-challenge-ink-detection-tutorial","metadata":{}},{"cell_type":"code","source":"# Declare configuaration variables \n\nnum_sub_images = 256 # number of images to process in each epoch\nbatch_size = 16 # number of images in each batch\nwindow_image_size = 101 # size of input images -- moving window on original fragment\n#window_image_size = 128 # size of input images -- moving window on original fragment\n#input_channels = 65\ninput_channels = 20 # number of input channels to read in -- depth of input \n#input_channels = 10 +1 # emp for testing with IR Channel\nsum_threshold = 5000 # threshold for the moving window, values between 0 and 10,000. Higher values should limit training to only windows that have text","metadata":{"execution":{"iopub.status.busy":"2023-03-19T00:43:49.744259Z","iopub.execute_input":"2023-03-19T00:43:49.744719Z","iopub.status.idle":"2023-03-19T00:43:49.751262Z","shell.execute_reply.started":"2023-03-19T00:43:49.744676Z","shell.execute_reply":"2023-03-19T00:43:49.750097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# These functions were generated using ChapGPT and then modified by me\n\n# read each input image into a different channel\ndef read_images(input_image_path, output_image_path, window_image_size):\n output_img = Image.open(output_image_path)\n output_arr = np.array(output_img)\n\n # Get the dimensions of the images\n width, height = output_img.size\n \n # Declare empty array to hold the image channels\n sublayers = np.array([]).reshape(height,width,0)\n\n # Calculate the maximum coordinates for the top-left corner of the sub-image window\n max_x = width - window_image_size\n max_y = height - window_image_size\n \n # Read in all the images into different channels of a numpy array\n for image_file in os.listdir(input_image_path)[0:input_channels]: # only do first 10 images for now to speed up testing --- TOM TODO fix this\n #for image_file in os.listdir(input_image_path): # loading all channels overloads memory\n image_path = os.path.join(input_image_path, image_file)\n\n # Load the image using PIL\n chan_img = Image.open(image_path)\n\n # Convert the image to a NumPy array and normalize its pixel values\n chan_array = np.array(chan_img) \n #chan_array = np.array(chan_img) / 255.0\n chan_array = np.array(chan_img) / 65535.0 # This TIFF file inputs seem to be 16 bits so the max is 65535 instead of 255\n # maybe we should be using Keras function, preprocess_input(), for this\n print(\"Reading image channel of shape = \",chan_array.shape, \" at \", image_file)\n\n # Add the sub-image to the list of sub-layers\n sublayers = np.dstack((sublayers, chan_array))\n \n # Temp test --- add IR image as a channel to test convergence. Remove this later since test data does not have IR\n #ir_img = Image.open(ir_image_file)\n #ir_arr = np.array(ir_img) / 255.0 # The IR image is only 8 bits so the max is 255\n #sublayers = np.dstack((sublayers, ir_arr))\n \n # Convert the list of sub-layers to a NumPy array\n sublayers_arr = np.array(sublayers)\n \n print(\"sublayers_arr shape = \",sublayers_arr.shape)\n print(\"output_arr shape = \",output_arr.shape)\n \n return sublayers_arr, output_arr\n\n# read each input image into a different channel\ndef get_matching_subimages(chan_arr, out_arr, window_image_size, sum_threshold):\n \n # Get the dimensions of the images\n height, width, channels = chan_arr.shape\n #print (\"width = \", width)\n #print (\"height = \",height)\n #print (\"channels = \",channels)\n\n # Calculate the maximum coordinates for the top-left corner of the sub-image window\n max_x = width - window_image_size\n max_y = height - window_image_size \n \n # Keep selecting random windows until you get one with good text\n while True:\n # Choose a random position for the top-left corner of the window\n start_x = random.randint(0, max_x)\n start_y = random.randint(0, max_y)\n # Extract the sub-image from the input image\n\n sub_out = out_arr[start_y:start_y+window_image_size, start_x:start_x+window_image_size]\n if sub_out.sum() > sum_threshold:\n #print(\"FINAL Ground truth window sum = \",sub_out.sum())\n break\n \n # Extract the sub-image from the input image\n sub_img = chan_arr[start_y:start_y+window_image_size, start_x:start_x+window_image_size, :]\n\n\n # Return the sub-images as NumPy arrays\n return sub_img, sub_out\n \nclass CustomImageGenerator(Sequence):\n def __init__(self, chan_arr, out_arr, batch_size, sub_image_size, sum_threshold):\n self.chan_arr = chan_arr\n self.out_arr = out_arr\n self.batch_size = batch_size\n self.sub_image_size = sub_image_size\n self.sum_threshold = sum_threshold\n\n def __len__(self):\n return num_sub_images\n\n def __getitem__(self, idx):\n batch_input = []\n batch_output = []\n for _ in range(self.batch_size):\n input_img, output_img = get_matching_subimages(chan_arr, out_arr, self.sub_image_size, sum_threshold)\n batch_input.append(input_img)\n batch_output.append(output_img)\n\n # Convert the input and output batches to NumPy arrays\n batch_input_arr = np.array(batch_input)\n batch_output_arr = np.array(batch_output)\n\n # Return the batch of input and output sub-images as a tuple\n return batch_input_arr, batch_output_arr\n \n","metadata":{"execution":{"iopub.status.busy":"2023-03-19T00:43:50.209121Z","iopub.execute_input":"2023-03-19T00:43:50.209881Z","iopub.status.idle":"2023-03-19T00:43:50.227053Z","shell.execute_reply.started":"2023-03-19T00:43:50.209836Z","shell.execute_reply":"2023-03-19T00:43:50.225628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read in first fragment and 20 Input channels \nchan_arr, out_arr = read_images(frag1, output_image_file, window_image_size)\n\n# Create a custom input image generator\nimage_generator = CustomImageGenerator(chan_arr, out_arr, batch_size, window_image_size, sum_threshold)","metadata":{"execution":{"iopub.status.busy":"2023-03-19T00:43:50.739839Z","iopub.execute_input":"2023-03-19T00:43:50.740281Z","iopub.status.idle":"2023-03-19T00:45:21.943446Z","shell.execute_reply.started":"2023-03-19T00:43:50.740239Z","shell.execute_reply":"2023-03-19T00:45:21.941885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"# Part 3 - U-net version","metadata":{}},{"cell_type":"markdown","source":"Some of this code is taken from the [TGS Salt Identification Challenge](https://www.kaggle.com/competitions/tgs-salt-identification-challenge/code) where I last used u-net models\n- https://www.kaggle.com/phoenigs/u-net-dropout-augmentation-stratification\n- https://www.kaggle.com/code/tgibbons/u-net-without-resizing-images\\\n- https://www.kaggle.com/code/shaojiaxin/u-net-with-simple-resnet-blocks-v2-new-loss","metadata":{}},{"cell_type":"code","source":"# An alternative model defination from https://pyimagesearch.com/2022/02/21/u-net-image-segmentation-in-keras/\n# ----- currently NOT using this version ------\n# This model has 34.5 million parameters\n\n# First, we create a function double_conv_block with layers Conv2D-ReLU-Conv2D-ReLU, which we will use in both the encoder (or the contracting path) and the bottleneck of the U-Net.\ndef double_conv_block(x, n_filters):\n # Conv2D then ReLU activation\n x = Conv2D(n_filters, 3, padding = \"same\", activation = \"relu\", kernel_initializer = \"he_normal\")(x)\n # Conv2D then ReLU activation\n x = Conv2D(n_filters, 3, padding = \"same\", activation = \"relu\", kernel_initializer = \"he_normal\")(x)\n return x\n\n# Then we define a downsample_block function for downsampling or feature extraction to be used in the encoder.\ndef downsample_block(x, n_filters):\n f = double_conv_block(x, n_filters)\n p = layers.MaxPool2D(2)(f)\n p = layers.Dropout(0.3)(p)\n\n return f, p\n\n# Finally, we define an upsampling function upsample_block for the decoder (or expanding path) of the U-Net.\ndef upsample_block(x, conv_features, n_filters):\n # upsample\n x = Conv2DTranspose(n_filters, 3, 2, padding=\"same\")(x)\n # concatenate\n x = concatenate([x, conv_features])\n # dropout\n x = Dropout(0.3)(x)\n # Conv2D twice with ReLU activation\n x = double_conv_block(x, n_filters)\n return x\n\n# First, we create a build_unet_model function, specify the inputs, encoder layers, bottleneck, decoder layers, and \n# finally the output layer with Conv2D with activation of softmax. Note the input image shape is 128x128x3. \n# The output has three channels corresponding to the three classes that the model will classify each pixel for: \n# background, foreground object, and object outline.\n\ndef build_unet_model_128(input_channels):\n # inputs\n inputs = layers.Input(shape=(128,128,input_channels))\n\n # encoder: contracting path - downsample\n # 1 - downsample\n f1, p1 = downsample_block(inputs, 64)\n # 2 - downsample\n f2, p2 = downsample_block(p1, 128)\n # 3 - downsample\n f3, p3 = downsample_block(p2, 256)\n # 4 - downsample\n f4, p4 = downsample_block(p3, 512)\n\n # 5 - bottleneck\n bottleneck = double_conv_block(p4, 1024)\n\n # decoder: expanding path - upsample\n # 6 - upsample\n u6 = upsample_block(bottleneck, f4, 512)\n # 7 - upsample\n u7 = upsample_block(u6, f3, 256)\n # 8 - upsample\n u8 = upsample_block(u7, f2, 128)\n # 9 - upsample\n u9 = upsample_block(u8, f1, 64)\n\n # outputs\n outputs = layers.Conv2D(1, 1, padding=\"same\", activation = \"softmax\")(u9)\n\n # unet model with Keras Functional API\n unet_model = Model(inputs, outputs, name=\"U-Net\")\n\n return unet_model\n","metadata":{"execution":{"iopub.status.busy":"2023-03-19T00:45:21.946049Z","iopub.execute_input":"2023-03-19T00:45:21.946565Z","iopub.status.idle":"2023-03-19T00:45:21.961708Z","shell.execute_reply.started":"2023-03-19T00:45:21.946511Z","shell.execute_reply":"2023-03-19T00:45:21.959988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---- Currently using this version ----\n# This model has 2 million parameters\n\ndef unet_model_101(input_layer, start_neurons):\n # standard size 128 -> 64 custome size 101 -> 50\n conv1 = Conv2D(start_neurons * 1, (3, 3), activation=\"relu\", padding=\"same\")(input_layer)\n conv1 = Conv2D(start_neurons * 1, (3, 3), activation=\"relu\", padding=\"same\")(conv1)\n pool1 = MaxPooling2D((2, 2))(conv1)\n #pool1 = Dropout(0.25)(pool1)\n \n # standard size 64 -> 32 custome size 50 -> 25\n conv2 = Conv2D(start_neurons * 2, (3, 3), activation=\"relu\", padding=\"same\")(pool1)\n conv2 = Conv2D(start_neurons * 2, (3, 3), activation=\"relu\", padding=\"same\")(conv2)\n pool2 = MaxPooling2D((2, 2))(conv2)\n #pool2 = Dropout(0.5)(pool2)\n\n # standard size 32 -> 16 custome size 25 -> 12\n conv3 = Conv2D(start_neurons * 4, (3, 3), activation=\"relu\", padding=\"same\")(pool2)\n conv3 = Conv2D(start_neurons * 4, (3, 3), activation=\"relu\", padding=\"same\")(conv3)\n pool3 = MaxPooling2D((2, 2))(conv3)\n #pool3 = Dropout(0.5)(pool3)\n \n # standard size 16 -> 8 custome size 12 -> 6\n conv4 = Conv2D(start_neurons * 8, (3, 3), activation=\"relu\", padding=\"same\")(pool3)\n conv4 = Conv2D(start_neurons * 8, (3, 3), activation=\"relu\", padding=\"same\")(conv4)\n pool4 = MaxPooling2D((2, 2))(conv4)\n #pool4 = Dropout(0.5)(pool4)\n \n # Middle\n convm = Conv2D(start_neurons * 16, (3, 3), activation=\"relu\", padding=\"same\")(pool4)\n convm = Conv2D(start_neurons * 16, (3, 3), activation=\"relu\", padding=\"same\")(convm)\n \n # standard size 8 -> 16 custome size 6-> 12\n deconv4 = Conv2DTranspose(start_neurons * 8, (3, 3), strides=(2, 2), padding=\"same\")(convm)\n uconv4 = concatenate([deconv4, conv4])\n #uconv4 = Dropout(0.5)(uconv4)\n uconv4 = Conv2D(start_neurons * 8, (3, 3), activation=\"relu\", padding=\"same\")(uconv4)\n uconv4 = Conv2D(start_neurons * 8, (3, 3), activation=\"relu\", padding=\"same\")(uconv4)\n \n # standard size 16 -> 32 custome size 12 -> 25\n # Changed padding from \"same\" to \"valid\" to round up image to next size\n #deconv3a = Conv2DTranspose(start_neurons * 4, (3, 3), strides=(2, 2), padding=\"same\")(uconv4)\n deconv3 = Conv2DTranspose(start_neurons * 4, (3, 3), strides=(2, 2), padding=\"valid\")(uconv4)\n uconv3 = concatenate([deconv3, conv3])\n #uconv3 = Dropout(0.5)(uconv3)\n uconv3 = Conv2D(start_neurons * 4, (3, 3), activation=\"relu\", padding=\"same\")(uconv3)\n uconv3 = Conv2D(start_neurons * 4, (3, 3), activation=\"relu\", padding=\"same\")(uconv3)\n \n # standard size 32 -> 64 custome size 25 -> 50\n deconv2 = Conv2DTranspose(start_neurons * 2, (3, 3), strides=(2, 2), padding=\"same\")(uconv3)\n uconv2 = concatenate([deconv2, conv2])\n #uconv2 = Dropout(0.5)(uconv2)\n uconv2 = Conv2D(start_neurons * 2, (3, 3), activation=\"relu\", padding=\"same\")(uconv2)\n uconv2 = Conv2D(start_neurons * 2, (3, 3), activation=\"relu\", padding=\"same\")(uconv2)\n\n # standard size 64 -> 128 custome size 50 -> 101\n # Changed padding from \"same\" to \"valid\" to round up image to next size\n #deconv1 = Conv2DTranspose(start_neurons * 1, (3, 3), strides=(2, 2), padding=\"same\")(uconv2)\n deconv1 = Conv2DTranspose(start_neurons * 1, (3, 3), strides=(2, 2), padding=\"valid\")(uconv2)\n uconv1 = concatenate([deconv1, conv1])\n #uconv1 = Dropout(0.5)(uconv1)\n uconv1 = Conv2D(start_neurons * 1, (3, 3), activation=\"relu\", padding=\"same\")(uconv1)\n uconv1 = Conv2D(start_neurons * 1, (3, 3), activation=\"relu\", padding=\"same\")(uconv1)\n\n #uconv1 = Dropout(0.5)(uconv1)\n output_layer = Conv2D(1, (1,1), padding=\"same\", activation=\"sigmoid\")(uconv1)\n \n return output_layer\n\ndef bulid_unet_model_101(start_neurons, input_channels):\n input_layer = Input((window_image_size, window_image_size, input_channels))\n output_layer = unet_model_101(input_layer, start_neurons)\n unet_model = Model(input_layer, output_layer)\n return unet_model","metadata":{"execution":{"iopub.status.busy":"2023-03-19T00:45:21.963646Z","iopub.execute_input":"2023-03-19T00:45:21.964558Z","iopub.status.idle":"2023-03-19T00:45:22.203424Z","shell.execute_reply.started":"2023-03-19T00:45:21.964516Z","shell.execute_reply":"2023-03-19T00:45:22.202154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nunet_model = bulid_unet_model_101(16, input_channels)\n#unet_model = build_unet_model_128(input_channels)\nprint (\"U-net model created and named 'unet_model'\")\n\n","metadata":{"execution":{"iopub.status.busy":"2023-03-19T00:45:22.206654Z","iopub.execute_input":"2023-03-19T00:45:22.207327Z","iopub.status.idle":"2023-03-19T00:45:22.835241Z","shell.execute_reply.started":"2023-03-19T00:45:22.207284Z","shell.execute_reply":"2023-03-19T00:45:22.833870Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# We evaluate how well your output image matches our reference image using the Sørensen–Dice coefficient, otherwise known as the F1 score.\n# Not 100% sure this works correctly...\n\ndef dice_coef(y_true, y_pred, smooth=1):\n y_true = tf.cast(y_true, dtype=tf.float32)\n y_pred = tf.cast(y_pred, dtype=tf.float32)\n intersection = K.sum(K.abs(y_true * y_pred), axis=-1)\n # Calculate the sums\n true_sum = K.sum(y_true, -1)\n pred_sum = K.sum(y_pred, -1)\n return (2.0 * intersection + smooth) / (true_sum + pred_sum + smooth)\n\ndef dice_loss(y_true, y_pred):\n return 1 - dice_coef(y_true, y_pred)\n\nfrom tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping, ModelCheckpoint\n\nlearning_rate_reduction = ReduceLROnPlateau(monitor='loss', \n patience=4, \n verbose=2, \n factor=0.5, \n min_lr=0.000001)\n\nearly_stops = EarlyStopping(monitor='loss', \n min_delta=0, \n patience=16, \n verbose=2, \n mode='auto')\n\ncheckpointer = ModelCheckpoint(filepath = 'cis3115.{epoch:02d}-{accuracy:.6f}.hdf5',\n verbose=2,\n save_best_only=True, \n save_weights_only = True)","metadata":{"execution":{"iopub.status.busy":"2023-03-19T00:45:22.836725Z","iopub.execute_input":"2023-03-19T00:45:22.837210Z","iopub.status.idle":"2023-03-19T00:45:22.850955Z","shell.execute_reply.started":"2023-03-19T00:45:22.837166Z","shell.execute_reply":"2023-03-19T00:45:22.849393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"unet_model.compile(loss=dice_loss, optimizer=\"adam\", metrics=[dice_coef])\n\n# Try a different loss function to see if that fixed the issue with null predictions\n#unet_model.compile(loss='binary_crossentropy', optimizer=\"adam\", metrics=[dice_coef])\n\nunet_model.summary()","metadata":{"execution":{"iopub.status.busy":"2023-03-19T01:11:57.315573Z","iopub.execute_input":"2023-03-19T01:11:57.316021Z","iopub.status.idle":"2023-03-19T01:11:57.434472Z","shell.execute_reply.started":"2023-03-19T01:11:57.315975Z","shell.execute_reply":"2023-03-19T01:11:57.433104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from https://www.kaggle.com/code/frlemarchand/high-res-samples-into-multi-input-cnn-keras\n\nhistory = unet_model.fit(\n image_generator,\n steps_per_epoch = int(num_sub_images/batch_size),\n callbacks=[learning_rate_reduction, early_stops],\n epochs=1)\n","metadata":{"execution":{"iopub.status.busy":"2023-03-19T01:12:06.568898Z","iopub.execute_input":"2023-03-19T01:12:06.569307Z","iopub.status.idle":"2023-03-19T01:12:28.372259Z","shell.execute_reply.started":"2023-03-19T01:12:06.569269Z","shell.execute_reply":"2023-03-19T01:12:28.371275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# We will display the loss and the accuracy of the model for each epoch\n# NOTE: this is a little fancy display than is shown in the textbook\ndef display_training_curves(training, validation, title, subplot):\n if subplot%10==1: # set up the subplots on the first call\n plt.subplots(figsize=(10,10), facecolor='#F0F0F0')\n plt.tight_layout()\n ax = plt.subplot(subplot)\n ax.set_facecolor('#F8F8F8')\n ax.plot(training)\n ax.plot(validation)\n ax.set_title('model '+ title)\n ax.set_ylabel(title)\n #ax.set_ylim(0.28,1.05)\n ax.set_xlabel('epoch')\n ax.legend(['train', 'valid.'])\n\n#display_training_curves(history.history['loss'], history.history['val_loss'], 'loss', 211)\n#display_training_curves(history.history['accuracy'], history.history['val_accuracy'], 'accuracy', 212)\ndisplay_training_curves(history.history['loss'], history.history['loss'], 'loss', 211)","metadata":{"execution":{"iopub.status.busy":"2023-03-19T00:47:27.121297Z","iopub.execute_input":"2023-03-19T00:47:27.121757Z","iopub.status.idle":"2023-03-19T00:47:27.534057Z","shell.execute_reply.started":"2023-03-19T00:47:27.121714Z","shell.execute_reply":"2023-03-19T00:47:27.532361Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Clean up RAM by deleting the image channel variable\n#del chan_arr","metadata":{"execution":{"iopub.status.busy":"2023-03-18T18:34:29.812108Z","iopub.execute_input":"2023-03-18T18:34:29.812473Z","iopub.status.idle":"2023-03-18T18:34:29.818292Z","shell.execute_reply.started":"2023-03-18T18:34:29.812439Z","shell.execute_reply":"2023-03-18T18:34:29.816774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Debugging code in cells below\nThe following two cells are only used for some debugging","metadata":{}},{"cell_type":"code","source":"# Declare configuaration variables \n\nnum_sub_images = 256 # number of images to process in each epoch\nbatch_size = 16 # number of images in each batch\nwindow_image_size = 101 # size of input images -- moving window on original fragment\n#window_image_size = 128 # size of input images -- moving window on original fragment\n#input_channels = 65\ninput_channels = 20 # number of input channels to read in -- depth of input \n#input_channels = 10 +1 # emp for testing with IR Channel\nsum_threshold = 5000 # threshold for the moving window, values between 0 and 10,000. Higher values should limit training to only windows that have text","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---- Debugging code \nimport random\nfrom PIL import Image\nimport numpy as np\nimport os\nfrom keras.utils import Sequence\n\nfrag1_layer1 = TRAIN_PATH / image_id/ \"surface_volume\" / \"00.tif\"\nimage_path = os.path.join(frag1_layer1)\n# Load the image using PIL\nchan_img = Image.open(image_path)\n\n# Convert the image to a NumPy array and normalize its pixel values\nchan_array = np.array(chan_img) \nprint(\"Before average of array = \", np.average(chan_array))\nchan_array2 = np.array(chan_img) / 255.0\n#img_array = img_array / 255.0\nprint(\"After average of array = \", np.average(chan_array2))\nprint(\"Reading image channel of size shape = \",chan_array.shape, \" at \", image_path)\n\nstart = random.randint(0, 5000)\nprint (chan_array[start:start+10,start:start+10] )\n","metadata":{"_kg_hide-input":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def debugDisplay():\n # ---- Debugging code --- select random window and display it\n # Get the dimensions of the images\n height, width, channels = chan_arr.shape\n # Declare empty array to hold the image channels\n sublayers = np.array([]).reshape(window_image_size,window_image_size,0)\n # Calculate the maximum coordinates for the top-left corner of the sub-image window\n max_x = width - window_image_size\n max_y = height - window_image_size \n # Convert the list of sub-layers to a NumPy array\n sublayers_arr = np.array(sublayers)\n\n # Keep selecting random windows until you get one with good text\n while True:\n # Choose a random position for the top-left corner of the window\n start_x = random.randint(0, max_x)\n start_y = random.randint(0, max_y)\n # Extract the sub-image from the input image\n\n sub_out = out_arr[start_y:start_y+window_image_size, start_x:start_x+window_image_size]\n if sub_out.sum() > sum_threshold:\n #print(\"FINAL Ground truth window sum = \",sub_out.sum())\n break\n\n print (\"start_x = \",start_x,\" and start_y = \", start_y)\n print(\"chan_arr shape = \",chan_arr.shape)\n print(\"out_arr shape = \",out_arr.shape)\n # Extract the sub-image from the input image\n sub_img = chan_arr[start_y:start_y+window_image_size, start_x:start_x+window_image_size, :]\n\n #print (chan_arr[2000:2005,2000:2005] )\n #print (out_arr[2000:2005,2000:2005] )\n print(\"sub_img sum = \",sub_img.sum())\n print(\"sub_out sum = \",sub_out.sum())\n #print (sub_img[50:55,50:55] )\n print (sub_out[20:30,50:80] )\n print (sub_out[50:60,50:80] )\n\n # Remember output values are in [0,1] and we want pixel values in [0,255]\n plt.imshow(sub_out*255)\n","metadata":{"execution":{"iopub.status.busy":"2023-03-19T01:13:51.764891Z","iopub.execute_input":"2023-03-19T01:13:51.765346Z","iopub.status.idle":"2023-03-19T01:13:51.778228Z","shell.execute_reply.started":"2023-03-19T01:13:51.765308Z","shell.execute_reply":"2023-03-19T01:13:51.776490Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"debugDisplay()","metadata":{"execution":{"iopub.status.busy":"2023-03-19T01:14:33.485134Z","iopub.execute_input":"2023-03-19T01:14:33.485573Z","iopub.status.idle":"2023-03-19T01:14:33.735858Z","shell.execute_reply.started":"2023-03-19T01:14:33.485535Z","shell.execute_reply":"2023-03-19T01:14:33.734034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n#resized_sub_img = sub_img.view(np.ndarray).reshape((sub_image_size, sub_image_size, -1))\nresized_sub_img = sub_img.reshape(( 1, window_image_size, window_image_size, -1))\n#print(\"resized_sub_img shape = \",resized_sub_img.shape)\npred_img = unet_model.predict(resized_sub_img, verbose=0)\n#print(\"resized_sub_img sum = \",resized_sub_img.sum(),\"pred_img sum = \",pred_img.sum())\n\n#print(\"output shape = \",output.shape)\n#print(\"pred_img shape before = \",pred_img.shape)\npred_img = pred_img.reshape(window_image_size, window_image_size)\nprint (pred_img[20:30,50:80] )\nprint (pred_img[50:60,50:80] )\nplt.imshow(pred_img*255)","metadata":{"execution":{"iopub.status.busy":"2023-03-19T01:14:37.442012Z","iopub.execute_input":"2023-03-19T01:14:37.442956Z","iopub.status.idle":"2023-03-19T01:14:37.823401Z","shell.execute_reply.started":"2023-03-19T01:14:37.442908Z","shell.execute_reply":"2023-03-19T01:14:37.822053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Part 4 - Predictions\n\nTrying to generate predictions. \n- Use model to generate prediction patches that cover the enter canvase\n- Stitch together all these prediction patches\n- Generate RLE of image\n\n\n## TODO \n- Add threshold to prediction\n- remove padding\n- generate RLE","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_id = \"a\"\nbase_path = Path(\"/kaggle/input/vesuvius-challenge-ink-detection/\")\nTEST_PATH = base_path / \"test\"\ntest_frag1_path = TEST_PATH / test_id / \"surface_volume\" \ntest_id_layer1 = TEST_PATH / test_id/ \"surface_volume\" / \"00.tif\"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Moving window prediction\n\nMove a window through the test image, generate a prediction for each window, then stitch together all the predictions\n\n\n\nImage from https://www.kaggle.com/code/jpposma/vesuvius-challenge-ink-detection-tutorial","metadata":{}},{"cell_type":"code","source":"\n# read each input image into a different channel\ndef read_input_only(input_image_path,input_image_sample):\n sample_img = Image.open(input_image_sample)\n sample_arr = np.array(input_image_sample)\n\n # Get the dimensions of the images\n width, height = sample_img.size\n \n # Declare empty array to hold the image channels\n sublayers = np.array([]).reshape(height,width,0)\n\n # Calculate the maximum coordinates for the top-left corner of the sub-image window\n max_x = width - sub_image_size\n max_y = height - sub_image_size\n \n # Read in all the images into different channels of a numpy array\n #for image_file in os.listdir(input_image_path)[30:40]: # only do 10 images for now to speed up testing \n for image_file in os.listdir(input_image_path)[0:input_channels]: # only do first 10 images for now to speed up testing --- TOM TODO fix this\n #for image_file in os.listdir(input_image_path): # loading all channels overloads memory\n image_path = os.path.join(input_image_path, image_file)\n\n # Load the image using PIL\n chan_img = Image.open(image_path)\n\n # Convert the image to a NumPy array and normalize its pixel values\n chan_array = np.array(chan_img) \n #chan_array = np.array(chan_img) / 255.0\n chan_array = np.array(chan_img) / 65535.0 # This TIFF file inputs seem to be 16 bits so the max is 65535 instead of 255\n print(\"Reading image channel of size shape = \",chan_array.shape, \" at \", image_file)\n\n # Add the sub-image to the list of sub-layers\n sublayers = np.dstack((sublayers, chan_array))\n \n # Convert the list of sub-layers to a NumPy array\n sublayers_arr = np.array(sublayers)\n \n print(\"sublayers_arr shape = \",sublayers_arr.shape)\n \n return sublayers_arr\n\n# Generated with the help of Chat GPT\ndef get_prediction_subimages(chan_arr, input_sample_path, sub_image_size=101, stride=101):\n # open up a sample input image for sizing\n sample_img = Image.open(input_sample_path)\n window_size = sub_image_size\n #height, width, channels = chan_arr.shape\n num_rows = chan_arr.shape[0]\n num_cols = chan_arr.shape[1]\n channels = chan_arr.shape[2]\n # Calculate the number of windows we need to traverse\n num_rows_padded = (num_rows // window_size) * window_size + window_size\n num_cols_padded = (num_cols // window_size) * window_size + window_size\n #print(\"num_rows_padded = \",num_rows_padded,\" and num_cols_padded = \",num_cols_padded)\n num_windows_row = (num_rows_padded - window_size) // stride + 1\n num_windows_col = (num_cols_padded - window_size) // stride + 1\n #print(\"num_windows_row = \",num_windows_row,\" and num_windows_col = \",num_windows_col)\n\n # Declare empty array to hold the image channels\n pre_windows = np.array([]).reshape(num_rows_padded,num_cols_padded,0)\n # Create a 2D numpy array the same size as the input image, with each pixel is from the prediction\n output = window_sums = np.zeros((num_rows_padded, num_cols_padded))\n \n # Traverse the image in window_size x window_size pixel windows and calculate the sum of each window\n for i in range(num_windows_row):\n print(\"Processing Row \",i,\" of \",num_windows_row)\n for j in range(num_windows_col):\n row_start = i * stride\n row_end = row_start + window_size\n col_start = j * stride\n col_end = col_start + window_size\n #print(\"row_start:\",row_start,\" row_end:\",row_end,\" col_start:\",col_start,\" col_end:\",col_end)\n # Extract the sub-image from the input image\n sub_img = chan_arr[row_start:row_end, col_start:col_end, :]\n # Pad the original image with zeros so that the moving window covers it all\n padded_image = np.zeros((sub_image_size, sub_image_size, channels))\n padded_image[:sub_img.shape[0], :sub_img.shape[1], :] = sub_img\n #print(\"sub_img shape = \",sub_img.shape)\n #print(\"padded_image shape = \",padded_image.shape)\n\n #resized_sub_img = sub_img.view(np.ndarray).reshape((sub_image_size, sub_image_size, -1))\n resized_sub_img = padded_image.reshape(( 1, sub_image_size, sub_image_size, -1))\n #print(\"resized_sub_img shape = \",resized_sub_img.shape)\n pred_img = unet_model.predict(resized_sub_img, verbose=0)\n #print(\"resized_sub_img sum = \",resized_sub_img.sum(),\"pred_img sum = \",pred_img.sum())\n\n #print(\"output shape = \",output.shape)\n #print(\"pred_img shape before = \",pred_img.shape)\n pred_img = pred_img.reshape(sub_image_size, sub_image_size)\n # Add the sub-image to the list of sub-layers\n #print(\"pred_img shape after = \",pred_img.shape)\n output[row_start:row_end, col_start:col_end] = pred_img\n \n return output\n\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nnum_sub_images = 256\nbatch_size = 16\nsub_image_size = 101\n\nchannel_arr = read_input_only(test_frag1_path,test_id_layer1)\nprediciton_image = get_prediction_subimages(channel_arr, test_id_layer1, sub_image_size)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print (prediciton_image.shape)\nprint(\"prediciton_image sum = \",prediciton_image.sum())\nprint (prediciton_image[50:55,50:55] )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Remember output values are in [0,1] and we want pixel values in [0,255]\nplt.imshow(prediciton_image*255)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Clear up RAM by deleting this variable\ndel channel_arr","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Running Length Encoding converter provided in instructions \n# From https://gist.github.com/janpaul123/ca3477c1db6de4346affca37e0e3d5b0\n\n#import numpy as np\n#from PIL import Image\n\n# Fast run length encoding, from https://www.kaggle.com/code/hackerpoet/even-faster-run-length-encoder/script\ndef rle (img):\n flat_img = img.flatten()\n flat_img = np.where(flat_img > 0.5, 1, 0).astype(np.uint8)\n\n starts = np.array((flat_img[:-1] == 0) & (flat_img[1:] == 1))\n ends = np.array((flat_img[:-1] == 1) & (flat_img[1:] == 0))\n starts_ix = np.where(starts)[0] + 2\n ends_ix = np.where(ends)[0] + 2\n lengths = ends_ix - starts_ix\n\n return starts_ix, lengths\n\n#inklabels = np.array(Image.open('inklabels.png'), dtype=np.uint8)\n#starts_ix, lengths = rle(inklabels)\n#inklabels_rle = \" \".join(map(str, sum(zip(starts_ix, lengths), ())))\n#print(\"Id,Predicted\\n1,\" + inklabels_rle, file=open('inklabels_rle.csv', 'w'))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# From https://www.kaggle.com/competitions/vesuvius-challenge-ink-detection/overview/evaluation\n# For a real-world example of what these files look like, see inklabels_rce.csv in the data directories, \n# which have been generated with this script. We also show how to output a file of this format in the Ink Detection tutorial.\n#import numpy as np\n#from PIL import Image\n\n# Fast run length encoding, from https://www.kaggle.com/code/hackerpoet/even-faster-run-length-encoder/script\ndef rle (img):\n flat_img = img.flatten()\n flat_img = np.where(flat_img > 0.5, 1, 0).astype(np.uint8)\n\n starts = np.array((flat_img[:-1] == 0) & (flat_img[1:] == 1))\n ends = np.array((flat_img[:-1] == 1) & (flat_img[1:] == 0))\n starts_ix = np.where(starts)[0] + 2\n ends_ix = np.where(ends)[0] + 2\n lengths = ends_ix - starts_ix\n\n return starts_ix, lengths\n\n#inklabels = np.array(Image.open('inklabels.png'), dtype=np.uint8)\n#starts_ix, lengths = rle(inklabels)\n#inklabels_rle = \" \".join(map(str, sum(zip(starts_ix, lengths), ())))\n#print(\"Id,Predicted\\n1,\" + inklabels_rle, file=open('inklabels_rle.csv', 'w'))","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}