{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 5.1 \n",
    "\n",
    "## Import Libraries\n",
    "Python requires importing libraries and functions you need to access specific tools like science (scipy), linear algebra (numpy), and graphics (matplotlib). These libraries can be installed using the ```pip``` command line tool. Alternatively you can install an python distribution like [Anaconda](https://www.continuum.io/downloads) or [Canopy](https://www.enthought.com/products/canopy/)  which have these and many other standard package pre-installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipywidgets as widgets  # add new widgets\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import os\n",
    "from IPython.display import display\n",
    "import numpy as np  # linear algebra / matrices\n",
    "from skimage.color import label2rgb\n",
    "from sklearn.metrics import roc_curve, auc  # roc curve tools\n",
    "from skimage.segmentation import mark_boundaries  # mark labels\n",
    "from skimage.io import imread  # read in images\n",
    "import matplotlib.pyplot as plt  # plotting\n",
    "%matplotlib inline\n",
    "# make the notebook interactive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_path = '../04-Segmentation/04-files'\n",
    "seg_path = os.path.join(base_path, 'DogVsMuffin_seg_bw.jpg')\n",
    "rgb_path = os.path.join(base_path, 'DogVsMuffin.jpg')\n",
    "face_path = os.path.join(base_path, 'DogVsMuffin_face.jpg')\n",
    "seg_img = imread(seg_path)[80:520:2, :450:2]\n",
    "rgb_img = imread(rgb_path)[80:520:2, :450:2, :]\n",
    "face_img = imread(face_path)\n",
    "print('RGB Size', rgb_img.shape, 'Seg Size',\n",
    "      seg_img.shape, 'Face Size', face_img.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 5))\n",
    "ax1.imshow(rgb_img)  # show the color image\n",
    "ax1.set_title(\"Color Image\")\n",
    "ax2.imshow(seg_img, cmap='gray')  # show the segments\n",
    "ax2.set_title(\"Ground Truth\")\n",
    "ax3.imshow(mark_boundaries(rgb_img, seg_img))\n",
    "ax3.set_title(\"Labeled Image\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a Simple ROC Curve\n",
    "We use the score function of taking the mean of the red green and blue channels\n",
    "$$ I = \\frac{R+G+B}{3} $$\n",
    "We then take the score by normalizing by the maximum value (since the image is 8bit this is 255)\n",
    "$$ s = \\frac{I}{255} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ground_truth_labels = seg_img.flatten() > 0\n",
    "score_value = 1-np.mean(rgb_img.astype(np.float32), 2).flatten()/255.0\n",
    "fpr, tpr, _ = roc_curve(ground_truth_labels, score_value)\n",
    "roc_auc = auc(fpr, tpr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(fpr, tpr, label='ROC curve (area = %0.2f)' % roc_auc)\n",
    "ax.plot([0, 1], [0, 1], 'k--')\n",
    "ax.set_xlim([0.0, 1.0])\n",
    "ax.set_ylim([0.0, 1.05])\n",
    "ax.set_xlabel('False Positive Rate')\n",
    "ax.set_ylabel('True Positive Rate')\n",
    "ax.set_title('Receiver operating characteristic example')\n",
    "ax.legend(loc=\"lower right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding Filters\n",
    "We can add a filter to this process by importing a ```uniform_filter``` and applying it before processing the image\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.ndimage.filters import uniform_filter\n",
    "%matplotlib inline\n",
    "filter_size = 45\n",
    "filtered_image = uniform_filter(np.mean(rgb_img, 2), filter_size)\n",
    "score_value = 1-filtered_image.astype(np.float32).flatten()/255.0\n",
    "fpr2, tpr2, _ = roc_curve(ground_truth_labels, score_value)\n",
    "roc_auc2 = auc(fpr2, tpr2)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(fpr, tpr, label='Raw ROC curve (area = %0.2f)' % roc_auc)\n",
    "ax.plot(fpr2, tpr2, label='Filtered ROC curve (area = %0.2f)' % roc_auc2)\n",
    "ax.plot([0, 1], [0, 1], 'k--')\n",
    "ax.set_xlim([0.0, 1.0])\n",
    "ax.set_ylim([0.0, 1.05])\n",
    "ax.set_xlabel('False Positive Rate')\n",
    "ax.set_ylabel('True Positive Rate')\n",
    "ax.set_title('Receiver operating characteristic example')\n",
    "ax.legend(loc=\"lower right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tasks\n",
    "1. How can you improve filtering in this analysis? \n",
    " - Which filter elements might improve the area under the ROC?\n",
    " - Try making workflows to test out a few different filters\n",
    "2. Where might morphological operations fit in?\n",
    " - How can you make them part of this workflow as well?\n",
    "3. (Challenge) Try and use the optimize toolbox of _scipy_ with the fmin function (```from scipy.optimize import fmin```) to find the optimum parmeters for the highers area (hint: fmin finds the minimum value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import fmin\n",
    "\n",
    "\n",
    "def calc_auc(rv, gv, bv, fsize):\n",
    "    filter_size = 45\n",
    "    gray_image = (rv*rgb_img[:, :, 0]+gv*rgb_img[:, :,\n",
    "                                                 1]+bv*rgb_img[:, :, 2])/(rv+gv+bv)\n",
    "    filtered_image = uniform_filter(gray_image, filter_size)\n",
    "    score_value = filtered_image.astype(np.float32).flatten()/255.0\n",
    "    fpr2, tpr2, _ = roc_curve(ground_truth_labels, score_value)\n",
    "    return {'fpr': fpr2, 'tpr': tpr2, 'auc': auc(fpr2, tpr2), 'gimg': gray_image, 'fimg': filtered_image}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test the function to make sure it works\n",
    "def min_func(args): return 1-calc_auc(*args)['auc']\n",
    "\n",
    "\n",
    "min_start = [1, 1, 1, 20]\n",
    "min_func(min_start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_res = fmin(min_func, min_start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_values = calc_auc(*opt_res)\n",
    "tprOpt = opt_values['tpr']\n",
    "fprOpt = opt_values['fpr']\n",
    "roc_aucOpt = opt_values['auc']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax_img, ax) = plt.subplots(1, 2, figsize=(20, 10))\n",
    "ax_img.imshow(opt_values['gimg'], cmap='gray')\n",
    "ax_img.set_title('Transformed Color Image')\n",
    "ax.plot(fpr, tpr, label='Raw ROC curve (area = %0.2f)' % roc_auc)\n",
    "ax.plot(fprOpt, tprOpt, label='Optimized ROC curve (area = %0.2f)' % roc_aucOpt)\n",
    "ax.plot([0, 1], [0, 1], 'k--')\n",
    "ax.set_xlim([0.0, 1.0])\n",
    "ax.set_ylim([0.0, 1.05])\n",
    "ax.set_xlabel('False Positive Rate')\n",
    "ax.set_ylabel('True Positive Rate')\n",
    "ax.set_title('Receiver operating characteristic example')\n",
    "ax.legend(loc=\"lower right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-linear optimization\n",
    "Here we use non-linear approaches to improve the quality of the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def relu(x):\n",
    "    return (x+np.abs(x))/2\n",
    "\n",
    "\n",
    "def calc_auc_nl(rv, rm, gv, gm, bv, bm):\n",
    "    filter_size = 45\n",
    "    gray_image = (rv*relu(rgb_img[:, :, 0]/255.0-rm)+gv*relu(rgb_img[:, :, 1]/255.0-gm) +\n",
    "                  bv*relu(rgb_img[:, :, 2]/255.0-bm))/(rv+gv+bv)\n",
    "    score_value = gray_image.astype(np.float32).flatten()\n",
    "    fpr2, tpr2, _ = roc_curve(ground_truth_labels, score_value)\n",
    "    return {'fpr': fpr2, 'tpr': tpr2, 'auc': auc(fpr2, tpr2), 'gimg': gray_image, 'fimg': filtered_image}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test the function to make sure it works\n",
    "def min_func(args): return 1-calc_auc_nl(*args)['auc']\n",
    "\n",
    "\n",
    "min_start = [1, 0, 1, 0, 1, 0]\n",
    "min_start[0] = opt_res[0]\n",
    "min_start[2] = opt_res[1]\n",
    "min_start[4] = opt_res[2]\n",
    "min_func(min_start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_res = fmin(min_func, min_start, maxiter=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_values_nl = calc_auc_nl(*opt_res)\n",
    "tprOpt_nl = opt_values_nl['tpr']\n",
    "fprOpt_nl = opt_values_nl['fpr']\n",
    "roc_aucOpt_nl = opt_values_nl['auc']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax_img, ax) = plt.subplots(1, 2, figsize=(20, 10))\n",
    "ax_img.imshow(opt_values_nl['gimg'], cmap='gray')\n",
    "ax_img.set_title('Transformed Color Image')\n",
    "ax.plot(fpr, tpr, label='Raw ROC curve (area = %0.2f)' % roc_auc)\n",
    "ax.plot(fprOpt, tprOpt, label='Optimized ROC curve (area = %0.2f)' % roc_aucOpt)\n",
    "ax.plot(fprOpt_nl, tprOpt_nl,\n",
    "        label='NL Optimized ROC curve (area = %0.2f)' % roc_aucOpt_nl)\n",
    "ax.plot([0, 1], [0, 1], 'k--')\n",
    "ax.set_xlim([0.0, 1.0])\n",
    "ax.set_ylim([0.0, 1.05])\n",
    "ax.set_xlabel('False Positive Rate')\n",
    "ax.set_ylabel('True Positive Rate')\n",
    "ax.set_title('Receiver operating characteristic example')\n",
    "ax.legend(loc=\"lower right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Next Steps\n",
    "Rather than simply adjusting basic parameters, we can adjust entire arrays of information. The example below is the a convolutional neural network with one two layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import fftconvolve\n",
    "\n",
    "def convolve(img1, img2): return fftconvolve(img1, img2, mode='same')\n",
    "\n",
    "\n",
    "CONV_SIZE = (10, 10, 1)\n",
    "grey_img = np.reshape(np.mean(rgb_img, 2)/255.0,\n",
    "                      (rgb_img.shape[0], rgb_img.shape[1], 1))\n",
    "\n",
    "def calc_auc_conv(rcoefs):\n",
    "    coefs = rcoefs.reshape(CONV_SIZE)/rcoefs.sum()\n",
    "    score_image = relu(convolve(grey_img, coefs))\n",
    "    score_value = score_image.flatten()\n",
    "    fpr2, tpr2, _ = roc_curve(ground_truth_labels, score_value)\n",
    "    return {'fpr': fpr2, 'tpr': tpr2, 'auc': auc(fpr2, tpr2), 'gimg': score_image}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a nice gaussian kernel\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(2019)\n",
    "from functools import reduce\n",
    "def gkern_nd(d=2, kernlen=21, nsigs=3, min_smooth_val=1e-2):\n",
    "    nsigs = [nsigs] * d\n",
    "    k_wid = (kernlen - 1) / 2\n",
    "    all_axs = [np.linspace(-k_wid, k_wid, kernlen)] * d\n",
    "    all_xxs = np.meshgrid(*all_axs)\n",
    "    all_dist = reduce(np.add, [\n",
    "        np.square(cur_xx) / (2 * np.square(np.clip(nsig, min_smooth_val,\n",
    "                                                   kernlen)))\n",
    "        for cur_xx, nsig in zip(all_xxs, nsigs)])\n",
    "    kernel_raw = np.exp(-all_dist)\n",
    "    return kernel_raw / kernel_raw.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test the function to make sure it works\n",
    "def min_func(rcoefs): return 1-calc_auc_conv(rcoefs)['auc']\n",
    "\n",
    "min_start = gkern_nd(2, CONV_SIZE[0]).ravel()\n",
    "min_func(min_start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_res_conv = min_start\n",
    "opt_res_conv = fmin(min_func,\n",
    "                    opt_res_conv,\n",
    "                    maxiter=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_values_conv = calc_auc_conv(opt_res_conv)\n",
    "tprOpt_conv = opt_values_conv['tpr']\n",
    "fprOpt_conv = opt_values_conv['fpr']\n",
    "roc_aucOpt_conv = opt_values_conv['auc']\n",
    "out_kernel = opt_res_conv.reshape(CONV_SIZE)/opt_res_conv.sum()\n",
    "fig, ax_all = plt.subplots(1, out_kernel.shape[2])\n",
    "for i, c_ax in enumerate(np.array(ax_all).flatten()):\n",
    "    c_ax.imshow(out_kernel[:, :, i])\n",
    "    c_ax.set_title(str(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax_img, ax) = plt.subplots(1, 2, figsize=(20, 10))\n",
    "ax_img.imshow(opt_values_conv['gimg'].squeeze(), cmap='gray')\n",
    "ax_img.set_title('Transformed Color Image')\n",
    "ax.plot(fpr, tpr, label='Raw ROC curve (area = %0.2f)' % roc_auc)\n",
    "ax.plot(fprOpt, tprOpt, label='Optimized ROC curve (area = %0.2f)' % roc_aucOpt)\n",
    "ax.plot(fprOpt_conv, tprOpt_conv,\n",
    "        label='CNN Optimized ROC curve (area = %0.2f)' % roc_aucOpt_conv)\n",
    "ax.plot([0, 1], [0, 1], 'k--')\n",
    "ax.set_xlim([0.0, 1.0])\n",
    "ax.set_ylim([0.0, 1.05])\n",
    "ax.set_xlabel('False Positive Rate')\n",
    "ax.set_ylabel('True Positive Rate')\n",
    "ax.set_title('Receiver operating characteristic example')\n",
    "ax.legend(loc=\"lower right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RGB CNN\n",
    "Using the RGB instead of the gray value for the CNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CONV_SIZE = (10, 10, 3)\n",
    "\n",
    "def calc_auc_conv2d(rcoefs):\n",
    "    coefs = rcoefs.reshape(CONV_SIZE)/rcoefs.sum()\n",
    "    score_image = relu(convolve(grey_img, coefs))\n",
    "    score_value = score_image.flatten()\n",
    "    fpr2, tpr2, _ = roc_curve(ground_truth_labels, score_value)\n",
    "    return {'fpr': fpr2, 'tpr': tpr2, 'auc': auc(fpr2, tpr2), 'gimg': score_image}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_func(rcoefs): return 1-calc_auc_conv2d(rcoefs)['auc']\n",
    "min_kernel = np.stack([gkern_nd(2, kernlen=CONV_SIZE[0])]*3, -1)\n",
    "min_start = min_kernel.ravel()\n",
    "for i in range(10):\n",
    "    min_func(min_start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_res_conv2d = fmin(min_func, min_start, maxfun=50, maxiter=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_values_conv = calc_auc_conv2d(opt_res_conv2d)\n",
    "tprOpt_conv = opt_values_conv['tpr']\n",
    "fprOpt_conv = opt_values_conv['fpr']\n",
    "roc_aucOpt_conv = opt_values_conv['auc']\n",
    "out_kernel = opt_res_conv2d.reshape(CONV_SIZE)/opt_res_conv.sum()\n",
    "fig, ax_all = plt.subplots(3, out_kernel.shape[2], figsize=(10, 10))\n",
    "for i, (c_ax, d_ax, cd_ax) in enumerate(ax_all.T):\n",
    "    c_ax.imshow(min_kernel[:, :, i])\n",
    "    c_ax.set_title('Initial {}'.format(i))\n",
    "    c_ax.axis('off')\n",
    "    d_ax.imshow(out_kernel[:, :, i])\n",
    "    d_ax.set_title('Updated {}'.format(i))\n",
    "    d_ax.axis('off')\n",
    "    cd_ax.imshow(out_kernel[:, :, i]-min_kernel[:, :, i],\n",
    "               cmap='RdBu', vmin=-1e-3, vmax=1e-3)\n",
    "    cd_ax.set_title('Difference {}'.format(i))\n",
    "    cd_ax.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax_img, ax) = plt.subplots(1, 2, figsize=(20, 10))\n",
    "ax_img.imshow(mark_boundaries(\n",
    "    opt_values_conv['gimg'].squeeze(), seg_img), cmap='gray')\n",
    "ax_img.set_title('Transformed Color Image')\n",
    "ax.plot(fpr, tpr, label='Raw ROC curve (area = %0.2f)' % roc_auc)\n",
    "ax.plot(fprOpt, tprOpt, label='Optimized ROC curve (area = %0.2f)' % roc_aucOpt)\n",
    "ax.plot(fprOpt_conv, tprOpt_conv,\n",
    "        label='RGBCNN Optimized ROC curve (area = %0.2f)' % roc_aucOpt_conv)\n",
    "ax.plot([0, 1], [0, 1], 'k--')\n",
    "ax.set_xlim([0.0, 1.0])\n",
    "ax.set_ylim([0.0, 1.05])\n",
    "ax.set_xlabel('False Positive Rate')\n",
    "ax.set_ylabel('True Positive Rate')\n",
    "ax.set_title('Receiver operating characteristic example')\n",
    "ax.legend(loc=\"lower right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Tasks\n",
    "1. How can you improve filtering in this analysis? \n",
    " - Which filter elements might improve the area under the ROC?\n",
    " - Try making workflows to test out a few different filters\n",
    "2. Where might morphological operations fit in?\n",
    " - How can you make them part of this workflow as well?\n",
    "3. (Challenge) How would you add multiple filter operations? Can you optimize all of the parameters? What problems do you run into as you make your model more complex?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}