{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "23356cce",
   "metadata": {},
   "source": [
    "2.4 Extension"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a89f5c2b",
   "metadata": {},
   "source": [
    "Further, notice how the bright regions in the Fourier-domain image are concentrated in a small set of frequency components. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c75bc83",
   "metadata": {},
   "outputs": [],
   "source": [
    "fourier_img = fourier_transform(img.flatten(), Q)\n",
    "fourier_img = M_fourier @ S_fourier @ fourier_img  # Feel free to experiment here!\n",
    "pixel_img = inv_fourier_transform(fourier_img, Q)\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.imshow(fourier_img.reshape(img.shape), cmap=\"gray\")\n",
    "plt.title(\"Fourier domain\")\n",
    "plt.axis(\"off\")\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.imshow(pixel_img.reshape(img.shape), cmap=\"gray\")\n",
    "plt.title(\"Pixel domain\")\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e9a9f2f",
   "metadata": {},
   "source": [
    "3.2 Extension"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ae38358",
   "metadata": {},
   "source": [
    "Implementing a numerically stable second-order Gaussian filter can be complicated. Compare your implementation, which should return a constant zero-valued second-derivative in the y-direction, with a package that implements Gaussian filters incorrectly... `scipy`!\n",
    "\n",
    "`scipy`'s implementation has errors of magnitude $10^3$, while yours should be < $10^{-13}$.\n",
    "\n",
    "Why does this occur? This discrepancy occurs because computing higher-order derivatives of Gaussian filters is numerically sensitive, especially when the filter is truncated to a finite length. In theory, the second derivative of a linear function is exactly zero. Since the test image is constructed so that its second derivative in the $y$-direction is zero, an ideal Gaussian second-derivative filter should produce values very close to zero everywhere. However, in practice, Gaussian filters must be sampled and truncated. If derivatives are computed by repeatedly applying finite-difference or discrete differentiation rules, small approximation errors are introduced. These errors accumulate and are amplified when taking higher-order derivatives, leading to large numerical inaccuracies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "270b3620",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct an image with a zero-valued second-derivative in the column-direction\n",
    "img = torch.zeros(1, 1, 31, 31, dtype=torch.float64)\n",
    "img[..., 7:24] += torch.arange(1, 18)\n",
    "img[..., 24:] += 17\n",
    "\n",
    "for sigma in [0.5]: # Modify as desired\n",
    "    # Apply your second-derivative filter\n",
    "    our_second_derivative = second_order_yy(img, sigma=sigma)\n",
    "\n",
    "    # Apply the second-derivative filter using scipy\n",
    "    scipy_second_derivative = gaussian_filter_scipy(\n",
    "        img.numpy().squeeze(), sigma=sigma, order=[2, 0]\n",
    "    )\n",
    "\n",
    "    # Plot the original image and the filtered images\n",
    "    plt.figure(figsize=(15, 3.75))\n",
    "    plt.subplot(131)\n",
    "    plt.imshow(img.squeeze(), cmap=\"gray\")\n",
    "    plt.title(\"Original Image\")\n",
    "    plt.colorbar()\n",
    "\n",
    "    plt.subplot(132)\n",
    "    plt.imshow(our_second_derivative.squeeze(), cmap=\"turbo\")\n",
    "    plt.title(\"Your Second Derivative\")\n",
    "    plt.colorbar()\n",
    "\n",
    "    plt.subplot(133)\n",
    "    plt.imshow(scipy_second_derivative, cmap=\"turbo\")\n",
    "    plt.title(\"Scipy Second Derivative\")\n",
    "    plt.colorbar()\n",
    "\n",
    "    plt.suptitle(f\"Sigma = {sigma}\")\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2aedb332",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 1\n",
    "fig, axs = imshow(\n",
    "    img := astronaut(),\n",
    "    zeroth_order(img, sigma),\n",
    "    first_order_x(img, sigma),\n",
    "    first_order_y(img, sigma),\n",
    "    second_order_xx(img, sigma),\n",
    "    second_order_yy(img, sigma),\n",
    ")\n",
    "axs[0].set_title(\"Original image\")\n",
    "axs[1].set_title(\"Smoothed\")\n",
    "axs[2].set_title(\"First x-derivative\")\n",
    "axs[3].set_title(\"First y-derivative\")\n",
    "axs[4].set_title(\"Second x-derivative\")\n",
    "axs[5].set_title(\"Second y-derivative\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7eb18823",
   "metadata": {},
   "source": [
    "3.3 Extension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7545ed2",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 1\n",
    "fig, axs = imshow(\n",
    "    img := astronaut(),\n",
    "    log(img, sigma),\n",
    ")\n",
    "axs[0].set_title(\"Original image\")\n",
    "axs[1].set_title(\"LoG filtered image\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e353a29",
   "metadata": {},
   "source": [
    "3.4 Extension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72f51b4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Verify your 2D filters are properly oriented\n",
    "fig, axs = plt.subplots(2, 4, figsize=(8, 4))\n",
    "for idx, ax in enumerate(axs.flatten()):\n",
    "    ax.imshow(oriented_filter(idx * torch.pi / 4, 2.5), cmap=\"bwr\")\n",
    "    ax.axis(\"off\")\n",
    "    ax.set_title(f\"$\\\\theta = {idx} \\\\pi / 4$\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".conda",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.8.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}