Adding bundle logic helper functions (#73)

* added to bundle adjust nb, up to setting up bundle iteration, known issues

* adding compute_residuals + test and generate_sensors

* Updated some docs

Co-authored-by: Jesse Mapel <jam826@nau.edu>

knoten/bundle.py

@@ -3,7 +3,51 @@ import pandas as pd
 import pvl
 import os
 import csmapi
+
+from pysis import isis
+from ale.drivers import loads
 from collections import OrderedDict
+from knoten.csm import create_csm
+
+def generate_sensors(cubes, directory=None, clean=False):
+    """
+    Generate a set of USGSCSM sensor models from a list of ISIS cube files
+
+    Parameters
+    ----------
+    cubes     : str
+              Directory/filename of a file containing ISIS cube file paths
+    directory : str
+              Output directory to save resulting json files. Defaults to the
+              same directory as cube list file
+    clean     : flag
+              Option to delete json file outputs
+
+    Returns
+    -------
+    sensors   : dictionary
+              Dictionary mapping ISIS serial numbers to USGSCSM sensor models
+    """
+    if directory is None:
+        directory = os.path.dirname(cubes)
+
+    isd_files = []
+    sensors = {}
+    for line in open(cubes):
+        basename = os.path.splitext(os.path.basename(line.strip()))[0]
+        isd = os.path.join(directory, basename+'.json')
+        isd_files.append(isd)
+        with open(isd, 'w+') as f:
+            f.write(loads(line.strip(), formatter='usgscsm'))
+
+        sn = isis.getsn(from_=line.strip()).strip().decode('utf-8')
+        sensors[sn] = create_csm(isd)
+
+    if clean:
+        for isd in isd_files:
+            os.remove(isd)
+
+    return sensors
 
 def closest_approach(points, direction):
     """
@@ -203,4 +247,36 @@ def compute_jacobian(network, sensors, parameters):
         jacobian[2*i : 2*i+2, image_column : image_column+len(params)] = compute_sensor_partials(sensor, params, ground_pt)
         jacobian[2*i : 2*i+2, point_column : point_column+3] = compute_ground_partials(sensor, ground_pt)
 
-    return jacobian
+    return jacobian, coefficient_columns
+
+def compute_residuals(network, sensors):
+    """
+    Compute the error in the observations by taking the difference between the
+    ground points groundToImage projections and measure values.
+
+    Parameters
+    ----------
+    network : DataFrame
+              The control network as a dataframe generated by plio
+    sensors : dict
+             A dictionary that maps ISIS serial numbers to CSM sensors
+
+    Returns
+    -------
+    V : np.array
+       The control network dataframe with updated ground points
+    """
+    num_meas = len(network)
+    V = np.zeros((num_meas, 2))
+
+    for i in range(num_meas):
+        row = network.iloc[i]
+        serial = row["serialnumber"]
+        ground_pt = row[["adjustedX", "adjustedY", "adjustedZ"]].values
+        ground_pt = csmapi.EcefCoord(ground_pt[0], ground_pt[1], ground_pt[2])
+        sensor = sensors[serial]
+        img_coord = sensor.groundToImage(ground_pt)
+        V[i,:] = [row['line'] - img_coord.line, row['sample'] - img_coord.samp]
+
+    V = V.reshape(num_meas*2)
+    return V

tests/test_bundle.py

@@ -119,7 +119,7 @@ def test_compute_jacobian(control_network, sensors):
     ground_partials = [-(i+1) * np.ones((2, 3)) for i in range(9)]
     with mock.patch('knoten.bundle.compute_sensor_partials', side_effect=sensor_partials) as sensor_par_mock, \
          mock.patch('knoten.bundle.compute_ground_partials', side_effect=ground_partials) as ground_par_mock:
-        J = bundle.compute_jacobian(control_network, sensors, parameters)
+        J, coefficient_columns = bundle.compute_jacobian(control_network, sensors, parameters)
 
     expected_J = [
     [1, 1, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, 0, 0, 0, 0, 0],
@@ -141,3 +141,17 @@ def test_compute_jacobian(control_network, sensors):
     [0, 0, 0, 0, 0, 0, 9, 9, 0, 0, 0, 0, 0, 0, -9, -9, -9],
     [0, 0, 0, 0, 0, 0, 9, 9, 0, 0, 0, 0, 0, 0, -9, -9, -9]]
     np.testing.assert_array_equal(J, expected_J)
+
+def test_compute_residuals(control_network, sensors):
+    # sensor.groundToImage.side_effect = [csmapi.ImageCoord(-0.1, 7.8), csmapi.ImageCoord(0.7, 6.6), csmapi.ImageCoord(1.5, 5.4),
+    #                                     csmapi.ImageCoord(2.3, 4.2), csmapi.ImageCoord(3.1, 4.9), csmapi.ImageCoord(5.8, 3.7),
+    #                                     csmapi.ImageCoord(6.6, 2.5), csmapi.ImageCoord(7.4, 1.3), csmapi.ImageCoord(8.2, 0.1)]
+
+    sensors['a'].groundToImage.side_effect = [csmapi.ImageCoord(-0.1, 7.8), csmapi.ImageCoord(5.8, 3.7)]
+    sensors['b'].groundToImage.side_effect = [csmapi.ImageCoord(0.7, 6.6), csmapi.ImageCoord(3.1, 4.9), csmapi.ImageCoord(6.6, 2.5)]
+    sensors['c'].groundToImage.side_effect = [csmapi.ImageCoord(1.5, 5.4), csmapi.ImageCoord(7.4, 1.3)]
+    sensors['d'].groundToImage.side_effect = [csmapi.ImageCoord(2.3, 4.2), csmapi.ImageCoord(8.2, 0.1)]
+
+    V = bundle.compute_residuals(control_network, sensors)
+    assert V.shape == (18,)
+    np.testing.assert_allclose(V, [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1])