{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Derivative Operation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Here we show you examples of derivative/laplacian matrices and how they work on the example image.\n",
    "\n",
    "Please see the [derivative theory page](../../user/theory/derivative.ipynb) if you want to know the\n",
    "details of the theory."
   ]
  },
  {
   "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.cm import ScalarMappable\n",
    "from matplotlib.colors import CenteredNorm\n",
    "from mpl_toolkits.axes_grid1 import ImageGrid\n",
    "from PIL import Image\n",
    "\n",
    "from cherab.inversion.derivative import derivative_matrix, laplacian_matrix"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize derivative matrix\n",
    "\n",
    "Let us compute the simple derivative matrix applying a 3 x 3 image array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly we show derivative matrices along to vertical direction (`axis=0` direction) with different schemes.\n",
    "The derivative matrices are shown by printing the matrix elements and plotting the matrix elements\n",
    "in 2-D image. The blue squares represent the matrix elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, dpi=150, tight_layout=True)\n",
    "for ax, scheme in zip(axes, [\"forward\", \"backward\", \"central\"]):  # noqa\n",
    "    # compute derivative matrix\n",
    "    dmat = derivative_matrix((3, 3), axis=0, scheme=scheme)\n",
    "\n",
    "    # print array\n",
    "    print(f\"Derivative matrix ({scheme}):\\n{dmat.toarray()}\\n\")\n",
    "\n",
    "    # plot derivative matrix as a sparse one\n",
    "    ax.spy(dmat)\n",
    "    ax.set_title(f\"{scheme}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also we show the horizontal derivative matrix (`axis=1` direction) case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, dpi=150, tight_layout=True)\n",
    "for ax, scheme in zip(axes, [\"forward\", \"backward\", \"central\"]):  # noqa\n",
    "    # compute derivative matrix\n",
    "    dmat = derivative_matrix((3, 3), axis=1, scheme=scheme)\n",
    "\n",
    "    # print array\n",
    "    print(f\"Derivative matrix ({scheme}):\\n{dmat.toarray()}\\n\")\n",
    "\n",
    "    # plot derivative matrix as a sparse one\n",
    "    ax.spy(dmat)\n",
    "    ax.set_title(f\"{scheme}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that some elements are zero due to the boundary condition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize laplacian matrix\n",
    "\n",
    "Let us compute the simple laplacian matrix applying a 3 x 3 image array and show it by printing the matrix elements and plotting the matrix elements in 2-D image.\n",
    "There are two types of laplacian matrix: difference between only orthogonal neighbors or all neighbors\n",
    "including diagonal neighbors. The step size along diagonal direction is $\\sqrt{2}h$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, dpi=150, tight_layout=True)\n",
    "for ax, diagonal in zip(axes, [False, True]):  # noqa\n",
    "    # compute laplacian matrix\n",
    "    lmat = laplacian_matrix((3, 3), diagonal=diagonal)\n",
    "\n",
    "    # plot laplacian matrix as a sparse one\n",
    "    ax.spy(lmat)\n",
    "    title = \"diagnoal\" if diagonal else \"no diagnoal\"\n",
    "    ax.set_title(f\"{title}\")\n",
    "\n",
    "    # print array\n",
    "    print(f\"Laplacian matrix with {title}:\\n{lmat.toarray()}\\n\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Apply the derivative/laplacian matrix to a sample image\n",
    "\n",
    "Next, let us to apply derivative/laplacian matrices to pixels of a sample image.\n",
    "We use sample image stored in `matplotlib` library.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us preview the sample image converted to monochrome image.\n"
   ]
  },
  {
   "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 to 1/4 of the original 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",
    "im = ax.imshow(arr_image, cmap=\"gray\")\n",
    "ax.set_xlabel(\"$x$\")\n",
    "ax.set_ylabel(\"$y$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then create the derivative/laplacian matrices.\n",
    "Note that the $x$, $y$ directions correspond to the 1, 0 axis of the image array, respectively, and\n",
    "the positive direction of the $y$ axis is downward."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute each derivative matrix along each axis with forward scheme\n",
    "dmat_x = derivative_matrix(arr_image.shape, axis=1, scheme=\"forward\")\n",
    "dmat_y = derivative_matrix(arr_image.shape, axis=0, scheme=\"forward\")\n",
    "\n",
    "# compute derivative matrices along y direction with different schemes\n",
    "dmat_y_backward = derivative_matrix(arr_image.shape, axis=0, scheme=\"backward\")\n",
    "dmat_y_central = derivative_matrix(arr_image.shape, axis=0, scheme=\"central\")\n",
    "\n",
    "# compute laplacian matrices with and without diagonal\n",
    "lmat_diag = laplacian_matrix(arr_image.shape, diagonal=True)\n",
    "lmat_no_diag = laplacian_matrix(arr_image.shape, diagonal=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check the difference effect of the derivative matrix between $x$ and $y$ directions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [],
   "source": [
    "# compute derivative of images\n",
    "arr_derivatives = []\n",
    "vmax = 0\n",
    "for dmat in [dmat_x, dmat_y]:\n",
    "    arr_derivative = dmat @ arr_image.flatten()\n",
    "    arr_derivative = arr_derivative.reshape(arr_image.shape)\n",
    "\n",
    "    # retrieve abosolute maximum value of derivative and choose the largest one\n",
    "    vmax = val if (val := np.abs(arr_derivative).max()) > vmax else vmax\n",
    "    arr_derivatives.append(arr_derivative)\n",
    "\n",
    "fig = plt.figure(dpi=150)\n",
    "grids = ImageGrid(fig, 111, nrows_ncols=(1, 2), cbar_mode=\"single\")\n",
    "norm = CenteredNorm(halfrange=vmax)\n",
    "for ax, arr, title in zip(grids, arr_derivatives, [\"$x$\", \"$y$\"]):  # noqa\n",
    "    # plot derivative\n",
    "    im = ax.imshow(arr, cmap=\"RdBu_r\", norm=norm)\n",
    "    ax.set_xlabel(\"$x$\")\n",
    "    ax.set_ylabel(\"$y$\")\n",
    "    ax.set_title(f\"{title} direction\")\n",
    "\n",
    "# plot colorbar\n",
    "mappable = ScalarMappable(norm=norm, cmap=\"RdBu_r\")\n",
    "cbar = plt.colorbar(mappable, cax=grids.cbar_axes[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that each derivative image varies depending on each image's gradient direction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check the difference derivative scheme (forward, backward, central)\n",
    "\n",
    "Let us check the difference between derivative schemes in the y direction.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# compute derivative of images\n",
    "arr_derivatives = []\n",
    "vmax = 0\n",
    "for dmat in [dmat_y, dmat_y_backward, dmat_y_central]:\n",
    "    arr_derivative = dmat @ arr_image.flatten()\n",
    "    arr_derivative = arr_derivative.reshape(arr_image.shape)\n",
    "\n",
    "    # retrieve abosolute maximum value of derivative and choose the largest one\n",
    "    vmax = val if (val := np.abs(arr_derivative).max()) > vmax else vmax\n",
    "    arr_derivatives.append(arr_derivative)\n",
    "\n",
    "fig = plt.figure(dpi=150)\n",
    "grids = ImageGrid(fig, 111, nrows_ncols=(1, 3), cbar_mode=\"single\")\n",
    "norm = CenteredNorm(halfrange=vmax)\n",
    "for ax, arr, title in zip(grids, arr_derivatives, [\"forward\", \"backward\", \"central\"]):  # noqa\n",
    "    # plot derivative\n",
    "    im = ax.imshow(arr, cmap=\"RdBu_r\", norm=norm)\n",
    "    ax.set_xlabel(\"$x$\")\n",
    "    ax.set_ylabel(\"$y$\")\n",
    "    ax.set_title(f\"{title}\")\n",
    "\n",
    "# plot colorbar\n",
    "mappable = ScalarMappable(norm=norm, cmap=\"RdBu_r\")\n",
    "cbar = plt.colorbar(mappable, cax=grids.cbar_axes[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no apparent difference between the derivative schemes.\n",
    "The difference seems be more apparent when the gradient is large or the area near the boundary because\n",
    "each scheme has different boundary condition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check the difference effect of the laplacian matrix w/ or w/o diagonal neighbors\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(dpi=150)\n",
    "grids = ImageGrid(fig, 111, nrows_ncols=(1, 2), cbar_mode=\"each\", axes_pad=0.6, cbar_pad=0.0)\n",
    "for ax, cax, lmat, title in zip(grids, grids.cbar_axes, [lmat_no_diag, lmat_diag], [\"w/o diagonal\", \"w/ diagonal\"]):  # noqa\n",
    "    # compute laplacian of image\n",
    "    arr_laplacian = lmat @ arr_image.flatten()\n",
    "    arr_laplacian = arr_laplacian.reshape(arr_image.shape)\n",
    "\n",
    "    # plot laplacian of image\n",
    "    im = ax.imshow(arr_laplacian, cmap=\"RdBu_r\", norm=CenteredNorm())\n",
    "    ax.set_xlabel(\"$x$\")\n",
    "    ax.set_ylabel(\"$y$\")\n",
    "    ax.set_title(f\"{title}\")\n",
    "\n",
    "    cbar = plt.colorbar(im, cax=cax)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both laplacian matrices can detect the edge of the image.\n",
    "The laplacian with diagonal neighbors generates larger values than that without diagonal neighbors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In image processing, the both derivative and laplacian matrices are used for edge detection.\n",
    "In the inversion problem, specific to the tomography, these matrices make a reconstructed image\n",
    "smoother by adding the derivative/laplacian regularization term to the objective function."
   ]
  }
 ],
 "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.11.7"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "2725905a4c02db19e04df9b8fdbbe5ec65a73ea52bebaf9474aa1cc98819834c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}