{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Derivative Operator\n",
    "===\n",
    "\n",
    "Here, we show you what the derivative operator and how it works on the example images"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filtering in an image means to perform some kind of processing on a pixel value $I(x, y)$\n",
    "using its neighboring pixel values as follows:\n",
    "$$\n",
    "I'(x, y) = \\sum_{i=-1}^{1}\\sum_{j=-1}^{1}K(i, j)I(x + i, y + j),\n",
    "$$\n",
    "where, $I'(x, y)$: performed pixel value, $K$: kernel matrix.\n",
    "\n",
    "The first derivative operator in each direction ($x$, $y$) is represented as following kernels:\n",
    "$$\n",
    "K_x =\n",
    "\\begin{bmatrix}\n",
    "    0 & 0 & 0 \\\\\n",
    "    -1 & 1 & 0 \\\\\n",
    "    0 & 0 & 0\n",
    "\\end{bmatrix},\n",
    "\\quad\n",
    "K_y =\n",
    "\\begin{bmatrix}\n",
    "    0 & -1 & 0 \\\\\n",
    "    0 & 1 & 0 \\\\\n",
    "    0 & 0 & 0\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Laplacian Operator is a also derivative operator which is used to find edges in an image and represented\n",
    "as following kernels:\n",
    "\n",
    "$$\n",
    "K_4 =\n",
    "\\begin{bmatrix}\n",
    "0 & 1 & 0\\\\\n",
    "1 & -4 & 1\\\\\n",
    "0 & 1 & 0\n",
    "\\end{bmatrix},\n",
    "\\quad\n",
    "K_8 =\n",
    "\\begin{bmatrix}\n",
    "1 & 1 & 1\\\\\n",
    "1 & -8 & 1\\\\\n",
    "1 & 1 & 1\n",
    "\\end{bmatrix},\n",
    "$$\n",
    "\n",
    "where $K_4$: a kernel that considers the contribution of 4 nearest neighbors\n",
    "(top, bottom, left, right) to the pixel of interest, $K_8$: a kernel that considers 8 nearest\n",
    "neighbors (top, bottom, left, right, diagonal) to the pixel of interest.\n",
    "\n",
    "To perform these operator to a image converted 1-D vector array, we generate derivative matrices.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.cbook import get_sample_data\n",
    "from matplotlib.colors import CenteredNorm\n",
    "from mpl_toolkits.axes_grid1 import ImageGrid\n",
    "from PIL import Image\n",
    "\n",
    "from cherab.phix.tools import compute_dmat"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize derivative matrix\n",
    "---\n",
    "\n",
    "Try to create the simple derivative matrix (10, 10).\n",
    "Firstly, create the mapping array denoting 2-D image shape and the element of which denotes a index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mapping array\n",
    "mapping_array = np.arange(0, 50, dtype=np.int32).reshape(10, 5)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot derivative matrix as a sparse matrix and compare $K_x$ and $K_4$ kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, dpi=150, tight_layout=True)\n",
    "for ax, kernel_type, kernel_name in zip(axes, [\"x\", \"laplacian4\"], [\"x\", \"4\"]):\n",
    "\n",
    "    # calculate derivative matrix\n",
    "    dmat = compute_dmat(mapping_array, kernel_type=kernel_type)\n",
    "\n",
    "    # plot sparse matrix\n",
    "    ax.spy(dmat, markersize=2)\n",
    "    ax.set_title(f\"derivative matrix $K_{kernel_name}$\", pad=25)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "show laplacian matrix $K_4$ in (10, 10) size as a numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmat.toarray()[0:10, 0:10]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply the derivative matrix to a sample image\n",
    "---\n",
    "Next, let us to apply derivative matrices to pixels of a sample image."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load sample image data from the matplotlib library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with get_sample_data(\"grace_hopper.jpg\") as file:\n",
    "    arr_image = plt.imread(file)\n",
    "\n",
    "# resize the image deu to the large size.\n",
    "with Image.fromarray(arr_image, mode=\"RGB\") as im:\n",
    "    (width, height) = (im.width // 4, im.height // 4)\n",
    "    arr_image = np.array(im.resize((width, height)))\n",
    "\n",
    "# convert RGB image to monotonic one\n",
    "arr_image = arr_image.mean(axis=2)\n",
    "\n",
    "# show image\n",
    "print(f\"image array shape: {arr_image.shape}\")\n",
    "fig, ax = plt.subplots(dpi=150)\n",
    "ax.imshow(arr_image, cmap=\"gray\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create mapping array\n",
    "image_map = np.arange(0, arr_image.size, dtype=np.int32).reshape(arr_image.shape)\n",
    "\n",
    "# create derivative matrix with Kx\n",
    "dmatx = compute_dmat(image_map, kernel_type=\"x\")\n",
    "\n",
    "# create derivative matrix with Ky\n",
    "dmaty = compute_dmat(image_map, kernel_type=\"y\")\n",
    "\n",
    "# create laplacian matrix with K8\n",
    "laplacian_mat = compute_dmat(image_map, kernel_type=\"laplacian8\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The filtered images are calculated by multiplying the image vector by each derivative matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_x = np.reshape(dmatx @ arr_image.ravel(), arr_image.shape)\n",
    "filtered_y = np.reshape(dmaty @ arr_image.ravel(), arr_image.shape)\n",
    "filtered_laplacian = np.reshape(laplacian_mat @ arr_image.ravel(), arr_image.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare $K_x$ and $K_y$ kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# extract max and min value\n",
    "vmax = max(filtered_x.max(), filtered_y.max())\n",
    "vmin = min(filtered_x.min(), filtered_y.min())\n",
    "half_range = max(abs(vmax), abs(vmin))\n",
    "norm = CenteredNorm(vcenter=0, halfrange=half_range)\n",
    "\n",
    "# show each image\n",
    "fig = plt.figure(dpi=150)\n",
    "grids = ImageGrid(fig, 111, nrows_ncols=(1, 2), axes_pad=0.1, cbar_mode=\"single\")\n",
    "for ax, filtered_image, title in zip(grids, [filtered_x, filtered_y], [\"$K_x$\", \"$K_y$\"]):\n",
    "    mappable = ax.imshow(filtered_image, cmap=\"seismic\", norm=norm)\n",
    "    ax.set_title(title)\n",
    "\n",
    "cbar = plt.colorbar(mappable, cax=grids.cbar_axes[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the laplacian filtered image $K_8$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(dpi=150)\n",
    "mappable = ax.imshow(filtered_laplacian, cmap=\"seismic\", norm=CenteredNorm())\n",
    "ax.set_title(\"$K_8$\")\n",
    "cbar = plt.colorbar(mappable)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results show that the edge of the image is emphasized clearly.\n",
    "So, we take advantage of this operator to smooth tomographic reconstructions."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cherab-phix-dev",
   "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.9.16"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "2725905a4c02db19e04df9b8fdbbe5ec65a73ea52bebaf9474aa1cc98819834c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}