{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dab9c79b",
   "metadata": {},
   "source": [
    "Illustrate separable 2D low-pass filtering of an image\n",
    "using the convolution matrix generated in HW1.\n",
    "- 2017-09-11, Jeff Fessler, University of Michigan\n",
    "- 2021-08-20 Julia 1.6.2\n",
    "- 2023-07-24 Julia 1.9.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23b0a6a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "using MIRTjim: jim # image display\n",
    "using ImagePhantoms: shepp_logan, SheppLoganEmis # test image\n",
    "using Random: seed!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fd9fe64",
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace with the path to your own function, or paste your function code here\n",
    "mydir = ENV[\"hw551test\"] # change path\n",
    "include(joinpath(mydir, \"convolution.jl\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58e61548",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a 2D image binary test image of a simple disk and display it\n",
    "Xtrue = shepp_logan(128, SheppLoganEmis())[15:114,:]\n",
    "@show (M,N) = size(Xtrue)\n",
    "seed!(0)\n",
    "X = Xtrue + 0.5 * randn(size(Xtrue)) # add some noise\n",
    "jim(X, xlabel=\"x\", ylabel=\"y\", \"original noisy image\") # jiffy image display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20784bb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1D moving-average filter represented as a matrix\n",
    "hx = ones(9)/9 # 9-pt moving average filter for x (along cols)\n",
    "hy = ones(5)/5 # 5-pt moving average filter for y (along rows)\n",
    "Hx,_ = convolution(hx, zeros(M))\n",
    "Hy,_ = convolution(hy, zeros(N));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8241c2b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# display one of the 1D convolution matrices\n",
    "# transpose to show it like a matrix H[i,j], rather than as an image f[x,y]\n",
    "jim(Hx', ylabel=\"i\", xlabel=\"j\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd15ad16",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = Hx * X * Hy'; # apply 1D filters via convolution to both cols and rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a3666d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# If the following line fails for you,\n",
    "# or if the result is unrelated to your name on Canvas,\n",
    "# then edit it to be `user =  \"(your name)\"`\n",
    "user = Sys.iswindows() ? ENV[\"USERNAME\"] : ENV[\"USER\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80e4428f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display filtered image - the noise is greatly reduced!\n",
    "# It is \"smeared out\" (filtered or smoothed) more along x than y.\n",
    "# Do you know why?\n",
    "jim(Y, xlabel=\"x\", ylabel=\"y\", \"filtered image by '$user'\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dca3e14",
   "metadata": {},
   "outputs": [],
   "source": [
    "#using Plots: savefig\n",
    "#savefig(\"hs003jf-fig.pdf\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.9.2",
   "language": "julia",
   "name": "julia-1.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}