{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coupled Decomposition - Part III"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import os\n",
    "import sys\n",
    "\n",
    "sys.path.insert(0, os.path.abspath('/data/autocnet'))\n",
    "\n",
    "import autocnet\n",
    "from autocnet import CandidateGraph\n",
    "\n",
    "# The GPU based extraction library that contains SIFT extraction and matching\n",
    "import cudasift as cs\n",
    "\n",
    "# A method to resize the images on the fly.\n",
    "from scipy.misc import imresize\n",
    "\n",
    "# Fundamental matrix computation\n",
    "from autocnet.transformation import fundamental_matrix as fm\n",
    "\n",
    "from autocnet.transformation.decompose import coupled_decomposition\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "%pylab inline\n",
    "figsize(16,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = 'AS15-P-0111_CENTER_LRG_CROPPED.png'\n",
    "b = 'AS15-P-0112_CENTER_LRG_CROPPED.png'\n",
    "\n",
    "adj = {a:[b],\n",
    "       b:[a]}\n",
    "\n",
    "cg = CandidateGraph.from_adjacency(adj)\n",
    "\n",
    "# Enable the GPU\n",
    "autocnet.cuda(enable=True, gpu=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Write a custom keypoint extraction function - this could get monkey patched onto the graph object...\n",
    "def extract(arr, downsample_amount=None, **kwargs):\n",
    "    total_size = arr.shape[0] * arr.shape[1]\n",
    "    if not downsample_amount:\n",
    "        downsample_amount = math.ceil(total_size / 12500**2)\n",
    "    shape = (int(arr.shape[0] / downsample_amount), int(arr.shape[1] / downsample_amount))\n",
    "    # Downsample\n",
    "    arr = imresize(arr, shape, interp='lanczos')\n",
    "    \n",
    "    npts = max(arr0.shape) / 3.5\n",
    "    sd = cs.PySiftData(npts)\n",
    "    cs.ExtractKeypoints(arr, sd, **kwargs)\n",
    "    kp, des = sd.to_data_frame()\n",
    "    kp = kp[['x', 'y', 'scale', 'sharpness', 'edgeness', 'orientation', 'score', 'ambiguity']]\n",
    "    kp['score'] = 0.0\n",
    "    kp['ambiguity'] = 0.0\n",
    "    \n",
    "    return kp, des, sd, downsample_amount, arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Write a generic decomposer\n",
    "def custom_decompose(arr0, arr1):\n",
    "    kp0, des0, sd0, downsample_amount0, arr0 = extract(arr0, thresh=1)\n",
    "    kp1, des1, sd1, downsample_amount1, arr1 = extract(arr1, thresh=1)\n",
    "\n",
    "    # Now apply matching, outlier detection, and compute a fundamental matrix\n",
    "    sd0 = cs.PySiftData.from_data_frame(kp0, des0)\n",
    "    sd1 = cs.PySiftData.from_data_frame(kp1, des1)\n",
    "\n",
    "    # Apply the matcher\n",
    "    cs.PyMatchSiftData(sd0, sd1)\n",
    "    matches, _ = sd0.to_data_frame()\n",
    "    # Generic decision about ambiguity and score based on quantiles\n",
    "    \n",
    "    # Apply outlier detection methods for the matches\n",
    "    ambiguity_threshold = matches.ambiguity.quantile(0.01)  # Grabbing the 1%s in this data set\n",
    "    score = matches.score.quantile(0.85)\n",
    "    submatches = matches.query('ambiguity <= {} and score >= {}'.format(ambiguity_threshold, score))\n",
    "    \n",
    "    # Compute a fundamental matrix\n",
    "    kpa = submatches[['x','y']]\n",
    "    kpb = submatches[['match_xpos', 'match_ypos']]\n",
    "    F, mask = fm.compute_fundamental_matrix(kpa, kpb, method='ransac', reproj_threshold=2.0)\n",
    "    F = fm.enforce_singularity_constraint(F)\n",
    "\n",
    "    # Grab the inliers\n",
    "    inliers = submatches[mask]\n",
    "    \n",
    "    # Prepare for coupled decomposition\n",
    "    midx = arr0.shape[1] / 2\n",
    "    midy = arr0.shape[0] / 2\n",
    "\n",
    "    mid = np.array([[midx, midy]])\n",
    "    dists = cdist(mid, inliers[['x', 'y']])\n",
    "    mid_correspondence = inliers.iloc[np.argmin(dists)]\n",
    "    mid_correspondence\n",
    "\n",
    "    # Decompose the images into quadrants\n",
    "    smembership, dmembership, = coupled_decomposition(arr0, arr1,\n",
    "                                                 sorigin=mid_correspondence[['x', 'y']],\n",
    "                                                 dorigin=mid_correspondence[['match_xpos', 'match_ypos']],\n",
    "                                                 theta=0)\n",
    "    \n",
    "    # Return the membership decisions\n",
    "    return smembership, dmembership"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0 11510 0 59720\n",
      "0 11568 0 59770\n",
      "0 0 59720 11510\n",
      "0 0 59770 11568\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "could not broadcast input array from shape (2302,11944) into shape (11510,59720)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-a10dc6618433>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     43\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     44\u001b[0m         \u001b[0;31m# Some fancy indexing to get the dmembership into the right place in the global membership\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 45\u001b[0;31m         \u001b[0mmembership0\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mminsy\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mmaxsy\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mminsx\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mmaxsx\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msmem\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     46\u001b[0m         \u001b[0mmembership1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdypart\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdxpart\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdmem\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0msdy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msdx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m         \u001b[0mpcounter\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m4\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: could not broadcast input array from shape (2302,11944) into shape (11510,59720)"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "gd0 = cg.node[0].geodata\n",
    "gd1 = cg.node[1].geodata\n",
    "membership0 = np.zeros(gd0.raster_size[::-1], dtype=np.int8)\n",
    "membership1 = np.zeros(gd1.raster_size[::-1], dtype=np.int8)\n",
    "\n",
    "# Recursively decompose twice.\n",
    "pcounter = 0\n",
    "for k in range(2):\n",
    "    print(k)\n",
    "    partitions = np.unique(membership0)\n",
    "    for p in partitions:\n",
    "        # Get the source extent as MBR\n",
    "        sypart, sxpart = np.where(membership0 == p)\n",
    "        minsy = np.min(sypart)\n",
    "        maxsy = np.max(sypart) + 1\n",
    "        minsx = np.min(sxpart)\n",
    "        maxsx = np.max(sxpart) + 1\n",
    "        \n",
    "        del sypart, sxpart\n",
    "        \n",
    "        # Get the destination extent as MBR      \n",
    "        dypart, dxpart = np.where(membership1 == p)\n",
    "        mindy = np.min(dypart)\n",
    "        maxdy = np.max(dypart) + 1\n",
    "        mindx = np.min(dxpart)\n",
    "        maxdx = np.max(dxpart) + 1\n",
    "        \n",
    "        # Offsets into classified array\n",
    "        sdy = dypart - min(dypart)\n",
    "        sdx = dxpart - min(dxpart)\n",
    "                \n",
    "        print(minsy, maxsy, minsx, maxsx)\n",
    "        print(mindy, maxdy, mindx, maxdx)\n",
    "        \n",
    "        arr0 = gd0.read_array(pixels=map(int, [minsx, minsy, maxsx-minsx, maxsy-minsy]))\n",
    "        arr1 = gd1.read_array(pixels=map(int,[mindx, mindy, maxdx-mindx, maxdy-mindy]))\n",
    "        \n",
    "        smem, dmem = custom_decompose(arr0, arr1)\n",
    "        smem += pcounter\n",
    "        dmem += pcounter\n",
    "        \n",
    "        # Some fancy indexing to get the dmembership into the right place in the global membership\n",
    "        membership0[minsy:maxsy,minsx:maxsx] = smem\n",
    "        membership1[dypart, dxpart] = dmem[sdy, sdx]\n",
    "        pcounter += 4\n",
    "        \n",
    "        # Force cleanup\n",
    "        arr0 = None\n",
    "        arr1 = None\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "AutoCNet",
   "language": "python",
   "name": "autocnet"
  },
  "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
