{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport glob\nfrom PIL import Image\nimport imageio\nimport scipy.ndimage as ndi\nimport matplotlib.pyplot as plt","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-04-12T00:50:08.552675Z","iopub.execute_input":"2023-04-12T00:50:08.553089Z","iopub.status.idle":"2023-04-12T00:50:08.559822Z","shell.execute_reply.started":"2023-04-12T00:50:08.553050Z","shell.execute_reply":"2023-04-12T00:50:08.558651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df_path = \"/kaggle/input/vesuvius-challenge-ink-detection/train\"\ntest_df_path  = \"/kaggle/input/vesuvius-challenge-ink-detection/test\"","metadata":{"execution":{"iopub.status.busy":"2023-04-12T00:50:08.564760Z","iopub.execute_input":"2023-04-12T00:50:08.565456Z","iopub.status.idle":"2023-04-12T00:50:08.572045Z","shell.execute_reply.started":"2023-04-12T00:50:08.565411Z","shell.execute_reply":"2023-04-12T00:50:08.570352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img_path = train_df_path+\"/1/surface_volume\"\nf_p = glob.glob(img_path+\"/*.tif\")","metadata":{"execution":{"iopub.status.busy":"2023-04-12T00:50:08.574195Z","iopub.execute_input":"2023-04-12T00:50:08.574569Z","iopub.status.idle":"2023-04-12T00:50:08.584349Z","shell.execute_reply.started":"2023-04-12T00:50:08.574536Z","shell.execute_reply":"2023-04-12T00:50:08.583549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load an image, I shall use this image for further demostrations\nim = np.array(Image.open(f_p[0]))\nim = im/im.max() # normalize the image\nim = im.astype(np.float32)","metadata":{"execution":{"iopub.status.busy":"2023-04-12T00:50:08.611866Z","iopub.execute_input":"2023-04-12T00:50:08.612186Z","iopub.status.idle":"2023-04-12T00:50:08.873741Z","shell.execute_reply.started":"2023-04-12T00:50:08.612151Z","shell.execute_reply":"2023-04-12T00:50:08.872199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#function to filter out nonzero pixel values\ndef nonzero_image(image):\n    pixels = image.ravel() #flatten\n    nonzero_pixels = pixels[np.nonzero(pixels)] # taking only non zero pixels coz thats what matter\n    return nonzero_pixels","metadata":{"execution":{"iopub.status.busy":"2023-04-12T01:20:40.478567Z","iopub.execute_input":"2023-04-12T01:20:40.478967Z","iopub.status.idle":"2023-04-12T01:20:40.484894Z","shell.execute_reply.started":"2023-04-12T01:20:40.478922Z","shell.execute_reply":"2023-04-12T01:20:40.483689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1. Histogram plot of pixel values (normalized)","metadata":{}},{"cell_type":"code","source":"plt.hist(nonzero_image(im),100000)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-12T01:21:00.961709Z","iopub.execute_input":"2023-04-12T01:21:00.962036Z","iopub.status.idle":"2023-04-12T01:23:27.071864Z","shell.execute_reply.started":"2023-04-12T01:21:00.962006Z","shell.execute_reply":"2023-04-12T01:23:27.070321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" The histogram shows which pixel values are the most frequent, the plot follows a kind of skewed normal (gaussian) distribution, with the highest density of pixel values being between the interval [0.2, 0.4]. This could be helpful to keep `overfitting` in check during training of the model.","metadata":{}},{"cell_type":"code","source":"#flatten the image, will require it later.\nflat_im = im.flatten()","metadata":{"execution":{"iopub.status.busy":"2023-04-12T00:50:08.904479Z","iopub.execute_input":"2023-04-12T00:50:08.904793Z","iopub.status.idle":"2023-04-12T00:50:08.945331Z","shell.execute_reply.started":"2023-04-12T00:50:08.904768Z","shell.execute_reply":"2023-04-12T00:50:08.943799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2. CDF Distribution of pixel values (normalized)","metadata":{}},{"cell_type":"code","source":"hist = ndi.histogram(im, min=0,max=1,bins=256) # no. of bins won't affect the shape/nature of CDF distribution\ncdf = hist.cumsum() / hist.sum()\nplt.plot(cdf)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-12T00:50:08.946523Z","iopub.execute_input":"2023-04-12T00:50:08.947980Z","iopub.status.idle":"2023-04-12T00:50:11.399036Z","shell.execute_reply.started":"2023-04-12T00:50:08.947915Z","shell.execute_reply":"2023-04-12T00:50:11.398056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The CDF follows the histogram. About 50% of the pixels reside within roughly 1/5th of the maxmimum pixel value, although one reason for this could be due to a fair amount of `plain black` pixels which have pixel value of 0.\n\nThe CDF contains of 256 intervals of pixels in incerasing order, with each interval of range [0,1].","metadata":{}},{"cell_type":"markdown","source":"##  3. Convolution using standard 3x3 matrix","metadata":{}},{"cell_type":"code","source":"im = imageio.v2.imread(f_p[0])\nweights = [[1, 1, 1],[ 0, 0, 0],[-1, -1, -1]]\nedges = ndi.convolve(im, weights)\nfig, ax = plt.subplots(figsize=(24,31))\nax.imshow(edges, cmap='seismic')","metadata":{"execution":{"iopub.status.busy":"2023-04-12T00:50:11.400138Z","iopub.execute_input":"2023-04-12T00:50:11.400583Z","iopub.status.idle":"2023-04-12T00:50:16.941713Z","shell.execute_reply.started":"2023-04-12T00:50:11.400547Z","shell.execute_reply":"2023-04-12T00:50:16.940358Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"When a standard 3x3 matrix is convolved over the input image, the output image contains visible gradients which can most probably be linked with extraction of primitive edge features during training of the model. This demonstrates that it should be possible for CNN models to detect edges of ink-labels given enough nodes and training time.","metadata":{}},{"cell_type":"code","source":"im = imageio.v2.imread(f_p[0])\nim = im/im.max()\nweights = [[1, 1, 1],[ 0, 0, 0],[-1, -1, -1]]\nedges = ndi.convolve(im, weights)\nfig, ax = plt.subplots(figsize=(24,31))\nax.imshow(edges, cmap='seismic')","metadata":{"execution":{"iopub.status.busy":"2023-04-12T01:18:16.709241Z","iopub.execute_input":"2023-04-12T01:18:16.709662Z","iopub.status.idle":"2023-04-12T01:18:22.324260Z","shell.execute_reply.started":"2023-04-12T01:18:16.709624Z","shell.execute_reply":"2023-04-12T01:18:22.323234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"When the same convolution matrix is convolved with the normalized version (pixel values between 0 an 1, inclusive) of the same input image, features although extracted are now harder to notice, but it shouldn't effect the training of the CNN model (if used) eitherway.","metadata":{}},{"cell_type":"markdown","source":"## 4. Sobel Filter","metadata":{}},{"cell_type":"code","source":"edges_0=ndi.sobel(im, axis=0)\nfig, ax = plt.subplots(figsize=(24,31))\nax.imshow(edges_0, cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2023-04-12T00:50:16.943193Z","iopub.execute_input":"2023-04-12T00:50:16.943587Z","iopub.status.idle":"2023-04-12T00:50:22.199478Z","shell.execute_reply.started":"2023-04-12T00:50:16.943548Z","shell.execute_reply":"2023-04-12T00:50:22.197436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"When Sobel Filter (horizontal) is passed over the input image, visible features are extracted from the image, going forward some of these features could be critical when performing edge detection over the pixels involving ink-labels during training.","metadata":{}},{"cell_type":"code","source":"edges_1=ndi.sobel(im, axis=1)\nfig, ax = plt.subplots(figsize=(24,31))\nax.imshow(edges_1, cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2023-04-12T00:50:22.201025Z","iopub.execute_input":"2023-04-12T00:50:22.201975Z","iopub.status.idle":"2023-04-12T00:50:27.442170Z","shell.execute_reply.started":"2023-04-12T00:50:22.201941Z","shell.execute_reply":"2023-04-12T00:50:27.441135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here the input image is passed with the vertical version of the sobel filter to give the output image. Features visible here are a result of that. ","metadata":{}},{"cell_type":"code","source":"im = im/im.max()\nedges_0=ndi.sobel(im, axis=0)\nedges_1=ndi.sobel(im, axis=1)\nedges=np.sqrt(np.square(edges_0) + np.square(edges_1))\nfig, ax = plt.subplots(figsize=(24,31))\nax.imshow(edges, cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2023-04-12T02:46:11.676829Z","iopub.execute_input":"2023-04-12T02:46:11.677439Z","iopub.status.idle":"2023-04-12T02:46:18.370603Z","shell.execute_reply.started":"2023-04-12T02:46:11.677381Z","shell.execute_reply":"2023-04-12T02:46:18.368894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is the output imageswhen both horizontal and vertical sobel filters are passed to the input image (normalized this time), features are harder to notice but for the model to pick features during training shouldn't be a problem. ","metadata":{}},{"cell_type":"markdown","source":"## 5. Feature extraction through error measurement","metadata":{}},{"cell_type":"code","source":"im_1=imageio.v2.imread(f_p[0])\nim_2=imageio.v2.imread(f_p[1])\nerr = im_1 - im_2\nfig, ax = plt.subplots(figsize=(24,31))\nax.imshow(err)","metadata":{"execution":{"iopub.status.busy":"2023-04-12T03:44:10.617226Z","iopub.execute_input":"2023-04-12T03:44:10.617634Z","iopub.status.idle":"2023-04-12T03:44:15.309115Z","shell.execute_reply.started":"2023-04-12T03:44:10.617598Z","shell.execute_reply":"2023-04-12T03:44:15.306567Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"im_1=imageio.v2.imread(f_p[0])\nim_2=imageio.v2.imread(f_p[1])\nim_1 = im_1/im_1.max()\nim_2 = im_2/im_2.max()\nerr = im_1 - im_2\nfig, ax = plt.subplots(figsize=(24,31))\nax.imshow(err)","metadata":{"execution":{"iopub.status.busy":"2023-04-12T03:42:52.304549Z","iopub.execute_input":"2023-04-12T03:42:52.304935Z","iopub.status.idle":"2023-04-12T03:42:58.157136Z","shell.execute_reply.started":"2023-04-12T03:42:52.304900Z","shell.execute_reply":"2023-04-12T03:42:58.155663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here the output images(non-normalized followed by normalized) are the `error` or raw pixel value differences between two images. In aspect of a 3D image, this shows that how features might be distributed across the image slices of a 3D image fragment.","metadata":{}},{"cell_type":"code","source":"im_1=imageio.v2.imread(f_p[0])\nim_2=imageio.v2.imread(f_p[1])\nerr = im_1 - im_2\nabs_err = np.abs(err)\nfig, ax = plt.subplots(figsize=(24,31))\nax.imshow(abs_err)","metadata":{"execution":{"iopub.status.busy":"2023-04-12T03:47:17.567745Z","iopub.execute_input":"2023-04-12T03:47:17.568117Z","iopub.status.idle":"2023-04-12T03:47:22.271838Z","shell.execute_reply.started":"2023-04-12T03:47:17.568086Z","shell.execute_reply":"2023-04-12T03:47:22.269212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"im_1=imageio.v2.imread(f_p[0])\nim_2=imageio.v2.imread(f_p[1])\nim_1 = im_1/im_1.max()\nim_2 = im_2/im_2.max()\nerr = im_1 - im_2\nabs_err = np.abs(err)\nfig, ax = plt.subplots(figsize=(24,31))\nax.imshow(abs_err)","metadata":{"execution":{"iopub.status.busy":"2023-04-12T03:48:53.732679Z","iopub.execute_input":"2023-04-12T03:48:53.733015Z","iopub.status.idle":"2023-04-12T03:48:59.390356Z","shell.execute_reply.started":"2023-04-12T03:48:53.732985Z","shell.execute_reply":"2023-04-12T03:48:59.389041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is the output image or `error` when absolute difference is taken into account. Different kind of features could be extracted other than that from the raw pixel difference could be extracted that might be useful during training. ","metadata":{}},{"cell_type":"code","source":"# normalized case\nmean_err = np.mean(abs_err)\nmean_err","metadata":{"execution":{"iopub.status.busy":"2023-04-12T03:43:39.827519Z","iopub.execute_input":"2023-04-12T03:43:39.827961Z","iopub.status.idle":"2023-04-12T03:43:39.858020Z","shell.execute_reply.started":"2023-04-12T03:43:39.827923Z","shell.execute_reply":"2023-04-12T03:43:39.857073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The mean error here indicates how different the pixel distribution on an average between two image slices are.","metadata":{}},{"cell_type":"markdown","source":"## 6. Distribution of RLE pixels.","metadata":{}},{"cell_type":"code","source":"#rle function\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","metadata":{"execution":{"iopub.status.busy":"2023-04-12T00:50:37.883552Z","iopub.execute_input":"2023-04-12T00:50:37.883829Z","iopub.status.idle":"2023-04-12T00:50:37.893241Z","shell.execute_reply.started":"2023-04-12T00:50:37.883802Z","shell.execute_reply":"2023-04-12T00:50:37.892441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x = rle(im)[0]\nx = x/x.max()\nplt.hist(x,100,density = True)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-12T00:50:38.158425Z","iopub.execute_input":"2023-04-12T00:50:38.158722Z","iopub.status.idle":"2023-04-12T00:50:38.431202Z","shell.execute_reply.started":"2023-04-12T00:50:38.158693Z","shell.execute_reply":"2023-04-12T00:50:38.429670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The Histogram here demonstrates how frequency might be related to RLE encoded pixels. The highest potential label pixel frequency seems somwhere around the 35% to 40% against the highest pixel index that was encoded via RLE. This could be useful to see which range of pixel values could be potentially close to actual ink-label pixels.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}