Commit 2c0738a1 authored by AustinSanders's avatar AustinSanders Committed by Jesse Mapel
Browse files

Added getIlluminationDirection for line scanners (#259) 



* Added getIlluminationDirection

* Added getSunVelocities

* Adjusted conditional statements to better meet internal conventions.

* Factored getSunPosition into its own function.

* Changed vector lengths from shorts to ints

Images could have tens of thousands of positions and/or velocities and shorts cannot hold that many, so they need to be ints.

Co-authored-by: default avatarJesse Mapel <jam826@nau.edu>
parent de568c63

include/usgscsm/UsgsAstroLsSensorModel.h

+13 −0
Original line number Diff line number Diff line
@@ -123,6 +123,10 @@ public: 
   std::vector<double> m_covariance;

   int          m_imageFlipFlag;



   std::vector<double> m_sunPosition;

   std::vector<double> m_sunVelocity;





   // Define logging pointer and file content

   std::string m_logFile;

   std::shared_ptr<spdlog::logger> m_logger;
@@ -908,6 +912,15 @@ public: 
       const std::vector<double>& adj,

       double attCorr[9]) const;



   virtual csm::EcefVector getSunPosition(

       const double imageTime) const;

    //> This method returns the position of the sun at the time the image point

    //  was recorded.  If multiple sun positions are available, the method uses

    //  lagrange interpolation.  If one sun position and at least one sun velocity

    //  are available, then the position is calculated using linear extrapolation.

    //  If only one sun position is available, then that value is returned.





