{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "041f659f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "#  Copyright 2021 - 2022 United Kingdom Research and Innovation\n",
    "#  Copyright 2021 - 2022 The University of Manchester\n",
    "#  Copyright 2021 - 2022 The Karlsruhe Institute of Technology\n",
    "#  Copyright 2021 - 2022 Technical University of Denmark \n",
    "#\n",
    "#  Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "#  you may not use this file except in compliance with the License.\n",
    "#  You may obtain a copy of the License at\n",
    "#\n",
    "#      http://www.apache.org/licenses/LICENSE-2.0\n",
    "#\n",
    "#  Unless required by applicable law or agreed to in writing, software\n",
    "#  distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "#  See the License for the specific language governing permissions and\n",
    "#  limitations under the License.\n",
    "#\n",
    "#   Authored by:    Evelina Ametova (KIT)\n",
    "#                   Jakob S. Jørgensen (DTU)\n",
    "#                   Gemma Fardell (UKRI-STFC)\n",
    "#                   Laura Murgatroyd (UKRI-STFC)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a297da2",
   "metadata": {},
   "source": [
    "# CT data preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bd28059",
   "metadata": {},
   "source": [
    "Sometimes we are lucky to have a (simulated) dataset which is ready to reconstruct. However sometimes we get raw data which we need to preprocess first to get sensible reconstruction.\n",
    "CIL provides a number of useful image manipulation tools - processors. \n",
    "In this notebook we will demonstrate some of them."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45e58aa7",
   "metadata": {},
   "source": [
    "## Learning objectives:\n",
    "- Read in and manipulate data\n",
    "- Compensate for centre-of-rotation offset\n",
    "- Slice and bin data\n",
    "- Remove hot/dead pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76fa5ad6",
   "metadata": {},
   "source": [
    "We start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04651e7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cil imports\n",
    "from cil.framework import ImageData, ImageGeometry\n",
    "from cil.framework import AcquisitionGeometry, AcquisitionData\n",
    "\n",
    "from cil.processors import CentreOfRotationCorrector, Slicer, \\\n",
    "    Binner, Masker, MaskGenerator, TransmissionAbsorptionConverter\n",
    "\n",
    "from cil.plugins.astra import FBP\n",
    "\n",
    "from cil.utilities import dataexample\n",
    "from cil.utilities.display import show2D, show_geometry\n",
    "\n",
    "\n",
    "# External imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import logging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1859bcfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set logging level for CIL processors:\n",
    "logging.basicConfig(level=logging.WARNING)\n",
    "cil_log_level = logging.getLogger('cil.processors')\n",
    "cil_log_level.setLevel(logging.INFO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8523d3cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up default colour map for visualisation\n",
    "cmap = \"gray\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97216f7d",
   "metadata": {},
   "source": [
    "## Loading dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b8edb1c",
   "metadata": {},
   "source": [
    "We use the steel-wire dataset from the Diamond Light Source (DLS). The dataset is included in CIL for demonstration purposes and can be loaded as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0a0d110",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up a reader object pointing to the Nexus data set\n",
    "data_raw = dataexample.SYNCHROTRON_PARALLEL_BEAM_DATA.get()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f13241a",
   "metadata": {},
   "source": [
    "We load not only the data array itself, but also its corresponding metadata, i.e. `AcquisitionGeometry`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36510dc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data_raw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75a357a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data_raw.geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f01b306",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualise data\n",
    "show2D(data_raw, slice_list=[('angle',0), ('angle', 30), ('angle',60)], \\\n",
    "        cmap=cmap, num_cols=1, size=(15,15), origin='upper-left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33fcc49c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# and show geometry\n",
    "show_geometry(data_raw.geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09aece51",
   "metadata": {},
   "source": [
    "## Transmission to absorption conversion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "399224d4",
   "metadata": {},
   "source": [
    "From the contrast and background values we can infer that the dataset has already been flat-field corrected. However the background values are less than 1. We simply rescale intensity values by taking the mean of an empty slice and dividing the data array by the mean value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6feba920",
   "metadata": {},
   "outputs": [],
   "source": [
    "background = data_raw.get_slice(vertical=20).mean()\n",
    "data_raw /= background"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0441155d",
   "metadata": {},
   "source": [
    "To convert from transmission contrast to absorption contrast we need to apply a negative logarithm.\n",
    "We have implemented `TransmissionAbsorptionConverter()` for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3063f16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert from transmission to attenuation \n",
    "data_exp = TransmissionAbsorptionConverter()(data_raw)\n",
    "# and plot again\n",
    "show2D(data_exp, slice_list=[('angle',0), ('angle', 30), ('angle',60)], \\\n",
    "        cmap=cmap, num_cols=1, size=(15,15), origin='upper-left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51e0cff7",
   "metadata": {},
   "source": [
    "The image contrast does not look right. The reason for this is a single dead pixel in the top left corner."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ebe56c3",
   "metadata": {},
   "source": [
    "## Removing bad pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c584e042",
   "metadata": {},
   "source": [
    "Dead/stuck/misperforming pixels are a very common problem. Sometimes there are only a few of them and they will be effectively filtered out during reconstruction. However sometimes the flat-field images look like the night sky. Misperforming pixels can significantly impair the reconstructed image quality and are best dealt with as a preprocessing step.\n",
    "\n",
    "CIL provides processors that can be used to correct these pixels. We'll step through them below but more advanced options will be discussed in the [advanced techniques](#advanced) section of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a28ba5c",
   "metadata": {},
   "source": [
    "For our wire dataset there are several ways to remove the bright pixel. The simplest is to crop the top slice. As there is only air in this slice anyway there will not be any loss of information.\n",
    "\n",
    "`Slicer()` is a processor used to slice the data, similar to numpy slicing. To crop the data pass the region of interest parameter `roi`. This is a dictionary where each element defines the behaviour along one dimension. \n",
    "\n",
    "To crop along an axis, pass a tuple containing the start index, the end index and step size. `roi={vertical: (index0, index1)}` will crop the data between `index0` and `index1` along the `vertical` dimension. Each of these values can be set to `None` to include all the data i.e., `{vertical: (index0, None)}` will crop the data only on one side. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f84635ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_crop = Slicer(roi={'vertical': (1, None)})(data_exp)\n",
    "\n",
    "show2D(data_crop, slice_list=[('angle',0), ('angle', 30), ('angle',60)], \\\n",
    "        cmap=cmap, num_cols=1, size=(15,15), origin='upper-left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16c489f9",
   "metadata": {},
   "source": [
    "These projections look much better now!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3122f51a",
   "metadata": {},
   "source": [
    "Alternatively we can use a processor for outlier detection. It is called `MaskGenerator()`. \n",
    "`MaskGenerator()` is a powerful tool to detect outliers, which was inspired by the MATLAB `rmoutliers`\n",
    "function. It supports a number of methods including simple threshold and quantiles along with statistical median, mean, moving median and moving mean methods."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c549914",
   "metadata": {},
   "source": [
    "In this case, a simple threshold is sufficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a54f4fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = MaskGenerator.threshold(max_val=10)(data_exp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3811a34",
   "metadata": {},
   "source": [
    "Now `mask` is a binary image which contains 0 where outliers were detected and 1 for other pixels. We use `Masker()` to mask out the detected outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "265145a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_masked = Masker.interpolate(mask=mask, method='nearest', axis='vertical')(data_exp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b07ad96b",
   "metadata": {},
   "source": [
    "Let's visualise the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71aa2922",
   "metadata": {},
   "outputs": [],
   "source": [
    "show2D(data_masked, slice_list=[('angle',0), ('angle', 30), ('angle',60)], \\\n",
    "        cmap=cmap, num_cols=1, size=(15,15), origin='upper-left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1eda7d93",
   "metadata": {},
   "source": [
    "Note that `data_crop` and `data_masked` will have different shapes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c6b58a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('data_crop shape: {}'.format(data_crop.shape))\n",
    "print('data_masked shape: {}'.format(data_masked.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fdb37f7",
   "metadata": {},
   "source": [
    "## FBP reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1a8c659",
   "metadata": {},
   "source": [
    "The next step is to reconstruct the dataset. In this notebook we will use a simple FBP reconstruction. More advanced methods will be discussed in the next notebooks. \n",
    "\n",
    "CIL supports different back-ends for which data order conventions may differ. Here we use the FBP algorithm from the ASTRA toolbox plugin.\n",
    "\n",
    "In 3D geometry the ASTRA toolbox requires the dataset in the form `['vertical','angle','horizontal']`, which doesn't match the DLS dataset. We can reorder the data in place, for use with the ASTRA plugin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84ca67c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('old dimension labels: {}'.format(data_crop.dimension_labels))\n",
    "data_crop.reorder(order='astra')\n",
    "print('new dimension labels: {}'.format(data_crop.dimension_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f4ffe4c",
   "metadata": {},
   "source": [
    "Now we are ready to run FBP reconstruction. Remember, reconstruction requires `ImageGeometry` and `AcquisitionGeometry`. `data_crop` contains the dataset itself along with all metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4c7e06b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get acquisition geometry\n",
    "ag = data_crop.geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1260c139",
   "metadata": {},
   "source": [
    "Here we use a default `ImageGeometry` calculated from the `AcquisitionGeometry`. `ImageGeometry` can be modified to reconstruct on a coarser/finer grid or perform ROI reconstruction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa2b5d5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create image geometry\n",
    "ig = ag.get_ImageGeometry()\n",
    "\n",
    "fbp_recon = FBP(ig, ag, device='gpu')(data_crop)\n",
    "\n",
    "# visualise reconstruction results\n",
    "show2D(fbp_recon, slice_list=[('vertical',80), ('vertical',100), ('horizontal_x',80)], \\\n",
    "        cmap=cmap, num_cols=1, size=(15,15), origin='upper-left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b057f51",
   "metadata": {},
   "source": [
    "We can see that the reconstructed slices do not look good. If you have ever looked at CT data you will probably recognise that there is an offset in the centre of rotation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d90dbc70",
   "metadata": {},
   "source": [
    "## Centre of Rotation correction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03557cd2",
   "metadata": {},
   "source": [
    "In a well aligned CT system, the axis of rotation is perpendicular to the X-ray beam and the rows of detector pixels. The centre of rotation is the projection of the axis of rotation on to the detector. The reconstruction assumes this is horizontally centred on the detector. An offset introduces blurring and artefacts in the reconstruction."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09e339e7",
   "metadata": {},
   "source": [
    "There are various ways to estimate the centre of rotation offset. For the parallel geometry case we can use cross-correlation between 0 and 180 degrees. CIL provides a processor which implements this method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b73523d",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_centred = CentreOfRotationCorrector.xcorrelation()(data_crop)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34905f17",
   "metadata": {},
   "source": [
    "Note, that `CentreOfRotationCorrector` doesn't modify the dataset but updates the corresponding geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b568842",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('data_crop rotation axis position: {}'.format(data_crop.geometry.config.system.rotation_axis.position))\n",
    "print('data_centred rotation axis position: {}'.format(data_centred.geometry.config.system.rotation_axis.position))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61b696cd",
   "metadata": {},
   "source": [
    "We use the `show_geometry` utility to illustrate `AcquisitionGeometry` before and after the correction. Note, after the correction the rotation axis position and the detector positon do not coincide anymore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cc875f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_geometry(data_crop.geometry)\n",
    "show_geometry(data_centred.geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6e47d99",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get acquisition geometry\n",
    "ag_centred = data_centred.geometry\n",
    "\n",
    "# FBP reconstruction\n",
    "fbp_recon_centred = FBP(ig, ag_centred, device='gpu')(data_centred)\n",
    "\n",
    "# visualize reconstruction results\n",
    "show2D(fbp_recon_centred,  slice_list = [('vertical',80), ('vertical', 100) ], \\\n",
    "        cmap=cmap, num_cols=1, size=(10,10), origin='upper-left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47c51f9b",
   "metadata": {},
   "source": [
    "## Slicing and binning data\n",
    "\n",
    "The data contains a redundant projection at 180 degrees, which can be discarded by keeping only the 90 angles.\n",
    "\n",
    "We could also crop both sides of the image to remove pixels that only see air. We want to keep only horizontal pixels from 20 to 140 out of 160.\n",
    "\n",
    "Both of these can be done in the same operation by using the `Slicer` Processor. This processor modifies both the data and the corersponding geometry and the trimmed data is printed, showing the horizontal dimension now reduced to 120. Note that `Slicer` supports negative indexing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9987bf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_centred = Slicer(roi={'angle':(0,90), 'horizontal':(20,-20,1)})(data_centred)\n",
    "print(data_centred)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c32ac7a2",
   "metadata": {},
   "source": [
    "Quite often we want to test methods with a lower number of projections. `Slicer` can be used to skip projections. Here we will use `Slicer` to generate new datasets with a lower number of projections. To speed up reconstrucion, we will work only with a single slice. Note, dimension labels refer to different dimensions therefore we can conviniently use the `get_slice` method to extract a single slice along the corresponding dimension."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83590378",
   "metadata": {},
   "source": [
    "We create a new 2D dataset from a single slice of the 3D data, and create a corresponding 2D `ImageGeometry`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b67ccbe3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get single slice\n",
    "data_slice = data_centred.get_slice(vertical=100)\n",
    "# and corresponding geometry\n",
    "ig_slice = data_slice.geometry.get_ImageGeometry()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c379f36",
   "metadata": {},
   "source": [
    "As you can see, the new acquisition geometry is already configured for us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3d5a7c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data_slice.geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb04c9d3",
   "metadata": {},
   "source": [
    "Here we will create 6 datasets each with fewer angles, and reconstruct each one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b9f42be",
   "metadata": {},
   "outputs": [],
   "source": [
    "step_list = [1,2,3,4,5,6]\n",
    "titles = []\n",
    "results = []\n",
    "\n",
    "for step in step_list:\n",
    "   \n",
    "    #slice acquisition data\n",
    "    data_sliced = Slicer(roi={'angle': (None, None, step)})(data_slice)\n",
    "\n",
    "    #Perform a fast reconstruction of the slice using FBP\n",
    "    FBP_output = FBP(ig_slice, data_sliced.geometry, device='gpu')(data_sliced)\n",
    "\n",
    "    #save the results\n",
    "    titles.append(\"# projections {}\".format(data_sliced.shape[0]))\n",
    "    results.append(FBP_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eba481b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot the results    \n",
    "show2D(results, titles, fix_range=True, cmap=cmap, origin='upper-left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04060bab",
   "metadata": {},
   "source": [
    "We might also want to bin the data. Instead of picking out values we can may want to average data together. CIL provides a `Binner` processor to do this. \n",
    "\n",
    "It is used in a similar way to `Slicer` but instead of skipping elements, it calculates their average. For instance, `Binner(roi={horizontal: (None, None,2)})` will calculate the average of every 2 elements along horizontal axis.\n",
    "\n",
    "This is demonstrated below, with the acquistion data with 2x binning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8112d8af",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_binned = Binner(roi={'horizontal': (None, None, 2)})(data_slice)\n",
    "show2D(data_slice, \"original sinogram\", fix_range=True, cmap=cmap)\n",
    "show2D(data_binned, \"binned sinogram\", fix_range=True, cmap=cmap)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4265e069",
   "metadata": {},
   "source": [
    "Now we will create 4 new datasets with increasing binning to see the effect binning the acquistion data has on the reconstructed volume.\n",
    "\n",
    "We bin the `AcquisitionData` but maintain the original `ImageData` with the original 120x120 resolution. This means the reconstructed slice contains the same number of voxels as before, but the spatial resolution of the reconstruction will degrade."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bbb57fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_list = [1,2,3,4]\n",
    "sinograms = []\n",
    "titles = []\n",
    "results = []\n",
    "\n",
    "for bin in bin_list:\n",
    "   \n",
    "    #slice acquisition data\n",
    "    data_binned = Binner(roi={'horizontal': (None, None, bin)})(data_slice)\n",
    "    sinograms.append(data_binned)\n",
    "    \n",
    "    #Perform a fast reconstruction of the slice using FBP\n",
    "    FBP_output = FBP(ig_slice, data_binned.geometry, device='gpu')(data_binned)\n",
    "\n",
    "    #save the results\n",
    "    titles.append(\"# pixels {}\".format(data_binned.shape[1]))\n",
    "    results.append(FBP_output)\n",
    "\n",
    "#plot the results    \n",
    "show2D(results, titles, fix_range=True, cmap=cmap, origin='upper-left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd86116b",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "In this notebook you learned how to use some basic processors provided by CIL. These processors support basic image manipulations and allow quick design of benchmark studies without manual modification of `AcquisitionGeometry`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "047b5dec",
   "metadata": {},
   "source": [
    "<a id='advanced'></a>\n",
    "## Advanced: working with bad pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "772d70cf",
   "metadata": {},
   "source": [
    "Often we have calibrated detector data with the bad pixels corrected before we get the acquistion data. When this isn't the case, or when calibration is not sufficient to remove bad pixels, we can use CIL's `MaskGenerator` and `Masker` to identify outliers and correct them."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffe231c2",
   "metadata": {},
   "source": [
    "We are going to add some bad pixels and columns of bad pixels to our dataset.\n",
    "\n",
    "The function below will return a new dataset, which is a copy of the input with the addition of `number_of_columns` corrupted columns and `number_of_hot_pix` hot pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae420cef",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_bad_pixels(data, number_of_columns, number_of_hot_pix, seed):\n",
    "\n",
    "    data_corrupted = data.copy()\n",
    "\n",
    "    # get intensity range\n",
    "    low = np.amin(data.as_array())\n",
    "    high = np.amax(data.as_array())\n",
    "\n",
    "    # we seed random number generator for repeatability\n",
    "    rng = np.random.RandomState(seed=seed) \n",
    "    # indices of bad columns\n",
    "    columns = rng.randint(0, data.shape[1], size=number_of_columns)\n",
    "    # indices of hot pixels\n",
    "    pix_row = rng.randint(0, data.shape[0], size=number_of_hot_pix)\n",
    "    pix_col = rng.randint(0, data.shape[1], size=number_of_hot_pix)\n",
    "    # values in hot pixels\n",
    "    pixel_values = rng.uniform(low=low, high=high, size=number_of_hot_pix)\n",
    "\n",
    "    for i in range(number_of_columns):\n",
    "        col_pattern = rng.uniform(low=low, high=high, size=data.shape[0])\n",
    "        data_corrupted.as_array()[:, columns[i]] = data.as_array()[:, columns[i]]+col_pattern\n",
    "\n",
    "    for i in range(number_of_hot_pix):\n",
    "        data_corrupted.as_array()[pix_row[i], pix_col[i]] = pixel_values[i]\n",
    "\n",
    "    return data_corrupted    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b78e2717",
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of 'bad' columns\n",
    "number_of_columns = 5\n",
    "# number of randomly located hot pixels\n",
    "number_of_hot_pix = 200\n",
    "# we seed random number generator for repeatability\n",
    "seed = 8392\n",
    "\n",
    "data_corrupted = add_bad_pixels(data_slice, number_of_columns, number_of_hot_pix, seed)\n",
    "\n",
    "show2D([data_slice, data_corrupted], \\\n",
    "        ['clean data', 'corrupted data'], \\\n",
    "        cmap=cmap, num_cols=2, size=(15,10))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "108c5248",
   "metadata": {},
   "source": [
    "Below we show the FBP reconstruction of `corrupted_data` and can see severe artifacts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0112041a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fbp = FBP(ig_slice, data_slice.geometry, device='gpu')\n",
    "fbp.set_input(data_slice)\n",
    "fbp_recon_clean = fbp.get_output()  \n",
    "\n",
    "fbp.set_input(data_corrupted)\n",
    "fbp_recon_corrupted = fbp.get_output()  \n",
    "\n",
    "show2D([fbp_recon_clean, fbp_recon_corrupted], \\\n",
    "        ['clean data', 'corrupted data'], \\\n",
    "        cmap=cmap, num_cols=2, size=(15,10), origin='upper-left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "645e8016",
   "metadata": {},
   "source": [
    "In this case, simple thresholding will not detect all bad pixels. We use `MaskGenerator` with the `median` method, and a moving window of 7 pixels, to detect outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f263bb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = MaskGenerator.median(threshold_factor=3, window=7)(data_corrupted)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be01610d",
   "metadata": {},
   "source": [
    "`MaskGenerator` returns a binary image which contains 0 where outliers were detected and 1 for other pixels. \n",
    "We can look at the generated mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37c2083c",
   "metadata": {},
   "outputs": [],
   "source": [
    "show2D(mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09aa99af",
   "metadata": {},
   "source": [
    "Now we have the mask we want to apply it by using `Masker`. We use linear interpolation in this case and pass it the mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3f20bf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_masked = Masker.interpolate(mask=mask, method='linear', axis='horizontal')(data_corrupted)\n",
    "\n",
    "# visualise corrected sinogram\n",
    "show2D([data_slice, data_masked, (data_slice-data_masked).abs()], \\\n",
    "        ['clean data', 'masked data', 'difference'], \\\n",
    "        cmap=cmap, num_cols=1, size=(15,15), origin='upper-left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6031e3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# reconstruct corrected data\n",
    "fbp.set_input(data_masked)\n",
    "fbp_recon_masked = fbp.get_output()  \n",
    "\n",
    "# visualise reconstructions\n",
    "show2D([fbp_recon_clean, fbp_recon_masked], \\\n",
    "        ['clean data', 'masked data'], \\\n",
    "        cmap=cmap, num_cols=2, size=(15,10), origin='upper-left')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.11",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  },
  "vscode": {
   "interpreter": {
    "hash": "cf07678abc5cc77bc6e1a7d19b1e87ab0c29b83e7ee41c2bc72506d16d80ed44"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}