{"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":"# COVID 19 - Image Enhancement\n\nIn this notebook, we will see the existing method for x-rays image enhancement.\nThe idea of this notebook, is do get a better visualization of the x-rays images, which could be useful for doctors to see possible disease or for deep learning techniques to preprocess the input images.\n\nIf you any remark, feel free to leave a comment. Good reading.","metadata":{}},{"cell_type":"code","source":"!conda install '/kaggle/input/pydicom-conda-helper/libjpeg-turbo-2.1.0-h7f98852_0.tar.bz2' -c conda-forge -y\n!conda install '/kaggle/input/pydicom-conda-helper/libgcc-ng-9.3.0-h2828fa1_19.tar.bz2' -c conda-forge -y\n!conda install '/kaggle/input/pydicom-conda-helper/gdcm-2.8.9-py37h500ead1_1.tar.bz2' -c conda-forge -y\n!conda install '/kaggle/input/pydicom-conda-helper/conda-4.10.1-py37h89c1867_0.tar.bz2' -c conda-forge -y\n!conda install '/kaggle/input/pydicom-conda-helper/certifi-2020.12.5-py37h89c1867_1.tar.bz2' -c conda-forge -y\n!conda install '/kaggle/input/pydicom-conda-helper/openssl-1.1.1k-h7f98852_0.tar.bz2' -c conda-forge -y","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-07-28T13:43:52.404272Z","iopub.execute_input":"2021-07-28T13:43:52.404637Z","iopub.status.idle":"2021-07-28T13:44:35.705403Z","shell.execute_reply.started":"2021-07-28T13:43:52.404603Z","shell.execute_reply":"2021-07-28T13:44:35.704258Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd \nimport os\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\n%matplotlib inline\nimport glob\nimport pydicom\n\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nfrom PIL import Image\n\nimport cv2 as cv\n\nimport random \n\nrandom.seed(42)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-28T13:44:35.707642Z","iopub.execute_input":"2021-07-28T13:44:35.708101Z","iopub.status.idle":"2021-07-28T13:44:35.717692Z","shell.execute_reply.started":"2021-07-28T13:44:35.708049Z","shell.execute_reply":"2021-07-28T13:44:35.716464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"TARGET_SIZE = 512\n\ndef dicom2array(path, voi_lut = True, fix_monochrome = True):\n    \"\"\"\n    Transform a dicom file to an array.\n    \n    - path : path of the dicom file\n    - voi_lut : Apply VOI LUT transformation\n    - fix_monochrome : Indicate if we fix the pixel value for specific files.\n    \n    VOI LUT (Value of Interest - Look Up Table) : The idea is to have a larger representation of the data.\n    Since, dicom files have larger pixel display range than usuall pictures. The idea is to keep a larger representation in order ot better see the subtle differences.\n    \n    Fix Monochrome : Some images have MONOCHROME1 interpretation. Which means that higher pixel values corresponding to the dark instead of the white.\n    \"\"\"\n    dicom = pydicom.read_file(path)\n    \n    # Apply the VOI LUT\n    if voi_lut:\n        data = apply_voi_lut(dicom.pixel_array, dicom)\n    else:\n        data = dicom.pixel_array\n        \n    # Fix the representation\n    if fix_monochrome and dicom.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n    \n    data = data - np.min(data)\n    data = data / np.max(data)\n    \n    data = data * 255\n    \n    data = data.astype(np.uint8)\n    \n    return data\n\ndef img_vizualisation(imgs, nb_samples = 5):\n    \n    fig, axes = plt.subplots(nrows=nb_samples // 5, ncols=min(5, nb_samples), figsize=(min(5, nb_samples) * 4, 4 * (nb_samples // 5)))\n    i = 0\n    for img in imgs:\n        axes[i // 5, i % 5].imshow(np.array(img), cmap=plt.cm.gray, aspect='auto')\n        axes[i // 5, i % 5].axis('off')\n        i += 1\n    fig.show()    \n","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:44:35.719821Z","iopub.execute_input":"2021-07-28T13:44:35.720277Z","iopub.status.idle":"2021-07-28T13:44:35.731154Z","shell.execute_reply.started":"2021-07-28T13:44:35.720233Z","shell.execute_reply":"2021-07-28T13:44:35.730139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Sample visualization\n\nTo see the different techniques, we will have a sample set that will serve as a reference.","metadata":{}},{"cell_type":"code","source":"TRAIN_PATH = \"../input/siim-covid19-detection/train/\"\n\npaths = glob.glob(TRAIN_PATH + \"*/*/*.dcm\")","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:44:35.732410Z","iopub.execute_input":"2021-07-28T13:44:35.732669Z","iopub.status.idle":"2021-07-28T13:44:40.509053Z","shell.execute_reply.started":"2021-07-28T13:44:35.732643Z","shell.execute_reply":"2021-07-28T13:44:40.507940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sampled_path = random.sample(paths, 10)\n\nsamples = []\nfor path in sampled_path:\n    img = dicom2array(path)\n    samples.append(img)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:44:40.510405Z","iopub.execute_input":"2021-07-28T13:44:40.510692Z","iopub.status.idle":"2021-07-28T13:44:41.329546Z","shell.execute_reply.started":"2021-07-28T13:44:40.510663Z","shell.execute_reply":"2021-07-28T13:44:41.328322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img_vizualisation(samples, 10)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:44:41.330930Z","iopub.execute_input":"2021-07-28T13:44:41.331240Z","iopub.status.idle":"2021-07-28T13:44:48.871268Z","shell.execute_reply.started":"2021-07-28T13:44:41.331210Z","shell.execute_reply":"2021-07-28T13:44:48.870172Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Histogram equalization\n\nHistogram equalization is a method that **increases the global contrast of images**, especially when the image is represented by a narrow range of intensity values. With this method, we can adjust the pixel intensity in order to obtain a better distribution on the histogram using the full range of intensities evenly. This allows for areas of lower local constrast to gain higher contrast. \n\nFor radiography, this could be usefull, to have a better view of bone structure in x-ray images. A disadvantge of this method is that it is indiscriminate. It may increase the contrast of background noise, while decreasing the usable signal.\n\n> https://en.wikipedia.org/wiki/Histogram_equalization","metadata":{}},{"cell_type":"markdown","source":"### Visualization of histogram equalization transformation\n\nTo have a better idea of the result of histogram equalization, we could take the top right picture in our reference set. With this picture, it's really difficult to see the content of the image. If we have a look at the histogram of value of this picture (see below), we see that the pixel values are squeezes through 0 and 75. With histogram equalization, we can resample the values, in order to use the full range of pixels intensities.","metadata":{}},{"cell_type":"code","source":"img_example = samples[4]\nimg_example_hist = cv.equalizeHist(img_example)\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))\n\nax1.imshow(img_example, cmap=plt.cm.gray)\nax1.axis('off')\nax1.set_title(\"Original image\")\n\nax2.imshow(img_example_hist, cmap=plt.cm.gray)\nax2.axis('off')\nax2.set_title(\"Histogram Equalization applied on the original image\")\n\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-28T13:44:48.872698Z","iopub.execute_input":"2021-07-28T13:44:48.873035Z","iopub.status.idle":"2021-07-28T13:44:50.008909Z","shell.execute_reply.started":"2021-07-28T13:44:48.873003Z","shell.execute_reply":"2021-07-28T13:44:50.007584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 4))\n\nax1.hist(img_example.flatten(), 256, [0, 256])\nax1.set_title(\"Original image\")\n\nax2.hist(img_example_hist.flatten(), 256, [0, 256])\nax2.set_title(\"Histogram Equalization applied on the original image\")\n\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-28T13:44:50.012443Z","iopub.execute_input":"2021-07-28T13:44:50.012952Z","iopub.status.idle":"2021-07-28T13:44:51.151905Z","shell.execute_reply.started":"2021-07-28T13:44:50.012859Z","shell.execute_reply":"2021-07-28T13:44:51.150845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Visualization on the sample data","metadata":{}},{"cell_type":"code","source":"import scipy.misc\n\nequalized_samples = []\nfor sample in samples:\n    img = cv.equalizeHist(sample)\n    equalized_samples.append(img)\n    \nimg_vizualisation(equalized_samples, 10)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:44:51.153817Z","iopub.execute_input":"2021-07-28T13:44:51.154175Z","iopub.status.idle":"2021-07-28T13:44:58.993555Z","shell.execute_reply.started":"2021-07-28T13:44:51.154141Z","shell.execute_reply":"2021-07-28T13:44:58.991863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CLAHE (Contrast Limited Adaptive Histogram Equalization) \n\nThe first histogram equalization we just saw, considers the global contrast of the image. In many cases, it is not a good idea. \nWe could have part of the image lost due to over-brightness.\n\nSo, to solve this problem, **adaptive histogram equalization** is used. To apply this method, image is divided into small blocks called \"tiles\". Then each of these blocks are histogram equalized as usual. So in a small area, histogram would confine to a small region (unless there is noise). If noise is there, it will be amplified. To avoid this, **contrast limiting** is applied. If any histogram bin is above the specified contrast limit, those pixels are clipped and distributed uniformly to other bins before applying histogram equalization. After equalization, to remove artifacts in tile borders, bilinear interpolation is applied.","metadata":{}},{"cell_type":"markdown","source":"### Visualization of CLAHE","metadata":{}},{"cell_type":"code","source":"clahe = cv.createCLAHE(clipLimit=40.0, tileGridSize=(8,8))\n\nimg_example = samples[9]\nimg_example_clahe = clahe.apply(img_example)\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))\n\nax1.imshow(img_example, cmap=plt.cm.gray)\nax1.axis('off')\nax1.set_title(\"Original image\")\n\nax2.imshow(img_example_clahe, cmap=plt.cm.gray)\nax2.axis('off')\nax2.set_title(\"CLAHE applied on the original image\")\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:44:58.996351Z","iopub.execute_input":"2021-07-28T13:44:58.996783Z","iopub.status.idle":"2021-07-28T13:45:01.489209Z","shell.execute_reply.started":"2021-07-28T13:44:58.996730Z","shell.execute_reply":"2021-07-28T13:45:01.488252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 4))\n\nax1.hist(img_example.flatten(), 256, [0, 256])\nax1.set_title(\"Original image\")\n\nax2.hist(img_example_clahe.flatten(), 256, [0, 256])\nax2.set_title(\"CLAHE applied on the original image\")\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:45:01.490655Z","iopub.execute_input":"2021-07-28T13:45:01.491247Z","iopub.status.idle":"2021-07-28T13:45:02.745384Z","shell.execute_reply.started":"2021-07-28T13:45:01.491203Z","shell.execute_reply":"2021-07-28T13:45:02.744505Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Visualization on the sample data","metadata":{}},{"cell_type":"code","source":"clahe_samples = []\n\nclahe = cv.createCLAHE(clipLimit=40.0, tileGridSize=(8,8))\n\nfor sample in samples:\n    img = clahe.apply(sample)\n    clahe_samples.append(img)\n    \nimg_vizualisation(clahe_samples, 10)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:45:02.746611Z","iopub.execute_input":"2021-07-28T13:45:02.746912Z","iopub.status.idle":"2021-07-28T13:45:10.344356Z","shell.execute_reply.started":"2021-07-28T13:45:02.746866Z","shell.execute_reply":"2021-07-28T13:45:10.341082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CLAHE with different parameters","metadata":{}},{"cell_type":"code","source":"clahe_samples_2 = []\n\nclahe = cv.createCLAHE(clipLimit=20.0, tileGridSize=(15, 15))\n\nfor sample in samples:\n    img = clahe.apply(sample)\n    clahe_samples_2.append(img)\n    \nimg_vizualisation(clahe_samples_2, 10)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:45:10.345502Z","iopub.execute_input":"2021-07-28T13:45:10.345931Z","iopub.status.idle":"2021-07-28T13:45:18.116408Z","shell.execute_reply.started":"2021-07-28T13:45:10.345898Z","shell.execute_reply":"2021-07-28T13:45:18.115683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Contrast Enhancement of Medical X-Ray ImageUsing Morphological Operators with OptimalStructuring Element\n\nPaper : Contrast Enhancement of Medical X-Ray ImageUsing Morphological Operators with OptimalStructuring Element\n> https://arxiv.org/pdf/1905.08545.pdf\n\nThe idea of this paper is to apply a combination of Top-hat and Bottom-hat Transform. \n\n- Top Hat : It is the difference between input image and Opening of the image. \n- Opening : Erosion followed by dilatation\n\n- Black Hat : It is the difference between the closing of the input image and input image. \n- Closing : Dilation followed by Erosion.\n\nI have made two implementations, one with skimage library using a disk as the footprint used for top and black hat. And another implementation with opencv using an ellipse. The opencv method is faster than the skimage. \n\nRegarding the result, I am not really satisfied by this method (I maybe miss something in the paper). I let you see the result.","metadata":{}},{"cell_type":"code","source":"import skimage.morphology\nfrom skimage.morphology import disk\n\nimg_example = samples[9]\n\nfootprint = disk(15)\n\ntop = skimage.morphology.white_tophat(img_example, footprint)\nbottom = skimage.morphology.black_tophat(img_example, footprint)\n\noutput = img_example + top - bottom\n\n\ncompare = np.concatenate((img_example, output), axis=1)\n\nplt.figure(figsize=(20,10))\nplt.imshow(compare, cmap=plt.cm.gray)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:45:18.117676Z","iopub.execute_input":"2021-07-28T13:45:18.118155Z","iopub.status.idle":"2021-07-28T13:46:36.621870Z","shell.execute_reply.started":"2021-07-28T13:45:18.118113Z","shell.execute_reply":"2021-07-28T13:46:36.620905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img_example = samples[9]\n\n# The size of the structured element has impact \nfilterSize = (15, 15)\nkernel = cv.getStructuringElement(cv.MORPH_RECT, filterSize) # MORPH_ELLIPSE\n\n\ntophat_img = cv.morphologyEx(img_example, cv.MORPH_TOPHAT, kernel)\nbothat_img = cv.morphologyEx(img_example, cv.MORPH_BLACKHAT, kernel) # Black --> Bottom\n\noutput = img_example + tophat_img - bothat_img\n\ncompare = np.concatenate((img_example, output), axis=1)\n\nplt.figure(figsize=(20,10))\nplt.imshow(compare, cmap=plt.cm.gray)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:46:36.623404Z","iopub.execute_input":"2021-07-28T13:46:36.623765Z","iopub.status.idle":"2021-07-28T13:46:38.799833Z","shell.execute_reply.started":"2021-07-28T13:46:36.623724Z","shell.execute_reply":"2021-07-28T13:46:38.798953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"hat_samples = []\n\nkernel = cv.getStructuringElement(cv.MORPH_RECT, (15, 15)) # MORPH_ELLIPSE\n\nfor sample in samples:\n\n    tophat = cv.morphologyEx(sample, cv.MORPH_TOPHAT, kernel)\n    bothat = cv.morphologyEx(sample, cv.MORPH_BLACKHAT, kernel)\n    img = sample + tophat - bothat\n\n    hat_samples.append(img)\n    \nimg_vizualisation(hat_samples, 10)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:46:38.801160Z","iopub.execute_input":"2021-07-28T13:46:38.801424Z","iopub.status.idle":"2021-07-28T13:46:46.323621Z","shell.execute_reply.started":"2021-07-28T13:46:38.801397Z","shell.execute_reply":"2021-07-28T13:46:46.322670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Change the approach\n\nI'm really not satisfy with the last approach. I think that maybe by applying histogram equalization first and then applying it the morphological operator could be interesting.","metadata":{}},{"cell_type":"code","source":"img_example = equalized_samples[9]\n\nkernel = cv.getStructuringElement(cv.MORPH_RECT, (15, 15)) # MORPH_ELLIPSE\n\ntophat_img = cv.morphologyEx(img_example, cv.MORPH_TOPHAT, kernel)\nbothat_img = cv.morphologyEx(img_example, cv.MORPH_BLACKHAT, kernel) # Black --> Bottom\n\noutput = img_example + tophat_img - bothat_img\n\ncompare = np.concatenate((img_example, output), axis=1)\n\nplt.figure(figsize=(20,10))\nplt.imshow(compare, cmap=plt.cm.gray)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:46:46.324851Z","iopub.execute_input":"2021-07-28T13:46:46.325154Z","iopub.status.idle":"2021-07-28T13:46:48.523949Z","shell.execute_reply.started":"2021-07-28T13:46:46.325126Z","shell.execute_reply":"2021-07-28T13:46:48.523186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 4))\n\nax1.hist(img_example.flatten(), 256, [0, 256])\nax1.set_title(\"Original image\")\n\nax2.hist(output.flatten(), 256, [0, 256])\nax2.set_title(\"This approach on the original image\")\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:50:34.607282Z","iopub.execute_input":"2021-07-28T13:50:34.607712Z","iopub.status.idle":"2021-07-28T13:50:36.099859Z","shell.execute_reply.started":"2021-07-28T13:50:34.607680Z","shell.execute_reply":"2021-07-28T13:50:36.098726Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"hat_samples = []\n\nkernel = cv.getStructuringElement(cv.MORPH_RECT, (15, 15)) # MORPH_ELLIPSE\n\nfor sample in equalized_samples:\n\n    tophat = cv.morphologyEx(sample, cv.MORPH_TOPHAT, kernel)\n    bothat = cv.morphologyEx(sample, cv.MORPH_BLACKHAT, kernel)\n    img = sample + tophat - bothat\n\n    hat_samples.append(img)\n    \nimg_vizualisation(hat_samples, 10)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:46:49.992921Z","iopub.execute_input":"2021-07-28T13:46:49.993192Z","iopub.status.idle":"2021-07-28T13:46:57.742303Z","shell.execute_reply.started":"2021-07-28T13:46:49.993166Z","shell.execute_reply":"2021-07-28T13:46:57.741318Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Noise reduction\n\nWhen we realize a x-rays images, we could had some noise during the process. Before doing image enhancement, it could be interesting to reduce/remove the noise present in our image. In this section, we are going to see two filter : the median filter and the DCT filter.\n\n## Median filter \n\nThe first method that we are going to use is the median filter. For each pixel value on our image, we will select the median value from the values next to the current pixel. In order to do that, we will pass over our image a box (for example : a box of 5x5). This will allow us to remove abnormal value.","metadata":{}},{"cell_type":"code","source":"img_example = samples[0]\n# Note : ksize need to be an odd number.\nnoised_samples_example = cv.medianBlur(img_example, 5)\n\ncompare = np.concatenate((img_example, noised_samples_example), axis=1)\n\nplt.figure(figsize=(20,10))\nplt.imshow(compare, cmap=plt.cm.gray)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:46:57.743466Z","iopub.execute_input":"2021-07-28T13:46:57.743739Z","iopub.status.idle":"2021-07-28T13:46:59.439231Z","shell.execute_reply.started":"2021-07-28T13:46:57.743712Z","shell.execute_reply":"2021-07-28T13:46:59.436560Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 4))\n\nax1.hist(img_example.flatten(), 256, [0, 256])\nax1.set_title(\"Original image\")\n\nax2.hist(noised_samples_example.flatten(), 256, [0, 256])\nax2.set_title(\"Noise reduction on the original image\")\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:46:59.440989Z","iopub.execute_input":"2021-07-28T13:46:59.441419Z","iopub.status.idle":"2021-07-28T13:47:00.700790Z","shell.execute_reply.started":"2021-07-28T13:46:59.441376Z","shell.execute_reply":"2021-07-28T13:47:00.699855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## DCT-based filter\n\nThe idea of DCT (Discrete Cosinus Transform) based filter is to decompose a signal as a sum of cosines function. The DCT is the same idea of Fourier transformation where we decompose a signal into frequencies. The DCT is well used for data compression. But, here we are going to use it for noise reduction.\n\nFirst, we need to understand that an image can be decomposed as signals from the two dimensions of our image. Signals from the x-axis and signals form the y-axis. Based on that, we can decompose our image by decomposing each signal of our image. Each decomposition tries to get the core of the signal. It tries to get the best continues signals possible.\n\nThus, if for a signal, we have an abnormal value for one given point (noise), by decomposing the signal we try to get the core of the signal. Then, this abnormal value will .\n\nFinally, with the decomposition of our image, we can reconstruct the image by applying the inverse function of DCT. By doing that, we get an image similar from the original one but with no abnormals values.  ","metadata":{}},{"cell_type":"code","source":"# https://stackoverflow.com/questions/7110899/how-do-i-apply-a-dct-to-an-image-in-python\nfrom scipy.fftpack import dct, idct\n\ndef dct2(a):\n    return dct(dct(a.T, norm='ortho').T, norm='ortho')\n\ndef idct2(a):\n    return idct(idct(a.T, norm='ortho').T, norm='ortho')  \n\ndef dtc_transform(img):\n    return idct2(dct2(img))\n\nimg_example = samples[0]\nimg_idct = dtc_transform(img_example)\n\nprint(\"Check if the image are similar :\", np.allclose(img_example, img_idct))\n\ncompare = np.concatenate((img_example, img_idct), axis=1)\n\nplt.figure(figsize=(20,10))\nplt.imshow(compare, cmap=plt.cm.gray)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:47:00.701975Z","iopub.execute_input":"2021-07-28T13:47:00.702238Z","iopub.status.idle":"2021-07-28T13:47:03.648065Z","shell.execute_reply.started":"2021-07-28T13:47:00.702212Z","shell.execute_reply":"2021-07-28T13:47:03.647051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 4))\n\nax1.hist(img_example.flatten(), 256, [0, 256])\nax1.set_title(\"Original image\")\n\nax2.hist(img_idct.flatten(), 256, [0, 256])\nax2.set_title(\"Noise reduction on the original image\")\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:47:03.650913Z","iopub.execute_input":"2021-07-28T13:47:03.651196Z","iopub.status.idle":"2021-07-28T13:47:05.138682Z","shell.execute_reply.started":"2021-07-28T13:47:03.651168Z","shell.execute_reply":"2021-07-28T13:47:05.137753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dct_samples = []\n\nfor sample in samples:\n    img = dtc_transform(sample)\n    dct_samples.append(img)\n    \nimg_vizualisation(dct_samples, 10)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:47:05.139997Z","iopub.execute_input":"2021-07-28T13:47:05.140251Z","iopub.status.idle":"2021-07-28T13:47:25.064569Z","shell.execute_reply.started":"2021-07-28T13:47:05.140226Z","shell.execute_reply":"2021-07-28T13:47:25.063582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Combine the different approaches\n\nAs we are dealing with grayscale image, it could be interesting to add different approaches and merge them in a single image in order to create a color image. The main advantage of this, is as we are usually using pretrained model, almost all of them use 3 channel. So, why not given additional information instead of just duplicate the channel.\n\nIn this proposition, I will keep the original image in the first channel and put in the two remaining channels transformed images from the original one.","metadata":{}},{"cell_type":"code","source":"channel_1 = samples[9]\nchannel_2 = clahe_samples_2[9]\nchannel_3 = hat_samples[9]\n\n# Merge all channel to create our image\noutput = np.dstack((channel_1, channel_2, channel_3))\n\nreference = np.dstack((channel_1, channel_1, channel_1))\n\n# See the difference\ncompare = np.concatenate((reference, output), axis=1)\n\nplt.figure(figsize=(20,10))\nplt.imshow(compare)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:47:25.066183Z","iopub.execute_input":"2021-07-28T13:47:25.066619Z","iopub.status.idle":"2021-07-28T13:47:27.637513Z","shell.execute_reply.started":"2021-07-28T13:47:25.066573Z","shell.execute_reply":"2021-07-28T13:47:27.636323Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"multi_channel_samples = []\n\nfor i in range(len(samples)):\n    channel_1 = samples[i]\n    channel_2 = clahe_samples_2[i]\n    # channel_2 = clahe_samples[i]\n    channel_3 = hat_samples[i]\n    \n    out_img = np.dstack((channel_1, channel_2, channel_3))\n    \n    multi_channel_samples.append(out_img)\n    \nimg_vizualisation(multi_channel_samples, 10)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T13:47:27.639053Z","iopub.execute_input":"2021-07-28T13:47:27.639420Z","iopub.status.idle":"2021-07-28T13:47:35.829500Z","shell.execute_reply.started":"2021-07-28T13:47:27.639383Z","shell.execute_reply":"2021-07-28T13:47:35.826159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Conclusion\n\nIn this notebook, we see some techniques that could be useful for image enhancement. \n\nThanks for your reading and your feedback. Don't hesitate if you have questions. \n\nIf you want to go further, there are other methods that I have not mentioned. Do not hesitate to go and see if you are interested.\n\n- MAHE (Median adaptive histogram equalization)\n- Multiple Morphological Gradient (mMG)\n- Wieners filter","metadata":{}},{"cell_type":"markdown","source":"## References \n\n- https://en.wikipedia.org/wiki/Histogram_equalization\n- https://docs.opencv.org/master/d5/daf/tutorial_py_histogram_equalization.html\n- https://www.ndt.net/article/icem2004/papers/64/64.htm\n- https://medium.com/@florestony5454/median-filtering-with-python-and-opencv-2bce390be0d1\n","metadata":{}}]}