Plane analysis filter

analysis/plane_filter.ipynb

@@ -6,18 +6,23 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# image processing imports\n",
     "from ipfml.processing.segmentation import divide_in_blocks\n",
     "from ipfml.processing import transform\n",
     "from ipfml import utils\n",
     "from PIL import Image\n",
     "from scipy import signal\n",
     "from skimage import color\n",
+    "import cv2\n",
     "import scipy.stats as stats\n",
+    "\n",
+    "# display imports\n",
     "import seaborn as sns\n",
-    "import cv2\n",
-    "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
-    "from numpy.linalg import svd\n",
+    "from mpl_toolkits.mplot3d import Axes3D\n",
+    "\n",
+    "# main imports\n",
+    "import numpy as np\n",
     "import os"
    ]
   },
@@ -83,59 +88,6 @@
     "    return zones_img"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def display_svd_reconstruction(interval, zones):\n",
-    "    \n",
-    "    output_images = []\n",
-    "    begin, end = interval\n",
-    "    for zone in zones:\n",
-    "        lab_img = metrics.get_LAB_L(zone)\n",
-    "        lab_img = np.array(lab_img, 'uint8')\n",
-    "        \n",
-    "        U, s, V = svd(lab_img, full_matrices=True)\n",
-    "        \n",
-    "        smat = np.zeros((end-begin, end-begin), dtype=complex)\n",
-    "        smat[:, :] = np.diag(s[begin:end])\n",
-    "        output_img = np.dot(U[:, begin:end],  np.dot(smat, V[begin:end, :]))\n",
-    "        \n",
-    "        print(output_img)\n",
-    "        print(np.allclose(lab_img, output_img))\n",
-    "        \n",
-    "        output_img = np.array(output_img, 'uint8')\n",
-    "        output_images.append(Image.fromarray(output_img))\n",
-    "        \n",
-    "    return output_images"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 6,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def display_images(dict_data, rec_images):\n",
-    "    \n",
-    "    indices = dict_data['indices']\n",
-    "    scene = dict_data['name']\n",
-    "    \n",
-    "    fig=plt.figure(figsize=(15, 8))\n",
-    "    columns = len(zones)\n",
-    "    rows = 1\n",
-    "    for i in range(1, columns*rows +1):\n",
-    "        index = i - 1\n",
-    "        fig.add_subplot(rows, columns, i)\n",
-    "        plt.imshow(rec_images[index], label=scene + '_' + str(indices[index]))\n",
-    "        img_path = 'tmp_images/' + dict_data['prefix'] + 'zone'+ str(current_dict['zone']) + '_reconstruct_' + str(indices[index]) + '.png'\n",
-    "        Image.fromarray(np.asarray(rec_images[index], 'uint8')).save(img_path)\n",
-    "    plt.show()\n",
-    "    "
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -145,7 +97,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -181,30 +133,9 @@
     "### Definition of parameters"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Here we define parameters for the rest of this study :\n",
