src-2/notebook_spat_filters.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assignment - 2\n",
    "### @Author: Honi, Selahaddin (honis@kaist.ac.kr)\n",
    "### @Date: 2022 April"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Program requires OpenCV for image I&O and matrix computations for NumPy. `Noise Addition` section contains Gaussian and S&P noise addition functions to apply these onto `lena` and `cameraman` images. Noisy images are stored in `A,B,C` lists. `Mean and Order-statistic Filters` functions are called inside a nested-loop for each neighborhood of each noisy images in `A,B,C` lists. Parameters (if exist) are determined in function arguments by default. `Adaptive Filters` have a bit different loop design such that `Local Noise Reduction Filter` calls `estimateNoiseVar()` function, firstly, to predict global image noise variance by calculating image histogram. `AdaptiveMedian()` filter function gets each neighborhood in the image and executes recursively. \n",
    "\n",
    "Neighborhood size `k=3` for all filters; however, adaptive median filter can change the neighborhood size if specified conditions are not satisfied. Therefore, maximum allowed neighborhood size is limited by `Smax=30`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Noise Addition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "def add_gaussian_noise(img,sigma):\n",
    "    gauss = np.random.normal(0,sigma,img.shape[:2])\n",
    "    img_noised = img + gauss\n",
    "    return img_noised\n",
    "\n",
    "def add_salt_and_pepper_noise(img,prob):\n",
    "    img_noised = img.copy()\n",
    "    probs = np.random.random(img.shape[:2])\n",
    "    img_noised[probs <     (prob / 2)] = 0\n",
    "    img_noised[probs > 1 - (prob / 2)] = 255\n",
    "    return img_noised   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = './img/'\n",
    "fnames = ['lena','cm']\n",
    "\n",
    "sigma = 20  # Gaussian noise\n",
    "prob = 0.25 # S&P noise\n",
    "\n",
    "O = [] # original images\n",
    "A = [] # S&P images\n",
    "B = [] # G images \n",
    "C = [] # S&P+G images\n",
    "\n",
    "for fname in fnames:\n",
    "    img = cv.imread(path+fname+'.png',0)\n",
    "\n",
    "    img_sp = add_salt_and_pepper_noise(img,prob=prob)\n",
    "    img_g = add_gaussian_noise(img,sigma=sigma)\n",
    "    img_spg = add_gaussian_noise(img_sp, sigma=sigma)\n",
    "\n",
    "    cv.imwrite(path+fname+'_sp'+'.png',img_sp)\n",
    "    cv.imwrite(path+fname+'_g'+'.png',img_g)\n",
    "    cv.imwrite(path+fname+'_spg'+'.png',img_spg)\n",
    "\n",
    "    O.append([img,fname])\n",
    "    A.append([img_sp,fname+'_sp'])\n",
    "    B.append([img_g,fname+'_g'])\n",
    "    C.append([img_spg,fname+'_spg'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of Mean and Order-statistic Filters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mean Filters \n",
    "# ================================================\n",
    "\n",
    "def arithmetic(neighborhood):\n",
    "    m,n = neighborhood.shape\n",
    "    val = 1/(m*n) * np.sum(neighborhood)\n",
    "    return val\n",
    "\n",
    "def geometric(neighborhood):\n",
    "    m,n = neighborhood.shape\n",
    "    prod = 1\n",
    "    for elem in neighborhood.flatten():\n",
    "        prod *= (elem + 1e-1)\n",
    "    val = prod ** (1/(m*n))\n",
    "    return val\n",
    "\n",
    "def harmonic(neighborhood):\n",
    "    m,n = neighborhood.shape\n",
    "    den = 1e-10 # prevent divide by zero \n",
    "    for elem in neighborhood.flatten():\n",
    "        den += 1/(elem+1e-2) # prevent divide by zero \n",
    "    val = (m*n) / den \n",
    "    return val\n",
    "\n",
    "def contraharmonic(neighborhood,Q=2):\n",
    "    nom = 0\n",
    "    den = 1e-10\n",
    "    for elem in neighborhood.flatten():\n",
    "        nom += elem**(Q+1)\n",
    "        den += elem**Q\n",
    "    val = nom/den\n",
    "    return val\n",
    "\n",
    "# Order-statistic Filters\n",
    "# ================================================\n",
    "\n",
    "def median(neighborhood):\n",
    "    m,n = neighborhood.shape\n",
    "    mid = m*n//2\n",
    "    sort = np.sort(neighborhood.flatten())\n",
    "    val = sort[mid]\n",
    "    return val\n",
    "\n",
    "def max(neighborhood):\n",
    "    val = neighborhood.max()\n",
    "    return val\n",
    "\n",
    "def min(neighborhood):\n",
    "    val = neighborhood.min()\n",
    "    return val\n",
    "\n",
    "def midpoint(neighborhood):\n",
    "    val = 0.5 * (neighborhood.max()+neighborhood.min())\n",
    "    return val\n",
    "\n",
    "def alpha_trimmed(neighborhood,d=2):\n",
    "    m,n = neighborhood.shape\n",
    "    sort = np.sort(neighborhood.flatten())\n",
    "    trim = sort[d//2:-d//2]\n",
    "    val = 1/(m*n-d) * np.sum(trim)\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3 # neigborhood size \n",
    "\n",
    "filters_mean = [arithmetic,geometric,harmonic,contraharmonic]\n",
    "filters_order_stat = [median, max, min, midpoint, alpha_trimmed]\n",
    "\n",
    "for imgs in [A,B,C]:\n",
    "    for img,fname in imgs:\n",
    "\n",
    "        out = img.copy()\n",
    "        H,W = img.shape\n",
    "\n",
    "        for applyFilter in filters_mean + filters_order_stat:\n",
    "\n",
    "            for y in range(H):\n",
    "                for x in range(W):\n",
    "\n",
    "                    # check neighborhood inside the image\n",
    "                    if y < k//2 or y >= H-(k//2) \\\n",
    "                        or x < k//2 or x >= W-(k//2):\n",
    "                        continue\n",
    "                    \n",
    "                    # y,x center point location\n",
    "                    # crop neighborhood area centered on (x,y) from original image\n",
    "                    neighborhood = img[y-(k//2):y+(k//2)+1, x-(k//2):x+(k//2)+1]\n",
    "                    # pass this neigborhood to specified filter\n",
    "                    out[y,x] = applyFilter(neighborhood)\n",
    "\n",
    "            out_path = path+fname+'_'+applyFilter.__name__+'_'+str(k)+'.png'\n",
    "            cv.imwrite(out_path,out)\n",
    "            print(fname+'_'+applyFilter.__name__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptive, Local Noise Reduction Filter "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimateNoiseVar(img,L=256):\n",
    "    # Use histogram to estimate noise variance\n",
    "    levels = np.arange(L)\n",
    "    p = []\n",
    "    m = 0\n",
    "\n",
    "    img1d = img.flatten()\n",
    "    pixel_total = len(img1d)\n",
    "\n",
    "    for r in levels:\n",
    "        idx = np.where(img1d==r)[0]\n",
    "        count = idx.shape[0]\n",
    "        p_of_r = count/pixel_total\n",
    "        p.append(p_of_r)\n",
    "        m += r*p_of_r\n",
    "\n",
    "    var_noise = 0\n",
    "    for r in range(L):\n",
    "        var_noise += ((levels[r]-m)**2) * p[r]\n",
    "    \n",
    "    return var_noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3 # neigborhood size \n",
    "\n",
    "for imgs in [A,B,C]:\n",
    "    for img,fname in imgs:\n",
    "\n",
    "        out = img.copy()\n",
    "        H,W = img.shape\n",
    "\n",
    "        # global image noise variance\n",
    "        var_noise = estimateNoiseVar(img)\n",
    "\n",
    "        for y in range(H):\n",
    "            for x in range(W):\n",
    "\n",
    "                # check neighborhood inside the image\n",
    "                if y < k//2 or y >= H-(k//2) \\\n",
    "                    or x < k//2 or x >= W-(k//2):\n",
    "                    continue\n",
    "                \n",
    "                # y,x center point location\n",
    "                # crop neighborhood area centered on (x,y) from original image\n",
    "                neighborhood = img[y-(k//2):y+(k//2)+1, x-(k//2):x+(k//2)+1]\n",
    "                \n",
    "                # local neighbothood variance \n",
    "                var_local = np.var(neighborhood.flatten())\n",
    "                # variance ratio check\n",
    "                ratio = 1 if var_noise>var_local else var_noise/var_local\n",
    "                # neighborhood intensity average\n",
    "                z_bar = np.mean(neighborhood.flatten())\n",
    "                # center intensity \n",
    "                g_of_xy = img[y,x]\n",
    "                \n",
    "                # adaptive restoration\n",
    "                out[y,x] = g_of_xy - ratio * (g_of_xy - z_bar)\n",
    "\n",
    "        out_path = path+fname+'_adaptive_reduction_'+str(k)+'.png'\n",
    "        cv.imwrite(out_path,out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptive, Median Filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AdaptiveMedian(z_ret,img,k,x,y,H,W,Smax=30):\n",
    "\n",
    "    # check neighborhood inside the image\n",
    "    if y < k//2 or y >= H-(k//2) \\\n",
    "        or x < k//2 or x >= W-(k//2):\n",
    "        return z_ret\n",
    "\n",
    "    # y,x center point location\n",
    "    # crop neighborhood area centered on (x,y) from original image\n",
    "    neighborhood = img[y-(k//2):y+(k//2)+1, x-(k//2):x+(k//2)+1]\n",
    "    \n",
    "    # min-max-median intensities\n",
    "    z_min = np.min(neighborhood)\n",
    "    z_max = np.max(neighborhood)\n",
    "    z_med = np.median(neighborhood)\n",
    "\n",
    "    # center intensity\n",
    "    z_xy = img[y,x]\n",
    "\n",
    "    # return intensity := z_med (by default)\n",
    "    z_ret = z_med\n",
    "\n",
    "    # Level-B\n",
    "    if z_min < z_med and z_med < z_max:\n",
    "        if z_min < z_xy and z_xy < z_max:\n",
    "            z_ret = z_xy\n",
    "        else:\n",
    "            z_ret = z_med\n",
    "    # Level-A\n",
    "    else:\n",
    "        k += 2 # increase neighborhood size \n",
    "        if k<=Smax:\n",
    "            z_ret = AdaptiveMedian(z_ret,img,k,x,y,H,W,Smax=Smax)\n",
    "        else:\n",
    "            z_ret = z_med\n",
    "\n",
    "    return z_ret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3 # neigborhood size \n",
    "\n",
    "for imgs in [A,B,C]:\n",
    "    for img,fname in imgs:\n",
    "\n",
    "        out = img.copy()\n",
    "        H,W = img.shape\n",
    "\n",
    "        for y in range(H):\n",
    "            for x in range(W):              \n",
    "\n",
    "                out[y,x] = AdaptiveMedian(1,img,k,x,y,H,W)\n",
    "\n",
    "        out_path = path+fname+'_adaptive_median_'+str(k)+'.png'\n",
    "        cv.imwrite(out_path,out)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "76ca30080b1287aa2231f5dfdce9b270326770c8d2601174cfcf63b351d9e30f"
  },
  "kernelspec": {
   "display_name": "Python 3.9.10 64-bit",
   "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.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}