{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from skimage.color import rgb2gray, gray2rgb\n",
    "from skimage.io import imread, imsave, imshow\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "# Increase these if figures appear small\n",
    "plt.rcParams[\"figure.figsize\"] = fx, fy = (18.40, 10.35)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Познакомимся с данными"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Загрузим и посмотрим на изображение"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!curl -O 'https://courses.cv-gml.ru/storage/seminars/basic-image-processing/msu.jpg'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = imread(\"msu.jpg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[img.dtype, img.shape]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[img.min(), img.mean(), img.max()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(img)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Что-то небо слишком синее. Давайте попробуем сделать его потемнее."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_dark_sky = img.copy()\n",
    "img_dark_sky[..., 2] -= 100\n",
    "# this doesn't work:\n",
    "# img_dark_sky = img_dark_sky.clip(0, None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(img_dark_sky)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ой! Деревья стали синими! Что пошло не так?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Давайте лучше будем работать в `float32` в диапазоне $\\left[0.0;\\;1.0\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_f32 = img.astype(np.float32) / 255"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[img_f32.dtype, img_f32.shape]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[img_f32.min(), img_f32.mean(), img_f32.max()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(img_f32)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "А теперь?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_dark_sky_f32 = img_f32.copy()\n",
    "img_dark_sky_f32[..., 2] -= 100 / 255\n",
    "img_dark_sky_f32 = img_dark_sky_f32.clip(0, None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(img_dark_sky_f32)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сконвертируем изображение в черно-белое"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gray_f32 = rgb2gray(img_f32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[gray_f32.dtype, gray_f32.shape]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(gray_f32)\n",
    "# plt.imshow(gray_f32, cmap=\"viridis\")\n",
    "# plt.imshow(gray_f32, cmap=\"hot\")\n",
    "# plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ой! Что пошло не так?\n",
    "\n",
    "Давайте попробуем сохранить картинку и посмотреть на неё без `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# msu_gray_u8 = gray_f32.astype(np.uint8)\n",
    "# msu_gray_u8 = (255 * gray_f32).astype(np.uint8)\n",
    "msu_gray_u8 = (255 * gray_f32).clip(0, 255).round().astype(np.uint8)\n",
    "\n",
    "imsave(\"msu_gray.png\", msu_gray_u8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Откроем сохраненный файл [`msu_gray.png`](./msu_gray.png)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://courses.cv-gml.ru/storage/seminars/basic-image-processing/black_hole.jpg\" width=480 >\n",
    "\n",
    "> In the image, the dark circle represents the \"shadow\" of the black hole and its boundary, created by the glowing material that surrounds it. However, the colors of the bright ring in the image aren't the actual hues of the gas; rather, they represent a color map chosen by EHT researchers to depict the brightness of the emissions, Fox explained.\n",
    "\n",
    "Исправим проблему"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(gray_f32, cmap=\"gray\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь все нормально?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Давайте сконвертируем в черно-белое изображение сами"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_manual_gray_f32 = img_f32.mean(axis=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[simple_manual_gray_f32.dtype, simple_manual_gray_f32.shape]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "И сравним с библиотечным"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(gray_f32, cmap=\"gray\")\n",
    "plt.title(\"gray_f32\")\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(simple_manual_gray_f32, cmap=\"gray\")\n",
    "plt.title(\"simple_manual_gray_f32\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(img)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Трава на нашей версии почему-то чуть темнее"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "abs(simple_manual_gray_f32 - gray_f32).max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Хм. Давайте посмотрим, что говорят в документации библиотечной версии."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(rgb2gray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Frequently Asked Questions about Color](https://www.poynton.ca/pdf/ColourFAQ.pdf) by Charles Poyton."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R, G, B = img_f32.transpose(2, 0, 1)\n",
    "rec709_luma_f32 = 0.2125 * R + 0.7154 * G + 0.0721 * B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[rec709_luma_f32.dtype, rec709_luma_f32.shape]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(gray_f32, cmap=\"gray\")\n",
    "plt.title(\"gray_f32\")\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(rec709_luma_f32, cmap=\"gray\")\n",
    "plt.title(\"rec709_luma_f32\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(img_f32)\n",
    "plt.title(\"img_f32\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "abs(rec709_luma_f32 - gray_f32).max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вот, другое дело"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Но все-таки `cmap=\"gray\"` - точно было правильно?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(gray_f32, cmap=\"gray\")\n",
    "# plt.imshow(gray_f32 - 100, cmap=\"gray\")\n",
    "# plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(gray2rgb(gray_f32))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "manual_gray_rgb_img = gray_f32[..., None].repeat(3, axis=-1)\n",
    "\n",
    "plt.imshow(manual_gray_rgb_img)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Работа с гистограммой изображения\n",
    "\n",
    "Гистограмма изображения позволяет оценить количество и разнообразие оттенков, а также общий уровень яркости. Гистограмма изображения не описывает художественных качеств, но с её помощью можно рассчитывать разные статистики в картинке. Например найти миниум, максимум, среднее, медиану."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_image_with_hist(image):\n",
    "    f, [ax0, ax1] = plt.subplots(1, 2, figsize=(fx, 0.5 * fy))\n",
    "\n",
    "    # flatten image cause in image's histogram\n",
    "    # there is no need to keep geometric information\n",
    "    values = image.flatten()\n",
    "\n",
    "    # plot the histogram\n",
    "    # ax0.hist(values)\n",
    "    ax0.hist(values, bins=256, range=(0, 1), density=True)\n",
    "\n",
    "    # show the image\n",
    "    ax1.imshow(gray2rgb(image))\n",
    "\n",
    "\n",
    "plot_image_with_hist(gray_f32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Гамма-коррекция"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Можно выделить общую форму для всех операторов преобразования изображений, как функцию от одной или более картинок, которая отображает аргументы в итоговую картинку.\n",
    "\n",
    "$$ y = g(x_{0}, x_{1}, ..., x_{n}) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Самым простым видом трансформаций цветовых компотент изображений являются точечные операторы, которые изменяют каждый пиксель независимо.\n",
    "\n",
    "Так например, гамма коррекция задается формулой\n",
    "\n",
    "$$ y = g(x) = x^{\\gamma} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Во времена черно-белых телевизоров фосфор в электронно-лучевых трубках нелинейно реагировал на входное напряжение. Отношение между входным напряжением и получающейся яркостью можно было описать гаммой $ \\gamma $, поскольку с грубой оценкой выполнялось\n",
    "\n",
    "$$\n",
    "brightness = voltage^{\\gamma}\n",
    "$$\n",
    "\n",
    "Где $ \\gamma = 2.2 $.\n",
    "Чтобы обратить этот эффект телевизоры перед подачей кадра обрабатывали его через\n",
    "\n",
    "$$\n",
    "Y' = Y^{\\frac{1}{\\gamma}},\\qquad\\dfrac{1}{\\gamma} = 0.45.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma = 2.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данные `uint8`, которые мы считываем с диска с помощью `imread(\"msu.jpg\")` были скорректированы по гамма-коррекции конвейером обработки камеры. Давайте обратим эту операцию, чтобы получить исходные линейные (с точки зрения энергии фотонов) значения яркости:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linear_gray_f32 = gray_f32 ** gamma\n",
    "\n",
    "plot_image_with_hist(linear_gray_f32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Давайте также применим гамма-коррекцию ещё раз к исходному изображению (исключительно в демонстрационных целях):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma_correction = gray_f32 ** (1 / gamma)\n",
    "\n",
    "plot_image_with_hist(gamma_correction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Frequently Asked Questions about Gamma](https://www.poynton.ca/pdf/GammaFAQ.pdf) by Charles Poyton."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Контраст изображений"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Самым простым видом трансформаций цветовых компотент изображений являются точечные операторы, которые изменяют каждый пиксель независимо. Например понизить яркость картинки можно поделив каждый пиксель на константу.\n",
    "\n",
    "$$ value = value\\ /\\ 3 $$\n",
    "\n",
    "Реализуем это преобразование и посмотрим на результат"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "darkened_f32 = gray_f32 / 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_image_with_hist(darkened_f32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Так. Что-то опять не получается. Давайте исправим."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "darkened_u8 = (255 * darkened_f32).clip(0, 255).round().astype(np.uint8)\n",
    "darkened_u8_f32 = darkened_u8.astype(np.float32) / 255\n",
    "brigthened_u8_f32 = 3 * darkened_u8_f32"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_image_with_hist(brigthened_u8_f32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_image_with_hist(3 * darkened_f32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Линейное контрастирование\n",
    "\n",
    "$$ y = g(x) = \\frac{x - x_{min}}{x_{max} - x_{min}} $$\n",
    "\n",
    "Используется чтобы повысить контрастность изображения. Один из самых примитивных и не очень хорошо работающих методов."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lin_stretch_darkened_img = darkened_f32 - darkened_f32.min()\n",
    "lin_stretch_darkened_img = lin_stretch_darkened_img / lin_stretch_darkened_img.max()\n",
    "\n",
    "plot_image_with_hist(lin_stretch_darkened_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "darkened_with_outliers = darkened_u8_f32.copy()\n",
    "darkened_with_outliers[42, 42] = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lin_stretch_darkened_with_outliers = darkened_with_outliers - darkened_with_outliers.min()\n",
    "lin_stretch_darkened_with_outliers = lin_stretch_darkened_with_outliers / lin_stretch_darkened_with_outliers.max()\n",
    "\n",
    "plot_image_with_hist(lin_stretch_darkened_with_outliers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как ещё мы можем повысить контрастность изображения? Чтобы сгладить резкие углы, используется эквивализация гистограммы. Алгоритм для выравнивания гистограммы цветов."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Histogram equalization\n",
    "\n",
    "Как найти такое преобразование? Трюк тот же, что и при сэмплировании случайных чисел для разнообразных распределений с использованием функции распределения $F(x)$. Ответом является [отображение значений пикселей в соответствующие им квантили](https://en.wikipedia.org/wiki/Inverse_transform_sampling).\n",
    "\n",
    "Почему?\n",
    "> После такого преобразования кумилятивная функция распределения итоговых значений $F(x)$ в идеальной картине мира без квантизации будет $F(x) = x$.\n",
    ">\n",
    "> Следовательно плотность распределения пикселей (гистограмма), будет равномерной."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Алгоритм:\n",
    "1. Посчитать гистограмму цветов h\n",
    "2. Посчитать префиксные суммы массива h, оно и есть CDF(cumulative distribution function) цветов картинки\n",
    "3. Чтобы получить картинку с линеаризованной CDF используем трансформацию пикселей по правилу: $ y = CDF(x) $, чтобы получить квантиль для исходного значения пикселя x.\n",
    "\n",
    "Подробнее об алгоритме можно почитать на wiki [Histogram_equalization](https://en.wikipedia.org/wiki/Histogram_equalization#Implementation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_cdf(gray_image):\n",
    "    hist, _ = np.histogram(gray_image, bins=256, range=(0, 1), density=True)\n",
    "    cdf = hist.cumsum() / 256\n",
    "    return cdf\n",
    "\n",
    "\n",
    "def hist_equalization(gray_image):\n",
    "    cdf = calculate_cdf(gray_image)\n",
    "\n",
    "    # map all values to cdf's values\n",
    "    # what we'd like to be able do: gray_image = cdf(gray_image)\n",
    "    gray_image = np.interp(gray_image, np.linspace(0, 1, 256), cdf)\n",
    "    # for uint8 gray_image we can do: gray_image = cdf[gray_image]\n",
    "\n",
    "    return gray_image.clip(0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!curl -O 'https://courses.cv-gml.ru/storage/seminars/basic-image-processing/tests.py'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Тесты\n",
    "from tests import test_he\n",
    "\n",
    "test_he(hist_equalization)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "he_f32 = hist_equalization(gray_f32)\n",
    "\n",
    "plot_image_with_hist(he_f32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_compare(images, names):\n",
    "    [nrows] = {len(images), len(names)}\n",
    "    [ncols] = {len(row) for arr in [images, names] for row in arr}\n",
    "    f, plots = plt.subplots(\n",
    "        nrows,\n",
    "        ncols,\n",
    "        sharex=\"all\",\n",
    "        sharey=\"all\",\n",
    "        figsize=(fx, fy * nrows / ncols),\n",
    "    )\n",
    "\n",
    "    for idx, (row_imgs, row_names) in enumerate(zip(images, names, strict=True)):\n",
    "        for jdx, (img, name) in enumerate(zip(row_imgs, row_names, strict=True)):\n",
    "            if img.ndim != 3:\n",
    "                img = gray2rgb(img)\n",
    "            ax = plots[jdx] if nrows == 1 else plots[idx, jdx]\n",
    "            ax.imshow(img)\n",
    "            ax.set_title(name)\n",
    "            ax.axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_compare(\n",
    "    [(gray_f32, he_f32)],\n",
    "    [(\"original\", \"histogram equalization\")],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "he_darkened_with_outliers_f32 = hist_equalization(darkened_with_outliers)\n",
    "\n",
    "plot_image_with_hist(he_darkened_with_outliers_f32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_compare(\n",
    "    [\n",
    "        (gray_f32, darkened_with_outliers),\n",
    "        (lin_stretch_darkened_with_outliers, he_darkened_with_outliers_f32),\n",
    "    ],\n",
    "    [\n",
    "        (\"original\", \"darkened with outliers\"),\n",
    "        (\"linear stretching\", \"histogram equalization\"),\n",
    "    ],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Locally adaptive histogram equalization\n",
    "\n",
    "Глобальная гистограмма цветов может быть очень полезной, но в то же время для некоторых картинок предпочтительней применить различного вида эквивализации к разным регионам. Добиться этого можно с помощью разбиения картинки на блоки размера $(M, M)$ и применить $HE$ независимо для каждого блока. Если блоки будут не пересекаться, то в картинке могут появиться артефакты."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def naive_lahe_filt(image, block_size=32):\n",
    "    # for simplicity, lets crop the image to a multiple of the block size\n",
    "    ylen = image.shape[0] // block_size\n",
    "    xlen = image.shape[1] // block_size\n",
    "    image = image[: ylen * block_size, : xlen * block_size]\n",
    "\n",
    "    new_image = np.zeros_like(image)\n",
    "\n",
    "    for i in range(ylen):\n",
    "        y = block_size * i\n",
    "        for j in range(xlen):\n",
    "            x = block_size * j\n",
    "\n",
    "            # Extract the block\n",
    "            block = image[y : y + block_size, x : x + block_size]\n",
    "\n",
    "            # Equalize histogram for this block\n",
    "            block = hist_equalization(block)\n",
    "\n",
    "            # Write the block back\n",
    "            new_image[y : y + block_size, x : x + block_size] = block\n",
    "\n",
    "    return new_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "naive_lahe_f32 = naive_lahe_filt(gray_f32)\n",
    "\n",
    "plot_image_with_hist(naive_lahe_f32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Чтобы избежать таких артефактов, можно использовать движущееся окно, а именно рассчитывать гистограммы для каждого блока $(M, M)$ и менять значение только в его центре.\n",
    "\n",
    "Давайте реализуем алгоритм через наш обобщенный алгоритм фильтрации."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Фильтрация изображений\n",
    "\n",
    "![](https://courses.cv-gml.ru/storage/seminars/basic-image-processing/filt.png)\n",
    "\n",
    "Фильтрацией изображений называют преобразование пикселей изображений с учетом региона определенной формы.\n",
    "\n",
    "$$ g(x)_{ij} = k\\left(x_{\\mathrm{region}\\left(i,\\, j\\right)}\\right)$$\n",
    "\n",
    "Функция $k$ применяется к каждому региону изображения и называется ядром фильтра.\n",
    "\n",
    "Давайте реализуем обобщенный алгоритм для фильтрации изображений"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def filt(image, kernel, *, kernel_size=3, padd=0, fill_mode=\"reflect\"):\n",
    "    \"\"\"\n",
    "    kernel - arbitrary function with one argument - window\n",
    "    fill_mode - 'constant', 'symmetric', 'reflect'\n",
    "    \"\"\"\n",
    "    if padd:\n",
    "        padding = [(padd, padd), (padd, padd)]\n",
    "        if image.ndim == 3:\n",
    "            padding += [(0, 0)]\n",
    "\n",
    "        # use np.pad for image padding\n",
    "        image = np.pad(image, padding, mode=fill_mode)\n",
    "\n",
    "    h_range = range(image.shape[0] - (kernel_size - 1))\n",
    "    w_range = range(image.shape[1] - (kernel_size - 1))\n",
    "    channels = image.shape[-1:] if image.ndim == 3 else ()\n",
    "\n",
    "    # define new image\n",
    "    new_image = np.empty((len(h_range), len(w_range)) + channels, dtype=image.dtype)\n",
    "\n",
    "    for i in tqdm(h_range, total=len(h_range)):\n",
    "        for j in w_range:\n",
    "            # define image slice - window\n",
    "            window = image[i : i + kernel_size, j : j + kernel_size]\n",
    "\n",
    "            # apply kernel to window\n",
    "            value = kernel(window)\n",
    "\n",
    "            # write out the result\n",
    "            new_image[i, j] = value\n",
    "\n",
    "    return new_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Тесты\n",
    "from tests import test_filt\n",
    "\n",
    "test_filt(filt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](https://courses.cv-gml.ru/storage/seminars/basic-image-processing/ahe3.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lahe_kernel(window):\n",
    "    h, w = window.shape[:2]\n",
    "\n",
    "    # calculate cdf of image's color histogram\n",
    "    cdf = calculate_cdf(window)\n",
    "\n",
    "    # extract the middle pixel from the window\n",
    "    center = window[h // 2, w // 2]\n",
    "\n",
    "    # map only the extracted value using the whole cdf\n",
    "    value = np.interp(center, np.linspace(0, 1, 256), cdf)\n",
    "\n",
    "    return value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Тесты\n",
    "from tests import test_lahe\n",
    "\n",
    "test_lahe(filt, lahe_kernel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lahe_f32 = filt(gray_f32, lahe_kernel, kernel_size=51, padd=25)\n",
    "plot_image_with_hist(lahe_f32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lahe_darkened_with_outliers_f32 = filt(darkened_with_outliers, lahe_kernel, kernel_size=51, padd=25)\n",
    "plot_image_with_hist(lahe_darkened_with_outliers_f32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_compare(\n",
    "    [\n",
    "        (he_f32, he_darkened_with_outliers_f32),\n",
    "        (lahe_f32, lahe_darkened_with_outliers_f32),\n",
    "    ],\n",
    "    [\n",
    "        (\"he(original)\", \"he(darkened_with_outliers)\"),\n",
    "        (\"lahe(original)\", \"lahe(darkened_with_outliers)\"),\n",
    "    ],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Шум"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Цифровой шум - это дефект изображения, вносимый фотосенсорами и электроникой устройств(как черно-белые телевизоры портили нам яркость) из-за несовершенства технологий или фотонной природы света. Некоторые виды шумов можно промоделировать и использовать для контроля качества алгоритмов понижения уровня шума на картинке. Например, гауссовским шумом на изображении называется колебание яркости распределенное по нормальному закону."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drop-out noise\n",
    "\n",
    "Или по-другому потеря информации. Соль и перец - случайные черные и белые пиксели, а только белые - импульсивный."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_idx(h, w, count=500):\n",
    "    y = np.random.randint(h, size=count)\n",
    "    x = np.random.randint(w, size=count)\n",
    "    return y, x\n",
    "\n",
    "\n",
    "def dropout_values(image, val, count=500):\n",
    "    image = image.copy()\n",
    "    random_coords = sample_idx(*image.shape[:2], count=count)\n",
    "    image[random_coords] = val\n",
    "    return image\n",
    "\n",
    "\n",
    "def salt_pepper(image, count=5000):\n",
    "    # random white and black pixels\n",
    "    image = dropout_values(image, 1, count=count)\n",
    "    image = dropout_values(image, 0, count=count)\n",
    "    return image\n",
    "\n",
    "\n",
    "# only white pixels\n",
    "impulsive_f32 = dropout_values(img_f32, 1, count=5000)\n",
    "# plus dark\n",
    "salt_peper_f32 = salt_pepper(img_f32)\n",
    "\n",
    "plot_compare(\n",
    "    [(img_f32, impulsive_f32, salt_peper_f32)],\n",
    "    [(\"original\", \"impulsive noise\", \"salt and peper noise\")],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Median filter\n",
    "\n",
    "Медианный фильтр - один из нелинейных фильтров, используемых для уменьшения уровня шума в изображении. Результатом преобразования одного региона является медиана его чисел."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def median_filter(window):\n",
    "    channels = window.shape[-1:] if window.ndim == 3 else ()\n",
    "    values = window.copy().reshape((-1,) + channels)\n",
    "\n",
    "    # sort values\n",
    "    # values.sort(axis=0)\n",
    "\n",
    "    # choose mid value\n",
    "    # middle, odd = divmod(values.shape[0], 2)\n",
    "    # median = values[middle - (not odd) : middle + 1].mean(axis=0)\n",
    "\n",
    "    median = np.median(values, axis=0)\n",
    "\n",
    "    return median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Тесты\n",
    "from tests import test_median\n",
    "\n",
    "test_median(filt, median_filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images = []\n",
    "names = []\n",
    "for name, noise in [\n",
    "    (\"impulsive noise\", impulsive_f32),\n",
    "    (\"salt and peper noise\", salt_peper_f32),\n",
    "]:\n",
    "    images += [(noise, filt(noise, median_filter, kernel_size=3, padd=1))]\n",
    "    names += [(name, f\"{name} + median filter\")]\n",
    "\n",
    "plot_compare(images, names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noisy_f32 = (img_f32 + np.random.randn(*img_f32.shape) * 0.1).clip(0, 1)\n",
    "plt.imshow(noisy_f32)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_compare(\n",
    "    [(noisy_f32, filt(noisy_f32, median_filter, kernel_size=3, padd=1))],\n",
    "    [(\"noisy\", \"noisy+median filter\")],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Границы объектов\n",
    "\n",
    "Основной целью обнаружения резких изменений яркости изображения является фиксация важных событий и изменений мира. Они могут отражать различные предположения о модели формирования изображения, изменения в яркости изображения могут указывать на:\n",
    "\n",
    "> изменения глубины;\n",
    "\n",
    "> изменения ориентации поверхностей;\n",
    "\n",
    "> изменения в свойствах материала;\n",
    "\n",
    "> различие в освещении сцены.\n",
    "\n",
    "В идеальном случае результатом выделения границ является набор связанных кривых, обозначающих границы объектов. Таким образом, применение фильтра выделения границ к изображению может существенно уменьшить количество обрабатываемых данных из-за того, что отфильтрованная часть изображения считается менее значимой, а наиболее важные структурные свойства изображения сохраняются.\n",
    "\n",
    "Один из самых простых методов выделения границ - оператор Собеля.\n",
    "\n",
    "![](https://courses.cv-gml.ru/storage/seminars/basic-image-processing/sobel.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](https://courses.cv-gml.ru/storage/seminars/basic-image-processing/grad_norm.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sobel_kernel = np.array(\n",
    "    [\n",
    "        [1.0, 2.0, 1.0],\n",
    "        [0.0, 0.0, 0.0],\n",
    "        [-1.0, -2.0, -1.0],\n",
    "    ]\n",
    ")\n",
    "\n",
    "\n",
    "def sobel_x(window):\n",
    "    return (window * sobel_kernel).sum()\n",
    "\n",
    "\n",
    "def sobel_y(window):\n",
    "    return (window * sobel_kernel.T).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sobel(gray_image):\n",
    "    der_x = filt(gray_image, sobel_x, kernel_size=3, padd=1)\n",
    "    der_y = filt(gray_image, sobel_y, kernel_size=3, padd=1)\n",
    "\n",
    "    edge = np.sqrt(der_x**2 + der_y**2)\n",
    "    return edge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Тесты\n",
    "from tests import test_sobel\n",
    "\n",
    "test_sobel(sobel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "edge = sobel(gray_f32)\n",
    "\n",
    "plt.imshow(edge)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Max filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def max_filter(window):\n",
    "    # return max-value among height and width\n",
    "    return window.max(axis=(0, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Тесты\n",
    "from tests import test_max\n",
    "\n",
    "test_max(filt, max_filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_f32 = (img_f32 + 0.5 * np.random.randn(*img_f32.shape)).clip(0, 1)\n",
    "median_f32 = filt(noise_f32, median_filter, kernel_size=5, padd=2)\n",
    "max_f32 = filt(median_f32, max_filter, kernel_size=5, padd=2)\n",
    "\n",
    "res = (edge[:, :, None] * max_f32).clip(0, 1)\n",
    "\n",
    "plt.imshow(res)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 4
}