-    "- the scene used\n",
-    "- the reconstructed interval (give reduced information from SVD decomposition) \n",
-    "- the displayed interval of SVD values"
-   ]
-  },
   {
    "cell_type": "code",
-   "execution_count": 8,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "current_dict = dict_appart\n",
-    "displayed_interval = (50, 200)\n",
-    "reconstructed_interval = (90, 200)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [
     {
@@ -242,7 +173,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -252,69 +183,111 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [],
    "source": [
-    "'''/Fonctionne : https://oomake.com/question/264689                                                                                                                  \n",
-    "import matplotlib.pyplot as plt\n",
-    "from mpl_toolkits.mplot3d import Axes3D\n",
-    "import numpy as np\n",
-    "N_POINTS = 10\n",
-    "TARGET_X_SLOPE = 2\n",
-    "TARGET_y_SLOPE = 3\n",
-    "TARGET_OFFSET  = 5\n",
-    "EXTENTS = 5\n",
-    "NOISE = 5\n",
-    "# create random data\n",
-    "xs = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]\n",
-    "ys = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]\n",
-    "zs = []\n",
-    "for i in range(N_POINTS):\n",
-    "    zs.append(xs[i]*TARGET_X_SLOPE + \\\n",
-    "              ys[i]*TARGET_y_SLOPE + \\\n",
-    "              TARGET_OFFSET + np.random.normal(scale=NOISE))\n",
-    "# plot raw data\n",
-    "plt.figure()\n",
-    "ax = plt.subplot(111, projection='3d')\n",
-    "ax.scatter(xs, ys, zs, color='b')\n",
-    "# do fit\n",
-    "tmp_A = []\n",
-    "tmp_b = []\n",
-    "for i in range(len(xs)):\n",
-    "    tmp_A.append([xs[i], ys[i], 1])\n",
-    "    tmp_b.append(zs[i])\n",
-    "b = np.matrix(tmp_b).T\n",
-    "A = np.matrix(tmp_A)\n",
-    "fit = (A.T * A).I * A.T * b\n",
-    "errors = b - A * fit\n",
-    "residual = np.linalg.norm(errors)\n",
-    "print \"solution:\"\n",
-    "print \"%f x + %f y + %f = z\" % (fit[0], fit[1], fit[2])\n",
-    "print \"errors:\"\n",
-    "print errors\n",
-    "print \"residual:\"\n",
-    "print residual\n",
-    "# plot plane\n",
-    "xlim = ax.get_xlim()\n",
-    "ylim = ax.get_ylim()\n",
-    "X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),\n",
-    "                  np.arange(ylim[0], ylim[1]))\n",
-    "Z = np.zeros(X.shape)\n",
-    "for r in range(X.shape[0]):\n",
-    "    for c in range(X.shape[1]):\n",
-    "        Z[r,c] = fit[0] * X[r,c] + fit[1] * Y[r,c] + fit[2]\n",
-    "ax.plot_wireframe(X,Y,Z, color='k')\n",
-    "ax.set_xlabel('x')\n",
-    "ax.set_ylabel('y')\n",
-    "ax.set_zlabel('z')\n",
-    "plt.show()\n",
-    "'''"
+    "# return residual information\n",
+    "def plane_kernel_filter(window):\n",
+    "    \n",
+    "    width, height = window.shape\n",
+    "\n",
+    "    # prepare data\n",
+    "    nb_elem = width * height\n",
+    "    xs = [int(i/height) for i in range(nb_elem)]\n",
+    "    ys = [i%height for i in range(nb_elem)]\n",
+    "    zs = np.array(window).flatten().tolist()\n",
+    "\n",
+    "    # get residual (error) from mean plane computed\n",
+    "    tmp_A = []\n",
+    "    tmp_b = []\n",
+    "\n",
+    "    for i in range(len(xs)):\n",
+    "        tmp_A.append([xs[i], ys[i], 1])\n",
+    "        tmp_b.append(zs[i])\n",
+    "\n",
+    "    b = np.matrix(tmp_b).T\n",
+    "    A = np.matrix(tmp_A)\n",
+    "\n",
+    "    fit = (A.T * A).I * A.T * b\n",
+    "\n",
+    "    errors = b - A * fit\n",
+    "    residual = np.linalg.norm(errors)\n",
+    "\n",
+    "    return residual\n",
+    "\n",
+    "\n",
+    "# return difference between min and max errors\n",
+    "def plane_kernel_filter_max_error(window):\n",
+    "    \n",
+    "    width, height = window.shape\n",
+    "\n",
+    "    # prepare data\n",
+    "    nb_elem = width * height\n",
+    "    xs = [int(i/height) for i in range(nb_elem)]\n",
+    "    ys = [i%height for i in range(nb_elem)]\n",
+    "    zs = np.array(window).flatten().tolist()\n",
+    "\n",
+    "    # get residual (error) from mean plane computed\n",
+    "    tmp_A = []\n",
+    "    tmp_b = []\n",
+    "\n",
+    "    for i in range(len(xs)):\n",
+    "        tmp_A.append([xs[i], ys[i], 1])\n",
+    "        tmp_b.append(zs[i])\n",
+    "\n",
+    "    b = np.matrix(tmp_b).T\n",
+    "    A = np.matrix(tmp_A)\n",
+    "\n",
+    "    fit = (A.T * A).I * A.T * b\n",
+    "\n",
+    "    errors = b - A * fit\n",
+    "    residual = np.linalg.norm(errors)\n",
+    "    \n",
+    "    errors = abs(np.array(errors))\n",
+    "\n",
+    "    return (errors.max() - errors.min())\n",
+    "\n",
+    "\n",
+    "def plane_custom_filter(img, kernel=(5, 5)):\n",
+    "    \n",
+    "    img = np.array(img)\n",
+    "    \n",
+    "    width, height = img.shape\n",
+    "    \n",
+    "    kernel_width, kernel_height = kernel\n",
+    "    \n",
+    "    if kernel_width % 2 == 0 or kernel_height % 2 == 0:\n",
+    "        raise ValueError(\"Invalid kernel size, need to be of odd size\")\n",
+    "        \n",
+    "    padding_height = (kernel_width - 1) / 2\n",
+    "    padding_width = (kernel_width - 1) / 2\n",
+    "    \n",
+    "    img_plane_error = []\n",
+    "    for i in range(width):\n",
+    "        \n",
+    "        if i >= padding_width and i < (width - padding_width):\n",
+    "            \n",
+    "            row_plane_error = []\n",
+    "            \n",
+    "            for j in range (height):\n",
+    "                \n",
+    "                if j >= padding_height and j < (height - padding_height):\n",
+    "                    \n",
+    "                    # pixel in the center of kernel window size, need to extract window from img\n",
+    "                    window = img[int(i-padding_width):int(i+padding_width + 1), int(j-padding_height):int(j+padding_height + 1)]\n",
+    "                    \n",
+    "                    diff = plane_kernel_filter(window)\n",
+    "                    row_plane_error.append(diff)\n",
+    "                    \n",
+    "            img_plane_error.append(row_plane_error)\n",
+    "        \n",
+    "    return np.array(img_plane_error)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 26,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -322,31 +295,45 @@
     "    \n",
     "    sub_zones = divide_in_blocks(zone, (20, 20))\n",
     "\n",
-    "    diff_list = []\n",
+    "    error_list = []\n",
     "\n",
     "    for sub_zone in sub_zones:\n",
     "        l_img = transform.get_LAB_L(sub_zone)\n",
-    "        diff = diff_custom_filter(utils.normalize_2D_arr(l_img), (5, 5), max)\n",
-    "        global_diff = np.std(diff)\n",
-    "        diff_list.append(global_diff)\n",
+    "        plane_error = plane_custom_filter(utils.normalize_2D_arr(l_img), (5, 5))\n",
+    "        global_diff = np.std(plane_error)\n",
+    "        error_list.append(plane_error)\n",
     "\n",
-    "    diff_list = np.array(diff_list)\n",
-    "    score = np.std(diff_list[0:int(len(sub_zones)/5)])\n",
-    "    print(score)"
+    "    error_list = np.array(error_list)\n",
+    "    score = np.std(error_list[0:int(len(sub_zones)/5)])\n",
+    "    print(score)\n",
+    "    \n",
+    "def apply_on_zone_plane_normed(zone, kernel=(5, 5)):\n",
+    "    \n",
+    "    l_img = transform.get_LAB_L(zone)\n",
+    "    \n",
+    "    plane_error = plane_custom_filter(utils.normalize_2D_arr(l_img), kernel)\n",
+    "    return np.mean(plane_error)\n",
+    "    \n",
+    "def apply_on_zone_plane(zone, kernel=(5, 5)):\n",
+    "    \n",
+    "    l_img = transform.get_LAB_L(zone)\n",
+    "    \n",
+    "    plane_error = plane_custom_filter(l_img, kernel)\n",
+    "    return np.mean(plane_error)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "<matplotlib.image.AxesImage at 0x7f63b9a36e48>"
+       "<matplotlib.image.AxesImage at 0x7fc3b51bf320>"
       ]
      },
-     "execution_count": 13,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -369,16 +356,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "<matplotlib.image.AxesImage at 0x7f63b61ccc88>"
+       "<matplotlib.image.AxesImage at 0x7fc3b194eeb8>"
       ]
      },
