RPC: Have rpc session working as valid option for stereo

RPC is largely going to have the same preprocessing filters as
Pinhole. We'll edit as we go.

I also made all the constants of the RPCModel go behind a method
call. This is to insure that the programmer is not modifying these
constants.

src/asp/Sessions/RPC/RPCModel.cc

@@ -14,6 +14,12 @@ using namespace asp;
 using namespace vw;
 
 void asp::RPCModel::initialize( DiskImageResourceGDAL* resource ) {
+  // Extract Datum (by means of GeoReference)
+  GeoReference georef;
+  read_georeference( georef, *resource );
+  m_datum = georef.datum();
+
+  // Extract RPC Info
   boost::shared_ptr<GDALDataset> dataset = resource->get_dataset_ptr();
   if ( !dataset )
     vw_throw( LogicErr() << "RPCModel: Could not read data. No file has been opened." );
@@ -24,19 +30,19 @@ void asp::RPCModel::initialize( DiskImageResourceGDAL* resource ) {
     vw_throw( LogicErr() << "RPCModel: GDAL resource appears not to have RPC metadata." );
 
   // Copy information over to our data structures.
-  lonlatheight_offset = Vector3(gdal_rpc.dfLONG_OFF,
-                                gdal_rpc.dfLAT_OFF,
-                                gdal_rpc.dfHEIGHT_OFF);
-  lonlatheight_scale = Vector3(gdal_rpc.dfLONG_SCALE,
-                               gdal_rpc.dfLAT_SCALE,
-                               gdal_rpc.dfHEIGHT_SCALE);
-  line_os = Vector2(gdal_rpc.dfLINE_OFF, gdal_rpc.dfLINE_SCALE);
-  sample_os = Vector2(gdal_rpc.dfSAMP_OFF, gdal_rpc.dfSAMP_SCALE);
-
-  m_line_num_coeff = Vector<double,20>(gdal_rpc.adfLINE_NUM_COEFF);
-  m_line_den_coeff = Vector<double,20>(gdal_rpc.adfLINE_DEN_COEFF);
-  m_sample_num_coeff = Vector<double,20>(gdal_rpc.adfSAMP_NUM_COEFF);
-  m_sample_den_coeff = Vector<double,20>(gdal_rpc.adfSAMP_DEN_COEFF);
+  m_lonlatheight_offset = Vector3(gdal_rpc.dfLONG_OFF,
+                                  gdal_rpc.dfLAT_OFF,
+                                  gdal_rpc.dfHEIGHT_OFF);
+  m_lonlatheight_scale = Vector3(gdal_rpc.dfLONG_SCALE,
+                                 gdal_rpc.dfLAT_SCALE,
+                                 gdal_rpc.dfHEIGHT_SCALE);
+  m_xy_offset = Vector2(gdal_rpc.dfSAMP_OFF,   gdal_rpc.dfLINE_OFF);
+  m_xy_scale  = Vector2(gdal_rpc.dfSAMP_SCALE, gdal_rpc.dfLINE_SCALE);
+
+  m_line_num_coeff =   CoeffVec(gdal_rpc.adfLINE_NUM_COEFF);
+  m_line_den_coeff =   CoeffVec(gdal_rpc.adfLINE_DEN_COEFF);
+  m_sample_num_coeff = CoeffVec(gdal_rpc.adfSAMP_NUM_COEFF);
+  m_sample_den_coeff = CoeffVec(gdal_rpc.adfSAMP_DEN_COEFF);
   m_lonlat_bbox = BBox2(Vector2(gdal_rpc.dfMIN_LONG, gdal_rpc.dfMIN_LAT),
                         Vector2(gdal_rpc.dfMAX_LONG, gdal_rpc.dfMAX_LAT));
 }
@@ -58,7 +64,7 @@ Vector2 asp::RPCModel::point_to_pixel( Vector3 const& point ) const {
 
   Vector3 lonlatrad = cartography::xyz_to_lon_lat_radius( point );
 
-  Vector<double,20> term =
+  CoeffVec term =
     calculate_terms( elem_quot(lonlatrad - lonlatheight_offset,
                                lonlatheight_scale) );
 

src/asp/Sessions/RPC/RPCModel.h

@@ -10,13 +10,28 @@
 #include <vw/FileIO/DiskImageResourceGDAL.h>
 #include <vw/Camera/CameraModel.h>
 
+namespace vw {
+  namespace cartography {
+    class Datum;
+  }
+}
+
 namespace asp {
 
   // TODO: Need to work out a different way to triangulate. Our
   //   standard midpoint method doesn't seem to work.
   // * Pull Datum from file. LLH measurement are geodetic.
 
   class RPCModel : public vw::camera::CameraModel {
+    Datum m_datum;
+
+    // Scaling parameters
+    vw::Vector<double,20> m_line_num_coeff, m_line_den_coeff,
+      m_sample_num_coeff, m_sample_den_coeff;
+    vw::Vector2 m_xy_offset;
+    vw::Vector2 m_xy_scale;
+    vw::Vector3 m_lonlatheight_offset;
+    vw::Vector3 m_lonlatheight_scale;
     vw::BBox2 m_lonlat_bbox;
 
     void initialize( vw::DiskImageResourceGDAL* resource );
@@ -37,9 +52,21 @@ namespace asp {
     void inverse_transform( vw::Vector2 const& pix, vw::Vector3& point,
                             vw::Vector3& direction ) const;
 
-    static vw::Vector<double,20>
-    calculate_terms( vw::Vector3 const& v ) {
-      vw::Vector<double,20> result;
+    // Access to constants
+    typedef vw::Vector<double,20> CoeffVec;
+    Datum const& datum() const { return m_datum; }
+    CoeffVec const& line_num_coeff() const   { return m_line_num_coeff; }
+    CoeffVec const& line_den_coeff() const   { return m_line_den_coeff; }
+    CoeffVec const& sample_num_coeff() const { return m_sample_num_coeff; }
+    CoeffVec const& sample_den_coeff() const { return m_sample_den_coeff; }
+    vw::Vector2 const& xy_offset() const     { return m_xy_offset; }
+    vw::Vector2 const& xy_scale() const      { return m_xy_scale; }
+    vw::Vector3 const& lonlatheight_offset() const { return m_lonlatheight_offset; }
+    vw::Vector3 const& lonlatheight_scale() const  { return m_lonlatheight_scale; }
+    vw::BBox2 const& lonlat_bbox() const     { return m_lonlat_bbox; }
+
+    static CoeffVec calculate_terms( vw::Vector3 const& v ) {
+      CoeffVec result;
       result[0] = 1.0;
       result[1] = v.x();
       result[2] = v.y();
@@ -64,13 +91,6 @@ namespace asp {
       return result;
     }
 
-    // Scaling parameters
-    vw::Vector<double,20> m_line_num_coeff, m_line_den_coeff,
-      m_sample_num_coeff, m_sample_den_coeff;
-    vw::Vector2 line_os; // [OFFSET, SCALE]
-    vw::Vector2 sample_os;
-    vw::Vector3 lonlatheight_offset;
-    vw::Vector3 lonlatheight_scale;
   };
 
 }

src/asp/Sessions/RPC/StereoSessionRPC.h

@@ -10,11 +10,11 @@
 #ifndef __STEREO_SESSION_RPC_H__
 #define __STEREO_SESSION_RPC_H__
 
-#include <asp/Sessions/StereoSession.h>
+#include <asp/Sessions/Pinhole/StereoSessionPinhole.h>
 
 namespace asp {
 
-  class StereoSessionRPC : public StereoSession {
+  class StereoSessionRPC : public StereoSessionPinhole {
 
   public:
     virtual ~StereoSessionRPC() {}

src/asp/Tools/stereo.cc

@@ -97,11 +97,18 @@ void vw::handle_arguments( int argc, char *argv[], Options& opt ) {
   // Some specialization here so that the user doesn't need to list
   // camera models on the command line for certain stereo session
   // types.  (e.g. isis).
+  //
+  // TODO: This modification of arguments should probably happen in
+  // initialization and not be dependent on Stereo knowing what
+  // session it is in.
   bool check_for_camera_models = true;
-  if ( opt.stereo_session_string == "isis" ) {
+  if ( opt.stereo_session_string == "isis" ||
+       opt.stereo_session_string == "rpc" ) {
     // Fix the ordering of the arguments if the user only supplies 3
-    if (opt.out_prefix.empty())
+    if (opt.out_prefix.empty()) {
       opt.out_prefix = opt.cam_file1;
+      opt.cam_file1.clear();
+    }
     check_for_camera_models = false;
   }
 

src/asp/Tools/stereo.in

@@ -96,7 +96,10 @@ if __name__ == '__main__':
         if ( opt.entry_point <= 3 ):
             run('stereo_fltr', args, msg='3: Filtering')
         if ( opt.entry_point <= 4 ):
-            run('stereo_tri', args,  msg='4: Triangulation')
+            tri = 'stereo_tri'
+            if ( opt.session.lowercase == "rpc" ):
+                tri += '_rpc'
+            run(tri,           args, msg='4: Triangulation')
     except Exception, e:
         if not opt.debug:
             die(e)

src/asp/Tools/stereo_tri_rpc.cc

@@ -45,6 +45,7 @@ int main( int argc, char* argv[] ) {
     // Triangulation with RPC is a little tricky. Here's the plan.
     // 1.) Solve for the mean XYZ for both cameras by looking at the
     // LLH boxes.
+
     // 2.) Solve for one ray's intersection at the mean height. This
     // will seed the rigoruos algorithm. I believe this is equation 24
     // in Grodecki 2004. We only need the first iteration's result.