{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os \n",
    "os.environ[\"ISISROOT\"] = \"/usgs/cpkgs/anaconda3_linux/envs/isis3.9.0\"\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from plio.io.io_gdal import GeoDataset \n",
    "import numpy as np\n",
    "\n",
    "from glob import glob\n",
    "\n",
    "from shapely import wkt\n",
    "from shapely.geometry import Point\n",
    "import geopandas as gpd\n",
    "\n",
    "from autocnet.examples import get_path\n",
    "from autocnet.graph.network import CandidateGraph\n",
    "from autocnet.graph.edge import Edge\n",
    "from autocnet.spatial.isis import point_info\n",
    "from autocnet.cg import change_detection as cd\n",
    "\n",
    "from plio.io.io_gdal import GeoDataset\n",
    "import kalasiris as isis\n",
    "from subprocess import CalledProcessError\n",
    "\n",
    "from IPython.display import display\n",
    "\n",
    "import shutil\n",
    "\n",
    "from autocnet.utils import hirise \n",
    "from autocnet.utils.utils import bytescale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x2ab10795da20>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 108936 characters]\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df1 = hirise.load_segments(\"/scratch/krodriguez/hirise/no_change/ESP_012024_1440/\")\n",
    "df2 = hirise.load_segments(\"/scratch/krodriguez/hirise/no_change/ESP_020661_1440/\")\n",
    "\n",
    "ax = df1.plot(figsize=(10,10), alpha=.3)\n",
    "df2.plot(figsize=(20,50), alpha=.3, ax=ax, color=\"orange\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12032 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = Point(129.393, -35.703)\n",
    "\n",
    "roi1 = df1[df1.intersects(p)].reindex()\n",
    "roi1_geoms = roi1[\"footprint\"]\n",
    "\n",
    "roi2 = df2[df2.intersects(roi1.footprint.unary_union)]\n",
    "roi2 = roi2[roi2.intersects(p)].iloc[0:1]\n",
    "\n",
    "\n",
    "ax = roi1.plot(alpha=.5)\n",
    "roi2.plot(ax=ax, alpha=.5, color=\"orange\")\n",
    "\n",
    "roi1 = roi1.iloc[0]\n",
    "roi2 = roi2.iloc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "new_dir = '/tmp/'\n",
    "\n",
    "roi1_path = os.path.join(new_dir, os.path.basename(roi1.path))\n",
    "roi2_path = os.path.join(new_dir, os.path.basename(roi2.path))\n",
    "\n",
    "shutil.copyfile(roi1.path, roi1_path)\n",
    "shutil.copyfile(roi2.path, roi2_path)\n",
    "\n",
    "with open(\"/tmp/test.lis\", \"w\") as f:\n",
    "    f.write(\"\\n\".join([roi1_path, roi2_path]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/krodriguez/.conda/envs/autocnetdev/lib/python3.6/site-packages/pandas/core/frame.py:4025: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  return super(DataFrame, self).rename(**kwargs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 33.9 s, sys: 716 ms, total: 34.6 s\n",
      "Wall time: 34.6 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "#Point to the adjacency Graph\n",
    "adjacency = {roi1_path: [roi2_path], roi2_path: [roi1_path]}\n",
    "\n",
    "cg = CandidateGraph.from_adjacency(adjacency)\n",
    "\n",
    "#Apply SIFT to extract features\n",
    "cg.extract_features(extractor_method='vlfeat', extractor_parameters={'nfeatures':100})\n",
    "cg.match()\n",
    "\n",
    "# Apply outlier detection\n",
    "cg.apply_func_to_edges(Edge.symmetry_check)\n",
    "cg.apply_func_to_edges(Edge.ratio_check)\n",
    "cg.minimum_spanning_tree()\n",
    "\n",
    "# Compute a homography and apply RANSAC\n",
    "cg.apply_func_to_edges(\"compute_fundamental_matrix\", clean_keys=['ratio', 'symmetry'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x2ab0bdc37978>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 182480 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cg[1][2][\"data\"].plot(clean_keys=[\"fundamental\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cg.generate_control_network(clean_keys=[\"fundamental\"])\n",
    "cg.to_isis(\"test.net\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/tmp/ESP_020661_1440_RED4.stitched.norm.30421_33069.cub\r\n",
      "/tmp/ESP_012024_1440_RED8.stitched.norm.22229_24877.cub\r\n"
     ]
    }
   ],
   "source": [
    "cat test.lis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Validating network...\n",
      "Validation complete!...\n",
      "starting iteration 1\n",
      "\n",
      "\n",
      "\n",
      "median deviation: 0.428527\n",
      "\n",
      "mad: 0.635334\n",
      "\n",
      "Rejection Limit: 2.73574\n",
      "\n",
      "\n",
      "Rejected Observations:0 (Rejection Limit:2.73574)\n",
      "\n",
      "Measures that came back: 0\n",
      "Iteration: 1 \n",
      "Sigma0:         0.7409386686 \n",
      "Observations: 636 \n",
      "Constrained Parameters:477 \n",
      "Unknowns: 489 \n",
      "Degrees of Freedom: 624 \n",
      "End of Iteration 1 \n",
      "Elapsed Time:         0.0200000000 \n",
      "Group = Iteration1\n",
      "  Sigma0                       = 0.74093866855463\n",
      "  Observations                 = 636\n",
      "  Constrained_Point_Parameters = 477\n",
      "  Constrained_Image_Parameters = 0\n",
      "  Unknown_Parameters           = 489\n",
      "  Degrees_of_Freedom           = 624\n",
      "  Rejected_Measures            = 0\n",
      "End_Group\n",
      "starting iteration 2\n",
      "\n",
      "\n",
      "\n",
      "median deviation: 0.46773\n",
      "\n",
      "mad: 0.693457\n",
      "\n",
      "Rejection Limit: 2.98239\n",
      "\n",
      "\n",
      "Rejected Observations:0 (Rejection Limit:2.98239)\n",
      "\n",
      "Measures that came back: 0\n",
      "Iteration: 2 \n",
      "Sigma0:         0.6506659496 \n",
      "Observations: 636 \n",
      "Constrained Parameters:477 \n",
      "Unknowns: 489 \n",
      "Degrees of Freedom: 624 \n",
      "End of Iteration 2 \n",
      "Elapsed Time:         0.0200000000 \n",
      "\n",
      "Group = Iteration2\n",
      "  Sigma0                       = 0.65066594956289\n",
      "  Observations                 = 636\n",
      "  Constrained_Point_Parameters = 477\n",
      "  Constrained_Image_Parameters = 0\n",
      "  Unknown_Parameters           = 489\n",
      "  Degrees_of_Freedom           = 624\n",
      "  Rejected_Measures            = 0\n",
      "End_Group\n",
      "starting iteration 3\n",
      "\n",
      "\n",
      "\n",
      "median deviation: 0.460487\n",
      "\n",
      "mad: 0.682719\n",
      "\n",
      "Rejection Limit: 2.9455\n",
      "\n",
      "\n",
      "Rejected Observations:0 (Rejection Limit:2.9455)\n",
      "\n",
      "Measures that came back: 0\n",
      "Iteration: 3 \n",
      "Sigma0:         0.6506067898 \n",
      "Observations: 636 \n",
      "Constrained Parameters:477 \n",
      "Unknowns: 489 \n",
      "Degrees of Freedom: 624 \n",
      "End of Iteration 3 \n",
      "Elapsed Time:         0.0300000000 \n",
      "\n",
      "Group = Iteration3\n",
      "  Sigma0                       = 0.65060678982728\n",
      "  Observations                 = 636\n",
      "  Constrained_Point_Parameters = 477\n",
      "  Constrained_Image_Parameters = 0\n",
      "  Unknown_Parameters           = 489\n",
      "  Degrees_of_Freedom           = 624\n",
      "  Rejected_Measures            = 0\n",
      "End_Group\n",
      "starting iteration 4\n",
      "\n",
      "\n",
      "\n",
      "median deviation: 0.46054\n",
      "\n",
      "mad: 0.682797\n",
      "\n",
      "Rejection Limit: 2.94579\n",
      "\n",
      "\n",
      "Rejected Observations:0 (Rejection Limit:2.94579)\n",
      "\n",
      "Measures that came back: 0\n",
      "Iteration: 4 \n",
      "Sigma0:         0.6506067571 \n",
      "Observations: 636 \n",
      "Constrained Parameters:477 \n",
      "Unknowns: 489 \n",
      "Degrees of Freedom: 624 \n",
      "End of Iteration 4 \n",
      "Elapsed Time:         0.0200000000 \n",
      "\n",
      "Group = Iteration4\n",
      "  Sigma0                       = 0.65060675706612\n",
      "  Observations                 = 636\n",
      "  Constrained_Point_Parameters = 477\n",
      "  Constrained_Image_Parameters = 0\n",
      "  Unknown_Parameters           = 489\n",
      "  Degrees_of_Freedom           = 624\n",
      "  Rejected_Measures            = 0\n",
      "End_Group\n",
      "starting iteration 5\n",
      "\n",
      "\n",
      "\n",
      "median deviation: 0.460549\n",
      "\n",
      "mad: 0.68281\n",
      "\n",
      "Rejection Limit: 2.94584\n",
      "\n",
      "\n",
      "Rejected Observations:0 (Rejection Limit:2.94584)\n",
      "\n",
      "Measures that came back: 0\n",
      "Iteration: 5 \n",
      "Sigma0:         0.6506067566 \n",
      "Observations: 636 \n",
      "Constrained Parameters:477 \n",
      "Unknowns: 489 \n",
      "Degrees of Freedom: 624 \n",
      "End of Iteration 5 \n",
      "Elapsed Time:         0.0200000000 \n",
      "\n",
      "Group = Iteration5\n",
      "  Sigma0                       = 0.65060675662698\n",
      "  Observations                 = 636\n",
      "  Constrained_Point_Parameters = 477\n",
      "  Constrained_Image_Parameters = 0\n",
      "  Unknown_Parameters           = 489\n",
      "  Degrees_of_Freedom           = 624\n",
      "  Rejected_Measures            = 0\n",
      "End_Group\n",
      "starting iteration 6\n",
      "\n",
      "\n",
      "\n",
      "median deviation: 0.46055\n",
      "\n",
      "mad: 0.682812\n",
      "\n",
      "Rejection Limit: 2.94585\n",
      "\n",
      "\n",
      "Rejected Observations:0 (Rejection Limit:2.94585)\n",
      "\n",
      "Measures that came back: 0\n",
      "Iteration: 6 \n",
      "Sigma0:         0.6506067566 \n",
      "Observations: 636 \n",
      "Constrained Parameters:477 \n",
      "Unknowns: 489 \n",
      "Degrees of Freedom: 624 \n",
      "Bundle has converged\n",
      "\n",
      "Starting Error Propagation     Time: 2020-01-30T14:05:54\n",
      "\n",
      "\r",
      "Error Propagation: Inverse Block 1 of 2; Point 1 of 159\r",
      "Error Propagation: Inverse Block 1 of 2; Point 101 of 159\r",
      "Error Propagation: Inverse Block 2 of 2; Point 1 of 159\r",
      "Error Propagation: Inverse Block 2 of 2; Point 101 of 159\n",
      "\n",
      "\r",
      "Filling point covariance matrices: Time 2020-01-30T14:05:54\n",
      "\n",
      "\r",
      "Error Propagation: Filling point covariance matrices 1 of 159\r",
      "\r",
      "Error Propagation: Filling point covariance matrices 101 of 159\r\n",
      "\n",
      "Error Propagation Complete\n",
      "\n",
      "Bundle Complete\n",
      "\n",
      "Group = \"Iteration6: Final\"\n",
      "  Sigma0                       = 0.65060675662956\n",
      "  Observations                 = 636\n",
      "  Constrained_Point_Parameters = 477\n",
      "  Constrained_Image_Parameters = 0\n",
      "  Unknown_Parameters           = 489\n",
      "  Degrees_of_Freedom           = 624\n",
      "  Rejected_Measures            = 0\n",
      "  Converged                    = TRUE\n",
      "  TotalElapsedTime             = 0.15\n",
      "  ErrorPropagationElapsedTime  = 0.0\n",
      "End_Group\n",
      "\n",
      "Generating report files\n",
      "\n",
      "\n",
      "Group = JigsawResults\n",
      "  Status = \"Camera pointing updated\"\n",
      "End_Group\n",
      "\n"
     ]
    }
   ],
   "source": [
    "try:\n",
    "    output = isis.jigsaw(fromlist='test.lis', cnet='test.net', onet='bundle_test.cnet', update='yes', radius='yes', errorpropagation='yes', outlier_rejection='yes', point_longitude_sigma=10, point_latitude_sigma=10, point_radius_sigma=2, maxits=20, camsolve='VELOCITIES', twist='yes', overexisting='yes', spsolve='NONE')\n",
    "    print(output.stdout)\n",
    "except CalledProcessError as e:\n",
    "    print(' '.join(e.cmd))\n",
    "    print(e.stdout\n",
    "    print(e.stderr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "roi1_proj = os.path.splitext(roi1_path)[0] + \".proj.cub\"\n",
    "roi2_proj = os.path.splitext(roi2_path)[0] + \".proj.cub\"\n",
    "\n",
    "try:\n",
    "    isis.cam2map(from_=roi1_path, to=roi1_proj, map='/usgs/cpkgs/isis3/data/base/templates/maps/equirectangular.map')\n",
    "    isis.cam2map(from_=roi2_path, to=roi2_proj, patchsize=8, map=roi1_proj, matchmap=True, warpalgorithm=\"REVERSEPATCH\")\n",
    "except CalledProcessError as e:\n",
    "    print(' '.join(e.cmd))\n",
    "    print(e.stdout\n",
    "    print(e.stderr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2ab115250748>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 81144 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "roi1_proj_geo = GeoDataset(roi1_proj)\n",
    "roi2_proj_geo = GeoDataset(roi2_proj)\n",
    "isis_null = isis.specialpixels.Real.Null\n",
    "\n",
    "proj_array1 = roi1_proj_geo.read_array()\n",
    "proj_array2 = roi2_proj_geo.read_array()\n",
    "\n",
    "proj_array1[proj_array1 == isis_null] = np.nan\n",
    "proj_array2[proj_array2 == isis_null] = np.nan\n",
    "\n",
    "plt.imshow(proj_array1, alpha=.4)\n",
    "plt.imshow(proj_array2, alpha=.4)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# apply change detection algo\n",
    "polys, weights, diffarr = cd.okubogar_detector(roi1_proj_geo, roi2_proj_geo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 79292 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib\n",
    "\n",
    "plt.imshow(proj_array1, cmap=\"Greys_r\")\n",
    "plt.imshow(proj_array2, cmap=\"Greys_r\")\n",
    "\n",
    "if weights:\n",
    "    cmap = matplotlib.cm.get_cmap('plasma')\n",
    "    colors = bytescale(np.asarray(weights))\n",
    "    ax = [plt.fill(*p.exterior.xy, c=cmap(colors[i]), alpha=.6) for i,p in enumerate(polys)]\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "autocnetdev",
   "language": "python",
   "name": "autocnetdev"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
