Commit e6360313 authored by Jeff Anderson's avatar Jeff Anderson
Browse files

Added Aerial photo software 

git-svn-id: http://subversion.wr.usgs.gov/repos/prog/isis3/trunk@5749 41f8697f-d340-4b68-9986-7bafba869bb8
parent 38e7b45d

isis/src/base/objs/CameraFocalPlaneMap/CameraFocalPlaneMap.cpp

+40 −0
Original line number Diff line number Diff line
@@ -23,7 +23,9 @@ 
#include "CameraFocalPlaneMap.h"



#include <cmath>

#include <QDebug>



#include "Affine.h"

#include "Camera.h"

#include "Spice.h"
@@ -49,6 +51,44 @@ namespace Isis { 
  CameraFocalPlaneMap::~CameraFocalPlaneMap(){

  }



/**

 * Added new method to allow programmer to pass in affine coefficients

 * 

 * @author janderson (3/25/2014)

 * 

 * @param parent Parent camera

 * @param affine Affine transform containing the coefficients for both 

 *               transforms to (samp,line) to (x,y) and reverse

 */

  CameraFocalPlaneMap::CameraFocalPlaneMap(Camera *parent, Affine &affine) {

    p_detectorSampleOrigin = 0.0;

    p_detectorLineOrigin = 0.0;

    p_detectorSampleOffset = 0.0;

    p_detectorLineOffset = 0.0;

    p_camera = parent;



    p_transx[0]  = affine.Coefficients(1)[2];

    p_transx[1]  = affine.Coefficients(1)[0];

    p_transx[2]  = affine.Coefficients(1)[1];



    p_transy[0]  = affine.Coefficients(2)[2];

    p_transy[1]  = affine.Coefficients(2)[0];

    p_transy[2]  = affine.Coefficients(2)[1];



    p_itranss[0]  = affine.InverseCoefficients(1)[2];

    p_itranss[1]  = affine.InverseCoefficients(1)[0];

    p_itranss[2]  = affine.InverseCoefficients(1)[1];



    p_itransl[0]  = affine.InverseCoefficients(2)[2];

    p_itransl[1]  = affine.InverseCoefficients(2)[0];

    p_itransl[2]  = affine.InverseCoefficients(2)[1];



    if (parent != 0) {

      p_camera->SetFocalPlaneMap(this);

    }

  }





  void CameraFocalPlaneMap::Init(Camera *parent, const int naifIkCode) {

    p_detectorSampleOrigin = 0.0;

    p_detectorLineOrigin = 0.0;

isis/src/base/objs/CameraFocalPlaneMap/CameraFocalPlaneMap.h

+7 −0
Original line number Diff line number Diff line
@@ -25,6 +25,7 @@ 




namespace Isis {

  class Affine;

  class Camera;

  /** Convert between distorted focal plane and detector coordinates

   *
@@ -87,12 +88,18 @@ namespace Isis { 
   *                           compliant with Isis standards. Added

   *                           documentation to member

   *                           variables.References #678.

   *   @history 2014-03-27 Jeff Anderson - Added constructor which takes an

   *                           Affine transform instead of only reading from 

   *                           NAIF kernels to obtain the transform.  This is

   *                           required for non-NAIF instruments such as

   *                           Aerial photos.

   *  

   */

  class CameraFocalPlaneMap {

    public:

      CameraFocalPlaneMap(Camera *parent, const int naifIkCode);

      CameraFocalPlaneMap(const int naifIkCode);

      CameraFocalPlaneMap(Camera *parent, Affine &affine);

      virtual ~CameraFocalPlaneMap();



      virtual bool SetDetector(const double sample, const double line);

isis/src/local/apps/apinit/Makefile

0 → 100644
+7 −0
Original line number Diff line number Diff line 
ifeq ($(ISISROOT), $(BLANK))

.SILENT:

error:

	echo "Please set ISISROOT";

else

	include $(ISISROOT)/make/isismake.apps

endif
 
 No newline at end of file

isis/src/local/apps/apinit/apinit.cpp

0 → 100644
+479 −0
Original line number Diff line number Diff line 
#include "Isis.h"



#include <QDate>

#include <QString>



#include <naif/SpiceUsr.h>

#include <naif/SpiceZfc.h>

#include <naif/SpiceZmc.h>



#include "Brick.h"

#include "Cube.h"

#include "FileName.h"

#include "History.h"

#include "IString.h"

#include "Projection.h"

#include "Pvl.h"

#include "PvlGroup.h"

#include "PvlKeyword.h"

#include "PvlObject.h"

#include "SpecialPixel.h"

#include "SpicePosition.h"

#include "SpiceRotation.h"

#include "Table.h"

#include "TProjection.h"



#include <vector>



using namespace Isis;

using namespace std;



PvlKeyword getArrayValues(const QString &parameter, const QString &keyword, int count);





void IsisMain() {

  // Fire up the user interface

  UserInterface &ui = Application::GetUserInterface();



  // Get the name of the cube to initialize

  Cube cube(ui.GetFileName("FROM"),"rw");

  

  // Create archive group with the following information:

  //   Roll Number

  //   Frame Number 

  //   Camera Identification

  //   Lens Identification

  //   Camera calibration date

  PvlGroup archiveGroup("Archive");



  PvlKeyword keyword;

  keyword.setName("RollNumber");

  keyword.setValue(ui.GetString("RollNumber") );

  archiveGroup.addKeyword(keyword);



  keyword.setName("FrameNumber");

  keyword.setValue(ui.GetString("FrameNumber") );

  archiveGroup.addKeyword(keyword);



  keyword.setName("CameraSerialNumber");

  keyword.setValue(ui.GetString("CameraSN") );

  archiveGroup.addKeyword(keyword);



  keyword.setName("LenseSerialNumber");

  keyword.setValue(ui.GetString("LenseSN") );

  archiveGroup.addKeyword(keyword);



  keyword.setName("CalibrationReport");

  keyword.setValue(ui.GetString("CalReport") );

  archiveGroup.addKeyword(keyword);



  cube.putGroup(archiveGroup);



  // Create instrument group with the following information:

  //   Mission/spacecraft name

  //   Target name

  //   Date and time of observation

  //   Ephemeris time

  //   Center longitude/latitude

  //   Altitude above mean radius

  //   Focal length

  //   Fidicual measurements in image

  //   Fidicual locations in focal plane

  //   Optical distortion information

  PvlGroup instrumentGroup("Instrument");



  keyword.setName("SpacecraftName");

  keyword.setValue("Aircraft");

  instrumentGroup.addKeyword(keyword);



  keyword.setName("InstrumentId");

  keyword.setValue("AerialPhoto");

  instrumentGroup.addKeyword(keyword);



  keyword.setName("TargetName");

  keyword.setValue("Earth");

  instrumentGroup.addKeyword(keyword);



  QString startTime = ui.GetString("DateTime");

  // Validate the observation date format.  See naif required reading for time format

  // http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/str2et_c.html

  ConstSpiceChar *spiceStartTime = startTime.toAscii().data();

  SpiceDouble junk;

  SpiceChar errmsg[500];

  tparse_c (spiceStartTime, 500, &junk, errmsg);

  if (errmsg[0] != 0) {

    QString msg = "Invalid date/time format [" + startTime + "].  See NAIF required reading "

                  "for acceptable formats at http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/time.html";

    msg += ".  Naif error message is [" + QString((char*) errmsg) + "]";

                        

    throw IException(IException::User, msg, _FILEINFO_);

  }

  keyword.setName("StartTime");

  keyword.setValue(startTime);

  instrumentGroup.addKeyword(keyword);



  FileName lsk("$base/kernels/lsk/naif0010.tls");

  furnsh_c(lsk.expanded().toAscii().constData());



  double et;

  spiceStartTime = startTime.toAscii().data();

  str2et_c(spiceStartTime, &et);



  keyword.setName("EphemerisTime");

  keyword.setValue(toString(et),"seconds");

  instrumentGroup.addKeyword(keyword);



  keyword.setName("EstimatedAltitude");

  keyword.setValue(toString(ui.GetDouble("Altitude")), "meters");

  instrumentGroup.addKeyword(keyword);



  keyword.setName("EstimatedCenterLatitude");

  keyword.setValue(toString(ui.GetDouble("CenterLat")), "degrees");

  instrumentGroup.addKeyword(keyword);



  keyword.setName("EstimatedCenterLongitude");

  keyword.setValue(toString(ui.GetDouble("CenterLon")), "degrees");

  instrumentGroup.addKeyword(keyword);



  keyword.setName("FocalLength");

  keyword.setValue(toString(ui.GetDouble("FocalLength")), "millimeters");

  instrumentGroup.addKeyword(keyword);



  keyword = getArrayValues("FiducialX", "FiducialX", 8);

  keyword.setUnits("millimeters");

  instrumentGroup.addKeyword(keyword);



  keyword = getArrayValues("FiducialY", "FiducialY", 8);

  keyword.setUnits("millimeters");

  instrumentGroup.addKeyword(keyword);



  keyword = getArrayValues("FiducialSample", "FiducialSample", 8);

  instrumentGroup.addKeyword(keyword);



  keyword = getArrayValues("FiducialLine", "FiducialLine", 8);

  instrumentGroup.addKeyword(keyword);



  keyword = getArrayValues("RadialDistortionCoefficients", "KCoefs", 5);

  instrumentGroup.addKeyword(keyword);



  keyword = getArrayValues("DecenteringDistortionCoefficients", "PCoefs", 4);

  instrumentGroup.addKeyword(keyword);



  keyword.setName("XCalibratedPrincipalPoint");

  keyword.setValue(toString(ui.GetDouble("XP")), "millimeters");

  instrumentGroup.addKeyword(keyword);



  keyword.setName("YCalibratedPrincipalPoint");

  keyword.setValue(toString(ui.GetDouble("YP")), "millimeters");

  instrumentGroup.addKeyword(keyword);



  keyword.setName("XIndicatedPrincipalPoint");

  keyword.setValue(toString(ui.GetDouble("XIPP")), "millimeters");

  instrumentGroup.addKeyword(keyword);



  keyword.setName("YIndicatedPrincipalPoint");

  keyword.setValue(toString(ui.GetDouble("YIPP")), "millimeters");

  instrumentGroup.addKeyword(keyword);



  cube.putGroup(instrumentGroup);



  // Create kernels group

  PvlGroup kernelsGroup("Kernels");



  keyword.setName("NaifFrameCode");

  keyword.setValue("-2000001");

  kernelsGroup.addKeyword(keyword);



  keyword.setName("LeapSecond");

  keyword.setValue("$base/kernels/lsk/naif0010.tls");

  kernelsGroup.addKeyword(keyword);



  keyword.setName("TargetAttitudeShape");

  keyword.setValue("$base/kernels/pck/pck00009.tpc");

  kernelsGroup.addKeyword(keyword);



  keyword.setName("TargetPosition");

  keyword.setValue("Table");

  keyword.addValue("$base/kernels/spk/de405.bsp");

  kernelsGroup.addKeyword(keyword);



  keyword.setName("InstrumentPointing");

  keyword.setValue("Table");

  kernelsGroup.addKeyword(keyword);



  keyword.setName("Instrument");

  keyword.clear();

  kernelsGroup.addKeyword(keyword);



  keyword.setName("SpacecraftClock");

  keyword.clear();

  kernelsGroup.addKeyword(keyword);



  keyword.setName("InstrumentPosition");

  keyword.setValue("Table");

  kernelsGroup.addKeyword(keyword);



  keyword.setName("InstrumentAddendum");

  keyword.clear();

  kernelsGroup.addKeyword(keyword);



  keyword.setName("ShapeModel");

  keyword.setValue(ui.GetAsString("ShapeModel"));

  kernelsGroup.addKeyword(keyword);



  keyword.setName("InstrumentPositionQuality");

  keyword.setValue("Predict");

  kernelsGroup.addKeyword(keyword);

  

  keyword.setName("InstrumentPointingQuality");

  keyword.setValue("Predict");

  kernelsGroup.addKeyword(keyword);



  keyword.setName("CameraVersion");

  keyword.setValue("1");

  kernelsGroup.addKeyword(keyword);



  cube.putGroup(kernelsGroup);



// Write out the naif keywords group.  This must be done because the camera model

// needs to know the earth's frame number and radii.  This is typically done by the spiceinit program

// but there are no SPICE kernels for aerial photos other than those that define the

// rotation matrix from j2000 to body-fixed for the earth.   Our group should look something

// like this in the labels

// 

// Object = NaifKeywords

//  BODY_FRAME_CODE                    = 10013

//  BODY399_RADII                      = (6378.14, 6378.14, 6356.75)

// End_Object



  PvlObject naifKeywords;

  naifKeywords.setName("NaifKeywords");

  

  keyword.setName("BODY_FRAME_CODE");

  keyword.setValue("10013");  

  naifKeywords.addKeyword(keyword);

    

  keyword.setName("BODY399_RADII");

  keyword.setValue("6378.14");  

  keyword.addValue("6378.14");  

  keyword.addValue("6356.75");

  keyword.setUnits("kilometers"); 

  naifKeywords.addKeyword(keyword);



  cube.label()->addObject(naifKeywords);



  // Create table with rotation matrices from bodyfixed to J2000 and write to the labels.  

  // Need to load the appropriatespk and pck kernels for this to work.  

  FileName pck("$base/kernels/pck/pck00009.tpc");

  furnsh_c(pck.expanded().toAscii().constData());

  FileName tspk("$base/kernels/spk/de405.bsp");

  furnsh_c(tspk.expanded().toAscii().constData());



  SpiceRotation bodyRotation(10013);

  bodyRotation.LoadCache(et-1.0, et+1.0, 2); // Get a small window about the observation time



  Table bodyTable = bodyRotation.Cache("BodyRotation");

  bodyTable.Label() += PvlKeyword("Description", "Created by apinit");

  bodyTable.Label() += PvlKeyword("Kernels");

  bodyTable.Label()["Kernels"].addValue("$base/kernels/spk/de405.bsp");

  bodyTable.Label()["Kernels"].addValue("$base/kernels/pck/pck00009.tpc");

  bodyTable.Label() += PvlKeyword("SolarLongitude","-9999"); // TODO: Need to compute the solar longitude



  cube.write(bodyTable);





  // Create table with position of sun relative to earth.  

  SpicePosition sunPosition(10, 399);

  sunPosition.LoadCache(et-1.0, et+1.0, 2); // Get a small window about the observation time



  Table sunTable = sunPosition.Cache("SunPosition");

  sunTable.Label() += PvlKeyword("Description", "Created by apinit");

  sunTable.Label() += PvlKeyword("Kernels");

  bodyTable.Label()["Kernels"].addValue("$base/kernels/spk/de405.bsp");



  cube.write(sunTable);



  // Create a table with apriori position of aircraft relative to earth in J2000.

  // This will take several steps.   First we need the radius for the aerial photo.   

  // The user supplied a lat/lon and DTM.   We will use the lat/lon to get a radius from

  // the DTM.  Those three values will be converted from polar coordinates (lat,lon,rad) to

  // P(x,y,z).  Then we will get a unit vector of P and multiply it by the

  // altitude to get Pa = (xhat*alt, yhat*alt, zhat*alt).   Then we will add P and Pa to 

  // estimate the aircraft position in body fixed.  Finally we will rotation Pa into J2000 

  // using the bodyRotation object (SpiceRotation class).  The velocity is unknown so we will not

  // write it out to our table.  We could improve the estimate of the spacecraft position by 

  // computing the surface normal at P off of an ellipsoid.  We won't go to this level of accuracy

  // because using jigsaw will be a requirement to improvement the pointing of the camera and the 

  // aircraft position should also be refined at that time.



  Cube dtm;

  dtm.open(ui.GetAsString("ShapeModel"));

  TProjection *mapproj = (TProjection *) dtm.projection();

  mapproj->SetGround(ui.GetDouble("CenterLat"), ui.GetDouble("CenterLon"));

  int line = mapproj->WorldY() + 0.5;

  int samp = mapproj->WorldX() + 0.5;

  Brick buf(1,1,1,dtm.pixelType());  

  buf.SetBasePosition(samp, line, 1);

  dtm.read(buf);



  if (IsSpecial(buf[0]) || buf[0] < 6300000) {

    QString msg = "DTM [" + ui.GetAsString("ShapeModel") + "] does not contain valid radius "

                  " at user specified latitude/longitude.  Read [" + toString(buf[0]) + "] at "

                  "sample/line [" + toString(samp) + "," + toString(line) + "]";

    throw IException(IException::Unknown, msg, _FILEINFO_);

  }



  SpiceDouble radius = buf[0] / 1000.0;

  SpiceDouble ulat = mapproj->UniversalLatitude() * Isis::PI / 180.0;

  SpiceDouble ulon = mapproj->UniversalLongitude() * Isis::PI / 180.0;

    

  SpiceDouble xyz[3];

  latrec_c(radius, ulon, ulat, xyz);

  

  SpiceDouble xyzhat[3];

  vhat_c(xyz, xyzhat);



  std::vector<double> apposB;  // Position of aircraft in body fixed

  apposB.push_back(xyz[0] + xyzhat[0] * ui.GetDouble("Altitude") / 1000.0); 

  apposB.push_back(xyz[1] + xyzhat[1] * ui.GetDouble("Altitude") / 1000.0); 

  apposB.push_back(xyz[2] + xyzhat[2] * ui.GetDouble("Altitude") / 1000.0); 



  bodyRotation.SetEphemerisTime(et);      

  std::vector<double> apposJ = bodyRotation.J2000Vector(apposB);



  TableRecord spkRecord;

  TableField x("J2000X", TableField::Double);

  TableField y("J2000Y", TableField::Double);

  TableField z("J2000Z", TableField::Double);

  TableField t("ET", TableField::Double);

  spkRecord += x;

  spkRecord += y;

  spkRecord += z;

  spkRecord += t;



  Table spkTable("InstrumentPosition", spkRecord);

  spkRecord[0] = apposJ[0];

  spkRecord[1] = apposJ[1];

  spkRecord[2] = apposJ[2];

  spkRecord[3] = et - 1.0;

  spkTable += spkRecord;



  spkRecord[3] = et + 1.0;

  spkTable += spkRecord;



  spkTable.Label() += PvlKeyword("Description", "Created by apinit");

  spkTable.Label() += PvlKeyword("Kernels");



  cube.write(spkTable);

  dtm.close();



  // Now we must write out the camera to J2000 rotation matrices.  Unfortunately

  // nothing exists that defines this information.   However, we can create the camera

  // pointing using photogramettric techniques by finding tie points between an aerial 

  // photo and a basemap that covers the same area of the photo.  

  // A basemap example is a DOQQ image that is map projected.  We won't do that work 

  // here since we already have programs that can help determine the tie points and pointing 

  // (findfeatures and jigsaw).  Instead we will just write out the identity rotation

  // such that we assume the aircraft is pointing straight down (roll and pitch of 0) which is

  // probably correct.   What we really lack is the yaw/twist to know if the aircraft is

  // flying north-south, south-north, or along some other track.   We will just assume 

  // the aircraft has no twist too.  It will be wrong but findfeatures and jigsaw will

  // figure it out.  (CJ is the J2000 to camera rotation matrix.  We'll set it to the 

  // identity and then convert to a quaternion so we can write out our table file.



  SpiceDouble CJ[9] = { 1.0, 0.0, 0.0,

                        0.0, 1.0, 0.0,

                        0.0, 0.0, 1.0 };

  SpiceDouble minusApposJ[3];

  //minusApposJ[0] = -apposJ[0];

  //minusApposJ[0] = -apposJ[0];

  //minusApposJ[1] = -apposJ[1];

  minusApposJ[2] = apposJ[2];

  minusApposJ[1] = apposJ[1];

  minusApposJ[2] = apposJ[2];



  SpiceDouble perpVec[3];

  SpiceDouble axis[3] = {0.0, 1.0, 0.0};

  vrotv_c(minusApposJ, axis, 90.0, perpVec);

  

  twovec_c (minusApposJ, 3, perpVec, 2, (SpiceDouble( *)[3]) &CJ[0]);



  SpiceDouble quat[4];

  m2q_c(CJ, quat);



  TableField q0("J2000Q0", TableField::Double);

  TableField q1("J2000Q1", TableField::Double);

  TableField q2("J2000Q2", TableField::Double);

  TableField q3("J2000Q3", TableField::Double);

  TableField qt("ET", TableField::Double);



  TableRecord ckRecord;

  ckRecord += q0;

  ckRecord += q1;

  ckRecord += q2;

  ckRecord += q3;

  ckRecord += qt;



  Table ckTable("InstrumentPointing", ckRecord);

  

  ckRecord[0] = quat[0];

  ckRecord[1] = quat[1];

  ckRecord[2] = quat[2];

  ckRecord[3] = quat[3];

  ckRecord[4] = et - 1.0;

  ckTable += ckRecord;



  ckTable.Label() += PvlKeyword("Description", "Created by apinit");

  ckTable.Label() += PvlKeyword("Kernels");



  cube.write(ckTable);  





  // Write out the history to the cube so we have a record of what the user ran

  History hist = History("IsisCube");

  try {

    cube.read(hist);

  }

  catch(IException &e) {

    // No history exist but don't throw an error. just add a history object

  }

  hist.AddEntry();

  cube.write(hist);



  // clean up

  cube.close();

}







PvlKeyword getArrayValues(const QString &parameter, const QString &uiKeyword, int count) {



  UserInterface &ui = Application::GetUserInterface();

  QString value = ui.GetString(uiKeyword);



  PvlKeyword keyword;

  keyword.setName(parameter);



  QStringList tokens = value.split(",");

  if (tokens.size() != count) {

    QString msg = "Invalid value for [" + uiKeyword.toUpper() + "].  Expecting exactly " +

                  QString::number(count) + " comma "

                  "separated values and found only [" + QString::number(tokens.size()) + "]";

    throw IException(IException::User, msg, _FILEINFO_);

  }



  for (int i=0; i<count; i++) {

    try {

      toDouble(tokens[i].trimmed()); 

    } 

    catch (IException &e) {

      QString msg = "Invalid value for [" + uiKeyword.toUpper() + "] at index [" + QString::number(i+1) + 

                    "].  Found [" + tokens[i] + "]";

      throw IException(e, IException::User, msg, _FILEINFO_);

    }

    keyword.addValue(tokens[i].trimmed());

  }



  return keyword;

}

isis/src/local/apps/apinit/apinit.xml

0 → 100644
+245 −0

File added.

Preview size limit exceeded, changes collapsed.