{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c59a354",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from skimage import data, img_as_float, exposure, filters, feature\n",
    "from scipy import ndimage\n",
    "\n",
    "image = img_as_float(data.camera())\n",
    "print(\"Shape:\", image.shape, \"| range:\", (image.min(), image.max()), \"| dtype:\", image.dtype)\n",
    "\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.imshow(image, cmap=\"gray\", vmin=0, vmax=1)\n",
    "plt.title(\"Original grayscale image\")\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01bde4bc",
   "metadata": {},
   "source": [
    "1. What does the image shape tell you?\n",
    "2. Is intensity 0.9 brighter or darker than 0.2?\n",
    "3. Why is it useful to keep the display range fixed at 0–1?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a50a1e8d",
   "metadata": {},
   "source": [
    "For normalized intensity \\(r\\), the negative is \\(s=1-r\\). Black becomes white and white becomes black."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33f64895",
   "metadata": {},
   "outputs": [],
   "source": [
    "negative = 1 - image\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
    "axes[0].imshow(image, cmap=\"gray\", vmin=0, vmax=1); axes[0].set_title(\"Original\")\n",
    "axes[1].imshow(negative, cmap=\"gray\", vmin=0, vmax=1); axes[1].set_title(\"Negative: 1 − image\")\n",
    "for ax in axes: ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()\n",
    "\n",
    "example = np.array([0.0, 0.2, 0.5, 0.8, 1.0])\n",
    "print(\"Input: \", example)\n",
    "print(\"Output:\", 1 - example)\n",
    "print(\"Applying negative twice, max error:\", np.max(np.abs((1-negative)-image)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d596eb2d",
   "metadata": {},
   "source": [
    "- Why does intensity 0.5 remain unchanged?\n",
    "- Apply the negative operation twice. Is the result exactly the same as the original?\n",
    "- Try it on another grayscale image."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f4ba60e",
   "metadata": {},
   "source": [
    "A log transform can expand differences among dark values while compressing differences among bright values. For normalized input we use \\(s=\\log(1+r)/\\log(2)\\), keeping output in [0, 1]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2493c61a",
   "metadata": {},
   "outputs": [],
   "source": [
    "log_image = np.log1p(image) / np.log(2)\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
    "axes[0].imshow(image, cmap=\"gray\", vmin=0, vmax=1); axes[0].set_title(\"Original\")\n",
    "axes[1].imshow(log_image, cmap=\"gray\", vmin=0, vmax=1); axes[1].set_title(\"Log transformed\")\n",
    "for ax in axes: ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()\n",
    "\n",
    "r = np.linspace(0, 1, 300)\n",
    "plt.figure(figsize=(6, 4))\n",
    "for c in [1, 5, 20]:\n",
    "    plt.plot(r, np.log1p(c*r)/np.log1p(c), label=f\"c={c}\")\n",
    "plt.plot(r, r, \"--\", label=\"Identity\")\n",
    "plt.xlabel(\"Input intensity\"); plt.ylabel(\"Output intensity\")\n",
    "plt.title(\"Log transformation curves\"); plt.grid(True); plt.legend(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fe79590",
   "metadata": {},
   "source": [
    "Change the `c` values to 0.5, 2, 10, and 50. Does the transform tend to brighten or darken the image? What does the curve tell you?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c679d0c2",
   "metadata": {},
   "source": [
    "For normalized intensities, gamma below 1 usually brightens; gamma above 1 usually darkens; gamma equal to 1 leaves the image unchanged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af957089",
   "metadata": {},
   "outputs": [],
   "source": [
    "gammas = [0.4, 1.0, 2.2]\n",
    "fig, axes = plt.subplots(1, len(gammas), figsize=(15, 5))\n",
    "for ax, gamma in zip(axes, gammas):\n",
    "    ax.imshow(image**gamma, cmap=\"gray\", vmin=0, vmax=1)\n",
    "    ax.set_title(f\"gamma={gamma}\"); ax.axis(\"off\")\n",
    "plt.suptitle(\"Gamma transformation\"); plt.tight_layout(); plt.show()\n",
    "\n",
    "r = np.linspace(0, 1, 300)\n",
    "plt.figure(figsize=(6, 4))\n",
    "for gamma in gammas: plt.plot(r, r**gamma, label=f\"gamma={gamma}\")\n",
    "plt.plot(r, r, \"--\", label=\"Identity\")\n",
    "plt.xlabel(\"Input intensity\"); plt.ylabel(\"Output intensity\")\n",
    "plt.title(\"Gamma curves\"); plt.grid(True); plt.legend(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12f75e9c",
   "metadata": {},
   "source": [
    "Add gamma values 0.7 and 1.5. Which makes dark regions easier to see? Why is gamma=1 a useful control?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b6bc139",
   "metadata": {},
   "source": [
    "A low-contrast image uses only a narrow part of the available intensity range. Min–max stretching maps its minimum to 0 and maximum to 1.\n",
    "This method can be sensitive to outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36a7425e",
   "metadata": {},
   "outputs": [],
   "source": [
    "low_contrast = 0.35 + 0.30*image\n",
    "lo, hi = low_contrast.min(), low_contrast.max()\n",
    "stretched = (low_contrast-lo)/(hi-lo)\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
    "for ax, im, title in zip(axes, [image, low_contrast, stretched],\n",
    "                         [\"Original\", \"Artificial low contrast\", \"Min–max stretched\"]):\n",
    "    ax.imshow(im, cmap=\"gray\", vmin=0, vmax=1); ax.set_title(title); ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()\n",
    "print(\"Low contrast range:\", (low_contrast.min(), low_contrast.max()))\n",
    "print(\"Stretched range:\", (stretched.min(), stretched.max()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d510cd2",
   "metadata": {},
   "source": [
    "Instead of using the absolute minimum and maximum, use percentiles to reduce the influence of extreme values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0052d890",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_lo, p_hi = np.percentile(low_contrast, (2, 98))\n",
    "percentile_stretched = exposure.rescale_intensity(low_contrast, in_range=(p_lo, p_hi), out_range=(0, 1))\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
    "axes[0].imshow(stretched, cmap=\"gray\", vmin=0, vmax=1); axes[0].set_title(\"Min–max stretch\")\n",
    "axes[1].imshow(percentile_stretched, cmap=\"gray\", vmin=0, vmax=1); axes[1].set_title(\"2nd–98th percentile stretch\")\n",
    "for ax in axes: ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()\n",
    "print(\"Percentile limits:\", p_lo, p_hi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31326627",
   "metadata": {},
   "source": [
    "Try percentile ranges `(1,99)`, `(5,95)`, and `(10,90)`. Look for stronger contrast, clipping, and lost details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28068e50",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change these percentiles and rerun\n",
    "p_lo, p_hi = np.percentile(low_contrast, (5, 95))\n",
    "custom_stretch = exposure.rescale_intensity(low_contrast, in_range=(p_lo, p_hi), out_range=(0, 1))\n",
    "plt.figure(figsize=(6, 5)); plt.imshow(custom_stretch, cmap=\"gray\", vmin=0, vmax=1)\n",
    "plt.title(f\"Percentile stretch ({p_lo:.2f}–{p_hi:.2f})\"); plt.axis(\"off\"); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "572549e9",
   "metadata": {},
   "source": [
    "Thresholding converts pixels to two groups: white when intensity is above threshold \\(T\\), otherwise black. It is a simple segmentation method, but can struggle with uneven illumination or overlapping foreground/background intensities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fff6b9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 0.5\n",
    "binary = image > threshold\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
    "axes[0].imshow(image, cmap=\"gray\", vmin=0, vmax=1); axes[0].set_title(\"Original\")\n",
    "axes[1].imshow(binary, cmap=\"gray\"); axes[1].set_title(f\"Threshold > {threshold}\")\n",
    "for ax in axes: ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()\n",
    "print(\"Fraction of pixels above threshold:\", binary.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1cc09dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholds = [0.25, 0.4, 0.6, 0.75]\n",
    "fig, axes = plt.subplots(1, len(thresholds), figsize=(16, 4))\n",
    "for ax, t in zip(axes, thresholds):\n",
    "    mask = image > t\n",
    "    ax.imshow(mask, cmap=\"gray\"); ax.set_title(f\"T={t}\\nwhite={mask.mean():.1%}\"); ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39a1c62c",
   "metadata": {},
   "source": [
    "Highlight only pixels in a selected intensity range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "534e1a10",
   "metadata": {},
   "outputs": [],
   "source": [
    "low, high = 0.35, 0.65\n",
    "slice_mask = (image >= low) & (image <= high)\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
    "axes[0].imshow(image, cmap=\"gray\", vmin=0, vmax=1); axes[0].set_title(\"Original\")\n",
    "axes[1].imshow(slice_mask, cmap=\"gray\"); axes[1].set_title(f\"Mask: {low}–{high}\")\n",
    "axes[2].imshow(slice_mask.astype(float), cmap=\"gray\"); axes[2].set_title(\"Selected range\")\n",
    "for ax in axes: ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9dad37a",
   "metadata": {},
   "source": [
    "Create a panel containing original, negative, log transform, two gamma transforms, contrast-stretched image, and threshold mask. Write 3–5 sentences: Which transformation best reveals the detail you care about, and what information might it suppress?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7e288c6",
   "metadata": {},
   "source": [
    "An edge is a rapid change in intensity. The gradient has horizontal and vertical components, \\(G_x\\) and \\(G_y\\), and magnitude \\(G=\\sqrt{G_x^2+G_y^2}\\)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffc3d2c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6, 6))\n",
    "plt.imshow(image[100:300, 100:300], cmap=\"gray\", vmin=0, vmax=1)\n",
    "plt.title(\"Image crop: inspect boundaries and texture\"); plt.axis(\"off\"); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faa230c7",
   "metadata": {},
   "source": [
    "Sobel estimates intensity changes in two directions. We combine responses into a gradient magnitude. A strong response means a strong local change in intensity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33356b2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sobel_h = filters.sobel_h(image)\n",
    "sobel_v = filters.sobel_v(image)\n",
    "sobel_mag = np.hypot(sobel_h, sobel_v)\n",
    "\n",
    "fig, axes = plt.subplots(1, 4, figsize=(18, 5))\n",
    "items = [(image, \"Original\"), (np.abs(sobel_h), \"Sobel horizontal response\"),\n",
    "         (np.abs(sobel_v), \"Sobel vertical response\"), (sobel_mag, \"Gradient magnitude\")]\n",
    "for ax, (im, title) in zip(axes, items):\n",
    "    ax.imshow(im, cmap=\"gray\"); ax.set_title(title); ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0c1ba00",
   "metadata": {},
   "source": [
    "`np.hypot(a,b)` computes `sqrt(a**2+b**2)` element by element.\n",
    "\n",
    "1. Compare the directional responses.\n",
    "2. Display a crop of `sobel_mag`.\n",
    "3. What information is hidden if you display only the absolute value of a signed response?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c0e97cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6, 5))\n",
    "plt.imshow(sobel_mag[100:300, 100:300], cmap=\"gray\")\n",
    "plt.title(\"Sobel magnitude crop\"); plt.axis(\"off\"); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfc06d2c",
   "metadata": {},
   "source": [
    "Prewitt has the same broad goal as Sobel but uses a different weighting pattern. Compare their results using the same image and display approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4be4dd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "prewitt_h = filters.prewitt_h(image)\n",
    "prewitt_v = filters.prewitt_v(image)\n",
    "prewitt_mag = np.hypot(prewitt_h, prewitt_v)\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
    "for ax, im, title in zip(axes, [sobel_mag, prewitt_mag, np.abs(sobel_mag-prewitt_mag)],\n",
    "                         [\"Sobel magnitude\", \"Prewitt magnitude\", \"Absolute difference\"]):\n",
    "    ax.imshow(im, cmap=\"gray\"); ax.set_title(title); ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f169da91",
   "metadata": {},
   "source": [
    "Which edges are visible in both outputs? Do not assume the visually brighter output is automatically better; display scaling affects appearance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "782e05ee",
   "metadata": {},
   "source": [
    "The Laplacian uses second-order derivatives and responds to rapid intensity transitions. It can be sensitive to noise. The signed response contains positive and negative values, so we visualize it with a diverging colour map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e1703d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "lap = ndimage.laplace(image)\n",
    "limit = np.percentile(np.abs(lap), 99)\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
    "axes[0].imshow(image, cmap=\"gray\", vmin=0, vmax=1); axes[0].set_title(\"Original\")\n",
    "axes[1].imshow(lap, cmap=\"seismic\", vmin=-limit, vmax=limit); axes[1].set_title(\"Signed Laplacian\")\n",
    "axes[2].imshow(np.abs(lap), cmap=\"gray\", vmin=0, vmax=limit); axes[2].set_title(\"Laplacian magnitude (clipped)\")\n",
    "for ax in axes: ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b46bc57",
   "metadata": {},
   "source": [
    "Why use a diverging colour map? Why can second derivatives respond strongly to noise? How does the Laplacian differ from gradient magnitude?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55c058cc",
   "metadata": {},
   "source": [
    "Canny combines smoothing, gradient estimation, non-maximum suppression, and hysteresis thresholding. It produces a binary edge map. The `sigma` argument controls internal Gaussian smoothing: larger values suppress fine texture but may remove small edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f91c83ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "canny = feature.canny(image, sigma=1.0)\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
    "axes[0].imshow(image, cmap=\"gray\", vmin=0, vmax=1); axes[0].set_title(\"Original\")\n",
    "axes[1].imshow(canny, cmap=\"gray\"); axes[1].set_title(\"Canny, sigma=1.0\")\n",
    "for ax in axes: ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()\n",
    "print(\"Fraction of pixels marked as edges:\", canny.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fb8ace1",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmas = [0.5, 1.0, 2.0, 3.0]\n",
    "fig, axes = plt.subplots(1, len(sigmas), figsize=(16, 4))\n",
    "for ax, sigma in zip(axes, sigmas):\n",
    "    edges = feature.canny(image, sigma=sigma)\n",
    "    ax.imshow(edges, cmap=\"gray\")\n",
    "    ax.set_title(f\"sigma={sigma}\\nedge pixels={edges.mean():.1%}\")\n",
    "    ax.axis(\"off\")\n",
    "plt.suptitle(\"Canny: effect of smoothing scale\"); plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d07f0d6",
   "metadata": {},
   "source": [
    "Change Canny's `sigma`. Compare fine texture, major boundaries, and isolated edge responses. Does increasing sigma always improve the result? Explain the trade-off."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a800b82",
   "metadata": {},
   "source": [
    "Noise creates rapid changes that can look like edges. Smooth first, then compute the Sobel magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4953ab7",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaussian_sigmas = [0.0, 0.8, 1.5, 2.5]\n",
    "fig, axes = plt.subplots(2, len(gaussian_sigmas), figsize=(16, 8))\n",
    "for col, sigma in enumerate(gaussian_sigmas):\n",
    "    smooth = image.copy() if sigma == 0 else ndimage.gaussian_filter(image, sigma=sigma)\n",
    "    grad = filters.sobel(smooth)\n",
    "    axes[0, col].imshow(smooth, cmap=\"gray\", vmin=0, vmax=1)\n",
    "    axes[0, col].set_title(f\"Smoothed, sigma={sigma}\")\n",
    "    axes[1, col].imshow(grad, cmap=\"gray\")\n",
    "    axes[1, col].set_title(\"Sobel magnitude\")\n",
    "    axes[0, col].axis(\"off\"); axes[1, col].axis(\"off\")\n",
    "plt.suptitle(\"Gaussian smoothing before edge detection\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e46255b0",
   "metadata": {},
   "source": [
    "Smoothing can remove noise and texture, but can also soften or erase small structures. The best sigma depends on the image and task."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51994e02",
   "metadata": {},
   "source": [
    "We add synthetic Gaussian noise, then compare edge responses before and after smoothing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9cc7e3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(7)\n",
    "noise_std = 0.08\n",
    "noisy = np.clip(image + rng.normal(0, noise_std, image.shape), 0, 1)\n",
    "edges_noisy = filters.sobel(noisy)\n",
    "denoised = ndimage.gaussian_filter(noisy, sigma=1.2)\n",
    "edges_denoised = filters.sobel(denoised)\n",
    "\n",
    "fig, axes = plt.subplots(1, 4, figsize=(18, 5))\n",
    "items = [(image, \"Original\"), (noisy, \"Noisy image\"),\n",
    "         (edges_noisy, \"Edges before smoothing\"), (edges_denoised, \"Edges after smoothing\")]\n",
    "for ax, (im, title) in zip(axes, items):\n",
    "    ax.imshow(im, cmap=\"gray\", vmin=0, vmax=1 if title in [\"Original\", \"Noisy image\"] else None)\n",
    "    ax.set_title(title); ax.axis(\"off\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6521b57",
   "metadata": {},
   "source": [
    "Create a comparison panel with original, Sobel magnitude, Prewitt magnitude, Laplacian magnitude, and Canny outputs at two sigma values."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "MMDP_Image_Labs_06_07_Intensity_and_Edge_Detection.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.x"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}