{"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":"
Understanding convolution layers
\n
Inputs, outputs and filters of a tensorflow Convolutional Neural Network
\n
for the Vesuvius Challenge - Ink Detection competition.
\n
\n \n### Notebook versions:\n\nV1: Initial\n\nV2: Included conv3d\n \n#

1. Introduction

\n\nThe Vesuvius Challenge invites us to build a model that can reconstruct inked text from 3-D X-ray scans of carbonized scrolls from Vesuvius, which are unable to be opened and read without destroying them. A convolutional neural network seems like the ideal model for such a task. However, I find the intricacies of convolution layers (e.g. `conv2d`, `conv3d`) confusing. In this notebook I hope to provide some insight into the following questions:\n\n* How can we construct a convolutional neural network with tensorflow/keras?\n* How do we align and reshape input matrices?\n* How do the input and output shapes work? Why do we need 5-D tensors for a 3-D convolution?\n* What is the difference between a 2-D and a 3-D convolution layer?\n\nThe ultimate aim is to have enough insight into convolution layers (and their friends such as pooling layers) to produce an appropriate Convolutional Neural Network (CNN) for the [Vesuvius Challenge competition](https://www.kaggle.com/competitions/vesuvius-challenge-ink-detection), as well as to be able to make informed choices on model architecture.\n\nI'm using the [tensorflow](https://www.tensorflow.org/) architecture, and in particular, I'll be looking to build a CNN using the [Keras](https://keras.io/) API. Typically, we'd use the Layers API, such as [`tf.keras.layers.Conv2d`](https://keras.io/api/layers/convolution_layers/convolution2d/) as part of a larger Keras model, but for now (in particular since we're not fitting a model to anything here, just exploring the convolution layer inputs and outputs), we're using lower-level tensorflow operations such as [`tf.nn.conv2d`](https://www.tensorflow.org/api_docs/python/tf/nn/conv2d). These allow us to compute the convolution with specified inputs and filters without having to fit a model. ","metadata":{"papermill":{"duration":0.022896,"end_time":"2023-04-16T09:23:29.441028","exception":false,"start_time":"2023-04-16T09:23:29.418132","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd\n\nimport os\n\nimport matplotlib.pyplot as plt\nimport tensorflow as tf\n\nimport matplotlib.image as mpimg\nimport matplotlib as mpl\n","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":10.824396,"end_time":"2023-04-16T09:23:40.285228","exception":false,"start_time":"2023-04-16T09:23:29.460832","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:37:41.240913Z","iopub.execute_input":"2023-04-24T13:37:41.241318Z","iopub.status.idle":"2023-04-24T13:37:51.645297Z","shell.execute_reply.started":"2023-04-24T13:37:41.241282Z","shell.execute_reply":"2023-04-24T13:37:51.643864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1.1 Input data\n\nThe input data for the Vesuvius Challenge competition is layered images representing X-ray scan depths, such as:","metadata":{"papermill":{"duration":0.019172,"end_time":"2023-04-16T09:23:40.324572","exception":false,"start_time":"2023-04-16T09:23:40.305400","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import glob\nimport torch\nimport PIL.Image as Image\n\nPREFIX = '/kaggle/input/vesuvius-challenge-ink-detection/train/1/'\nZ_START = 27 # First slice in the z direction to use\nZ_DIM = 10 # Number of slices in the z direction\nDEVICE = torch.device(\"cpu\")\n\n# Load the 3d x-ray scan, one slice at a time\nimages = [np.array(Image.open(filename), dtype=np.float32)/65535.0 for filename in sorted(glob.glob(PREFIX+\"surface_volume/*.tif\"))[Z_START:Z_START+Z_DIM]]\nimage_stack = torch.stack([torch.from_numpy(image) for image in images], dim=0).to(DEVICE)\n\nfig, axes = plt.subplots(1, len(images), figsize=(15, 3))\nfor image, ax in zip(images, axes):\n ax.imshow(np.array(Image.fromarray(image).resize((image.shape[1]//20, image.shape[0]//20)), dtype=np.float32), cmap='gray')\n ax.set_xticks([]); ax.set_yticks([])\nfig.tight_layout()\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-24T13:37:51.647489Z","iopub.execute_input":"2023-04-24T13:37:51.648243Z","iopub.status.idle":"2023-04-24T13:38:16.989527Z","shell.execute_reply.started":"2023-04-24T13:37:51.648196Z","shell.execute_reply":"2023-04-24T13:38:16.988576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"(Code for this image courtesy of [Vesuvius Challenge Tutorial Notebook](https://www.kaggle.com/code/jpposma/vesuvius-challenge-ink-detection-tutorial).) \n\nFor now, let's take a simple RGB image to play with. We can think of the RGB channels as being similar to the scan layers from the competition.","metadata":{}},{"cell_type":"code","source":"BW_fname = '/kaggle/input/protea/GS_Australian_Flower_Pink.jpg'\ncol_fname = '/kaggle/input/protea/Australian_Flower_Pink.jpg'\nBW_img = mpimg.imread(BW_fname)\ncol_img = mpimg.imread(col_fname)","metadata":{"papermill":{"duration":0.055694,"end_time":"2023-04-16T09:23:40.399912","exception":false,"start_time":"2023-04-16T09:23:40.344218","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:16.990955Z","iopub.execute_input":"2023-04-24T13:38:16.991730Z","iopub.status.idle":"2023-04-24T13:38:17.022746Z","shell.execute_reply.started":"2023-04-24T13:38:16.991685Z","shell.execute_reply":"2023-04-24T13:38:17.021666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from matplotlib.gridspec import GridSpec\nfig = plt.figure(figsize=(13*7/3/2, 17/2))\n\n_=plt.suptitle('Protea image decomposed into RGB channels')\n\ngs2 = GridSpec(1,3, width_ratios=[1, 1/3,1])\nBWax=fig.add_subplot(gs2[0])\n_=BWax.imshow(BW_img, cmap='gray')\n\ncolax=fig.add_subplot(gs2[2])\n_=colax.imshow(col_img)\n\nnorm_scale = mpl.colors.Normalize(vmin=0, vmax=400) # Just renormalize the plotting\n # colours otherwise very intense color intensities come out as black\n\ngs = GridSpec(3,3, width_ratios=[1, 1/3,1], height_ratios=[1, 1,1])\nRax = fig.add_subplot(gs[1])\n_=Rax.imshow(col_img[:,:,0],\n cmap = 'Reds',\n norm = norm_scale)\n \nGax = fig.add_subplot(gs[4])\n_=Gax.imshow(col_img[:,:,1],\n cmap = 'Greens',\n norm = norm_scale)\n\nBax = fig.add_subplot(gs[7])\n_=Bax.imshow(col_img[:,:,2],\n cmap = 'Blues',\n norm = norm_scale)\n","metadata":{"papermill":{"duration":0.957646,"end_time":"2023-04-16T09:23:41.378490","exception":false,"start_time":"2023-04-16T09:23:40.420844","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:17.025805Z","iopub.execute_input":"2023-04-24T13:38:17.026633Z","iopub.status.idle":"2023-04-24T13:38:17.880160Z","shell.execute_reply.started":"2023-04-24T13:38:17.026578Z","shell.execute_reply":"2023-04-24T13:38:17.878966Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our sample image is a small photo of a protea, taken from [Wikimedia Commons](https://commons.wikimedia.org/wiki/File:Australian_Flower_Pink.jpg) (Tammcd7, CC0, via Wikimedia Commons). The above figure shows the grayscale vesion of the figure in the left panel, the colour (RGB) version on the right, and the RGB channels in the middle. Blue is mostly used here to provide the whitish pixels on the edges of the petals and sepals of the flower. Let's start with the grayscale image, and consider multi-channel (RGB) images later.\n\n#

2. Two-dimensional convolution

\n\nA convolution can be thought of as a fancy (convoluted?) matrix multiplication. We align filters and inputs, and then the convolution operation itself is just element-wise matrix multiplication and summing. Actually, it might help to [think of convolution](https://stackoverflow.com/questions/43086557/convolve2d-just-by-using-numpy) as essentially a particular [`np.einsum`](https://numpy.org/doc/stable/reference/generated/numpy.einsum.html) operation. Then again, it might help to not think of it this way as well 🤔.\n\nWe begin by defining filters. In normal usage, the filters are generated automatically when we fit the CNN, but for illustrative purposes, we can define four filters that pick out vertical lines, horizontal lines and diagonal lines.","metadata":{"papermill":{"duration":0.030287,"end_time":"2023-04-16T09:23:41.439684","exception":false,"start_time":"2023-04-16T09:23:41.409397","status":"completed"},"tags":[]}},{"cell_type":"code","source":"filter_size=8\ndef create_filters(nchannels, filter_size):\n filters = np.zeros(shape=(filter_size,) * 2 + (nchannels,4), dtype=np.float32)\n filters[:,2,:,0]=1\n filters[:,filter_size-3,:,0]=-1\n filters[2,:,:,1] = 1\n filters[filter_size-3,:,:,1] = -1\n diag_offset = 5\n A = np.diagflat([-1,]*(filter_size-diag_offset), diag_offset) \n B = np.diagflat([-1,]*(filter_size-diag_offset), -diag_offset)\n filters[:,:,:,2] = np.reshape(np.eye(filter_size) + A + B, (filter_size,)*2 + (1,))\n filters[:,:,:,3] = np.reshape(np.fliplr(np.eye(filter_size)) + \n np.fliplr(A) + np.fliplr(B), \n (filter_size,)*2 + (1,))\n return filters\nfilters = create_filters(nchannels = 1, filter_size=filter_size)\n\nplt.figure(figsize=(12,8/4))\nfor i in range(4):\n plt.subplot(1,4,i+1)\n plt.imshow(np.reshape(filters[:,:,:,i],(filter_size,)*2),cmap='gray')\n","metadata":{"papermill":{"duration":0.43312,"end_time":"2023-04-16T09:23:41.903443","exception":false,"start_time":"2023-04-16T09:23:41.470323","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:17.881629Z","iopub.execute_input":"2023-04-24T13:38:17.882009Z","iopub.status.idle":"2023-04-24T13:38:18.353090Z","shell.execute_reply.started":"2023-04-24T13:38:17.881973Z","shell.execute_reply":"2023-04-24T13:38:18.351810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Intuitively, a 2-D convolution operation just slides a filter around in two dimensions (hence the name), and creates a new cell in the output that is a convolution (multiply and sum) of the filter with sections of the image:","metadata":{"papermill":{"duration":0.03002,"end_time":"2023-04-16T09:23:41.964031","exception":false,"start_time":"2023-04-16T09:23:41.934011","status":"completed"},"tags":[]}},{"cell_type":"code","source":"_=plt.imshow(mpimg.imread('/kaggle/input/protea/conv2d_with_padding.png'))\n_=plt.gca().axis('off')","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.40468,"end_time":"2023-04-16T09:23:42.399545","exception":false,"start_time":"2023-04-16T09:23:41.994865","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:18.354578Z","iopub.execute_input":"2023-04-24T13:38:18.355076Z","iopub.status.idle":"2023-04-24T13:38:18.716611Z","shell.execute_reply.started":"2023-04-24T13:38:18.355037Z","shell.execute_reply":"2023-04-24T13:38:18.715489Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the figure above, the image is padded so the output shape is the same as the input shape, and stride is equal to one (i.e. the filter is moved around with overlap).\n\nThe first filter will light up when a 8x8 part of the image contains a vertical line, and so on.\n\nThe [tensorflow docs](https://www.tensorflow.org/api_docs/python/tf/nn/conv2d) describe the required dimensions of the image and filters as follows:\n\n>The input tensor may have rank 4 or higher, where shape dimensions [:-3] are considered batch dimensions (batch_shape).\n>\n>Given an input tensor of shape `batch_shape + [in_height, in_width, in_channels`] and a filter / kernel tensor of shape `[filter_height, filter_width, in_channels, out_channels]`, this op performs the following:\n>\n> 1. Flattens the filter to a 2-D matrix with shape `[filter_height * filter_width * in_channels, output_channels]`.\n>\n> 2. Extracts image patches from the input tensor to form a virtual tensor of shape `[batch, out_height, out_width, filter_height * filter_width * in_channels]`.\n>\n> 3. For each patch, right-multiplies the filter matrix and the image patch vector.\n\nTensorflow likes float32 `dtypes` rather than int8, which is the native `dtype` for images read in by `imgread`, so convert the image `dtype`, and reshape the image into a 4-D tensor\n\n# 2.1 `conv2d` on a grayscale image","metadata":{"papermill":{"duration":0.035486,"end_time":"2023-04-16T09:23:42.472645","exception":false,"start_time":"2023-04-16T09:23:42.437159","status":"completed"},"tags":[]}},{"cell_type":"code","source":"BW_img=np.array(BW_img, dtype=np.float32)\nBW_img_tensor = np.reshape(BW_img,(1,) + BW_img.shape + (1,))\n\nBW_img_tensor.shape","metadata":{"papermill":{"duration":0.048282,"end_time":"2023-04-16T09:23:42.556825","exception":false,"start_time":"2023-04-16T09:23:42.508543","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:18.718480Z","iopub.execute_input":"2023-04-24T13:38:18.719293Z","iopub.status.idle":"2023-04-24T13:38:18.728282Z","shell.execute_reply.started":"2023-04-24T13:38:18.719239Z","shell.execute_reply":"2023-04-24T13:38:18.727059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The shape of the image tensor is `(BATCH_SIZE, HEIGHT, WIDTH, NCHANNELS)`. As mentioned above, we can drop the `BATCH_SIZE` dimension, but it is useful to think of the inputs as 4-D tensors since, when we get to building models, the batch size (usually greater than one) is more important. The number of channels in this case is 1 as we have a grayscale image.\n\nAs specified in the tensorflow docs, the shape of the filter tensor is `(FILTER_HEIGHT, FILTER_WIDTH, IN_CHANNELS, OUT_CHANNELS)`. `IN_CHANNELS` is the grayscale channel, and `OUT_CHANNELS` is really the number of filters (which we've set as four).","metadata":{"execution":{"iopub.execute_input":"2023-04-16T00:56:30.462574Z","iopub.status.busy":"2023-04-16T00:56:30.461376Z","iopub.status.idle":"2023-04-16T00:56:30.488911Z","shell.execute_reply":"2023-04-16T00:56:30.487218Z","shell.execute_reply.started":"2023-04-16T00:56:30.462524Z"},"papermill":{"duration":0.03586,"end_time":"2023-04-16T09:23:42.628733","exception":false,"start_time":"2023-04-16T09:23:42.592873","status":"completed"},"tags":[]}},{"cell_type":"code","source":"filters.shape","metadata":{"papermill":{"duration":0.046691,"end_time":"2023-04-16T09:23:42.711027","exception":false,"start_time":"2023-04-16T09:23:42.664336","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:18.730144Z","iopub.execute_input":"2023-04-24T13:38:18.730791Z","iopub.status.idle":"2023-04-24T13:38:18.746303Z","shell.execute_reply.started":"2023-04-24T13:38:18.730734Z","shell.execute_reply":"2023-04-24T13:38:18.744909Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, we are ready to perform the 2-D convolution operation:","metadata":{"papermill":{"duration":0.035698,"end_time":"2023-04-16T09:23:42.782373","exception":false,"start_time":"2023-04-16T09:23:42.746675","status":"completed"},"tags":[]}},{"cell_type":"code","source":"output = tf.nn.conv2d(BW_img_tensor, \n filters, \n strides=1,\n padding='SAME')\noutput.shape","metadata":{"papermill":{"duration":0.253316,"end_time":"2023-04-16T09:23:43.072252","exception":false,"start_time":"2023-04-16T09:23:42.818936","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:18.748217Z","iopub.execute_input":"2023-04-24T13:38:18.748566Z","iopub.status.idle":"2023-04-24T13:38:18.947783Z","shell.execute_reply.started":"2023-04-24T13:38:18.748532Z","shell.execute_reply":"2023-04-24T13:38:18.946842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The dimensions of the output are `(BATCH_SIZE, HEIGHT, WIDTH, NFILTERS)`. Altering `strides` and/or `padding` will change the height and width of the output (see below). Let's look at the output from the convolution.","metadata":{"papermill":{"duration":0.035825,"end_time":"2023-04-16T09:23:43.143690","exception":false,"start_time":"2023-04-16T09:23:43.107865","status":"completed"},"tags":[]}},{"cell_type":"code","source":"fig=plt.figure(figsize=(12,2*8/4))\n\ngs3 = GridSpec(2,4, \n height_ratios=[1/3,1])\nfor i in range(4):\n _=fig.add_subplot(gs3[i]).imshow(np.reshape(filters[:,:,:,i],(filter_size,)*2),cmap='gray')\n _=fig.add_subplot(gs3[4+i]).imshow(output[0,:,:,i], cmap='gray')\n","metadata":{"papermill":{"duration":1.060336,"end_time":"2023-04-16T09:23:44.240032","exception":false,"start_time":"2023-04-16T09:23:43.179696","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:18.953238Z","iopub.execute_input":"2023-04-24T13:38:18.953615Z","iopub.status.idle":"2023-04-24T13:38:19.886638Z","shell.execute_reply.started":"2023-04-24T13:38:18.953577Z","shell.execute_reply":"2023-04-24T13:38:19.885158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The vertical filter picks out the vertical leaves, the horizontal filter mostly picks out the edges of the leaves, and the first diagonal filter mostly light ups in the top-left corner of the image. The difference between the two diagonal filters can clearly be seen here.\n\n\n## 2.2 Changing `padding` and `strides`\n\nIf we look closely, we can see artifacts on the edges of the image due to the padding. This is the price we pay to get the convolution output shape equal to the image shape. Removing the padding would result in a smaller output shape. For example, if we take the output from the first filter using \"VALID\" padding, we remove the artifacts on the vertical sides of the image, but get a smaller output shape:","metadata":{"papermill":{"duration":0.039458,"end_time":"2023-04-16T09:23:44.320010","exception":false,"start_time":"2023-04-16T09:23:44.280552","status":"completed"},"tags":[]}},{"cell_type":"code","source":"plt.figure(figsize=(17/3,13/3))\nvalid_padding_output = tf.nn.conv2d(BW_img_tensor, \n filters, \n strides=1,\n padding='VALID')\n_=plt.imshow(valid_padding_output[0,:,:,0],cmap='gray')\nvalid_padding_output.shape","metadata":{"papermill":{"duration":0.330312,"end_time":"2023-04-16T09:23:44.692260","exception":false,"start_time":"2023-04-16T09:23:44.361948","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:19.888269Z","iopub.execute_input":"2023-04-24T13:38:19.888624Z","iopub.status.idle":"2023-04-24T13:38:20.174646Z","shell.execute_reply.started":"2023-04-24T13:38:19.888589Z","shell.execute_reply":"2023-04-24T13:38:20.172651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Increasing the stride length decreases (or eliminates) the overlap of the filter, and likewise produces a smaller output. This can be useful for dimensionality reduction.","metadata":{"papermill":{"duration":0.0417,"end_time":"2023-04-16T09:23:44.778451","exception":false,"start_time":"2023-04-16T09:23:44.736751","status":"completed"},"tags":[]}},{"cell_type":"code","source":"plt.figure(figsize=(17/3.5,12/3.5))\nstride_3_output = tf.nn.conv2d(BW_img_tensor, \n filters, \n strides=3,\n padding='SAME')\n_=plt.imshow(stride_3_output[0,:,:,0],cmap='gray')\nstride_3_output.shape","metadata":{"papermill":{"duration":0.257163,"end_time":"2023-04-16T09:23:45.077191","exception":false,"start_time":"2023-04-16T09:23:44.820028","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:20.176808Z","iopub.execute_input":"2023-04-24T13:38:20.178726Z","iopub.status.idle":"2023-04-24T13:38:20.397093Z","shell.execute_reply.started":"2023-04-24T13:38:20.178665Z","shell.execute_reply":"2023-04-24T13:38:20.395810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.3 Max pooling\n\nLet's have a look now at max pooling layers. Often, convolution layers are followed by pooling layers. Pooling applies a filter over the output layers, and extracts the maximum value within the filter. As an example, consider a 2x2 pooling layer, with stride equal to one (`None`):","metadata":{"papermill":{"duration":0.042801,"end_time":"2023-04-16T09:23:45.162562","exception":false,"start_time":"2023-04-16T09:23:45.119761","status":"completed"},"tags":[]}},{"cell_type":"code","source":"pooled_output = tf.nn.max_pool2d(output,\n ksize=5,strides=None,padding='SAME')\npooled_output.shape","metadata":{"papermill":{"duration":0.069081,"end_time":"2023-04-16T09:23:45.274372","exception":false,"start_time":"2023-04-16T09:23:45.205291","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:20.398504Z","iopub.execute_input":"2023-04-24T13:38:20.399827Z","iopub.status.idle":"2023-04-24T13:38:20.420474Z","shell.execute_reply.started":"2023-04-24T13:38:20.399749Z","shell.execute_reply":"2023-04-24T13:38:20.419100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(17/3,13/3))\n_=plt.imshow(pooled_output[0,:,:,0],cmap='gray')","metadata":{"papermill":{"duration":0.305317,"end_time":"2023-04-16T09:23:45.622389","exception":false,"start_time":"2023-04-16T09:23:45.317072","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:20.422168Z","iopub.execute_input":"2023-04-24T13:38:20.423123Z","iopub.status.idle":"2023-04-24T13:38:20.683997Z","shell.execute_reply.started":"2023-04-24T13:38:20.423074Z","shell.execute_reply":"2023-04-24T13:38:20.682727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This looks like a blurred version of the output from the vertical filter. Typically, the stride is increased in order to reduce the shape of the output layers, and to subsequently reduce the number of parameters in subsequent layers:","metadata":{"papermill":{"duration":0.043043,"end_time":"2023-04-16T09:23:45.711076","exception":false,"start_time":"2023-04-16T09:23:45.668033","status":"completed"},"tags":[]}},{"cell_type":"code","source":"_=plt.figure(figsize=(17/3.5,12/3.5))\nstride_5_pooled_output = tf.nn.max_pool2d(output,\n ksize=5,strides=5,padding='SAME')\n_=plt.imshow(stride_5_pooled_output[0,:,:,0],\n cmap='gray')","metadata":{"papermill":{"duration":0.255275,"end_time":"2023-04-16T09:23:46.009916","exception":false,"start_time":"2023-04-16T09:23:45.754641","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:20.685894Z","iopub.execute_input":"2023-04-24T13:38:20.686672Z","iopub.status.idle":"2023-04-24T13:38:20.900613Z","shell.execute_reply.started":"2023-04-24T13:38:20.686619Z","shell.execute_reply":"2023-04-24T13:38:20.899177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.4 `conv2d` on a color image\n\nLet's now turn to the color image which differs from a grayscale in having an extra dimension (3 RGB channels in addition to the height and width of the image).\n\n\n### What is a channel?\n\nThe natural example of channels on an image are the RGB channels. In subsequent layers, the number of channels is equal to the number of filters specified in previous layers. In this way, we can think of the 3 RGB channels as the output of a layer that just picks out the intensities of red, green, and blue pixels from an image. \n\nFor 3-dimensional stacks of images, we can consider this as having a height, width and depth of the stack, and 1 channel; or alternatively, the images having height and width, a depth equal to one, and a number of channels. Of course, we can have 3-dimensional stacks of images (height, width and depth) with multiple channels as well. The canonical example of this is a colour \"movie\", where the depth dimension is the time, and the RGB channels are represented as channels.\n\nThe main difference between depth and channel dimensions is that channels are summed after convolution see [this article](https://towardsdatascience.com/a-comprehensive-introduction-to-different-types-of-convolutions-in-deep-learning-669281e58215) (or [this stack overflow question](https://stackoverflow.com/q/37095783)) for a nice illustration of this.","metadata":{"papermill":{"duration":0.044392,"end_time":"2023-04-16T09:23:46.099276","exception":false,"start_time":"2023-04-16T09:23:46.054884","status":"completed"},"tags":[]}},{"cell_type":"code","source":"fig = plt.figure(figsize=(13*7/3/2, 17/2))\n\n_=plt.suptitle('Protea image decomposed into RGB channels')\n\ngs2 = GridSpec(1,3, width_ratios=[1, 1/3,1])\nBWax=fig.add_subplot(gs2[0])\n_=BWax.imshow(BW_img, cmap='gray')\n\ncolax=fig.add_subplot(gs2[2])\n_=colax.imshow(col_img)\n\nnorm_scale = mpl.colors.Normalize(vmin=0, vmax=400) # Just renormalize the plotting\n # colours otherwise very intense color intensities come out as black\n\ngs = GridSpec(3,3, width_ratios=[1, 1/3,1], height_ratios=[1, 1,1])\nRax = fig.add_subplot(gs[1])\n_=Rax.imshow(col_img[:,:,0],\n cmap = 'Reds',\n norm = norm_scale)\n \nGax = fig.add_subplot(gs[4])\n_=Gax.imshow(col_img[:,:,1],\n cmap = 'Greens',\n norm = norm_scale)\n\nBax = fig.add_subplot(gs[7])\n_=Bax.imshow(col_img[:,:,2],\n cmap = 'Blues',\n norm = norm_scale)\n","metadata":{"papermill":{"duration":0.897381,"end_time":"2023-04-16T09:23:47.040647","exception":false,"start_time":"2023-04-16T09:23:46.143266","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:20.902254Z","iopub.execute_input":"2023-04-24T13:38:20.902749Z","iopub.status.idle":"2023-04-24T13:38:22.018375Z","shell.execute_reply.started":"2023-04-24T13:38:20.902698Z","shell.execute_reply":"2023-04-24T13:38:22.017254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Once again, convert the image `dtype`:","metadata":{"papermill":{"duration":0.052173,"end_time":"2023-04-16T09:23:47.148213","exception":false,"start_time":"2023-04-16T09:23:47.096040","status":"completed"},"tags":[]}},{"cell_type":"code","source":"col_img = np.array(col_img, dtype=np.float32)\ncol_img.shape","metadata":{"papermill":{"duration":0.067494,"end_time":"2023-04-16T09:23:47.269123","exception":false,"start_time":"2023-04-16T09:23:47.201629","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:22.019802Z","iopub.execute_input":"2023-04-24T13:38:22.020947Z","iopub.status.idle":"2023-04-24T13:38:22.029153Z","shell.execute_reply.started":"2023-04-24T13:38:22.020904Z","shell.execute_reply":"2023-04-24T13:38:22.027950Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"filters_3channels = create_filters(nchannels = 3, filter_size=filter_size)\nfilters_3channels.shape","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:22.030706Z","iopub.execute_input":"2023-04-24T13:38:22.031174Z","iopub.status.idle":"2023-04-24T13:38:22.040647Z","shell.execute_reply.started":"2023-04-24T13:38:22.031136Z","shell.execute_reply":"2023-04-24T13:38:22.039384Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we have 12 \"filters\", being the four 2-D (8x8) filters for each of the three RGB channels.","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[]}},{"cell_type":"code","source":"col_img.shape","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:22.042192Z","iopub.execute_input":"2023-04-24T13:38:22.042504Z","iopub.status.idle":"2023-04-24T13:38:22.050773Z","shell.execute_reply.started":"2023-04-24T13:38:22.042474Z","shell.execute_reply":"2023-04-24T13:38:22.049688Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"col_img_tensor = np.reshape(col_img, \n (1,) + col_img.shape)\ncol_img_tensor.shape","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:22.052343Z","iopub.execute_input":"2023-04-24T13:38:22.053150Z","iopub.status.idle":"2023-04-24T13:38:22.063075Z","shell.execute_reply.started":"2023-04-24T13:38:22.053100Z","shell.execute_reply":"2023-04-24T13:38:22.062105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"output2 = tf.nn.conv2d(col_img_tensor,\n filters_3channels, \n strides=1,\n padding='SAME')\noutput2.shape","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:22.064412Z","iopub.execute_input":"2023-04-24T13:38:22.064732Z","iopub.status.idle":"2023-04-24T13:38:22.085814Z","shell.execute_reply.started":"2023-04-24T13:38:22.064703Z","shell.execute_reply":"2023-04-24T13:38:22.084576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Notice here that the shape of the output is (BATCH_SIZE, HEIGHT, WIDTH, NFILTERS), where NFILTERS or equivalently NCHANNELS. The input channels get summed in the conv2d operation, as described above.\n","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(12,2*8/4))\nfor i in range(4):\n plt.subplot(1,4,i+1)\n plt.imshow(output2[0,:,:,i], cmap='gray')","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:22.087220Z","iopub.execute_input":"2023-04-24T13:38:22.087661Z","iopub.status.idle":"2023-04-24T13:38:22.774344Z","shell.execute_reply.started":"2023-04-24T13:38:22.087625Z","shell.execute_reply":"2023-04-24T13:38:22.772901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This looks very similar to the grayscale convolution we performed before. Actually, it is identical, since we used the same filters for each RGB channel, and the 2-D convolution just sums up over the channels. We can get a different result if we change the filters. Let's change the vertical green channel filter:","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[]}},{"cell_type":"code","source":"modfilters_3channels = filters_3channels.copy()\nmodfilters_3channels[:,:,1,0] = np.reshape(np.random.normal(size=8*8)*5,(8,8))\n\nmodoutput2 = tf.nn.conv2d(col_img_tensor,\n modfilters_3channels, \n strides=1,\n padding='SAME')\nmodoutput2.shape","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:22.775898Z","iopub.execute_input":"2023-04-24T13:38:22.776700Z","iopub.status.idle":"2023-04-24T13:38:22.792688Z","shell.execute_reply.started":"2023-04-24T13:38:22.776654Z","shell.execute_reply":"2023-04-24T13:38:22.791453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12,2*8/4))\nfor i in range(4):\n plt.subplot(1,4,i+1)\n plt.imshow(modoutput2[0,:,:,i], cmap='gray')","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:22.794258Z","iopub.execute_input":"2023-04-24T13:38:22.795261Z","iopub.status.idle":"2023-04-24T13:38:23.477725Z","shell.execute_reply.started":"2023-04-24T13:38:22.795220Z","shell.execute_reply":"2023-04-24T13:38:23.476523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we can see that the first filter is corrupted by the random green filter layer, with the result that the output from the vertical filter looks more like the original image.","metadata":{}},{"cell_type":"markdown","source":"#

3. Three-dimensional convolution

\n\n\nWhereas a 2-D convolution takes 2-D filters and moves them around in two dimensions. A 3-D convolution layer moves 3-D filters in three dimensions. We need to add one extra dimension (depth) to the image and filters. \n\nFrom the [Keras documentation](https://keras.io/api/layers/convolution_layers/convolution3d/), `conv3d` takes arguments as follows:\n\n`input`: A Tensor. Shape: `[batch, in_depth, in_height, in_width, in_channels]`.\n\n`filters`: A Tensor. Must have the same type as input. Shape: `[filter_depth, filter_height, filter_width, in_channels, out_channels]`. in_channels must match between input and filters. \n\n## 3.1 Adding depth to the image stack\n\nTo demonstrate, let's reorder the input image so that the RGB channels are considered as the depth dimension.\n","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[]}},{"cell_type":"code","source":"col_img_tensor.shape","metadata":{"execution":{"iopub.status.busy":"2023-04-24T13:38:23.479497Z","iopub.execute_input":"2023-04-24T13:38:23.480275Z","iopub.status.idle":"2023-04-24T13:38:23.488045Z","shell.execute_reply.started":"2023-04-24T13:38:23.480220Z","shell.execute_reply":"2023-04-24T13:38:23.486773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We need to use `tf.transpose` here to [reorder dimensions](https://www.tensorflow.org/api_docs/python/tf/transpose). The previous `np.reshape` operations have just been to insert additional dimensions (such as `BATCH_SIZE=1`), whereas using `reshape` to reorder dimensions with shape greater than one would not work as the values would get mixed up.","metadata":{}},{"cell_type":"code","source":"tr_col_img_tensor = tf.transpose(col_img_tensor, (0,3,1,2))\ntr_col_img_tensor.shape","metadata":{"execution":{"iopub.status.busy":"2023-04-24T13:38:23.489355Z","iopub.execute_input":"2023-04-24T13:38:23.489697Z","iopub.status.idle":"2023-04-24T13:38:23.515231Z","shell.execute_reply.started":"2023-04-24T13:38:23.489663Z","shell.execute_reply":"2023-04-24T13:38:23.513652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next add in a channels dimension (just equal to one in this example):","metadata":{}},{"cell_type":"code","source":"image_3d = np.array(tr_col_img_tensor).reshape((1, # BATCH_SIZE\n 3, # Image depth - was number of channels (RGB)\n 171, # Image height\n 128, # Image width\n 1)) # Number of channels\nimage_3d.shape","metadata":{"execution":{"iopub.status.busy":"2023-04-24T13:38:23.516611Z","iopub.execute_input":"2023-04-24T13:38:23.516990Z","iopub.status.idle":"2023-04-24T13:38:23.526164Z","shell.execute_reply.started":"2023-04-24T13:38:23.516954Z","shell.execute_reply":"2023-04-24T13:38:23.524690Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.2 Adding depth to the filter stack\n\nFor demonstration purposes, let's take the horizontal, vertical and diagonal filters from before but set the depth of the filter stack equal to two. ","metadata":{}},{"cell_type":"code","source":"f2c = filters_3channels[:,:,:2,:]\nf2c.shape","metadata":{"execution":{"iopub.status.busy":"2023-04-24T13:38:23.532487Z","iopub.execute_input":"2023-04-24T13:38:23.532904Z","iopub.status.idle":"2023-04-24T13:38:23.542812Z","shell.execute_reply.started":"2023-04-24T13:38:23.532864Z","shell.execute_reply":"2023-04-24T13:38:23.541647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As before, `transpose` this to reorder the dimensions:","metadata":{}},{"cell_type":"code","source":"f2c_tr = tf.transpose(f2c,[2,0,1,3])","metadata":{"execution":{"iopub.status.busy":"2023-04-24T13:38:23.544482Z","iopub.execute_input":"2023-04-24T13:38:23.545005Z","iopub.status.idle":"2023-04-24T13:38:23.552698Z","shell.execute_reply.started":"2023-04-24T13:38:23.544954Z","shell.execute_reply":"2023-04-24T13:38:23.551470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f2c_3d = np.array(f2c_tr).reshape((2, # filter depth\n 8, # filter height\n 8, # filter width,\n 1, # in channels, now 1 as RGB is in the depth dimension\n 4)) # out channels , i.e. number of filters\nf2c_3d.shape","metadata":{"execution":{"iopub.status.busy":"2023-04-24T13:38:23.554295Z","iopub.execute_input":"2023-04-24T13:38:23.555426Z","iopub.status.idle":"2023-04-24T13:38:23.568548Z","shell.execute_reply.started":"2023-04-24T13:38:23.555387Z","shell.execute_reply":"2023-04-24T13:38:23.567386Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.3 `conv3d` output","metadata":{}},{"cell_type":"code","source":"output3 = tf.nn.conv3d(image_3d, \n f2c_3d,\n strides=[1,]*5,\n padding='VALID')\nprint(output3.shape)\n","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:23.569882Z","iopub.execute_input":"2023-04-24T13:38:23.572059Z","iopub.status.idle":"2023-04-24T13:38:23.614782Z","shell.execute_reply.started":"2023-04-24T13:38:23.572020Z","shell.execute_reply":"2023-04-24T13:38:23.613924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The output shape here is `(BATCH_SIZE, DEPTH, HEIGHT, WIDTH, CHANNELS_OUT)` where `CHANNELS_OUT`, the number of channels in the output, is equal to the number of filters used in the convolution.\n\nHere we had no padding (`padding='VALID'`), and so the output depth is 2 since a depth-2 filter can take two positions over a depth-3 image stack. In this example, the depth-2 filter convolves over red/green and green/blue channels. (If we had `'SAME'` padding, the output depth would be 3).","metadata":{}},{"cell_type":"code","source":"fig=plt.figure(figsize=(12,2*12/4))\ngs3d = GridSpec(4,4, \n height_ratios=[1/3,1,1,1])\n\nimport matplotlib.colors as clr\nyellow_cmap = clr.LinearSegmentedColormap.from_list('yellow (R+G)', ['#000000','#FFFF00'], N=256)\ncyan_cmap = clr.LinearSegmentedColormap.from_list('cyan (G+B)', ['#000000','#00FFFF'], N=256)\nfor i in range(4):\n fig.add_subplot(gs3d[i]).imshow(f2c_3d[:,:,:,0,i].sum(0), cmap='gray')\n ax=fig.add_subplot(gs3d[i+4])\n ax.imshow(output3[0,0,:,:,i], cmap=yellow_cmap)\n if i == 0:\n ax.set_ylabel('Filter over\\nR+G channels')\n ax=fig.add_subplot(gs3d[i+8])\n ax.imshow(output3[0,1,:,:,i], cmap=cyan_cmap)\n if i == 0:\n ax.set_ylabel('Filter over\\nG+B channels')\n ax=fig.add_subplot(gs3d[i+12])\n ax.imshow(output3[0,0,:,:,i]-output3[0,1,:,:,i], cmap='gray')\n if i == 0:\n ax.set_ylabel('Difference')","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[],"execution":{"iopub.status.busy":"2023-04-24T13:38:23.615689Z","iopub.execute_input":"2023-04-24T13:38:23.616016Z","iopub.status.idle":"2023-04-24T13:38:25.209963Z","shell.execute_reply.started":"2023-04-24T13:38:23.615984Z","shell.execute_reply":"2023-04-24T13:38:25.208691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In contrast to the 2d convolution over the image stack with RGB layers as channels, applying the depth-2 3-d convolution over the RGB layers as depths allows us to separate out the R+G and G+B channels rather than summing all three after convolution. The difference between these two output channels is shown in the last row.\n\nAnother benefit of 3d convolution is that we can potentially have different filters by depth in the depth-2 filter, whereas the 2d convolution uses the same filter over each channel.\n\n","metadata":{}},{"cell_type":"markdown","source":"#

4. Application to Vesuvius competition

\n\n
\nUnder construction
\n\n#

5. Conclusions

\n\n
\nUnder construction
\n","metadata":{}},{"cell_type":"markdown","source":"#

References

\n\nGeron, A. (2019) Hands-On Machine Learning with Sckit-Learn, Keras & Tensorflow, 2nd ed. O'Reilly.\n\nhttps://towardsdatascience.com/a-comprehensive-introduction-to-different-types-of-convolutions-in-deep-learning-669281e58215\n\nAllosteric via Stack Overflow (2017) Convolve2d just by using Numpy. (https://stackoverflow.com/questions/43086557/convolve2d-just-by-using-numpy)\n\nAragorn via Stack Overflow (2016) How is a convolution calculated on an image with three (RGB) channels? https://stackoverflow.com/q/37095783\n\nKeras API Documentation (https://keras.io)","metadata":{"papermill":{"duration":null,"end_time":null,"exception":null,"start_time":null,"status":"pending"},"tags":[]}}]}