diff --git a/CMakeLists.txt b/CMakeLists.txt index c8ebdffa0be2827abc5a497f82b1a071389f7d2a..58f2ab12c649f8202f2486524f7760c19f90681c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,8 +26,8 @@ if(USGSCSM_EXTERNAL_DEPS) PATH_SUFFIXES "csm" PATHS $ENV{CONDA_PREFIX}/include/) find_library(CSM_LIBRARY csmapi PATHS $ENV{CONDA_PREFIX}/lib) - message("--Found external CSM Library: ${CSM_LIBRARY}") - message("--Found external CSM Include Directory: ${CSM_INCLUDE_DIR}") + message("-- Found external CSM Library: ${CSM_LIBRARY}") + message("-- Found external CSM Include Directory: ${CSM_INCLUDE_DIR}") # Nlohmann JSON find_package(nlohmann_json REQUIRED) @@ -61,11 +61,16 @@ else() set(EIGEN3_INCLUDE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/ale/eigen) endif(USGSCSM_EXTERNAL_DEPS) +# Proj +find_package(PROJ REQUIRED CONFIG) +set(PROJ_TARGET PROJ::proj) + add_library(usgscsm SHARED src/UsgsAstroPlugin.cpp src/UsgsAstroFrameSensorModel.cpp src/UsgsAstroPushFrameSensorModel.cpp src/UsgsAstroLsSensorModel.cpp + src/UsgsAstroProjectedLsSensorModel.cpp src/UsgsAstroSarSensorModel.cpp src/Distortion.cpp src/Utilities.cpp @@ -75,11 +80,11 @@ set_target_properties(usgscsm PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION 1 ) - set(USGSCSM_INCLUDE_DIRS "${CMAKE_CURRENT_SOURCE_DIR}/include/usgscsm" "${CMAKE_CURRENT_SOURCE_DIR}/include" "${CMAKE_CURRENT_SOURCE_DIR}" - "${EIGEN3_INCLUDE_DIR}") + "${EIGEN3_INCLUDE_DIR}" + "${PROJ_INCLUDE_DIR}") # These will be copied upon installation to ${CMAKE_INSTALL_INCLUDEDIR}/include set(USGSCSM_INSTALL_INCLUDE_DIRS "${CMAKE_CURRENT_SOURCE_DIR}/include/usgscsm" @@ -94,6 +99,7 @@ target_include_directories(usgscsm target_link_libraries(usgscsm ${CSM_LIBRARY} ${ALE_TARGET} + ${PROJ_TARGET} nlohmann_json::nlohmann_json) add_executable(usgscsm_cam_test bin/usgscsm_cam_test.cc) diff --git a/environment.yml b/environment.yml index 38195c3da0baa9bc173aef32eeb44ab21b4e3fe7..159b2f1d8b13ae26cd50258af81dd5b501fed62d 100644 --- a/environment.yml +++ b/environment.yml @@ -4,8 +4,9 @@ channels: - default dependencies: - - cmake>=3.15 - ale>=0.8.8 + - cmake>=3.15 - csm - - nlohmann_json - eigen + - nlohmann_json + - proj diff --git a/include/usgscsm/UsgsAstroProjectedLsSensorModel.h b/include/usgscsm/UsgsAstroProjectedLsSensorModel.h new file mode 100644 index 0000000000000000000000000000000000000000..56b6fbb0df050c7400cada2f2c160e3554f81200 --- /dev/null +++ b/include/usgscsm/UsgsAstroProjectedLsSensorModel.h @@ -0,0 +1,266 @@ +/** Copyright © 2017-2022 BAE Systems Information and Electronic Systems Integration Inc. + +Redistribution and use in source and binary forms, with or without modification, are permitted +provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this list of conditions +and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, this list of +conditions and the following disclaimer in the documentation and/or other materials provided +with the distribution. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to +endorse or promote products derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. **/ + +#ifndef INCLUDE_USGSCSM_USGSASTROPROJECTEDLSSENSORMODEL_H_ +#define INCLUDE_USGSCSM_USGSASTROPROJECTEDLSSENSORMODEL_H_ + +#include <RasterGM.h> +#include <SettableEllipsoid.h> + +#include<utility> +#include<memory> +#include<string> +#include<vector> + +#include "ale/Orientations.h" +#include "ale/States.h" + +#include "spdlog/spdlog.h" + +#include "UsgsAstroLsSensorModel.h" + +class UsgsAstroProjectedLsSensorModel : public UsgsAstroLsSensorModel { + public: + // Initializes the class from state data as formatted + // in a string by the toString() method + void setState(const std::string& state); + + virtual void replaceModelState(const std::string& stateString); + //> This method attempts to initialize the current model with the state + // given by argState. The argState argument can be a string previously + // retrieved from the getModelState method. + // + // If argState contains a valid state for the current model, + // the internal state of the model is updated. + // + // If the model cannot be updated to the given state, a csm::Error is + // thrown and the internal state of the model is undefined. + // + // If the argument state string is empty, the model remains unchanged. + //< + + // This method checks to see if the model name is recognized + // in the input state string. + static std::string getModelNameFromModelState(const std::string& model_state); + + std::string constructStateFromIsd(const std::string imageSupportData, + csm::WarningList* list); + + // State data elements; + std::vector<double> m_geoTransform; + std::string m_projString; + + // Define logging pointer and file content + std::shared_ptr<spdlog::logger> m_logger = spdlog::get("usgscsm_logger"); + + // Hardcoded + static const std::string _SENSOR_MODEL_NAME; // state date element 0 + + static const std::string _STATE_KEYWORD[]; + + // Set to default values + void reset(); + + //-------------------------------------------------------------- + // Constructors/Destructor + //-------------------------------------------------------------- + + UsgsAstroProjectedLsSensorModel(); + ~UsgsAstroProjectedLsSensorModel(); + + virtual std::string getModelState() const; + + // Set the sensor model based on the input state data + void set(const std::string& state_data); + + //---------------------------------------------------------------- + // The following public methods are implementations of + // the methods inherited from RasterGM and SettableEllipsoid. + // These are defined in the CSM API. + //---------------------------------------------------------------- + + //--- + // Core Photogrammetry + //--- + virtual csm::ImageCoord groundToImage( + const csm::EcefCoord& groundPt, double desiredPrecision = 0.001, + double* achievedPrecision = NULL, + csm::WarningList* warnings = NULL) const; + + //> This method converts the given groundPt (x,y,z in ECEF meters) to a + // returned image coordinate (line, sample in full image space pixels). + // + // Iterative algorithms will use desiredPrecision, in meters, as the + // convergence criterion, otherwise it will be ignored. + // + // If a non-NULL achievedPrecision argument is received, it will be + // populated with the actual precision, in meters, achieved by iterative + // algorithms and 0.0 for deterministic algorithms. + // + // If a non-NULL warnings argument is received, it will be populated + // as applicable. + //< + + virtual csm::ImageCoordCovar groundToImage( + const csm::EcefCoordCovar& groundPt, double desiredPrecision = 0.001, + double* achievedPrecision = NULL, + csm::WarningList* warnings = NULL) const; + //> This method converts the given groundPt (x,y,z in ECEF meters and + // corresponding 3x3 covariance in ECEF meters squared) to a returned + // image coordinate with covariance (line, sample in full image space + // pixels and corresponding 2x2 covariance in pixels squared). + // + // Iterative algorithms will use desiredPrecision, in meters, as the + // convergence criterion, otherwise it will be ignored. + // + // If a non-NULL achievedPrecision argument is received, it will be + // populated with the actual precision, in meters, achieved by iterative + // algorithms and 0.0 for deterministic algorithms. + // + // If a non-NULL warnings argument is received, it will be populated + // as applicable. + //< + + virtual csm::EcefCoord imageToGround(const csm::ImageCoord& imagePt, + double height, + double desiredPrecision = 0.001, + double* achievedPrecision = NULL, + csm::WarningList* warnings = NULL) const; + //> This method converts the given imagePt (line,sample in full image + // space pixels) and given height (in meters relative to the WGS-84 + // ellipsoid) to a returned ground coordinate (x,y,z in ECEF meters). + // + // Iterative algorithms will use desiredPrecision, in meters, as the + // convergence criterion, otherwise it will be ignored. + // + // If a non-NULL achievedPrecision argument is received, it will be + // populated with the actual precision, in meters, achieved by iterative + // algorithms and 0.0 for deterministic algorithms. + // + // If a non-NULL warnings argument is received, it will be populated + // as applicable. + //< + + virtual csm::EcefCoordCovar imageToGround( + const csm::ImageCoordCovar& imagePt, double height, double heightVariance, + double desiredPrecision = 0.001, double* achievedPrecision = NULL, + csm::WarningList* warnings = NULL) const; + //> This method converts the given imagePt (line, sample in full image + // space pixels and corresponding 2x2 covariance in pixels squared) + // and given height (in meters relative to the WGS-84 ellipsoid) and + // corresponding heightVariance (in meters) to a returned ground + // coordinate with covariance (x,y,z in ECEF meters and corresponding + // 3x3 covariance in ECEF meters squared). + // + // Iterative algorithms will use desiredPrecision, in meters, as the + // convergence criterion, otherwise it will be ignored. + // + // If a non-NULL achievedPrecision argument is received, it will be + // populated with the actual precision, in meters, achieved by iterative + // algorithms and 0.0 for deterministic algorithms. + // + // If a non-NULL warnings argument is received, it will be populated + // as applicable. + //< + + virtual csm::EcefLocus imageToProximateImagingLocus( + const csm::ImageCoord& imagePt, const csm::EcefCoord& groundPt, + double desiredPrecision = 0.001, double* achievedPrecision = NULL, + csm::WarningList* warnings = NULL) const; + //> This method, for the given imagePt (line, sample in full image space + // pixels), returns the position and direction of the imaging locus + // nearest the given groundPt (x,y,z in ECEF meters). + // + // Note that there are two opposite directions possible. Both are + // valid, so either can be returned; the calling application can convert + // to the other as necessary. + // + // Iterative algorithms will use desiredPrecision, in meters, as the + // convergence criterion for the locus position, otherwise it will be + // ignored. + // + // If a non-NULL achievedPrecision argument is received, it will be + // populated with the actual precision, in meters, achieved by iterative + // algorithms and 0.0 for deterministic algorithms. + // + // If a non-NULL warnings argument is received, it will be populated + // as applicable. + //< + + virtual csm::EcefLocus imageToRemoteImagingLocus( + const csm::ImageCoord& imagePt, double desiredPrecision = 0.001, + double* achievedPrecision = NULL, + csm::WarningList* warnings = NULL) const; + //> This method, for the given imagePt (line, sample in full image space + // pixels), returns the position and direction of the imaging locus + // at the sensor. + // + // Note that there are two opposite directions possible. Both are + // valid, so either can be returned; the calling application can convert + // to the other as necessary. + // + // Iterative algorithms will use desiredPrecision, in meters, as the + // convergence criterion for the locus position, otherwise it will be + // ignored. + // + // If a non-NULL achievedPrecision argument is received, it will be + // populated with the actual precision, in meters, achieved by iterative + // algorithms and 0.0 for deterministic algorithms. + // + // If a non-NULL warnings argument is received, it will be populated + // as applicable. + // + // Notes: + // + // The remote imaging locus is only well-defined for optical sensors. + // It is undefined for SAR sensors and might not be available for + // polynomial and other non-physical models. The + // imageToProximateImagingLocus method should be used instead where + // possible. + //< + + //--- + // Error Correction + //--- + + virtual csm::Version getVersion() const; + //> This method returns the version of the model code. The Version + // object can be compared to other Version objects with its comparison + // operators. Not to be confused with the CSM API version. + //< + + virtual std::string getModelName() const; + //> This method returns a string identifying the name of the model. + //< + + virtual std::string getPedigree() const; + //> This method returns a string that identifies the sensor, + // the model type, its mode of acquisition and processing path. + // For example, an optical sensor model or a cubic rational polynomial + // model created from the same sensor's support data would produce + // different pedigrees for each case. + //< +}; + +#endif // INCLUDE_USGSCSM_USGSASTROPROJECTEDLSSENSORMODEL_H_ diff --git a/include/usgscsm/Utilities.h b/include/usgscsm/Utilities.h index d847e1e1d960979db643134165538c2025ad10dd..f1106299a48f361967f5e81c213e2f2c37b5a229 100644 --- a/include/usgscsm/Utilities.h +++ b/include/usgscsm/Utilities.h @@ -199,4 +199,8 @@ void applyRotationTranslationToXyzVec(ale::Rotation const& r, ale::Vec3d const& // to a calendar time string, such as 2000-01-01T00:16:40Z. std::string ephemTimeToCalendarTime(double ephemTime); +std::vector<double> pixelToMeter(double line, double sample, std::vector<double> geoTransform); + +std::vector<double> meterToPixel(double meter_x, double meter_y, std::vector<double> geoTransform); + #endif // INCLUDE_USGSCSM_UTILITIES_H_ diff --git a/src/UsgsAstroLsSensorModel.cpp b/src/UsgsAstroLsSensorModel.cpp index e5d3b2dcf8c009d421b225d84d258346d8292593..8d1d9857a3a78732a8a4462f92112c63a9380beb 100644 --- a/src/UsgsAstroLsSensorModel.cpp +++ b/src/UsgsAstroLsSensorModel.cpp @@ -702,7 +702,7 @@ csm::ImageCoord UsgsAstroLsSensorModel::groundToImage( ground_pt.x, ground_pt.y, ground_pt.z, desired_precision); // The public interface invokes the private interface with no adjustments. - csm::ImageCoord imagePt = groundToImage( + csm::ImageCoord imagePt = UsgsAstroLsSensorModel::groundToImage( ground_pt, _no_adjustment, desired_precision, achieved_precision, warnings); MESSAGE_LOG( spdlog::level::info, @@ -837,7 +837,7 @@ csm::ImageCoordCovar UsgsAstroLsSensorModel::groundToImage( gp.z = groundPt.z; csm::ImageCoord ip = - groundToImage(gp, desired_precision, achieved_precision, warnings); + UsgsAstroLsSensorModel::groundToImage(gp, desired_precision, achieved_precision, warnings); csm::ImageCoordCovar result(ip.line, ip.samp); // Compute partials ls wrt XYZ @@ -968,20 +968,20 @@ csm::EcefCoordCovar UsgsAstroLsSensorModel::imageToGround( const double DELTA_GROUND = m_gsd; csm::ImageCoord ip(image_pt.line, image_pt.samp); - csm::EcefCoord gp = imageToGround(ip, height, desired_precision, - achieved_precision, warnings); + csm::EcefCoord gp = UsgsAstroLsSensorModel::imageToGround(ip, height, desired_precision, + achieved_precision, warnings); // Compute numerical partials xyz wrt to lsh ip.line = image_pt.line + DELTA_IMAGE; ip.samp = image_pt.samp; - csm::EcefCoord gpl = imageToGround(ip, height, desired_precision); + csm::EcefCoord gpl = UsgsAstroLsSensorModel::imageToGround(ip, height, desired_precision); double xpl = (gpl.x - gp.x) / DELTA_IMAGE; double ypl = (gpl.y - gp.y) / DELTA_IMAGE; double zpl = (gpl.z - gp.z) / DELTA_IMAGE; ip.line = image_pt.line; ip.samp = image_pt.samp + DELTA_IMAGE; - csm::EcefCoord gps = imageToGround(ip, height, desired_precision); + csm::EcefCoord gps = UsgsAstroLsSensorModel::imageToGround(ip, height, desired_precision); double xps = (gps.x - gp.x) / DELTA_IMAGE; double yps = (gps.y - gp.y) / DELTA_IMAGE; double zps = (gps.z - gp.z) / DELTA_IMAGE; @@ -989,7 +989,7 @@ csm::EcefCoordCovar UsgsAstroLsSensorModel::imageToGround( ip.line = image_pt.line; ip.samp = image_pt.samp; csm::EcefCoord gph = - imageToGround(ip, height + DELTA_GROUND, desired_precision); + UsgsAstroLsSensorModel::imageToGround(ip, height + DELTA_GROUND, desired_precision); double xph = (gph.x - gp.x) / DELTA_GROUND; double yph = (gph.y - gp.y) / DELTA_GROUND; double zph = (gph.z - gp.z) / DELTA_GROUND; @@ -1064,7 +1064,7 @@ csm::EcefLocus UsgsAstroLsSensorModel::imageToProximateImagingLocus( // Ground point on object ray with the same elevation csm::EcefCoord gp1 = - imageToGround(image_pt, height, desired_precision, achieved_precision); + UsgsAstroLsSensorModel::imageToGround(image_pt, height, desired_precision, achieved_precision); // Vector between 2 ground points above double dx1 = x - gp1.x; @@ -1072,8 +1072,8 @@ csm::EcefLocus UsgsAstroLsSensorModel::imageToProximateImagingLocus( double dz1 = z - gp1.z; // Unit vector along object ray - csm::EcefCoord gp2 = imageToGround(image_pt, height - DELTA_GROUND, - desired_precision, achieved_precision); + csm::EcefCoord gp2 = UsgsAstroLsSensorModel::imageToGround(image_pt, height - DELTA_GROUND, + desired_precision, achieved_precision); double dx2 = gp2.x - gp1.x; double dy2 = gp2.y - gp1.y; double dz2 = gp2.z - gp1.z; @@ -1093,8 +1093,8 @@ csm::EcefLocus UsgsAstroLsSensorModel::imageToProximateImagingLocus( double hLocus = computeEllipsoidElevation(gp2.x, gp2.y, gp2.z, m_majorAxis, m_minorAxis, desired_precision); - locus.point = imageToGround(image_pt, hLocus, desired_precision, - achieved_precision, warnings); + locus.point = UsgsAstroLsSensorModel::imageToGround(image_pt, hLocus, desired_precision, + achieved_precision, warnings); locus.direction.x = dx2; locus.direction.y = dy2; @@ -2735,7 +2735,7 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd( } throw csm::Error(csm::Error::SENSOR_MODEL_NOT_CONSTRUCTIBLE, "ISD is invalid for creating the sensor model.", - "UsgsAstroFrameSensorModel::constructStateFromIsd"); + "UsgsAstroLsSensorModel::constructStateFromIsd"); } // The state data will still be updated when a sensor model is created since diff --git a/src/UsgsAstroPlugin.cpp b/src/UsgsAstroPlugin.cpp index 8c5395eb0e08f90ea85fc61d147249f9c0fc13b8..ba24543c7b35b086583b0b4d0de821a464b7e43c 100644 --- a/src/UsgsAstroPlugin.cpp +++ b/src/UsgsAstroPlugin.cpp @@ -26,6 +26,7 @@ OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "UsgsAstroFrameSensorModel.h" #include "UsgsAstroLsSensorModel.h" +#include "UsgsAstroProjectedLsSensorModel.h" #include "UsgsAstroPushFrameSensorModel.h" #include "UsgsAstroSarSensorModel.h" @@ -133,6 +134,7 @@ std::string UsgsAstroPlugin::getModelName(size_t modelIndex) const { std::vector<std::string> supportedModelNames = { UsgsAstroFrameSensorModel::_SENSOR_MODEL_NAME, UsgsAstroLsSensorModel::_SENSOR_MODEL_NAME, + UsgsAstroProjectedLsSensorModel::_SENSOR_MODEL_NAME, UsgsAstroSarSensorModel::_SENSOR_MODEL_NAME, UsgsAstroPushFrameSensorModel::_SENSOR_MODEL_NAME}; MESSAGE_LOG(spdlog::level::debug, "Get Model Name: {}. Used index: {}", @@ -186,7 +188,7 @@ bool UsgsAstroPlugin::canModelBeConstructedFromState( if (warnings) { warnings->push_back(csm::Warning (csm::Warning::UNKNOWN_WARNING, fullMsg, - "UsgsAstroFrameSensorModel::canModelBeConstructedFromState()")); + "UsgsAstroPlugin::canModelBeConstructedFromState()")); } return false; @@ -210,7 +212,7 @@ bool UsgsAstroPlugin::canModelBeConstructedFromISD( MESSAGE_LOG(spdlog::level::warn, msg); warnings->push_back(csm::Warning( csm::Warning::UNKNOWN_WARNING, msg, - "UsgsAstroFrameSensorModel::canModelBeConstructedFromISD()")); + "UsgsAstroPlugin::canModelBeConstructedFromISD()")); } } catch (...) { if (warnings) { @@ -220,7 +222,7 @@ bool UsgsAstroPlugin::canModelBeConstructedFromISD( MESSAGE_LOG(spdlog::level::warn, msg); warnings->push_back(csm::Warning( csm::Warning::UNKNOWN_WARNING, msg, - "UsgsAstroFrameSensorModel::canModelBeConstructedFromISD()")); + "UsgsAstroPlugin::canModelBeConstructedFromISD()")); } } return false; @@ -294,7 +296,7 @@ bool UsgsAstroPlugin::canISDBeConvertedToModelState( MESSAGE_LOG(spdlog::level::warn, msg); warnings->push_back(csm::Warning( csm::Warning::UNKNOWN_WARNING, msg, - "UsgsAstroFrameSensorModel::canISDBeConvertedToModelState()")); + "UsgsAstroPlugin::canISDBeConvertedToModelState()")); } return false; } @@ -367,6 +369,26 @@ csm::Model *UsgsAstroPlugin::constructModelFromISD( throw csm::Error(aErrorType, aMessage, aFunction); } return model; + } else if (modelName == UsgsAstroProjectedLsSensorModel::_SENSOR_MODEL_NAME) { + UsgsAstroProjectedLsSensorModel *model = new UsgsAstroProjectedLsSensorModel(); + try { + MESSAGE_LOG(spdlog::level::debug, "Trying to construct a UsgsAstroProjectedLsSensorModel"); + model->replaceModelState( + model->constructStateFromIsd(stringIsd, warnings)); + } catch (std::exception &e) { + delete model; + csm::Error::ErrorType aErrorType = + csm::Error::SENSOR_MODEL_NOT_CONSTRUCTIBLE; + std::string aMessage = "Could not construct model ["; + aMessage += modelName; + aMessage += "] with error ["; + aMessage += e.what(); + aMessage += "]"; + std::string aFunction = "UsgsAstroPlugin::constructModelFromISD()"; + MESSAGE_LOG(spdlog::level::err, aMessage); + throw csm::Error(aErrorType, aMessage, aFunction); + } + return model; } else if (modelName == UsgsAstroSarSensorModel::_SENSOR_MODEL_NAME) { UsgsAstroSarSensorModel *model = new UsgsAstroSarSensorModel(); MESSAGE_LOG(spdlog::level::debug, "Trying to construct a UsgsAstroSarSensorModel"); @@ -433,7 +455,12 @@ csm::Model *UsgsAstroPlugin::constructModelFromState( UsgsAstroLsSensorModel *model = new UsgsAstroLsSensorModel(); model->replaceModelState(modelState); return model; - } else if (modelName == UsgsAstroSarSensorModel::_SENSOR_MODEL_NAME) { + } else if (modelName == UsgsAstroProjectedLsSensorModel::_SENSOR_MODEL_NAME) { + MESSAGE_LOG(spdlog::level::debug, "Constructing a UsgsAstroProjectedLsSensorModel"); + UsgsAstroProjectedLsSensorModel *model = new UsgsAstroProjectedLsSensorModel(); + model->replaceModelState(modelState); + return model; + }else if (modelName == UsgsAstroSarSensorModel::_SENSOR_MODEL_NAME) { MESSAGE_LOG(spdlog::level::debug, "Constructing a UsgsAstroSarSensorModel"); UsgsAstroSarSensorModel *model = new UsgsAstroSarSensorModel(); model->replaceModelState(modelState); diff --git a/src/UsgsAstroProjectedLsSensorModel.cpp b/src/UsgsAstroProjectedLsSensorModel.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d2eb6395dbacb0a42c7354834852adea82fb6d8c --- /dev/null +++ b/src/UsgsAstroProjectedLsSensorModel.cpp @@ -0,0 +1,398 @@ +/** Copyright © 2017-2022 BAE Systems Information and Electronic Systems Integration Inc. + +Redistribution and use in source and binary forms, with or without modification, are permitted +provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this list of conditions +and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, this list of +conditions and the following disclaimer in the documentation and/or other materials provided +with the distribution. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to +endorse or promote products derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. **/ + +#include "UsgsAstroProjectedLsSensorModel.h" +#include "Utilities.h" + +#include <proj.h> + +#include <Error.h> +#include <nlohmann/json.hpp> + +#include "ale/Util.h" + +#define MESSAGE_LOG(...) \ + if (m_logger) { \ + m_logger->log(__VA_ARGS__); \ + } + +using json = nlohmann::json; + +const std::string UsgsAstroProjectedLsSensorModel::_SENSOR_MODEL_NAME = + "USGS_ASTRO_PROJECTED_LINE_SCANNER_SENSOR_MODEL"; + +const std::string UsgsAstroProjectedLsSensorModel::_STATE_KEYWORD[] = { + "m_modelName", + "m_imageIdentifier", + "m_sensorName", + "m_nLines", + "m_nSamples", + "m_platformFlag", + "m_intTimeLines", + "m_intTimeStartTimes", + "m_intTimes", + "m_startingEphemerisTime", + "m_centerEphemerisTime", + "m_detectorSampleSumming", + "m_detectorSampleSumming", + "m_startingDetectorSample", + "m_startingDetectorLine", + "m_ikCode", + "m_focalLength", + "m_zDirection", + "m_distortionType", + "m_opticalDistCoeffs", + "m_iTransS", + "m_iTransL", + "m_detectorSampleOrigin", + "m_detectorLineOrigin", + "m_majorAxis", + "m_minorAxis", + "m_platformIdentifier", + "m_sensorIdentifier", + "m_minElevation", + "m_maxElevation", + "m_dtEphem", + "m_t0Ephem", + "m_dtQuat", + "m_t0Quat", + "m_numPositions", + "m_numQuaternions", + "m_positions", + "m_velocities", + "m_quaternions", + "m_currentParameterValue", + "m_parameterType", + "m_referencePointXyz", + "m_sunPosition", + "m_sunVelocity", + "m_gsd", + "m_flyingHeight", + "m_halfSwath", + "m_halfTime", + "m_covariance", + "m_geoTransform", + "m_projString", +}; + +//*************************************************************************** +// UsgsAstroLineScannerSensorModel::replaceModelState +//*************************************************************************** +void UsgsAstroProjectedLsSensorModel::replaceModelState(const std::string& stateString) { + reset(); + + auto j = stateAsJson(stateString); + m_geoTransform = j["m_geoTransform"].get<std::vector<double>>(); + m_projString = j["m_projString"]; + MESSAGE_LOG( + spdlog::level::trace, + "m_geoTransform: {} " + "m_projString: {} ", + j["m_geoTransform"].dump(), j["m_projString"].dump()); + UsgsAstroLsSensorModel::replaceModelState(stateString); +} + +//*************************************************************************** +// UsgsAstroLineScannerSensorModel::getModelNameFromModelState +//*************************************************************************** +std::string UsgsAstroProjectedLsSensorModel::getModelNameFromModelState( + const std::string& model_state) { + // Parse the string to JSON + auto j = stateAsJson(model_state); + // If model name cannot be determined, return a blank string + std::string model_name; + + if (j.find("m_modelName") != j.end()) { + model_name = j["m_modelName"]; + } else { + csm::Error::ErrorType aErrorType = csm::Error::INVALID_SENSOR_MODEL_STATE; + std::string aMessage = "No 'm_modelName' key in the model state object."; + std::string aFunction = "UsgsAstroProjectedLsPlugin::getModelNameFromModelState"; + csm::Error csmErr(aErrorType, aMessage, aFunction); + throw(csmErr); + } + if (model_name != _SENSOR_MODEL_NAME) { + csm::Error::ErrorType aErrorType = csm::Error::SENSOR_MODEL_NOT_SUPPORTED; + std::string aMessage = "Sensor model not supported."; + std::string aFunction = "UsgsAstroProjectedLsPlugin::getModelNameFromModelState()"; + csm::Error csmErr(aErrorType, aMessage, aFunction); + throw(csmErr); + } + return model_name; +} + +//*************************************************************************** +// UsgsAstroLineScannerSensorModel::getModelState +//*************************************************************************** +std::string UsgsAstroProjectedLsSensorModel::getModelState() const { + auto state = stateAsJson(UsgsAstroLsSensorModel::getModelState()); + state["m_geoTransform"] = m_geoTransform; + state["m_projString"] = m_projString; + MESSAGE_LOG( + spdlog::level::trace, + "m_geoTransform: {}, {}, {}, {}, {}, {} " + "m_projString: {} ", + m_geoTransform[0], + m_geoTransform[1], + m_geoTransform[2], + m_geoTransform[3], + m_geoTransform[4], + m_geoTransform[5], + m_projString); + // Use dump(2) to avoid creating the model string as a single long line + std::string stateString = getModelName() + "\n" + state.dump(2); + return stateString; +} + +//*************************************************************************** +// UsgsAstroLineScannerSensorModel::reset +//*************************************************************************** +void UsgsAstroProjectedLsSensorModel::reset() { + MESSAGE_LOG(spdlog::level::debug, "Running reset()"); + + m_geoTransform = std::vector<double>(6, 0.0); + m_projString = ""; +} + +//***************************************************************************** +// UsgsAstroProjectedLsSensorModel Constructor +//***************************************************************************** +UsgsAstroProjectedLsSensorModel::UsgsAstroProjectedLsSensorModel() : UsgsAstroLsSensorModel() {} + +//***************************************************************************** +// UsgsAstroProjectedLsSensorModel Destructor +//***************************************************************************** +UsgsAstroProjectedLsSensorModel::~UsgsAstroProjectedLsSensorModel() {} + +//--------------------------------------------------------------------------- +// Core Photogrammetry +//--------------------------------------------------------------------------- + +//*************************************************************************** +// UsgsAstroLsSensorModel::groundToImage +//*************************************************************************** +csm::ImageCoord UsgsAstroProjectedLsSensorModel::groundToImage( + const csm::EcefCoord &ground_pt, double desired_precision, + double *achieved_precision, csm::WarningList *warnings) const { + + PJ_CONTEXT *C = proj_context_create(); + + /* Create a projection. */ + PJ *isdProj = proj_create(C, (m_projString + " +type=crs").c_str()); + if (0 == isdProj) { + MESSAGE_LOG( + spdlog::level::debug, + "Failed to create isd transformation object"); + return csm::ImageCoord(0, 0); + } + + /* Create the geocentric projection for our target */ + std::string radius_a = "+a=" + std::to_string(m_majorAxis); + std::string radius_b = "+b=" + std::to_string(m_minorAxis); + std::string projString = "+proj=geocent " + radius_a + " " + radius_b + " +type=crs"; + PJ *ecefProj = proj_create(C, projString.c_str()); + if (0 == ecefProj) { + MESSAGE_LOG( + spdlog::level::debug, + "Failed to create geocent transformation object"); + return csm::ImageCoord(0, 0); + } + + // Compute the transformation from our ISIS projection to ecef + PJ *isdProj2ecefProj = proj_create_crs_to_crs_from_pj(C, isdProj, ecefProj, 0, 0); + PJ_COORD c_in; + c_in.xyz.x = ground_pt.x; + c_in.xyz.y = ground_pt.y; + c_in.xyz.z = ground_pt.z; + MESSAGE_LOG( + spdlog::level::info, + "Ground point {}, {}, {}", + c_in.xyz.x, c_in.xyz.y, c_in.xyz.z); + PJ_COORD c_out = proj_trans(isdProj2ecefProj, PJ_INV, c_in); + MESSAGE_LOG( + spdlog::level::info, + "Meters {}, {}", + c_out.xyz.x, c_out.xyz.y); + std::vector<double> lineSampleCoord = meterToPixel(c_out.xyz.x, c_out.xyz.y, m_geoTransform); + csm::ImageCoord imagePt(lineSampleCoord[0], lineSampleCoord[1]); + MESSAGE_LOG( + spdlog::level::info, + "groundToImage result of ({}, {})", + imagePt.line, imagePt.samp); + return imagePt; +} + +//*************************************************************************** +// UsgsAstroLsSensorModel::groundToImage +//*************************************************************************** +csm::ImageCoordCovar UsgsAstroProjectedLsSensorModel::groundToImage( + const csm::EcefCoordCovar &groundPt, double desired_precision, + double *achieved_precision, csm::WarningList *warnings) const { + csm::ImageCoordCovar imageCoordCovar = UsgsAstroLsSensorModel::groundToImage(groundPt, desired_precision, achieved_precision, warnings); + csm::ImageCoord projImagePt = groundToImage(groundPt); + + imageCoordCovar.line = projImagePt.line; + imageCoordCovar.samp = projImagePt.samp; + return imageCoordCovar; +} + +//*************************************************************************** +// UsgsAstroProjectedLsSensorModel::imageToGround +//*************************************************************************** +csm::EcefCoord UsgsAstroProjectedLsSensorModel::imageToGround( + const csm::ImageCoord& image_pt, double height, double desired_precision, + double* achieved_precision, csm::WarningList* warnings) const { + MESSAGE_LOG( + spdlog::level::info, + "Computing imageToGround for {}, {}, {}, with desired precision {}", + image_pt.line, image_pt.samp, height, desired_precision); + + double x = 0, y = 0, z = 0; + double meterLine, meterSamp; + std::vector<double> meterCoord = pixelToMeter(image_pt.line, image_pt.samp, m_geoTransform); + meterLine = meterCoord[0]; + meterSamp = meterCoord[1]; + PJ_CONTEXT *C = proj_context_create(); + + /* Create a projection. */ + PJ *isdProj = proj_create(C, (m_projString + " +type=crs").c_str()); + if (0 == isdProj) { + MESSAGE_LOG( + spdlog::level::debug, + "Failed to create isd transformation object"); + return csm::EcefCoord(x, y, z); + } + + /* Create the geocentric projection for our target */ + std::string radius_a = "+a=" + std::to_string(m_majorAxis); + std::string radius_b = "+b=" + std::to_string(m_minorAxis); + std::string projString = "+proj=geocent " + radius_a + " " + radius_b + " +type=crs"; + PJ *ecefProj = proj_create(C, projString.c_str()); + if (0 == ecefProj) { + MESSAGE_LOG( + spdlog::level::debug, + "Failed to create geocent transformation object"); + return csm::EcefCoord(x, y, z); + } + + // Compute the transformation from our ISIS projection to ecef + PJ *isdProj2ecefProj = proj_create_crs_to_crs_from_pj(C, isdProj, ecefProj, 0, 0); + PJ_COORD c_in; + c_in.xy.x = meterSamp; + c_in.xy.y = meterLine; + PJ_COORD c_out = proj_trans(isdProj2ecefProj, PJ_FWD, c_in); + x = c_out.xyz.x, y = c_out.xyz.y, z = c_out.xyz.z; + MESSAGE_LOG( + spdlog::level::info, + "imageToGround result {} {} {}", + x, y, z); + return csm::EcefCoord(x, y, z); +} + +//*************************************************************************** +// UsgsAstroProjectedLsSensorModel::imageToGround +//*************************************************************************** +csm::EcefCoordCovar UsgsAstroProjectedLsSensorModel::imageToGround( + const csm::ImageCoordCovar& image_pt, double height, double heightVariance, + double desired_precision, double* achieved_precision, + csm::WarningList* warnings) const { + csm::EcefCoord groundCoord = imageToGround(image_pt, height); + csm::ImageCoord cameraImagePt = UsgsAstroLsSensorModel::groundToImage(groundCoord); + csm::ImageCoordCovar cameraImagePtCovar(cameraImagePt.line, cameraImagePt.samp); + + return UsgsAstroLsSensorModel::imageToGround(cameraImagePtCovar, height, heightVariance, desired_precision, achieved_precision, warnings); +} + +//*************************************************************************** +// UsgsAstroProjectedLsSensorModel::imageToProximateImagingLocus +//*************************************************************************** +csm::EcefLocus UsgsAstroProjectedLsSensorModel::imageToProximateImagingLocus( + const csm::ImageCoord& image_pt, const csm::EcefCoord& ground_pt, + double desired_precision, double* achieved_precision, + csm::WarningList* warnings) const { + csm::EcefCoord projGround = imageToGround(image_pt, 0); + csm::ImageCoord cameraImagePt = UsgsAstroLsSensorModel::groundToImage(projGround); + + return UsgsAstroLsSensorModel::imageToProximateImagingLocus(cameraImagePt, ground_pt, desired_precision, achieved_precision, warnings); +} + +//*************************************************************************** +// UsgsAstroProjectedLsSensorModel::imageToRemoteImagingLocus +//*************************************************************************** +csm::EcefLocus UsgsAstroProjectedLsSensorModel::imageToRemoteImagingLocus( + const csm::ImageCoord& image_pt, double desired_precision, + double* achieved_precision, csm::WarningList* warnings) const { + // Go from proj x, y to latlon then ground to image + // Convert imagept to camera imagept + csm::EcefCoord groundCoord = imageToGround(image_pt, 0); + csm::ImageCoord cameraImagePt = UsgsAstroLsSensorModel::groundToImage(groundCoord); + + return UsgsAstroLsSensorModel::imageToRemoteImagingLocus(cameraImagePt, desired_precision, achieved_precision, warnings); +} + +//--------------------------------------------------------------------------- +// Sensor Model Information +//--------------------------------------------------------------------------- + +//*************************************************************************** +// UsgsAstroProjectedLsSensorModel::getPedigree +//*************************************************************************** +std::string UsgsAstroProjectedLsSensorModel::getPedigree() const { + return "USGS_PROJECTED_LINE_SCANNER"; +} + +//*************************************************************************** +// UsgsAstroProjectedLsSensorModel::getSensorModelName +//*************************************************************************** +std::string UsgsAstroProjectedLsSensorModel::getModelName() const { + return UsgsAstroProjectedLsSensorModel::_SENSOR_MODEL_NAME; +} + +//*************************************************************************** +// UsgsAstroProjectedLsSensorModel::getVersion +//*************************************************************************** +csm::Version UsgsAstroProjectedLsSensorModel::getVersion() const { + return csm::Version(1, 0, 0); +} + +//*************************************************************************** +// UsgsAstroLineScannerSensorModel::constructStateFromIsd +//*************************************************************************** +std::string UsgsAstroProjectedLsSensorModel::constructStateFromIsd( + const std::string imageSupportData, csm::WarningList* warnings) { + // return UsgsAstroLsSensorModel::constructStateFromIsd(imageSupportData, warnings); + json lsState = json::parse(UsgsAstroLsSensorModel::constructStateFromIsd(imageSupportData, warnings)); + json state = json::parse(imageSupportData); + + lsState["m_geoTransform"] = ale::getGeoTransform(state); + lsState["m_projString"] = ale::getProjection(state); + MESSAGE_LOG( + spdlog::level::trace, + "m_geoTransform: {} " + "m_projString: {} ", + lsState["m_geoTransform"].dump(), lsState["m_projString"].dump()); + // The state data will still be updated when a sensor model is created since + // some state data is not in the ISD and requires a SM to compute them. + return lsState.dump(); +} diff --git a/src/Utilities.cpp b/src/Utilities.cpp index 14ed72ee316e19b8583db27d2de3826557616ce2..5599c252995bd5b394eb387afc4f130a3d2ec690 100644 --- a/src/Utilities.cpp +++ b/src/Utilities.cpp @@ -1528,3 +1528,23 @@ std::string ephemTimeToCalendarTime(double ephemTime) { buffer[21] = '\0'; return buffer; } + +std::vector<double> pixelToMeter(double line, double sample, std::vector<double> geoTransform) { + double meter_x = (sample * geoTransform[1]) + geoTransform[0]; + double meter_y = (line * geoTransform[5]) + geoTransform[3]; + + meter_x += geoTransform[1] * 0.5; + meter_y += geoTransform[5] * 0.5; + + return {meter_y, meter_x}; +} + +std::vector<double> meterToPixel(double meter_x, double meter_y, std::vector<double> geoTransform) { + meter_x -= geoTransform[1] * 0.5; + meter_y -= geoTransform[5] * 0.5; + + double sample = (meter_x - geoTransform[0]) / geoTransform[1]; + double line = (meter_y - geoTransform[3]) / geoTransform[5]; + + return {line, sample}; +} \ No newline at end of file diff --git a/tests/UtilitiesTests.cpp b/tests/UtilitiesTests.cpp index 75797876143af358019e68117374b62434dc96bb..36a2294b6d3a0b29b2f95f63b26f7dcdc374ed2a 100644 --- a/tests/UtilitiesTests.cpp +++ b/tests/UtilitiesTests.cpp @@ -636,3 +636,17 @@ TEST(UtilitiesTests, sanitizeTest) { sanitize(input); EXPECT_STREQ(input.c_str(), "\nHello World\n"); } + +TEST(UtilitiesTests, pixelToMeterTest) { + std::vector<double> geoTransform = {-104607.78948592, 4.980137561815, 0.0, -570564.40018202, 0.0, -4.980137561815}; + std::vector<double> ret = pixelToMeter(13.1191, 4981.08, geoTransform); + EXPECT_NEAR(ret[0], -570632.225173488, 1e-4); + EXPECT_NEAR(ret[1], -79798.83581073358, 1e-4); +} + +TEST(UtilitiesTests, meterToPixelTest) { + std::vector<double> geoTransform = {-104607.78948592, 4.980137561815, 0.0, -570564.40018202, 0.0, -4.980137561815}; + std::vector<double> ret = meterToPixel(-79798.83581073358, -570632.225173488, geoTransform); + EXPECT_NEAR(ret[0], 13.1191, 1e-4); + EXPECT_NEAR(ret[1], 4981.0800000000099, 1e-4); +} \ No newline at end of file