-     "execution_count": 14,
+     "execution_count": 11,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -401,105 +388,99 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 18,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "0.019411785375902203\n",
-      "0.022010065900329633\n",
-      "0.023538288175525304\n",
-      "0.024255978068903696\n",
-      "0.024555364941822696\n",
-      "0.025371094905613262\n",
-      "0.0257215863991888\n",
-      "0.025645858204273037\n",
-      "0.02586192220176574\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
-    "for zone in zones_appart:\n",
-    "    apply_on_zone(zone)"
+    "def display_computed_data(zones, dict_scene):\n",
+    "    \n",
+    "    errors = []\n",
+    "    errors_normed = []\n",
+    "    \n",
+    "    print(\"---------------------------------------------------------------------------------------\")\n",
+    "    print(\"Compute error on \" + dict_scene[\"name\"] + \" scene (zone \" + str(dict_scene[\"zone\"]) + \")\") \n",
+    "    \n",
+    "    for index, zone in enumerate(zones):\n",
+    "        plane_mean_error = apply_on_zone_plane(zone)\n",
+    "        plane_mean_error_normed = apply_on_zone_plane_normed(zone)\n",
+    "    \n",
+    "        errors.append(plane_mean_error)\n",
+    "        errors_normed.append(plane_mean_error_normed)\n",
+    "        \n",
+    "        print(dict_scene[\"prefix\"] +dict_scene[\"indices\"][index], \"=> score\",\"{0:.8f}\".format(plane_mean_error),\"| normed :\",\"{0:.8f}\".format(plane_mean_error_normed))\n",
+    "        \n",
+    "    return errors, errors_normed"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": 24,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "0.016907859437032855\n",
-      "0.022510883215429146\n",
-      "0.02774428742446508\n",
-      "0.030267113383204484\n",
-      "0.032330051361439856\n",
-      "0.032666832076238675\n",
-      "0.03325133529876735\n",
-      "0.033614614352306754\n",
-      "0.03435194845496702\n",
-      "0.035516845913060015\n",
-      "0.035844373588709205\n",
-      "0.03626505901706533\n"
+      "---------------------------------------------------------------------------------------\n",
+      "Compute error on Appart1opt02 scene (zone 9)\n",
+      "appartAopt_00050 => score 13.30254427 | normed : 0.15702255\n",
+      "appartAopt_00100 => score 11.09155512 | normed : 0.13641112\n",
+      "appartAopt_00200 => score 9.53841914 | normed : 0.11535671\n",
+      "appartAopt_00300 => score 8.93721088 | normed : 0.10914040\n",
+      "appartAopt_00400 => score 8.57067712 | normed : 0.10599665\n",
+      "appartAopt_00600 => score 8.16948645 | normed : 0.10148266\n",
+      "appartAopt_00700 => score 8.04385333 | normed : 0.09958658\n",
+      "appartAopt_00800 => score 7.96053962 | normed : 0.09888709\n",
+      "appartAopt_00900 => score 7.89593901 | normed : 0.09808461\n",
+      "---------------------------------------------------------------------------------------\n",
+      "Compute error on Cuisine01 scene (zone 6)\n",
+      "cuisine01_00050 => score 14.26744032 | normed : 0.15361529\n",
+      "cuisine01_00100 => score 11.91046744 | normed : 0.12842483\n",
+      "cuisine01_00200 => score 10.20911653 | normed : 0.11054688\n",
+      "cuisine01_00300 => score 9.45622725 | normed : 0.10291345\n",
+      "cuisine01_00400 => score 9.01093786 | normed : 0.09767920\n",
+      "cuisine01_00600 => score 8.50241410 | normed : 0.09253297\n",
+      "cuisine01_00700 => score 8.33067727 | normed : 0.09066393\n",
+      "cuisine01_00800 => score 8.21224043 | normed : 0.08937497\n",
+      "cuisine01_00900 => score 8.10144871 | normed : 0.08782028\n",
+      "cuisine01_01000 => score 8.02322747 | normed : 0.08731791\n",
+      "cuisine01_01100 => score 7.94861294 | normed : 0.08650587\n",
+      "cuisine01_01200 => score 7.89148254 | normed : 0.08592211\n"
      ]
     }
    ],
    "source": [
-    "for zone in zones_cuisine:\n",
-    "    apply_on_zone(zone)"
+    "appart_errors = display_computed_data(zones_appart, dict_appart)\n",
+    "cuisine_errors = display_computed_data(zones_cuisine, dict_cuisine)"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": 160,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "1.55 s ± 15.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
-     ]
-    },
-    {
-     "data": {
-      "text/plain": [
-       "(194, 194)"
-      ]
-     },
-     "execution_count": 160,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
    "source": [
-    "%timeit diff_img = diff_custom_filter(l_img, kernel=(3, 3), func=max)\n",
-    "diff_img.shape"
+    "## Performance indication"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 123,
+   "execution_count": 23,
    "metadata": {},
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "1.6049964585730372"
-      ]
-     },
-     "execution_count": 123,
-     "metadata": {},
-     "output_type": "execute_result"
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "8.48 s ± 267 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
+      "8.66 s ± 293 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
+      "9.56 s ± 679 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
+     ]
     }
    ],
    "source": [
-    "np.mean(diff_img)"
+    "first_zone = zones_appart[0]\n",
+    "%timeit plane_error_img = apply_on_zone_plane(first_zone, kernel=(3, 3))\n",
+    "%timeit plane_error_img = apply_on_zone_plane(first_zone, kernel=(5, 5))\n",
+    "%timeit plane_error_img = apply_on_zone_plane(first_zone, kernel=(7, 7))"
    ]
   }
  ],