private:



   void determineSensorCovarianceInImageSpace(

include/usgscsm/Utilities.h

+1 −0
Original line number Diff line number Diff line
@@ -121,6 +121,7 @@ DistortionType getDistortionModel(nlohmann::json isd, csm::WarningList *list=nul 
std::vector<double> getDistortionCoeffs(nlohmann::json isd, csm::WarningList *list=nullptr);

std::vector<double> getRadialDistortion(nlohmann::json isd, csm::WarningList *list=nullptr);

std::vector<double> getSunPositions(nlohmann::json isd, csm::WarningList *list=nullptr);

std::vector<double> getSunVelocities(nlohmann::json isd, csm::WarningList *list=nullptr);

std::vector<double> getSensorPositions(nlohmann::json isd, csm::WarningList *list=nullptr);

std::vector<double> getSensorVelocities(nlohmann::json isd, csm::WarningList *list=nullptr);

std::vector<double> getSensorOrientations(nlohmann::json isd, csm::WarningList *list=nullptr);

src/UsgsAstroLsSensorModel.cpp

+72 −7
Original line number Diff line number Diff line
@@ -102,6 +102,8 @@ const std::string UsgsAstroLsSensorModel::_STATE_KEYWORD[] = 
   "m_currentParameterValue",

   "m_parameterType",

   "m_referencePointXyz",

   "m_sunPosition",

   "m_sunVelocity",

   "m_gsd",

   "m_flyingHeight",

   "m_halfSwath",
@@ -255,6 +257,8 @@ void UsgsAstroLsSensorModel::replaceModelState(const std::string &stateString ) 
   m_quaternions = j["m_quaternions"].get<std::vector<double>>();

   m_currentParameterValue = j["m_currentParameterValue"].get<std::vector<double>>();

   m_covariance = j["m_covariance"].get<std::vector<double>>();

   m_sunPosition = j["m_sunPosition"].get<std::vector<double>>();

   m_sunVelocity = j["m_sunVelocity"].get<std::vector<double>>();



   m_logFile = j["m_logFile"].get<std::string>();

   if (m_logFile.empty()) {
@@ -443,6 +447,12 @@ std::string UsgsAstroLsSensorModel::getModelState() const { 
                             m_referencePointXyz.x, m_referencePointXyz.y,

                             m_referencePointXyz.z)



      state["m_sunPosition"] = m_sunPosition;

      MESSAGE_LOG(m_logger, "num sun positions: {} ", m_sunPosition.size())



      state["m_sunVelocity"] = m_sunVelocity;

      MESSAGE_LOG(m_logger, "num sun velocities: {} ", m_sunVelocity.size())



      state["m_logFile"] = m_logFile;



      return state.dump();
@@ -518,6 +528,10 @@ void UsgsAstroLsSensorModel::reset() 
  m_referencePointXyz.x = 0.0;

  m_referencePointXyz.y = 0.0;

  m_referencePointXyz.z = 0.0;



  m_sunPosition = std::vector<double>(3, 0.0);

  m_sunVelocity = std::vector<double>(3, 0.0);



  m_gsd = 1.0;

  m_flyingHeight = 1000.0;

  m_halfSwath = 1000.0;
@@ -1668,14 +1682,23 @@ UsgsAstroLsSensorModel::getValidImageRange() const 
csm::EcefVector UsgsAstroLsSensorModel::getIlluminationDirection(

   const csm::EcefCoord& groundPt) const

{

  MESSAGE_LOG(m_logger, "Accessing illimination direction of ground point"

              "{} {} {}."

              "Illimination direction is not supported, throwing exception",

  MESSAGE_LOG(m_logger, "Accessing illumination direction of ground point"

              "{} {} {}.",

              groundPt.x, groundPt.y, groundPt.z);

   throw csm::Error(

      csm::Error::UNSUPPORTED_FUNCTION,

      "Unsupported function",

      "UsgsAstroLsSensorModel::getIlluminationDirection");



  csm::EcefVector sunPosition = getSunPosition(getImageTime(groundToImage(groundPt)));

  csm::EcefVector illuminationDirection = csm::EcefVector(groundPt.x - sunPosition.x,

                                                          groundPt.y - sunPosition.y,

                                                          groundPt.z - sunPosition.z);



  double scale = sqrt(illuminationDirection.x * illuminationDirection.x +

                      illuminationDirection.y * illuminationDirection.y +

                      illuminationDirection.z * illuminationDirection.z);



  illuminationDirection.x /= scale;

  illuminationDirection.y /= scale;

  illuminationDirection.z /= scale;

  return illuminationDirection;

}



//---------------------------------------------------------------------------
@@ -2755,6 +2778,10 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imag 
  MESSAGE_LOG(m_logger, "m_referencePointXyz: {} ",

                        state["m_referencePointXyz"].dump())



  // sun_position and velocity are required for getIlluminationDirection

  state["m_sunPosition"]= getSunPositions(isd, parsingWarnings);

  state["m_sunVelocity"]= getSunVelocities(isd, parsingWarnings);



  // leave these be for now.

  state["m_gsd"] = 1.0;

  state["m_flyingHeight"] = 1000.0;
@@ -2800,6 +2827,7 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imag 
                        state["m_detectorSampleOrigin"].dump(),

                        state["m_detectorLineOrigin"].dump())





  // These are exlusive to LineScanners, leave them here for now.

  try {

    state["m_dtEphem"] = isd.at("dt_ephemeris");
@@ -2942,3 +2970,40 @@ std::shared_ptr<spdlog::logger> UsgsAstroLsSensorModel::getLogger() { 
void UsgsAstroLsSensorModel::setLogger(std::shared_ptr<spdlog::logger> logger) {

  m_logger = logger;

}





csm::EcefVector UsgsAstroLsSensorModel::getSunPosition(

  const double imageTime) const

{



  int numSunPositions = m_sunPosition.size();

  int numSunVelocities = m_sunVelocity.size();

  csm::EcefVector sunPosition = csm::EcefVector();



  // If there are multiple positions, use Lagrange interpolation

  if ((numSunPositions/3) > 1) {

    double sunPos[3];

    double endTime = m_t0Ephem + (m_dtEphem * ((m_numPositions/3)));

    double sun_dtEphem = (endTime - m_t0Ephem) / (numSunPositions/3);

    lagrangeInterp(numSunPositions/3, &m_sunPosition[0], m_t0Ephem, sun_dtEphem,

                   imageTime, 3, 8, sunPos);

    sunPosition.x = sunPos[0];

    sunPosition.y = sunPos[1];

    sunPosition.z = sunPos[2];

  }

  else if ((numSunVelocities/3) >= 1){

    // If there is one position triple with at least one velocity triple

    //  then the illumination direction is calculated via linear extrapolation.

      sunPosition.x = (imageTime * m_sunVelocity[0] + m_sunPosition[0]);

      sunPosition.y = (imageTime * m_sunVelocity[1] + m_sunPosition[1]);

      sunPosition.z = (imageTime * m_sunVelocity[2] + m_sunPosition[2]);

  }

  else {

    // If there is one position triple with no velocity triple, then the

    //  illumination direction is the difference of the original vectors.

      sunPosition.x = m_sunPosition[0];

      sunPosition.y = m_sunPosition[1];

      sunPosition.z = m_sunPosition[2];

  }

  return sunPosition;

}

src/Utilities.cpp

+24 −0
Original line number Diff line number Diff line
@@ -1039,6 +1039,30 @@ std::vector<double> getSunPositions(json isd, csm::WarningList *list) { 
  return positions;

}



std::vector<double> getSunVelocities(json isd, csm::WarningList *list) {

  std::vector<double> velocities;

  try {

    json jayson = isd.at("sun_position");

    json unit = jayson.at("unit");

    for (auto& location : jayson.at("velocities")) {

      velocities.push_back(metric_conversion(location[0].get<double>(), unit.get<std::string>()));

      velocities.push_back(metric_conversion(location[1].get<double>(), unit.get<std::string>()));

      velocities.push_back(metric_conversion(location[2].get<double>(), unit.get<std::string>()));

    }

  }

  catch (...) {

    std::cout<<"Could not parse sun velocity"<<std::endl;

    if (list) {

      list->push_back(

        csm::Warning(

          csm::Warning::DATA_NOT_AVAILABLE,

          "Could not parse the sun velocities.",

          "Utilities::getSunVelocities()"));

    }

  }

  return velocities;

}



std::vector<double> getSensorPositions(json isd, csm::WarningList *list) {

  std::vector<double> positions;

  try {

tests/LineScanCameraTests.cpp

+69 −0
Original line number Diff line number Diff line
@@ -122,6 +122,75 @@ TEST_F(ConstVelocityLineScanSensorModel, calculateAttitudeCorrection) { 
  EXPECT_NEAR(attCorr[8], 0, 1e-8);

}





TEST_F(OrbitalLineScanSensorModel, getIlluminationDirectionStationary) {

  // Get state information, replace sun position / velocity to hit third case:

  //  One position, no velocity.

  std::string state = sensorModel->getModelState();

  json jState = json::parse(state);

  jState["m_sunPosition"] = std::vector<double>{100.0,100.0,100.0};

  jState["m_sunVelocity"] = std::vector<double>{};

  sensorModel->replaceModelState(jState.dump());



  csm::ImageCoord imagePt(8.5,8);

  csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0);

  csm::EcefVector direction = sensorModel->getIlluminationDirection(groundPt);



  // Calculate expected sun direction

  // These are the ground point coordinates minus constant sun positions.

  double expected_x = 999899.680000017;

  double expected_y = -100;

  double expected_z = -899.99991466668735;



  //normalize

  double scale = sqrt((expected_x * expected_x) + (expected_y * expected_y) + (expected_z * expected_z));

  expected_x /= scale;

  expected_y /= scale;

  expected_z /= scale;



  EXPECT_DOUBLE_EQ(direction.x, expected_x);

  EXPECT_DOUBLE_EQ(direction.y, expected_y);

  EXPECT_DOUBLE_EQ(direction.z, expected_z);

}



TEST_F(OrbitalLineScanSensorModel, getSunPositionLagrange){

  std::cout<<sensorModel->m_t0Ephem<<std::endl;

  csm::EcefVector sunPosition = sensorModel->getSunPosition(-.6);

  EXPECT_DOUBLE_EQ(sunPosition.x, 125);

  EXPECT_DOUBLE_EQ(sunPosition.y, 125);

  EXPECT_DOUBLE_EQ(sunPosition.z, 125);

}



TEST_F(OrbitalLineScanSensorModel, getSunPositionLinear){

  // Get state information, replace sun position / velocity to hit third case:

  //  One position, no velocity.

  std::string state = sensorModel->getModelState();

  json jState = json::parse(state);

  jState["m_sunPosition"] = std::vector<double>{100.0,100.0,100.0};

  jState["m_sunVelocity"] = std::vector<double>{50.0, 50.0, 50.0};

  sensorModel->replaceModelState(jState.dump());



  csm::EcefVector sunPosition = sensorModel->getSunPosition(.5);

  EXPECT_DOUBLE_EQ(sunPosition.x, 125);

  EXPECT_DOUBLE_EQ(sunPosition.y, 125);

  EXPECT_DOUBLE_EQ(sunPosition.z, 125);

}



TEST_F(OrbitalLineScanSensorModel, getSunPositionStationary){

  // Get state information, replace sun position / velocity to hit third case:

  //  One position, no velocity.

  std::string state = sensorModel->getModelState();

  json jState = json::parse(state);

  jState["m_sunPosition"] = std::vector<double>{100.0,100.0,100.0};

  jState["m_sunVelocity"] = std::vector<double>{};

  sensorModel->replaceModelState(jState.dump());



  csm::EcefVector sunPosition = sensorModel->getSunPosition(1);

  EXPECT_DOUBLE_EQ(sunPosition.x, 100);

  EXPECT_DOUBLE_EQ(sunPosition.y, 100);

  EXPECT_DOUBLE_EQ(sunPosition.z, 100);

}



TEST_F(OrbitalLineScanSensorModel, Center) {

  csm::ImageCoord imagePt(8.5, 8.0);

  csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0);