{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Imports" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from pysis import isis\n", "\n", "from plio.io import io_controlnetwork\n", "from knoten.csm import create_csm\n", "from scipy import sparse\n", "import ale\n", "import csmapi\n", "import numpy as np\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "from knoten.bundle import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load Network and Generate Sensors\n", "- Generate a set of USGSCSM sensor models from a list of ISIS cube files\n", "- Generate a plio dataframe from an ISIS control network\n", "- Compute a priori ground points for all of the free points in a control network." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cubes = 'data/bundle-adjust/cubelist.lis'\n", "sensors = generate_sensors(cubes, only_isis_spice=True, verbose=True)\n", "\n", "network = 'data/bundle-adjust/controlnet.net'\n", "cnet = io_controlnetwork.from_isis(network)\n", "cnet = compute_apriori_ground_points(cnet, sensors) # autoseed did not generate ground points, calculate and repopulate the data frame" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Determine Which Sensor Parameters to Solve For\n", "Get a set of the CSM parameters for each CSM sensor in the sensors set" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Image: MRO/CTX/1016816309:030\n", " IT Pos. Bias | 0 | 0.0\n", " CT Pos. Bias | 1 | 0.0\n", " Rad Pos. Bias | 2 | 0.0\n", " IT Vel. Bias | 3 | 0.0\n", " CT Vel. Bias | 4 | 0.0\n", " Rad Vel. Bias | 5 | 0.0\n", " Omega Bias | 6 | 0.0\n", " Phi Bias | 7 | 0.0\n", " Kappa Bias | 8 | 0.0\n", " Omega Rate | 9 | 0.0\n", " Phi Rate | 10 | 0.0\n", " Kappa Rate | 11 | 0.0\n", " Omega Accl | 12 | 0.0\n", " Phi Accl | 13 | 0.0\n", " Kappa Accl | 14 | 0.0\n", " Focal Bias | 15 | 0.0\n", "Image: MRO/CTX/0884541979:200\n", " IT Pos. Bias | 0 | 0.0\n", " CT Pos. Bias | 1 | 0.0\n", " Rad Pos. Bias | 2 | 0.0\n", " IT Vel. Bias | 3 | 0.0\n", " CT Vel. Bias | 4 | 0.0\n", " Rad Vel. Bias | 5 | 0.0\n", " Omega Bias | 6 | 0.0\n", " Phi Bias | 7 | 0.0\n", " Kappa Bias | 8 | 0.0\n", " Omega Rate | 9 | 0.0\n", " Phi Rate | 10 | 0.0\n", " Kappa Rate | 11 | 0.0\n", " Omega Accl | 12 | 0.0\n", " Phi Accl | 13 | 0.0\n", " Kappa Accl | 14 | 0.0\n", " Focal Bias | 15 | 0.0\n", "Image: MRO/CTX/0891288961:044\n", " IT Pos. Bias | 0 | 0.0\n", " CT Pos. Bias | 1 | 0.0\n", " Rad Pos. Bias | 2 | 0.0\n", " IT Vel. Bias | 3 | 0.0\n", " CT Vel. Bias | 4 | 0.0\n", " Rad Vel. Bias | 5 | 0.0\n", " Omega Bias | 6 | 0.0\n", " Phi Bias | 7 | 0.0\n", " Kappa Bias | 8 | 0.0\n", " Omega Rate | 9 | 0.0\n", " Phi Rate | 10 | 0.0\n", " Kappa Rate | 11 | 0.0\n", " Omega Accl | 12 | 0.0\n", " Phi Accl | 13 | 0.0\n", " Kappa Accl | 14 | 0.0\n", " Focal Bias | 15 | 0.0\n" ] } ], "source": [ "all_parameters = {sn: get_sensor_parameters(sensor) for sn, sensor in sensors.items()}\n", "for sn, parameters in all_parameters.items():\n", " print(f\"Image: {sn}\")\n", " for param in parameters:\n", " print(f\" {param.name} | {param.index} | {param.value}\")" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Solve for angles and angular rates\n", "solve_parameters = {sn: params[6:12] for sn, params in all_parameters.items()}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute the Column Indices for Parameters\n", "Obtain dictionary that maps serial numbers and point IDs to the column range their parameters are in the Jacobian matrix." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "column_dict = compute_coefficient_columns(cnet, sensors, solve_parameters)\n", "# num_parameters = max(col_range[1] for col_range in column_dict.values())" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "source": [ "## Compute the Weight Matrix\n", "- According to the weighted Normal equation (J.TWJ), W needs to be a square matrix the size of (# of measures)x2. So it is the weight of the observations. In ISIS, the weight of the observations are an inverted function of the size of the pixels on the focal plane (resolution). However, in csm we do not have access to that information. \n", "- For the time being, since we are working exclusively with CTX images we are going to set the weight matrix equal to the identity matrix -> all observations have the same weight." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "num_observations = 2 * len(cnet)\n", "W_observations = np.eye(num_observations) # this is a place holder until Jesse adds his calculations\n", "W_params = compute_parameter_weights(cnet, sensors, solve_parameters, column_dict)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate Initial Sigma0\n", "Compute the resulting standard deviation of the residuals for the current state of the bundle network." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.5961263976217976\n" ] } ], "source": [ "V = compute_residuals(cnet, sensors)\n", "dX = np.zeros(W_params.shape[0])\n", "sigma0 = compute_sigma0(V, dX, W_params, W_observations)\n", "\n", "print((sigma0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Populate Jacobian\n", "Compute the Jacobian matrix with controlnet, set of sensors, solve parameters, and coefficient columns" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "J = compute_jacobian(cnet, sensors, solve_parameters, column_dict)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bundle Iteration" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": true }, "outputs": [], "source": [ "def bundle_iteration(J, V, W_parameters, W_observations):\n", " \"\"\"\n", " Parameters\n", " ----------\n", " J : ndarray\n", " The Jacobian matrix\n", " V : np.array\n", " An array of residuals of the difference between registered measure \n", " and back projected ground points in image space.\n", " W_parameters : ndarray \n", " The parameter weight matrix (i.e.: sensor parameters and point weights)\n", " W_observations : ndarray\n", " The observation weight matrix (i.e.: measure weights)\n", " \n", " Returns\n", " -------\n", " N : np.ndarray\n", " Normal equation matrix \n", " \n", " dX : np.ndarray\n", " An array of updated parameter values\n", " \"\"\"\n", " \n", " N = J.T.dot(W_observations).dot(J) + W_parameters\n", " C = J.T.dot(W_observations).dot(V)\n", " dX = np.linalg.inv(N).dot(C)\n", " return N, dX" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1365,)\n" ] } ], "source": [ "N, dX = bundle_iteration(J, V, W_params, W_observations)\n", "print(dX.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate Updated Sigma0" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.676667701883667\n" ] } ], "source": [ "dof = W_observations.shape[0] - W_params.shape[0]\n", "updated_sigma0 = np.sqrt((V.dot(W_observations).dot(V) - dX.dot(J.T).dot(W_observations).dot(V))/dof)\n", "print(updated_sigma0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Redundancy Number" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Minimum redundancy: 0.9236070227160993\n", "Maximum redundancy: 0.9887066971286086\n" ] }, { "data": { "text/plain": [ "{'whiskers': [,\n", " ],\n", " 'caps': [,\n", " ],\n", " 'boxes': [],\n", " 'medians': [],\n", " 'fliers': [],\n", " 'means': []}" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAp3klEQVR4nO3df1DV153/8RfeXuAagZighET0oqjcWWwSrxaFMCOTKQ4VK1V3cF3ddUdt3LgzMabdDYl+tzETGSeFtZsIFZWJP3aKHWWS1NBsmck6wWD3Lte4lSpgmiCuYimsARMU2Mvn+4df7ndvQOsl1nu4Ph8zn1HOPZ8P709mMvfl+ZzPORGWZVkCAAAw2JhQFwAAAPDHEFgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMb7RqgLuFsGBgZ0+fJlxcTEKCIiItTlAACAO2BZlq5du6ZHH31UY8bcehwlbALL5cuXlZSUFOoyAADACFy8eFGTJk265edhE1hiYmIk3bzh2NjYEFcDAADuRHd3t5KSkvzf47cSNoFl8DFQbGwsgQUAgFHmj03nYNItAAAwHoEFAAAYj8ACAACMN6LAUlpaquTkZEVHR8vtdqu2tva2/Xft2iWXyyWHw6GZM2fqwIEDAZ/39/dr27ZtmjZtmqKjo/X444/r/fffH0lpAAAgDAUdWA4fPqxNmzbp5Zdf1scff6ysrCzl5uaqtbV12P5lZWUqLCzUj370I/32t7/VK6+8oo0bN+oXv/iFv8+WLVu0e/duvfHGGzp79qw2bNig733ve/r4449HfmcAACBsRFiWZQVzQnp6umbPnq2ysjJ/m8vlUn5+voqKiob0z8jIUGZmpl5//XV/26ZNm1RfX68TJ05Ikh599FG9/PLL2rhxo79Pfn6+xo0bp0OHDt1RXd3d3YqLi1NXVxdvCQEAMErc6fd3UCMsfX198nq9ysnJCWjPyclRXV3dsOf09vYqOjo6oM3hcMjj8ai/v/+2fQYDza2u293dHXAAAIDwFFRg6ejokM/nU0JCQkB7QkKCrly5Muw5Cxcu1N69e+X1emVZlurr61VRUaH+/n51dHT4+5SUlOj8+fMaGBhQTU2N3nnnHbW1td2ylqKiIsXFxfkPVrkFACB8jWjS7VcXd7Es65YLvmzdulW5ubmaN2+e7Ha7lixZojVr1kiSbDabJOknP/mJpk+frtTUVEVGRurv/u7v9Dd/8zf+z4dTWFiorq4u/3Hx4sWR3AoAABgFggos8fHxstlsQ0ZT2tvbh4y6DHI4HKqoqFBPT49aWlrU2toqp9OpmJgYxcfHS5ImTJigt99+W19++aUuXLigxsZGjRs3TsnJybesJSoqyr+qLavbAgAQ3oIKLJGRkXK73aqpqQlor6mpUUZGxm3PtdvtmjRpkmw2myorK5WXlzdkV8bo6Gg99thj+p//+R8dPXpUS5YsCaY8AAAQpoLeS2jz5s1avXq15syZo/nz56u8vFytra3asGGDpJuPai5duuRfa6W5uVkej0fp6em6evWqSkpK1NDQoP379/uv+e///u+6dOmSnnjiCV26dEk/+tGPNDAwoL//+7+/S7cJAABGs6ADS0FBgTo7O7Vt2za1tbUpLS1N1dXVmjJliiSpra0tYE0Wn8+n4uJiNTU1yW63Kzs7W3V1dXI6nf4+N27c0JYtW/Tpp59q3Lhx+s53vqODBw/qwQcf/No3CCA0enp61NjYeFeudf36dbW0tMjpdMrhcHzt66Wmpmrs2LF3oTIA90rQ67CYinVYALOcOnVKbrc71GUMy+v1avbs2aEuA4Du/Ps76BEWALgTqamp8nq9d+Va586d06pVq3To0CG5XK6vfb3U1NS7UBWAe4nAAuBPYuzYsXd9FMPlcjEyAtyn2K0ZAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgjCiylpaVKTk5WdHS03G63amtrb9t/165dcrlccjgcmjlzpg4cODCkz86dOzVz5kw5HA4lJSXp+eef140bN0ZSHgAACDPfCPaEw4cPa9OmTSotLVVmZqZ2796t3NxcnT17VpMnTx7Sv6ysTIWFhdqzZ4/mzp0rj8ej9evXa/z48Vq8eLEk6V/+5V/04osvqqKiQhkZGWpubtaaNWskSf/0T//09e4QAACMekGPsJSUlGjt2rVat26dXC6Xdu7cqaSkJJWVlQ3b/+DBg3rmmWdUUFCgqVOnasWKFVq7dq127Njh73Py5EllZmZq5cqVcjqdysnJ0V/8xV+ovr5+5HcGAADCRlCBpa+vT16vVzk5OQHtOTk5qqurG/ac3t5eRUdHB7Q5HA55PB719/dLkp566il5vV55PB5J0qeffqrq6motWrTolrX09vaqu7s74AAAAOEpqMDS0dEhn8+nhISEgPaEhARduXJl2HMWLlyovXv3yuv1yrIs1dfXq6KiQv39/ero6JAkrVixQq+++qqeeuop2e12TZs2TdnZ2XrxxRdvWUtRUZHi4uL8R1JSUjC3AgAARpERTbqNiIgI+NmyrCFtg7Zu3arc3FzNmzdPdrtdS5Ys8c9PsdlskqTjx4/rtddeU2lpqU6dOqWqqiodO3ZMr7766i1rKCwsVFdXl/+4ePHiSG4FAACMAkEFlvj4eNlstiGjKe3t7UNGXQY5HA5VVFSop6dHLS0tam1tldPpVExMjOLj4yXdDDWrV6/WunXrNGvWLH3ve9/T9u3bVVRUpIGBgWGvGxUVpdjY2IADAACEp6ACS2RkpNxut2pqagLaa2pqlJGRcdtz7Xa7Jk2aJJvNpsrKSuXl5WnMmJu/vqenx//3QTabTZZlybKsYEoEAABhKOjXmjdv3qzVq1drzpw5mj9/vsrLy9Xa2qoNGzZIuvmo5tKlS/61Vpqbm+XxeJSenq6rV6+qpKREDQ0N2r9/v/+aixcvVklJiZ588kmlp6frk08+0datW/Xd737X/9gIAADcv4IOLAUFBers7NS2bdvU1tamtLQ0VVdXa8qUKZKktrY2tba2+vv7fD4VFxerqalJdrtd2dnZqqurk9Pp9PfZsmWLIiIitGXLFl26dEkTJkzQ4sWL9dprr339OwQAAKNehBUmz1y6u7sVFxenrq4u5rMAYebUqVNyu93yer2aPXt2qMsBcBfd6fc3ewkBAADjBf1ICED4O3/+vK5duxbqMvzOnTsX8KcpYmJiNH369FCXAdwXCCwAApw/f14zZswIdRnDWrVqVahLGKK5uZnQAtwDBBYAAQZHVg4dOiSXyxXiam66fv26Wlpa5HQ65XA4Ql2OpJujPatWrTJqJAoIZwQWAMNyuVxGTXDNzMwMdQkAQohJtwAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYb0SBpbS0VMnJyYqOjpbb7VZtbe1t++/atUsul0sOh0MzZ87UgQMHAj5fsGCBIiIihhyLFi0aSXkAACDMfCPYEw4fPqxNmzaptLRUmZmZ2r17t3Jzc3X27FlNnjx5SP+ysjIVFhZqz549mjt3rjwej9avX6/x48dr8eLFkqSqqir19fX5z+ns7NTjjz+uP//zP/8atwYAAMJF0CMsJSUlWrt2rdatWyeXy6WdO3cqKSlJZWVlw/Y/ePCgnnnmGRUUFGjq1KlasWKF1q5dqx07dvj7PPTQQ3rkkUf8R01NjcaOHUtgAQAAkoIMLH19ffJ6vcrJyQloz8nJUV1d3bDn9Pb2Kjo6OqDN4XDI4/Gov79/2HP27dunFStW6IEHHrhlLb29veru7g44AABAeAoqsHR0dMjn8ykhISGgPSEhQVeuXBn2nIULF2rv3r3yer2yLEv19fWqqKhQf3+/Ojo6hvT3eDxqaGjQunXrbltLUVGR4uLi/EdSUlIwtwIAAEaREU26jYiICPjZsqwhbYO2bt2q3NxczZs3T3a7XUuWLNGaNWskSTabbUj/ffv2KS0tTd/61rduW0NhYaG6urr8x8WLF0dyKwAAYBQIKrDEx8fLZrMNGU1pb28fMuoyyOFwqKKiQj09PWppaVFra6ucTqdiYmIUHx8f0Lenp0eVlZV/dHRFkqKiohQbGxtwAACA8BRUYImMjJTb7VZNTU1Ae01NjTIyMm57rt1u16RJk2Sz2VRZWam8vDyNGRP463/+85+rt7dXq1atCqYsAAAQ5oJ+rXnz5s1avXq15syZo/nz56u8vFytra3asGGDpJuPai5duuRfa6W5uVkej0fp6em6evWqSkpK1NDQoP379w+59r59+5Sfn6+HH374a94WAAAIJ0EHloKCAnV2dmrbtm1qa2tTWlqaqqurNWXKFElSW1ubWltb/f19Pp+Ki4vV1NQku92u7Oxs1dXVyel0Bly3ublZJ06c0K9+9auvd0cAACDsBB1YJOnZZ5/Vs88+O+xnb731VsDPLpdLH3/88R+95owZM2RZ1kjKAQAAYY69hAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeCNamh9AeHtkXIQcnzdLl/k3za04Pm/WI+MiQl0GcN8gsAAY4hl3pFwfPiN9GOpKzOXSzf9OAO4NAguAIXZ7+1Twf96SKzU11KUY61xjo3YXr9R3Q10IcJ8gsAAY4soXlq4/OEN69IlQl2Ks61cGdOULdpgH7hUeUAMAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHrs1AwjQ09MjSTp16lSIK/n/rl+/rpaWFjmdTjkcjlCXI0k6d+5cqEsA7isEFgABGhsbJUnr168PcSWjQ0xMTKhLAO4LBBYAAfLz8yVJqampGjt2bGiL+X/OnTunVatW6dChQ3K5XKEuxy8mJkbTp08PdRnAfYHAAiBAfHy81q1bF+oyhuVyuTR79uxQlwEgBJh0CwAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMN6LAUlpaquTkZEVHR8vtdqu2tva2/Xft2iWXyyWHw6GZM2fqwIEDQ/p8/vnn2rhxoxITExUdHS2Xy6Xq6uqRlAcAAMJM0EvzHz58WJs2bVJpaakyMzO1e/du5ebm6uzZs5o8efKQ/mVlZSosLNSePXs0d+5ceTwerV+/XuPHj9fixYslSX19ffr2t7+tiRMn6siRI5o0aZIuXrzIpmIAAECSFGFZlhXMCenp6Zo9e7bKysr8bS6XS/n5+SoqKhrSPyMjQ5mZmXr99df9bZs2bVJ9fb1OnDghSfrpT3+q119/XY2NjbLb7SO6ke7ubsXFxamrq0uxsbEjugYAM506dUput1ter5e9hIAwc6ff30E9Eurr65PX61VOTk5Ae05Ojurq6oY9p7e3V9HR0QFtDodDHo9H/f39kqR3331X8+fP18aNG5WQkKC0tDRt375dPp8vmPIAAECYCiqwdHR0yOfzKSEhIaA9ISFBV65cGfachQsXau/evfJ6vbIsS/X19aqoqFB/f786OjokSZ9++qmOHDkin8+n6upqbdmyRcXFxXrttdduWUtvb6+6u7sDDgAAEJ5GNOk2IiIi4GfLsoa0Ddq6datyc3M1b9482e12LVmyRGvWrJEk2Ww2SdLAwIAmTpyo8vJyud1urVixQi+//HLAY6evKioqUlxcnP9ISkoaya0AAIBRIKjAEh8fL5vNNmQ0pb29fcioyyCHw6GKigr19PSopaVFra2tcjqdiomJUXx8vCQpMTFRM2bM8AcY6ea8mCtXrqivr2/Y6xYWFqqrq8t/XLx4MZhbAQAAo0hQgSUyMlJut1s1NTUB7TU1NcrIyLjtuXa7XZMmTZLNZlNlZaXy8vI0ZszNX5+ZmalPPvlEAwMD/v7Nzc1KTExUZGTksNeLiopSbGxswAEAAMJT0I+ENm/erL1796qiokLnzp3T888/r9bWVm3YsEHSzZGPv/qrv/L3b25u1qFDh3T+/Hl5PB6tWLFCDQ0N2r59u7/P3/7t36qzs1PPPfecmpub9d5772n79u3auHHjXbhFAAAw2gW9DktBQYE6Ozu1bds2tbW1KS0tTdXV1ZoyZYokqa2tTa2trf7+Pp9PxcXFampqkt1uV3Z2turq6uR0Ov19kpKS9Ktf/UrPP/+8vvnNb+qxxx7Tc889p3/4h3/4+ncIAABGvaDXYTEV67AA4Yt1WIDw9SdZhwUAACAUCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8UYUWEpLS5WcnKzo6Gi53W7V1tbetv+uXbvkcrnkcDg0c+ZMHThwIODzt956SxEREUOOGzdujKQ8AAAQZr4R7AmHDx/Wpk2bVFpaqszMTO3evVu5ubk6e/asJk+ePKR/WVmZCgsLtWfPHs2dO1cej0fr16/X+PHjtXjxYn+/2NhYNTU1BZwbHR09glsCAADhJujAUlJSorVr12rdunWSpJ07d+pf//VfVVZWpqKioiH9Dx48qGeeeUYFBQWSpKlTp+rXv/61duzYERBYIiIi9Mgjj4z0PgAAQBgL6pFQX1+fvF6vcnJyAtpzcnJUV1c37Dm9vb1DRkocDoc8Ho/6+/v9bV988YWmTJmiSZMmKS8vTx9//PFta+nt7VV3d3fAAQAAwlNQgaWjo0M+n08JCQkB7QkJCbpy5cqw5yxcuFB79+6V1+uVZVmqr69XRUWF+vv71dHRIUlKTU3VW2+9pXfffVc/+9nPFB0drczMTJ0/f/6WtRQVFSkuLs5/JCUlBXMrAABgFBnRpNuIiIiAny3LGtI2aOvWrcrNzdW8efNkt9u1ZMkSrVmzRpJks9kkSfPmzdOqVav0+OOPKysrSz//+c81Y8YMvfHGG7esobCwUF1dXf7j4sWLI7kVAAAwCgQVWOLj42Wz2YaMprS3tw8ZdRnkcDhUUVGhnp4etbS0qLW1VU6nUzExMYqPjx++qDFjNHfu3NuOsERFRSk2NjbgAAAA4SmowBIZGSm3262ampqA9pqaGmVkZNz2XLvdrkmTJslms6myslJ5eXkaM2b4X29Zlk6fPq3ExMRgygMAAGEq6LeENm/erNWrV2vOnDmaP3++ysvL1draqg0bNki6+ajm0qVL/rVWmpub5fF4lJ6erqtXr6qkpEQNDQ3av3+//5qvvPKK5s2bp+nTp6u7u1v//M//rNOnT2vXrl136TYBAMBoFnRgKSgoUGdnp7Zt26a2tjalpaWpurpaU6ZMkSS1tbWptbXV39/n86m4uFhNTU2y2+3Kzs5WXV2dnE6nv8/nn3+u73//+7py5Yri4uL05JNP6sMPP9S3vvWtr3+HAABg1IuwLMsKdRF3Q3d3t+Li4tTV1cV8FiDMnDp1Sm63W16vV7Nnzw51OQDuojv9/mYvIQAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8b4R6gIAhKeenh41NjbelWudO3cu4M+vKzU1VWPHjr0r1wJwbxBYAPxJNDY2yu1239Vrrlq16q5cx+v1avbs2XflWgDuDQILgD+J1NRUeb3eu3Kt69evq6WlRU6nUw6H42tfLzU19S5UBeBeirAsywp1EXdDd3e34uLi1NXVpdjY2FCXAwAA7sCdfn8z6RYAABiPR0IAjObz+VRbW6u2tjYlJiYqKytLNpst1GUBuMcYYQFgrKqqKqWkpCg7O1srV65Udna2UlJSVFVVFerSANxjjLAAMFJVVZWWL1+uRYsW6Yc//KEcDoeuX7+uX/7yl1q+fLmOHDmipUuXhrpMAPcIk24BGMfn8yklJUXx8fHq6OhQS0uL/zOn06n4+Hh1dnbq/PnzPB4CRrk7/f5mhAWAcWpra9XS0qILFy4oOjo64LPf//73unDhgizLUm1trRYsWBCaIgHcUyOaw1JaWqrk5GRFR0fL7Xartrb2tv137doll8slh8OhmTNn6sCBA7fsW1lZqYiICOXn54+kNABh4NKlS5Iky7L09NNP6+TJk7p27ZpOnjypp59+WoMDw4P9AIS/oEdYDh8+rE2bNqm0tFSZmZnavXu3cnNzdfbsWU2ePHlI/7KyMhUWFmrPnj2aO3euPB6P1q9fr/Hjx2vx4sUBfS9cuKAf/OAHysrKGvkdARj1rly5Ikn65je/qXfeeUdjxtz8t9W8efP0zjvv6IknntCZM2f8/QCEv6BHWEpKSrR27VqtW7dOLpdLO3fuVFJSksrKyobtf/DgQT3zzDMqKCjQ1KlTtWLFCq1du1Y7duwI6Ofz+fSXf/mXeuWVVzR16tSR3Q2AsPDf//3fkqQHHnhg2M8H2wf7AQh/QQWWvr4+eb1e5eTkBLTn5OSorq5u2HN6e3uHPIN2OBzyeDzq7+/3t23btk0TJkzQ2rVr76iW3t5edXd3BxwAwsPgiMrJkyeVn58f8EgoPz9fv/71rwP6AQh/Qf3f3tHRIZ/Pp4SEhID2hISEWw7NLly4UHv37pXX65VlWaqvr1dFRYX6+/vV0dEhSfroo4+0b98+7dmz545rKSoqUlxcnP9ISkoK5lYAGGxwIq3L5dJvfvMbZWRkKDY2VhkZGTpz5ox/LyAm3AL3jxH98yQiIiLgZ8uyhrQN2rp1q3JzczVv3jzZ7XYtWbJEa9askSTZbDZdu3ZNq1at0p49exQfH3/HNRQWFqqrq8t/XLx4cSS3AsBACxYs0IQJE3Tu3DmlpaXpzTff1L59+/Tmm2/qz/7sz9TY2KiJEycSWID7SFCTbuPj42Wz2YaMprS3tw8ZdRnkcDhUUVGh3bt36/e//70SExNVXl6umJgYxcfH6ze/+Y1aWloCJuAODAzcLO4b31BTU5OmTZs25LpRUVGKiooKpnwAo4TNZtNPf/pTLVu2TB988IHee+89/2djx46VdHNCP2uwAPePoEZYIiMj5Xa7VVNTE9BeU1OjjIyM255rt9s1adIk2Ww2VVZWKi8vT2PGjFFqaqrOnDmj06dP+4/vfve7ys7O1unTp3nUA9ynli5dqqNHj2rixIkB7RMnTtTRo0dZ5Ra4zwT9WvPmzZu1evVqzZkzR/Pnz1d5eblaW1u1YcMGSTcf1Vy6dMm/1kpzc7M8Ho/S09N19epVlZSUqKGhQfv375ckRUdHKy0tLeB3PPjgg5I0pB3A/WXp0qVasmQJmx8CCD6wFBQUqLOzU9u2bVNbW5vS0tJUXV2tKVOmSJLa2trU2trq7+/z+VRcXKympibZ7XZlZ2errq5OTqfzrt0EgPBls9mYqwKAvYQAAEDo3On3N4sYAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwXtCvNQPAveTz+ViHBQAjLADMVVVVpWnTpik7O1srV65Udna2pk2bpqqqqlCXBuAeI7AAMFJVVZWWLVum9vb2gPb29nYtW7aM0ALcZ1g4DoBxfD6fEhMT9Yc//EGLFi3Sd77zHTkcDl2/fl3V1dV67733NHHiRF2+fJnHQ8Aod6ff38xhAWCc48eP6w9/+INSU1P129/+NmC3ZqfTqdTUVDU2Nur48eN6+umnQ1gpgHuFwALAOMePH5ckNTU1adGiRfrhD3/oH2H55S9/6Q8wBBbg/kFgAWCcgYEBSdL06dN15swZHTt2zP/ZlClTNH36dDU3N/v7AQh/TLoFYJyHHnpIktTc3KxZs2bp5MmTunbtmk6ePKlZs2apubk5oB+A8EdgAWCciRMn+v9uWdaQY7h+AMIbj4QAGKezs9P/9w8++CBg0u3YsWOH7QcgvDHCAsA4EyZMkCQ9+eSTQ0ZRJk6cqCeffDKgH4DwxwgLAOM89thjkqTTp08PeUvo/fff94+4DPYDEP5YOA6AcXw+n1JSUhQfH6+Ojg61tLT4P0tOTtbDDz+szs5OnT9/noXjgFGOheMAjFo2m03FxcVavny5Fi1apB/84AdDRliOHDlCWAHuIwQWAEZaunSpjhw5ohdeeCFgHZbk5GQdOXJES5cuDWF1AO41HgkBMJrP51Ntba3a2tqUmJiorKwsRlaAMMIjIQBhwWazacGCBaEuA0CI8VozAAAwHoEFAAAYj8ACAACMR2ABAADGY9ItAKPxlhAAiREWAAarqqpSSkqKsrOztXLlSmVnZyslJUVVVVWhLg3APUZgAWCkqqoqLV++XLNmzdLJkyd17do1nTx5UrNmzdLy5csJLcB9hoXjABhncC+hWbNm6ejRo/roo4/8j4QyMzO1bNkyNTQ0sJcQEAbu9PubERYAxqmtrVVLS4syMjI0Y8aMgEdCM2bM0Pz58/XZZ5+ptrY21KUCuEcILACM09bWJkkqLCwc9pHQSy+9FNAPQPgjsAAwzsSJEyVJTz31lI4ePaobN27oF7/4hW7cuKGjR48qMzMzoB+A8MdrzQCM1dHRoRkzZqilpcXf5nQ6FR0dHbqiAIQEIywAjNPe3i5Jamxs1PXr11VeXq7Lly+rvLxc169fV2NjY0A/AOGPERYAxhl81ONyudTT06Pvf//7/s+cTqdSU1PV2NjIIyHgPkJgAWCshx9+WP/5n/855LXm7OzsUJcG4B7jkRAA4ww+6jlx4oSWLVumqKgo5eXlKSoqSsuWLdNHH30U0A9A+COwADBOYmKiJKmoqEhnzpxRRkaGYmNjlZGRoYaGBm3fvj2gH4DwxyMhAMbJysqS0+lUXV2dmpubh13pNjk5WVlZWaEuFcA9MqIRltLSUiUnJys6Olput/uPrja5a9cuuVwuORwOzZw5UwcOHAj4vKqqSnPmzNGDDz6oBx54QE888YQOHjw4ktIAhAGbzabi4mIdO3Zs2EdCx44d049//GOW5QfuJ1aQKisrLbvdbu3Zs8c6e/as9dxzz1kPPPCAdeHChWH7l5aWWjExMVZlZaX1u9/9zvrZz35mjRs3znr33Xf9ff7t3/7Nqqqqss6ePWt98skn1s6dOy2bzWa9//77d1xXV1eXJcnq6uoK9pYAGOro0aOW0+m0JPmP5ORk6+jRo6EuDcBdcqff30Fvfpienq7Zs2errKzM3+ZyuZSfn6+ioqIh/TMyMpSZmanXX3/d37Zp0ybV19frxIkTt/w9s2fP1qJFi/Tqq6/eUV1sfgiEJ5/Pp9raWv8joaysLEZWgDByp9/fQc1h6evrk9fr1YsvvhjQnpOTo7q6umHP6e3tHbIqpcPhkMfjUX9/v+x2e8BnlmXpgw8+UFNTk3bs2HHLWnp7e9Xb2+v/ubu7O5hbATBK2Gw2LViwINRlAAixoOawdHR0yOfzKSEhIaA9ISFBV65cGfachQsXau/evfJ6vbIsS/X19aqoqFB/f786Ojr8/bq6ujRu3DhFRkZq0aJFeuONN/Ttb3/7lrUUFRUpLi7OfyQlJQVzKwAAYBQZ0aTbiIiIgJ8tyxrSNmjr1q3Kzc3VvHnzZLfbtWTJEq1Zs0aSAoZ1Y2JidPr0af3Hf/yHXnvtNW3evFnHjx+/ZQ2FhYXq6uryHxcvXhzJrQAAgFEgqMASHx8vm802ZDSlvb19yKjLIIfDoYqKCvX09KilpUWtra1yOp2KiYlRfHz8/y9kzBilpKToiSee0AsvvKDly5cPOydmUFRUlGJjYwMOAAAQnoIKLJGRkXK73aqpqQlor6mpUUZGxm3PtdvtmjRpkmw2myorK5WXl6cxY2796y3LCpijAgAA7l9BLxy3efNmrV69WnPmzNH8+fNVXl6u1tZWbdiwQdLNRzWXLl3yr7XS3Nwsj8ej9PR0Xb16VSUlJWpoaND+/fv91ywqKtKcOXM0bdo09fX1qbq6WgcOHAh4EwnA/Ym3hABIIwgsBQUF6uzs1LZt29TW1qa0tDRVV1drypQpkqS2tja1trb6+/t8PhUXF6upqUl2u13Z2dmqq6uT0+n09/nyyy/17LPP6r/+67/kcDiUmpqqQ4cOqaCg4OvfIYBRq6qqSi+88IJaWlr8bU6nU8XFxVq6dGnoCgNwzwW9DoupWIcFCC9VVVVavny58vLy9NJLLyktLc2/j9CxY8d05MgRQgsQBu70+5vAAsA4Pp9PKSkpmjVrlt5+++2A+W4DAwPKz89XQ0ODzp8/z+MhYJS70+9vdmsGYJza2lq1tLTopZdeGjI5f8yYMSosLNRnn332R/cxAxA+CCwAjNPW1iZJSktLG/bzwfbBfgDCH4EFgHESExMlSQ0NDcN+Ptg+2A9A+COwADBOVlaWnE6ntm/froGBgYDPBgYGVFRUpOTkZGVlZYWoQgD3GoEFgHFsNpuKi4t17Ngx5efn6+TJk7p27ZpOnjyp/Px8HTt2TD/+8Y+ZcAvcR4JehwUA7oWlS5fqyJEjeuGFFwJW0k5OTuaVZuA+xGvNAIzW19en0tJS/e53v9O0adP07LPPKjIyMtRlAbhL7vT7mxEWAMYabqXbn/zkJ6x0C9yHmMMCwEiDK93OmjUrYA7LrFmztHz5clVVVYW6RAD3EI+EABiHlW6B+wcr3QIYtVjpFsBXEVgAGIeVbgF8FYEFgHFY6RbAVxFYABiHlW4BfBWBBYBxWOkWwFexDgsAI7HSLYD/jdeaARjN5/OptrZWbW1tSkxMVFZWFiMrQBhhpVsAYcFms2nBggWhLgNAiDGHBQAAGI/AAgAAjMcjIQBGYw4LAIkRFgAGq6qqUkpKirKzs7Vy5UplZ2crJSWFjQ+B+xCBBYCR2K0ZwP/Ga80AjMNuzcD9g92aAYxa7NYM4KsILACMw27NAL6KwALAOOzWDOCrCCwAjMNuzQC+isACwDjs1gzgq1g4DoCR2K0ZwP/Ga80AjMZKt0B4Y7dmAGGB3ZoBSMxhAQAAowCBBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwXtisdDu4w0B3d3eIKwEAAHdq8Hv7j+0UFDaB5dq1a5KkpKSkEFcCAACCde3aNcXFxd3y87DZ/HBgYECXL19WTEyMIiIiQl0OgLuou7tbSUlJunjxIpubAmHGsixdu3ZNjz76qMaMufVMlbAJLADCF7uxA2DSLQAAMB6BBQAAGI/AAsB4UVFR+sd//EdFRUWFuhQAIcIcFgAAYDxGWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBYCxPvzwQy1evFiPPvqoIiIi9Pbbb4e6JAAhQmABYKwvv/xSjz/+uN58881QlwIgxMJm80MA4Sc3N1e5ubmhLgOAARhhAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPN4SAmCsL774Qp988on/588++0ynT5/WQw89pMmTJ4ewMgD3Grs1AzDW8ePHlZ2dPaT9r//6r/XWW2/d+4IAhAyBBQAAGI85LAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAY7/8CUdJEexAmRcIAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# redundancy for every measure\n", "# vector will hold same order as the measures in the cnet df\n", "# def compute_measure_redundancy\n", "def compute_redundancy(N, W_observations, J):\n", " Qxx = np.linalg.inv(N)\n", " Qvv = np.linalg.inv(W_observations) - J.dot(Qxx).dot(J.T)\n", " r = np.diagonal(Qvv.dot(W_observations))\n", " \n", " return r\n", "\n", "r = compute_redundancy(N, W_observations, J)\n", "print(f'Minimum redundancy: {min(r)}')\n", "print(f'Maximum redundancy: {max(r)}')\n", "plt.boxplot(r)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Whole bundle process in a loop" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "iteration 0: sigma0 = 2.5728422824613784\n", "\n", "corrections: mean = 0.653346079031689 min = -15.919884218441833 max = 10.327444239483894\n", "iteration 1: sigma0 = 0.9391071877144452\n", "\n", "corrections: mean = 0.007578983174291483 min = -0.6175605679370911 max = 2.362208151398094\n", "iteration 2: sigma0 = 0.9219212955967535\n", "\n", "corrections: mean = 0.00047562127535878755 min = -6.720164033293535e-05 max = 0.34999571717706907\n", "iteration 3: sigma0 = 0.9219206529217796\n", "\n", "corrections: mean = -0.00043949558479830525 min = -0.35037971701721626 max = 0.0013604163271051927\n", "iteration 4: sigma0 = 0.9219208109877487\n", "\n", "corrections: mean = 0.00044131364530114236 min = -5.7892344840921925e-06 max = 0.35036852412427094\n", "iteration 5: sigma0 = 0.9219206528104202\n", "\n", "corrections: mean = -0.00044126496471851365 min = -0.35037040745737597 max = 1.4554057900294883e-05\n", "iteration 6: sigma0 = 0.9219208110024477\n", "\n", "corrections: mean = 0.00044126488421801513 min = -5.783747480588058e-06 max = 0.3503829032363816\n", "iteration 7: sigma0 = 0.921920652818743\n", "\n", "corrections: mean = -0.0004411037049932667 min = -0.3503829421081344 max = 1.2348549615825461e-05\n", "iteration 8: sigma0 = 0.9219208109707007\n", "\n", "corrections: mean = 0.0004409523138932092 min = -0.00011674460187552498 max = 0.35038181894977977\n", "iteration 9: sigma0 = 0.9219206528067627\n", "\n", "corrections: mean = -0.0004411144651418224 min = -0.3503813294832801 max = 6.1212020049696e-06\n", "iteration 10: sigma0 = 0.9219208109789235\n", "\n" ] } ], "source": [ "sensors = generate_sensors(cubes) # generate sensors\n", "cnet = io_controlnetwork.from_isis(network) # load in network\n", "cnet = compute_apriori_ground_points(cnet, sensors) # calculate ground points\n", "\n", "### INPUTS ###\n", "all_parameters = {sn: get_sensor_parameters(sensor) for sn, sensor in sensors.items()} #all parameters\n", "parameters = {sn: parameter[:3] for sn, parameter in all_parameters.items()} #just solving for camera angles and angle velocity\n", "##############\n", "\n", "column_dict = compute_coefficient_columns(cnet, sensors, parameters)\n", "num_parameters = max(col_range[1] for col_range in column_dict.values())\n", "num_observations = 2 * len(cnet)\n", "W_observations = np.eye(num_observations)\n", "W_params = compute_parameter_weights(cnet, sensors, parameters, column_dict)\n", "\n", "iteration = 0\n", "V = compute_residuals(cnet, sensors)\n", "dX = np.zeros(W_params.shape[0]) #initialize for sigma calculatioN\n", "sigma0 = compute_sigma0(V, dX, W_params, W_observations)\n", "print(f'iteration {iteration}: sigma0 = {sigma0}\\n')\n", "\n", "max_iterations = 10\n", "tol = 1e-10\n", "total_correction = np.zeros(num_parameters)\n", "for i in range(max_iterations): \n", " iteration += 1\n", " old_sigma0 = sigma0\n", " \n", " J = compute_jacobian(cnet, sensors, parameters, column_dict) \n", " N = J.T.dot(W_observations).dot(J) + W_params # calculate the normal equation\n", " C = J.T.dot(W_observations).dot(V) - W_params.dot(total_correction)\n", " dX = np.linalg.inv(N).dot(C) #calculate change in camera parameters and ground points\n", " total_correction += dX\n", " print(f'corrections: mean = {dX.mean()} min = {dX.min()} max = {dX.max()}')\n", " \n", " update_parameters(sensors, parameters, cnet, dX, column_dict)\n", " \n", " V = compute_residuals(cnet, sensors)\n", " sigma0 = compute_sigma0(V, dX, W_params, W_observations)\n", " sigma0 = np.sqrt((V.dot(W_observations).dot(V) + dX.dot(W_params).dot(dX))/dof)\n", " print(f'iteration {iteration}: sigma0 = {sigma0}\\n')\n", " \n", " if (abs(sigma0 - old_sigma0) < tol):\n", " print(f'change in sigma0 of {abs(sigma0 - old_sigma0)} converged!')\n", " break\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.18" } }, "nbformat": 4, "nbformat_minor": 4 }