From a3d433c0fdda95b89b2348af9db9458039786a6d Mon Sep 17 00:00:00 2001
From: Jesse Mapel <jam826@nau.edu>
Date: Fri, 22 Feb 2019 11:47:36 -0700
Subject: [PATCH] Merge dev into master for 1.1 release (#190)
* Updated ISD and tests (#87)
* updated tests + RIP TestyMcTestFace
* keys updated
* Probably should not hard code .cpps in the CMakeLists.txt for GTest
* updated specs to mimic synthetic json
* appveyor assumed broken from standup, disabling until fixed
* eulers are now quats
* merge from jay's changes
* merge markers
* renamed quats as per Jessie's reccomendation
* only appends sensor_orientation to missing params
* stupid markers
* updated test ISD
* changes as per some comments
* updated as per comments
* defaults are all 0
* #pragma once to C style header guard
* stragler change
* missed some params
* missed another
* moving things around
* patched segfault error + added detector center as defined in json schema (#94)
* patched segfault error + added detector center as defined in the swagger schema
* verbose travi
* Fixed optical distortion and pixel to focal transforms. (#115)
* Update to use quaternions in test to match new ISD
* Fixes #121
* Fix 1 remaining failing omega rotation test
* Updates token to be in quotes
* Still trying to get the builds working.
* Updates to force upload
* Adds newline
* Update for formatting
* Trying to force travis to deploy with labels
I am editing this on the repo in order to get Travis to do its thing - I can not do this via a PR.
* Update .travis.yml
* Update meta.yaml
to try and force a dev tag on the builds to get labels to work.
* Update meta.yaml
* Update .travis.yml
* test x scaling for reset model state
* test FL500 conversion
* test FL500 conversion
* Add tests to include
* commit to force run of tests on mac with debug output
* switch to using tolerances from dev
* 0.5 pixel offset
* try moving set of filename into setup
* initialize differently test
* more debug output
* Add matrix to debug output
* Even more debug output from matrix calculation
* Add missing quaternion component to state string
* cleanup debug output now that problem is fixed
* fix added spacing
* Adds in check for PR for upload
* echo to check travis env
* fixes equality check for PR
* Clean up some spacing (non-functional) in UsgsAstroFrameModel
* Fixed canModelBeConstructedFromISD throwing an exception when the metadata file doesn't exist. (#134)
* fixed can be converted error
* Added not constructible test for frame plugin
* Removed old LS ISD header (#137)
* Update some of the tests that set a new state value to use a function in the fixture, rather than repeating code
* Windows Build (#139)
* adds win build
* Windows build
* Updates submodules
* Refactoring to move to one plugin and removed unused classes. (#138)
* merge
* changed varrs
* reset test
* First iteration
* it compiles, although renaming this is still neccessary
* added source
* cleaned up, tests still need to pass
* post merge clean up
* validation updated, validation tests are now passing
* Addressed comments from Jesse
* last Jesse comment, convertISDToModelState doesn't check if pointer is invalid anymore as the sensor model should except
* model_name -> m_modelName for frame getState
* copy paste error in FrameSensorModel
* Resolve merge conflicts with FrameIsdTest vs. SimpleFrameIsdTest names
* Fix hopefully last merge problem
* Conda build on win (#143)
* Conda build on win
* Trying a straight build first
* Unneeded csmapi
* Now trying to upload build
* Trying a build and upload
* Trying syntax change
* trying ps
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Update meta.yaml
* Update bld.bat
This is a test to see where appveyor is failing.
* Update bld.bat
* Update bld.bat
* Update meta.yaml
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Update .appveyor.yml
* Adds if/else into appveyor (#146)
* Adds if/else into appveyor
* Update .appveyor.yml
* Updates to use gcc7
* Update .travis.yml
* Update .appveyor.yml
* Fixed Line Scan construction and condensed plugin tests (#145)
* First pass at test LS ISD
* Initial LS plugin test base
* fixed a test in frame plugin tests
* Initial LS plugin test suite
* Moved to single plugin test file
* Added some new plugin tests
* Fixed LS construction
* Re-updated submodules
* Reverted gtest
* removed debuf prints left in and made getModelName functional (#148)
* Changed line scan sensor model to new ISD spec (#149)
* Changed line scan sensor model to new ISD spec
* Updated LS to the proper spec
* Changed model_name to name_model
* Updated tests to use new name_ format in ISD
* Updated LS test data to new spec
* Fixed typo in ls test ISD
* Moved framer to new ISD format. Added bad debug statements.
* Updated LS to new spec
* Fixed focal length epsilon name
* Fixes 142 for framer
* Removes unused state keyword counter
* Updates Ls
* Updates for parametercov and ref point
* Adds sensor and platform names
* Updates with collection identifier
* updates for vector parameters in state
* Added model state test for framer (#152)
* Adds current param value and referencePt
* unsupported pushed in
* fixed nan issue (#153)
* fixed nan issue
* left in for loop
* Removed switch statement for ls distortion: (#150)
* Updated build recipe to use release and sha hash
* Copied distortion inversion from ISIS3 Distortion Map (#161)
* Fixed ground to image and image to ground not giving inverse results. (#162)
* Copied distortion inversion from ISIS3 Distortion Map
* mid-testing commit to fix LS i2g and g2i inconsistency
* Fixed ground to image taking excessive iterations. Added simple line scan test and test stubs
* Added more comments and replaced some hard coded tolerances
* Added debug output to ground to image error message. (#168)
* Initial steps to modularize losToEcf in UsgsAstroLsSensorModel
* Try changing line endings
* Re-converted line endings to dos
* Modularize some of computeViewingPixel
* Cleanup
* Reconvert to dos line endings
* Reformat member variables in header and switch more to use const pass-by-reference parameters to be consistent with rest of code
* Updates in response to code review (part 1)
* Update calculateRotationMatrixFromQuaternions -- removed invert parameter. Now, will always return the matrix that rotates the in the same direction as the quaternions.
* Updates the version number in the meta. (#159)
* Updates the version number in the meta.
* Update meta.yaml
* Build to hash
This should be using the short hash now.
* Update meta.yaml
* Update meta.yaml
* Finish Modularizing computeViewingPixel (#170)
* Finishing modularization of UsgsAstroLsSensorModel::computeViewingPixel
* Fixes per request
* Removing ISIS references in variable names while I'm here
* Example for LS testing (#171)
* Removed odd height gradient descent search from image to ground when not intersecting.
* Added first off body test.
* Added angular velocity line scan test
* Added failing, copy paste tests
* Commented out line scan pan tests until 0 velocity bug is addressed
* Removed extraneous focal epsilon from LS test ISDs
* Extracts distortion models from the Line Scanner and Framer cameras into its own function set (#172)
* Moved transverse distortion from line scanner into new file
* Updated fixtures for distortion tests
* Forgot a semicolon
* Extracted radial destortion, and built associated tests
* Removed invertDistortion commented code from Linescanner
* Reconstructed distortion functions return results
* Incorperated invertDistortion function
* Corrected spelling, and removed commented out functions
* Renamed all dpoint variables to distortionPoint
* Removed computeUndistortedFocalPlaneCoordinates function
* Renamed distortionPoint to undistortedPoint and distortedPoint where applicable
* Fixed typo
* Removed doc string parameters
* Further extracted transverse destortion out of the Framer
* Pulled losToEcf shared functions out into Utilities class (#174)
* Initial steps to pull common camera model or math utility functions into a separate Utilities file.
* dos2unix line endings again
* ../include/usgscsm/UsgsAstroLsSensorModel.h
* updates to modularized functions
* dos2unix UsgsAstroFrameSensorModel
* Forgot to actually add the Utilities class
* Keep needed to unix2dos to re-fix line ending caracters
* We have some dos line ending and some unix line endings in our code
* Copied unit conversion from framer to ls (#177)
* Makes GTest fully optional (#175)
* removes unused setup.py (#178)
* Refactored ISD parsing into utilities and added warnings (#179)
* Added utility ISD parsing functions
* ISD Parsing Test stub
* Added ISD parsing method tests
* Added ISD Parsing utilities to Framer
* Replaced framer ISD parsing
* added check for warning output to bad fram ISD test
* Added a bit more testing
* Removed ground to image 10m errormessage (#180)
* Fixed state being written out and read in differently (#183)
* Fixed state being written out and read in differently
* Fixed canModelBeConstructedFromState
* Updated README (#154)
* Updated README
* Added BAE license and removed open source license
* Gtests for Utilities.cpp (#181)
* Finishing modularization of UsgsAstroLsSensorModel::computeViewingPixel
* Fixes per request
* Removing ISIS references in variable names while I'm here
* Beginning to implement tests for modularization changes in LS sensor model class
* Adding tests for Utilities as well as a few minor modifications
* Making suggested changes and adding test for caculateAttitudeCorrection
* Fixing vector size
* Updated LS ISD parsing code to match Framer (#186)
* updated LS ISD parsing to match style in framer
* reverted the 1 to 0
* Fixed attitude test, removed all cerr lines
* removed image flip flag
* added tests
* Fixed computeDistortedFocalPlaneCoordinates test (#187)
* Distortion Integration (#184)
* Moved transverse distortion from line scanner into new file
* Extracted radial destortion, and built associated tests
* Reconstructed distortion functions return results
* Corrected spelling, and removed commented out functions
* Renamed all dpoint variables to distortionPoint
* Removed computeUndistortedFocalPlaneCoordinates function
* Renamed distortionPoint to undistortedPoint and distortedPoint where applicable
* Removed doc string parameters
* Added combined new distortion model
* Half way through redefining parameters
* Major update to distortion and coefficient storage
* More debugging
* More debugging
* More debugging
* debugging
* debugging
* Removed debugging messages and set the output value for remove distortion
* Brought changes inline with dev
* Updated isd parsing to handle either radial or transverse
* Added Distortion Parse function
* Split apply and remove distortion, and added the distortion type to each model state
* Removed cmath
* Removing prints, and other fixes
* Updated how coefficients are parsed
* Reverted notebook
* Function namespace and parameter clean up
* More function and parameter clean up
* Reverted old code and updated distortion defaults
* Small changes for distortion
* Updated utilities and fixed failing test
* Final LS tests and release prep (#188)
* Added image locus tests
* Pan LS tests now pass
* updated release date in plugin
* Forgot to update plugin release date test
---
CMakeLists.txt | 20 +-
LICENSE | 24 -
README.md | 154 +--
include/usgscsm/Distortion.h | 30 +
include/usgscsm/UsgsAstroFrameSensorModel.h | 22 +-
include/usgscsm/UsgsAstroLsSensorModel.h | 125 ++-
include/usgscsm/UsgsAstroPlugin.h | 4 -
include/usgscsm/Utilities.h | 88 ++
license.txt | 11 +
setup.py | 54 -
src/Distortion.cpp | 237 ++++
src/UsgsAstroFrameSensorModel.cpp | 623 ++++-------
src/UsgsAstroLsSensorModel.cpp | 1055 +++++++-----------
src/UsgsAstroPlugin.cpp | 17 +-
src/Utilities.cpp | 729 ++++++++++++
tests/CMakeLists.txt | 6 +-
tests/DistortionTests.cpp | 151 +++
tests/Fixtures.h | 61 +
tests/FrameCameraTests.cpp | 164 +--
tests/ISDParsingTests.cpp | 313 ++++++
tests/LineScanCameraTests.cpp | 127 +++
tests/PluginTests.cpp | 26 +-
tests/UtilitiesTests.cpp | 69 ++
tests/data/constAngularVelocityLineScan.json | 92 ++
tests/data/constVelocityLineScan.json | 13 +-
tests/data/empty.json | 1 +
tests/data/simpleFramerISD.json | 3 +
27 files changed, 2755 insertions(+), 1464 deletions(-)
delete mode 100644 LICENSE
create mode 100644 include/usgscsm/Distortion.h
create mode 100644 include/usgscsm/Utilities.h
create mode 100644 license.txt
delete mode 100644 setup.py
create mode 100644 src/Distortion.cpp
create mode 100644 src/Utilities.cpp
create mode 100644 tests/DistortionTests.cpp
create mode 100644 tests/ISDParsingTests.cpp
create mode 100644 tests/LineScanCameraTests.cpp
create mode 100644 tests/UtilitiesTests.cpp
create mode 100644 tests/data/constAngularVelocityLineScan.json
create mode 100644 tests/data/empty.json
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 76ab4a6..625de6a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,9 +1,9 @@
cmake_minimum_required(VERSION 3.10)
project(usgscsm VERSION 0.0.1 DESCRIPTION "usgscsm library")
+include(GoogleTest)
include(cmake/gtest.cmake)
include(GNUInstallDirs)
-include(GoogleTest)
set(CMAKE_CXX_STANDARD 11)
@@ -18,7 +18,7 @@ else()
PATH_SUFFIXES "csm"
PATHS $ENV{CONDA_PREFIX}/include/)
find_library(CSM_LIBRARY csmapi PATHS $ENV{CONDA_PREFIX}/lib)
-
+
message("--Found CSM Library: ${CSM_LIBRARY}")
message("--Found CSM Include Directory: ${CSM_INCLUDE_DIR}")
endif(BUILD_CSM)
@@ -26,7 +26,9 @@ endif(BUILD_CSM)
add_library(usgscsm SHARED
src/UsgsAstroPlugin.cpp
src/UsgsAstroFrameSensorModel.cpp
- src/UsgsAstroLsSensorModel.cpp)
+ src/UsgsAstroLsSensorModel.cpp
+ src/Distortion.cpp
+ src/Utilities.cpp)
set_target_properties(usgscsm PROPERTIES
VERSION ${PROJECT_VERSION}
@@ -42,12 +44,8 @@ target_include_directories(usgscsm
${CSM_INCLUDE_DIR}
)
-# Setup for GoogleTest
-find_package (Threads)
-
target_link_libraries(usgscsm
- ${CSM_LIBRARY}
- gtest ${CMAKE_THREAD_LIBS_INIT})
+ ${CSM_LIBRARY})
if(WIN32)
install(TARGETS usgscsm RUNTIME DESTINATION ${CMAKE_INSTALL_LIBDIR})
@@ -60,6 +58,12 @@ install(DIRECTORY ${USGSCSM_INCLUDE_DIRS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR
# Optional build or link against CSM
option (BUILD_TESTS "Build tests" ON)
if(BUILD_TESTS)
+
+ # Setup for GoogleTest
+ find_package (Threads)
+
+ target_link_libraries(usgscsm
+ gtest ${CMAKE_THREAD_LIBS_INIT})
include(CTest)
enable_testing()
add_subdirectory(tests)
diff --git a/LICENSE b/LICENSE
deleted file mode 100644
index cf1ab25..0000000
--- a/LICENSE
+++ /dev/null
@@ -1,24 +0,0 @@
-This is free and unencumbered software released into the public domain.
-
-Anyone is free to copy, modify, publish, use, compile, sell, or
-distribute this software, either in source code form or as a compiled
-binary, for any purpose, commercial or non-commercial, and by any
-means.
-
-In jurisdictions that recognize copyright laws, the author or authors
-of this software dedicate any and all copyright interest in the
-software to the public domain. We make this dedication for the benefit
-of the public at large and to the detriment of our heirs and
-successors. We intend this dedication to be an overt act of
-relinquishment in perpetuity of all present and future rights to this
-software under copyright law.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
-OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
-ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
-OTHER DEALINGS IN THE SOFTWARE.
-
-For more information, please refer to <http://unlicense.org>
diff --git a/README.md b/README.md
index 9f4d5c6..473f62d 100644
--- a/README.md
+++ b/README.md
@@ -1,100 +1,66 @@
# CSM-CameraModel
-ISD Specification: https://github.com/USGS-Astrogeology/pfeffernusse/blob/master/swagger.yaml
-
-## Setting up dependencies with conda (RECOMMENDED)
-
-Install conda if you do not already have it.
-```bash
-wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
-bash miniconda.sh -b
-```
-> You can add a `-p <install-prefix>` to choose where to install miniconda. By default, it will install it to `$HOME/miniconda3`.
-
-### Setting up conda for bash
-Copy and paste the following into a terminal running the `bash` shell:
-```bash
-echo -e "\n\n# Adding miniconda3 to PATH" >> $HOME/.bashrc && \
-echo -e "export PATH=$HOME/miniconda3/bin:\$PATH" >> $HOME/.bashrc && \
-source $HOME/.bashrc && \
-which conda
-```
-> *For more information: [bash installation](https://conda.io/docs/user-guide/install/linux.html "Reference to bash conda install")*
-
-### Setting up conda for tcsh
-Copy and paste the following into a terminal running the `tcsh` shell:
-```tcsh
-echo "\n\n# Setting up miniconda3 for tcsh" >> $HOME/.cshrc && \
-echo "source $HOME/miniconda3/etc/profile.d/conda.csh > /dev/null" >> $HOME/.cshrc && \
-source $HOME/.cshrc && \
-which conda
-```
-> *For more information: [tcsh installation](https://github.com/ESMValGroup/ESMValTool/issues/301 "Reference to tcsh conda install")*
-
-### Creating an isolated conda environment
-Run the following commands to create a self-contained dev environment for CSM-CameraModel (type `y` to confirm creation):
-```bash
-conda create -n csmdev -c usgs-astrogeology cmake libcsm
-```
-> *For more information: [conda environments](https://conda.io/docs/user-guide/tasks/manage-environments.html)*
-
-### Activating the environment
-After creating the `csmdev` environment and installing cmake and libcsm into it, we need to activate it. Right now, cmake, libcsm, and their dependencies are isolated to a conda environment and we need to tell conda that we want to use it. The activation command depends on your shell.
-* **bash**: `source activate csmdev`
-* **tcsh**: `conda activate csmdev`
-> *You can add these to the end of your $HOME/.bashrc or $HOME/.cshrc if you want the `csmdev` environment to be active in every new terminal.*
+Community Sensor Model (CSM) compliant sensor models created by USGS Astrogeology
+Science Center.
-## Building CSM-CameraModel
-After you've set up conda, you can build CSM-CameraModel:
-
-1. Fork `USGS-Astrogeology/CSM-CameraModel` if you don't already have a fork.
-2. Clone your fork of `CSM-CameraModel` *with `--recursive` option to get the gtest submodule*.
-```bash
-git clone --recursive git@github.com:<your-username>/CSM-CameraModel.git
-cd CSM-CameraModel
-git remote add upstream git@github.com:USGS-Astrogeology/CSM-CameraModel.git
-```
-3. Sync your fork with `upstream` and ensure the gtest submodule is init'd if your fork is old.
-```bash
-git pull upstream master
-git submodule update --init --recursive
-git push -u origin master
-```
-4. `mkdir build && cd build`
-5. `cmake .. && make`
-6. `ctest`
+CSM-CameraModel contains two different sensor models. The first, is a generic
+framing camera model written from scratch. The second is a generic line scan
+camera model based on code from BAE Systems Information and Electronic Systems
+Integration Inc.
+
+## Using CSM-CameraModel
+
+This library is a CSM plugin library that is intended to be dynamically loaded
+at run time along side the
+[CSM API library](https://github.com/USGS-Astrogeology/csm).
+
+Once, the library is loaded, it can be accessed through the CSM Plugin interface.
+For an example of how to do through the CSM c++ interface see the SensorModelFactory
+class in [SensorUtils](https://github.com/USGS-Astrogeology/SensorUtils).
+For an example of how to do this through the CSM Python bindings see this
+[notebook](http://nbviewer.jupyter.org/gist/thareUSGS/4c0eb72799edc33ff4816b2587027148).
+
+From the CSM Plugin interface, a generic framing camera model
+(USGS_ASTRO_FRAME_SENSOR_MODEL) or generic line scan camera model
+(USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL) can be instantiated from suitable Image
+Support Data (ISD). Under the CSM standard, each plugin library can define its
+own ISD format. This library uses an auxiliary JSON formatted file that must be
+next to the image file passed to the CSM::ISD class. We provide an OpenAPI
+server for generating these,
+[pfeffernusse](https://github.com/USGS-Astrogeology/pfeffernusse). The swagger
+specification is located on
+[swaggerhub](https://app.swaggerhub.com/apis/USGS-Astro/pfeffernusse2/0.1.4-oas3).
---
-## Building without a package manager
-To build:
-
-1. Install [libcsmapi](https://github.com/sminster/csm "CSM API")
- > You can install this with an `INSTDIR` of your choice, or let it default (see [libcsmapi README](https://github.com/sminster/csm/blob/master/README))
-```bash
-mkdir $HOME/csmenv
-cd $HOME
-git clone git@github.com:sminster/csm.git
-cd csm
-make -f Makefile.linux64 all install clean INSTDIR="$csmenv"
-```
-2. Install cmake >= 3.10
-```bash
-cd $HOME
-wget https://cmake.org/files/v3.12/cmake-3.12.0-Linux-x86_64.tar.gz
-tar xzf cmake-3.12.0-Linux-x86_64.tar.gz
-rsync -azv cmake-3.12.0-Linux-x86_64/ $HOME/csmenv/
-echo -e "\n#Prepending csm env to path\nsetenv PATH "$HOME/csmenv/bin:$PATH" >> $HOME/.cshrc
-source $HOME/.cshrc
-```
-3. Fork and clone down this repo and its submodules (gtest)
-```bash
-git clone --recursive git@github.com:<username>/CSM-CameraModel.git
-cd CSM-CameraModel
-git remote add upstream git@github.com:USGS-Astrogeology/CSM-CameraModel.git
-git pull upstream master
-git submodule update --init --recursive
-git push -u origin master
-```
-4. `mkdir build` && `cd build`
-5. `cmake -DCSM_INCLUDE_DIR="${csmenv}/include/csm -DCSM_LIBRARY="${csmenv}/lib/libcsmapi.so .. && make`
+## Build Requirements
+
+* cmake 3.10 or newer
+* GNU-compatible Make
+* a c++11 compliant compiler
+
+This repository has all of its external c++ dependencies included in it. The
+excellent header-only JSON library
+[JSON for Modern C++](https://github.com/nlohmann/json) is included directly in
+the source code. The other two dependencies, the CSM API library, and gtest
+are included as git submodules. When you clone this library make sure you add
+the `--recursive` flag to your `git clone` command. Alaterntively, you can run
+`git submodule update --init --recursive` after cloning. The library can also be
+compiled against an installed versions of the CSM API by setting the BUILD_CSM
+flag to OFF during cmake configuration.
+
+## Building CSM-CameraModel
+
+CSM-CameraModel uses a standard cmake build system. To compile the library, and
+tests use the following commands:
+
+1. `mkdir build && cd build`
+2. `cmake .. && cmake --build .`
+
+## Testing CSM-CameraModel
+
+All of the tests for CSM-CameraModel are written in the googletests framework
+and are run via ctest. To run all of the tests simply run `ctest` in the build.
+
+All of the tests are purposefully written to use generic data that values have
+been hand validated for. This data can be found under `tests/data`.
diff --git a/include/usgscsm/Distortion.h b/include/usgscsm/Distortion.h
new file mode 100644
index 0000000..574a970
--- /dev/null
+++ b/include/usgscsm/Distortion.h
@@ -0,0 +1,30 @@
+#ifndef Distortion_h
+#define Distortion_h
+
+#include <vector>
+#include <math.h>
+#include <tuple>
+#include <iostream>
+
+enum DistortionType {
+ RADIAL,
+ TRANSVERSE
+};
+
+// Transverse Distortion
+void distortionJacobian(double x, double y, double *jacobian,
+ const std::vector<double> opticalDistCoeffs);
+
+void computeTransverseDistortion(double ux, double uy, double &dx, double &dy,
+ const std::vector<double> opticalDistCoeffs);
+
+void removeDistortion(double dx, double dy, double &ux, double &uy,
+ const std::vector<double> opticalDistCoeffs,
+ DistortionType distortionType,
+ const double tolerance = 1.0E-6);
+
+void applyDistortion(double ux, double uy, double &dx, double &dy,
+ const std::vector<double> opticalDistCoeffs,
+ DistortionType distortionType,
+ const double desiredPrecision = 0, const double tolerance = 1.0E-6);
+#endif
diff --git a/include/usgscsm/UsgsAstroFrameSensorModel.h b/include/usgscsm/UsgsAstroFrameSensorModel.h
index 7da5473..7f3b2f1 100644
--- a/include/usgscsm/UsgsAstroFrameSensorModel.h
+++ b/include/usgscsm/UsgsAstroFrameSensorModel.h
@@ -8,6 +8,12 @@
#include <gtest/gtest.h>
#include "RasterGM.h"
#include "CorrelationModel.h"
+#include "Distortion.h"
+#include "Utilities.h"
+
+#include <json.hpp>
+using json = nlohmann::json;
+
#include <json.hpp>
using json = nlohmann::json;
@@ -320,12 +326,6 @@ protected:
FRIEND_TEST(FrameSensorModel, setFocalPlane_AlternatingOnes);
FRIEND_TEST(FrameSensorModel, distortMe_AlternatingOnes);
- virtual bool setFocalPlane(double dx,double dy,double &undistortedX,double &undistortedY) const;
- virtual void distortionFunction(double ux, double uy, double &dx, double &dy) const;
- virtual void distortionJacobian(double x, double y, double &Jxx,
- double &Jxy, double &Jyx, double &Jyy) const;
-
-
private:
// Input parameters
@@ -335,8 +335,8 @@ protected:
std::vector<double> m_currentParameterCovariance;
std::vector<csm::param::Type> m_parameterType;
std::vector<double> m_noAdjustments;
- std::vector<double> m_odtX;
- std::vector<double> m_odtY;
+ DistortionType m_distortionType;
+ std::vector<double> m_opticalDistCoeffs;
std::vector<double> m_transX;
std::vector<double> m_transY;
std::vector<double> m_spacecraftVelocity;
@@ -356,6 +356,10 @@ protected:
double m_startingDetectorLine;
std::string m_targetName;
std::string m_modelName;
+ std::string m_sensorName;
+ std::string m_platformName;
+ std::string m_imageIdentifier;
+ std::string m_collectionIdentifier;
double m_ifov;
std::string m_instrumentID;
double m_focalLengthEpsilon;
@@ -369,6 +373,8 @@ protected:
int m_nSamples;
int m_nParameters;
+ csm::EcefCoord m_referencePointXyz;
+
json _state;
static const int _NUM_STATE_KEYWORDS;
static const int NUM_PARAMETERS;
diff --git a/include/usgscsm/UsgsAstroLsSensorModel.h b/include/usgscsm/UsgsAstroLsSensorModel.h
index 2fc6d7e..5329ede 100644
--- a/include/usgscsm/UsgsAstroLsSensorModel.h
+++ b/include/usgscsm/UsgsAstroLsSensorModel.h
@@ -30,7 +30,7 @@
#include <RasterGM.h>
#include <SettableEllipsoid.h>
#include <CorrelationModel.h>
-
+#include "Distortion.h"
class UsgsAstroLsSensorModel : public csm::RasterGM, virtual public csm::SettableEllipsoid
{
@@ -64,60 +64,56 @@ public:
std::string constructStateFromIsd(const std::string imageSupportData, csm::WarningList *list) const;
// State data elements;
- std::string m_imageIdentifier; // 1
- std::string m_sensorType; // 2
- int m_totalLines; // 3
- int m_totalSamples; // 4
- double m_offsetLines; // 5
- double m_offsetSamples; // 6
- int m_platformFlag; // 7
- int m_aberrFlag; // 8
- int m_atmRefFlag; // 9
+ std::string m_imageIdentifier;
+ std::string m_sensorName;
+ int m_nLines;
+ int m_nSamples;
+ int m_platformFlag;
std::vector<double> m_intTimeLines;
std::vector<double> m_intTimeStartTimes;
std::vector<double> m_intTimes;
- double m_startingEphemerisTime; // 11
- double m_centerEphemerisTime; // 12
- double m_detectorSampleSumming; // 13
- double m_startingSample; // 14
- int m_ikCode; // 15
- double m_focal; // 16
- double m_isisZDirection; // 17
- double m_opticalDistCoef[3]; // 18
- double m_iTransS[3]; // 19
- double m_iTransL[3]; // 20
- double m_detectorSampleOrigin; // 21
- double m_detectorLineOrigin; // 22
- double m_detectorLineOffset; // 23
- double m_mountingMatrix[9]; // 24
- double m_semiMajorAxis; // 25
- double m_semiMinorAxis; // 26
- std::string m_referenceDateAndTime; // 27
- std::string m_platformIdentifier; // 28
- std::string m_sensorIdentifier; // 29
- std::string m_trajectoryIdentifier; // 30
- std::string m_collectionIdentifier; // 31
- double m_refElevation; // 32
- double m_minElevation; // 33
- double m_maxElevation; // 34
- double m_dtEphem; // 35
- double m_t0Ephem; // 36
- double m_dtQuat; // 37
- double m_t0Quat; // 38
- int m_numEphem; // 39
- int m_numQuaternions; // 40
- std::vector<double> m_ephemPts; // 41
- std::vector<double> m_ephemRates; // 42
- std::vector<double> m_quaternions; // 43
- std::vector<double> m_parameterVals; // 44
- std::vector<csm::param::Type> m_parameterType; // 45
- csm::EcefCoord m_referencePointXyz; // 46
- double m_gsd; // 47
- double m_flyingHeight; // 48
- double m_halfSwath; // 49
- double m_halfTime; // 50
- std::vector<double> m_covariance; // 51
- int m_imageFlipFlag; // 52
+ double m_startingEphemerisTime;
+ double m_centerEphemerisTime;
+ double m_detectorSampleSumming;
+ double m_startingSample;
+ int m_ikCode;
+ double m_focalLength;
+ double m_zDirection;
+ DistortionType m_distortionType;
+ std::vector<double> m_opticalDistCoeffs;
+ double m_iTransS[3];
+ double m_iTransL[3];
+ double m_detectorSampleOrigin;
+ double m_detectorLineOrigin;
+ double m_mountingMatrix[9];
+ double m_majorAxis;
+ double m_minorAxis;
+ std::string m_referenceDateAndTime;
+ std::string m_platformIdentifier;
+ std::string m_sensorIdentifier;
+ std::string m_trajectoryIdentifier;
+ std::string m_collectionIdentifier;
+ double m_refElevation;
+ double m_minElevation;
+ double m_maxElevation;
+ double m_dtEphem;
+ double m_t0Ephem;
+ double m_dtQuat;
+ double m_t0Quat;
+ int m_numPositions;
+ int m_numQuaternions;
+ std::vector<double> m_positions;
+ std::vector<double> m_velocities;
+ std::vector<double> m_quaternions;
+ std::vector<double> m_currentParameterValue;
+ std::vector<csm::param::Type> m_parameterType;
+ csm::EcefCoord m_referencePointXyz;
+ double m_gsd;
+ double m_flyingHeight;
+ double m_halfSwath;
+ double m_halfTime;
+ std::vector<double> m_covariance;
+ int m_imageFlipFlag;
// Hardcoded
static const std::string _SENSOR_MODEL_NAME; // state date element 0
@@ -892,6 +888,11 @@ public:
//> This method sets the planetary ellipsoid.
//<
+ void calculateAttitudeCorrection(
+ const double& time,
+ const std::vector<double>& adj,
+ double attCorr[9]) const;
+
private:
void determineSensorCovarianceInImageSpace(
@@ -916,7 +917,16 @@ private:
double* achievedPrecision = NULL,
csm::WarningList* warnings = NULL) const;
- // This method computes the imaging locus.
+ void reconstructSensorDistortion(
+ double& focalX,
+ double& focalY,
+ const double& desiredPrecision) const;
+
+ void getQuaternions(const double& time,
+ double quaternion[4]) const;
+
+// This method computes the imaging locus.
+// imaging locus : set of ground points associated with an image pixel.
void losToEcf(
const double& line, // CSM image convention
const double& sample, // UL pixel center == (0.5, 0.5)
@@ -927,9 +937,9 @@ private:
double& vx, // output sensor x velocity
double& vy, // output sensor y velocity
double& vz, // output sensor z cvelocity
- double& xl, // output line-of-sight x coordinate
- double& yl, // output line-of-sight y coordinate
- double& zl ) const;
+ double& bodyFixedX, // output line-of-sight x coordinate
+ double& bodyFixedY, // output line-of-sight y coordinate
+ double& bodyFixedZ ) const;
// Computes the LOS correction due to light aberration
void lightAberrationCorr(
@@ -1019,7 +1029,8 @@ private:
csm::ImageCoord computeViewingPixel(
const double& time, // The time to use the EO at
const csm::EcefCoord& groundPoint, // The ground coordinate
- const std::vector<double>& adj // Parameter Adjustments for partials
+ const std::vector<double>& adj, // Parameter Adjustments for partials
+ const double& desiredPrecision // Desired precision for distortion inversion
) const;
// The linear approximation for the sensor model is used as the starting point
diff --git a/include/usgscsm/UsgsAstroPlugin.h b/include/usgscsm/UsgsAstroPlugin.h
index a3a3c76..8789106 100644
--- a/include/usgscsm/UsgsAstroPlugin.h
+++ b/include/usgscsm/UsgsAstroPlugin.h
@@ -55,11 +55,7 @@ private:
static const std::string _MANUFACTURER_NAME;
static const std::string _RELEASE_DATE;
static const int _N_SENSOR_MODELS;
- static const int _NUM_ISD_KEYWORDS;
static const std::string _ISD_KEYWORD[];
- static const int _NUM_STATE_KEYWORDS;
- static const std::string _STATE_KEYWORD[];
- static const json MODEL_KEYWORDS;
typedef csm::Model* (*sensorConstructor)(void);
static std::map<std::string, sensorConstructor> MODELS;
diff --git a/include/usgscsm/Utilities.h b/include/usgscsm/Utilities.h
new file mode 100644
index 0000000..719a724
--- /dev/null
+++ b/include/usgscsm/Utilities.h
@@ -0,0 +1,88 @@
+#ifndef Utilities_h
+#define Utilities_h
+
+#include "Distortion.h"
+
+#include <vector>
+#include <math.h>
+#include <tuple>
+#include <string>
+
+#include <json.hpp>
+
+#include <Warning.h>
+
+// methods pulled out of los2ecf and computeViewingPixel
+
+// for now, put everything in here.
+// TODO: later, consider if it makes sense to pull sample/line offsets out
+// Compute distorted focalPlane coordinates in mm
+void computeDistortedFocalPlaneCoordinates(
+ const double& line,
+ const double& sample,
+ const double& sampleOrigin,
+ const double& lineOrigin,
+ const double& sampleSumming,
+ const double& startingSample,
+ const double iTransS[],
+ const double iTransL[],
+ std::tuple<double, double>& natFocalPlane);
+
+void calculateRotationMatrixFromQuaternions(
+ double quaternions[4],
+ double cameraToBody[9]);
+
+void calculateRotationMatrixFromEuler(
+ double euler[],
+ double rotationMatrix[]);
+
+void createCameraLookVector(
+ const double& undistortedFocalPlaneX,
+ const double& undistortedFocalPlaneY,
+ const double& zDirection,
+ const double& focalLength,
+ double cameraLook[]);
+
+//void calculateAttitudeCorrection(
+// const double& time,
+//
+// double attCorr[9]);
+//
+
+// Methods for checking/accessing the ISD
+
+double metric_conversion(double val, std::string from, std::string to="m");
+std::string getSensorModelName(nlohmann::json isd, csm::WarningList *list=nullptr);
+std::string getImageId(nlohmann::json isd, csm::WarningList *list=nullptr);
+std::string getSensorName(nlohmann::json isd, csm::WarningList *list=nullptr);
+std::string getPlatformName(nlohmann::json isd, csm::WarningList *list=nullptr);
+int getTotalLines(nlohmann::json isd, csm::WarningList *list=nullptr);
+int getTotalSamples(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getStartingTime(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getCenterTime(nlohmann::json isd, csm::WarningList *list=nullptr);
+std::vector<double> getIntegrationStartLines(nlohmann::json isd, csm::WarningList *list=nullptr);
+std::vector<double> getIntegrationStartTimes(nlohmann::json isd, csm::WarningList *list=nullptr);
+std::vector<double> getIntegrationTimes(nlohmann::json isd, csm::WarningList *list=nullptr);
+int getSampleSumming(nlohmann::json isd, csm::WarningList *list=nullptr);
+int getLineSumming(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getFocalLength(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getFocalLengthEpsilon(nlohmann::json isd, csm::WarningList *list=nullptr);
+std::vector<double> getFocal2PixelLines(nlohmann::json isd, csm::WarningList *list=nullptr);
+std::vector<double> getFocal2PixelSamples(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getDetectorCenterLine(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getDetectorCenterSample(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getDetectorStartingLine(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getDetectorStartingSample(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getMinHeight(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getMaxHeight(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getSemiMajorRadius(nlohmann::json isd, csm::WarningList *list=nullptr);
+double getSemiMinorRadius(nlohmann::json isd, csm::WarningList *list=nullptr);
+DistortionType getDistortionModel(nlohmann::json isd, csm::WarningList *list=nullptr);
+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> 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);
+
+#endif
diff --git a/license.txt b/license.txt
new file mode 100644
index 0000000..756adf6
--- /dev/null
+++ b/license.txt
@@ -0,0 +1,11 @@
+Copyright © 2017 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.
diff --git a/setup.py b/setup.py
deleted file mode 100644
index be02bd9..0000000
--- a/setup.py
+++ /dev/null
@@ -1,54 +0,0 @@
-import os
-import pkg_resources
-import sys
-import sysconfig
-from setuptools import setup, Extension, find_packages
-from Cython.Build import cythonize
-
-# Look for the csmapi headers in the standard location
-incdir = os.path.dirname(sysconfig.get_path('include'))
-
-INCLUDE_DIRS = ['include/json', 'include/genericframe', 'include/orex',
- 'include/genericls',
- 'include/', incdir, os.path.join(incdir, 'csm')]
-LIBRARY_DIRS = [] # This assumes that libcsmapi is installed in a standard place
-LIBRARIES = ['csmapi']
-COMPILE_ARGS = ['-g', '-std=c++11']#, '-stdlib=libc++']
-
-if sys.platform == 'darwin':
- COMPILE_ARGS.append('-mmacosx-version-min=10.9')
-elif sys.platform.startswith("win"):
- try:
- INCLUDE_DIRS.append(os.path.join(os.environ['LIBRARY_INC'], 'csm'))
- except: pass
- COMPILE_ARGS = []
-
-def generate_extension(path_name, sources):
- return Extension(path_name,
- sources=sources,
- extra_compile_args=COMPILE_ARGS,
- language='c++',
- include_dirs=INCLUDE_DIRS,
- runtime_library_dirs=LIBRARY_DIRS,
- libraries=LIBRARIES)
-
-# Create the extensions
-extensions = [generate_extension('usgscam.genericframe', ['usgscam/genericframe.pyx',
- 'src/UsgsAstroFramePlugin.cpp',
- 'src/UsgsAstroFrameSensorModel.cpp']),
- generate_extension('usgscam.orex', ['usgscam/orex.pyx',
- 'src/ORexPlugin.cpp',
- 'src/ORexSensorModel.cpp']),
- generate_extension('usgscam.genericls', ['usgscam/genericls.pyx',
- 'src/UsgsAstroLsPlugin.cpp',
- 'src/UsgsAstroLsSensorModel.cpp',
- 'src/UsgsAstroLsStateData.cpp'])]
-
-setup(
- name='usgscam',
- version='0.1.1',
- ext_modules=cythonize(extensions),
- description='Cython wrapper to the USGS MDIS Camera Model',
- author='Jay Laura',
- packages = find_packages(),
- install_requires=['cycsm'])
diff --git a/src/Distortion.cpp b/src/Distortion.cpp
new file mode 100644
index 0000000..7a43ceb
--- /dev/null
+++ b/src/Distortion.cpp
@@ -0,0 +1,237 @@
+#include "Distortion.h"
+
+/**
+ * @description Jacobian of the distortion function. The Jacobian was computed
+ * algebraically from the function described in the distortionFunction
+ * method.
+ *
+ * @param x
+ * @param y
+ * @param odtX opticalDistCoef In X
+ * @param odtY opticalDistCoef In Y
+ *
+ * @returns jacobian a jacobian vector of vectors as
+ [0][0]: xx, [0][1]: xy
+ [1][0]: yx, [1][1]: yy
+ */
+
+void distortionJacobian(double x, double y, double *jacobian,
+ const std::vector<double> opticalDistCoeffs) {
+
+ double d_dx[10];
+ d_dx[0] = 0;
+ d_dx[1] = 1;
+ d_dx[2] = 0;
+ d_dx[3] = 2 * x;
+ d_dx[4] = y;
+ d_dx[5] = 0;
+ d_dx[6] = 3 * x * x;
+ d_dx[7] = 2 * x * y;
+ d_dx[8] = y * y;
+ d_dx[9] = 0;
+ double d_dy[10];
+ d_dy[0] = 0;
+ d_dy[1] = 0;
+ d_dy[2] = 1;
+ d_dy[3] = 0;
+ d_dy[4] = x;
+ d_dy[5] = 2 * y;
+ d_dy[6] = 0;
+ d_dy[7] = x * x;
+ d_dy[8] = 2 * x * y;
+ d_dy[9] = 3 * y * y;
+
+ jacobian[0] = 0; // xx
+ jacobian[1] = 0; // xy
+ jacobian[2] = 0; // yx
+ jacobian[3] = 0; // yy
+
+ int xPointer = 0;
+ int yPointer = opticalDistCoeffs.size() / 2;
+
+ for (int i = 0; i < 10; xPointer++, yPointer++, i++) {
+ jacobian[0] = jacobian[0] + d_dx[i] * opticalDistCoeffs[xPointer];
+ jacobian[1] = jacobian[1] + d_dy[i] * opticalDistCoeffs[xPointer];
+ jacobian[2] = jacobian[2] + d_dx[i] * opticalDistCoeffs[yPointer];
+ jacobian[3] = jacobian[3] + d_dy[i] * opticalDistCoeffs[yPointer];
+ }
+}
+
+/**
+ * @description Compute distorted focal plane (dx,dy) coordinate given an undistorted focal
+ * plane (ux,uy) coordinate. This uses the third order Taylor approximation to the
+ * distortion model.
+ *
+ * @param ux Undistored x
+ * @param uy Undistored y
+ * @param opticalDistCoeffs For both X and Y coefficients
+ *
+ * @returns distortedPoint Newly adjusted focal plane coordinates as an x, y tuple
+ */
+void computeTransverseDistortion(double ux, double uy, double &dx, double &dy,
+ const std::vector<double> opticalDistCoeffs) {
+
+ double f[10];
+ f[0] = 1;
+ f[1] = ux;
+ f[2] = uy;
+ f[3] = ux * ux;
+ f[4] = ux * uy;
+ f[5] = uy * uy;
+ f[6] = ux * ux * ux;
+ f[7] = ux * ux * uy;
+ f[8] = ux * uy * uy;
+ f[9] = uy * uy * uy;
+
+ int xPointer = 0;
+ int yPointer = opticalDistCoeffs.size() / 2;
+
+ dx = 0.0;
+ dy = 0.0;
+
+ for (int i = 0; i < 10; xPointer++, yPointer++, i++) {
+ dx = dx + f[i] * opticalDistCoeffs[xPointer];
+ dy = dy + f[i] * opticalDistCoeffs[yPointer];
+ }
+}
+
+void removeDistortion(double dx, double dy, double &ux, double &uy,
+ const std::vector<double> opticalDistCoeffs,
+ DistortionType distortionType,
+ const double tolerance) {
+ ux = dx;
+ uy = dy;
+
+ switch (distortionType) {
+ // Compute undistorted focal plane coordinate given a distorted
+ // coordinate set and the distortion coefficients
+ case RADIAL: {
+ double rr = dx * dx + dy * dy;
+
+ if (rr > tolerance)
+ {
+ double dr = opticalDistCoeffs[0] + (rr * (opticalDistCoeffs[1] + rr * opticalDistCoeffs[2]));
+ ux = dx * (1.0 - dr);
+ uy = dy * (1.0 - dr);
+ }
+ }
+ break;
+ // Computes undistorted focal plane (x,y) coordinates given a distorted focal plane (x,y)
+ // coordinate. The undistorted coordinates are solved for using the Newton-Raphson
+ // method for root-finding if the distortionFunction method is invoked.
+ case TRANSVERSE: {
+ // Solve the distortion equation using the Newton-Raphson method.
+ // Set the error tolerance to about one millionth of a NAC pixel.
+ // The maximum number of iterations of the Newton-Raphson method.
+ const int maxTries = 20;
+
+ double x;
+ double y;
+ double fx;
+ double fy;
+ double jacobian[4];
+
+ // Initial guess at the root
+ x = dx;
+ y = dy;
+
+ computeTransverseDistortion(x, y, fx, fy, opticalDistCoeffs);
+
+ for (int count = 1; ((fabs(fx) +fabs(fy)) > tolerance) && (count < maxTries); count++) {
+
+ computeTransverseDistortion(x, y, fx, fy, opticalDistCoeffs);
+
+ fx = dx - fx;
+ fy = dy - fy;
+
+ distortionJacobian(x, y, jacobian, opticalDistCoeffs);
+
+ // Jxx * Jyy - Jxy * Jyx
+ double determinant = jacobian[0] * jacobian[3] - jacobian[1] * jacobian[2];
+ if (fabs(determinant) < 1E-6) {
+ ux = x;
+ uy = y;
+ //
+ // Near-zero determinant. Add error handling here.
+ //
+ //-- Just break out and return with no convergence
+ return;
+ }
+
+ x = x + (jacobian[3] * fx - jacobian[1] * fy) / determinant;
+ y = y + (jacobian[0] * fy - jacobian[2] * fx) / determinant;
+ }
+
+ if ((fabs(fx) + fabs(fy)) <= tolerance) {
+ // The method converged to a root.
+ ux = x;
+ uy = y;
+
+ return;
+ }
+ // Otherwise method did not converge to a root within the maximum
+ // number of iterations
+ }
+ break;
+ }
+}
+
+void applyDistortion(double ux, double uy, double &dx, double &dy,
+ const std::vector<double> opticalDistCoeffs,
+ DistortionType distortionType,
+ const double desiredPrecision, const double tolerance)
+{
+ dx = ux;
+ dy = uy;
+
+ switch (distortionType) {
+ // Compute undistorted focal plane coordinate given a distorted
+ // focal plane coordinate. This case works by iteratively adding distortion
+ // until the new distorted point, r, undistorts to within a tolerance of the
+ // original point, rp.
+ case RADIAL: {
+ double rp2 = (ux * ux) + (uy * uy);
+
+ if (rp2 > tolerance) {
+ double rp = sqrt(rp2);
+ // Compute first fractional distortion using rp
+ double drOverR = opticalDistCoeffs[0]
+ + (rp2 * (opticalDistCoeffs[1] + (rp2 * opticalDistCoeffs[2])));
+ // Compute first distorted point estimate, r
+ double r = rp + (drOverR * rp);
+ double r_prev, r2_prev;
+ int iteration = 0;
+ do {
+ // Don't get in an end-less loop. This algorithm should
+ // converge quickly. If not then we are probably way outside
+ // of the focal plane. Just set the distorted position to the
+ // undistorted position. Also, make sure the focal plane is less
+ // than 1km, it is unreasonable for it to grow larger than that.
+ if (iteration >= 15 || r > 1E9) {
+ drOverR = 0.0;
+ break;
+ }
+
+ r_prev = r;
+ r2_prev = r * r;
+
+ // Compute new fractional distortion:
+ drOverR = opticalDistCoeffs[0]
+ + (r2_prev * (opticalDistCoeffs[1] + (r2_prev * opticalDistCoeffs[2])));
+
+ // Compute new estimate of r
+ r = rp + (drOverR * r_prev);
+ iteration++;
+ }
+ while (fabs(r * (1 - drOverR) - rp) > desiredPrecision);
+ dx = ux / (1.0 - drOverR);
+ dy = uy / (1.0 - drOverR);
+ }
+ }
+ break;
+ case TRANSVERSE: {
+ computeTransverseDistortion(ux, uy, dx, dy, opticalDistCoeffs);
+ }
+ break;
+ }
+}
diff --git a/src/UsgsAstroFrameSensorModel.cpp b/src/UsgsAstroFrameSensorModel.cpp
index 2ef4aea..a54c2c4 100644
--- a/src/UsgsAstroFrameSensorModel.cpp
+++ b/src/UsgsAstroFrameSensorModel.cpp
@@ -26,14 +26,16 @@ const std::string UsgsAstroFrameSensorModel::m_parameterName[] = {
"v3" // 6
};
-const int UsgsAstroFrameSensorModel::_NUM_STATE_KEYWORDS = 32;
-
UsgsAstroFrameSensorModel::UsgsAstroFrameSensorModel() {
reset();
}
void UsgsAstroFrameSensorModel::reset() {
m_modelName = _SENSOR_MODEL_NAME;
+ m_platformName = "";
+ m_sensorName = "";
+ m_imageIdentifier = "";
+ m_collectionIdentifier = "";
m_majorAxis = 0.0;
m_minorAxis = 0.0;
m_focalLength = 0.0;
@@ -59,14 +61,17 @@ void UsgsAstroFrameSensorModel::reset() {
m_ccdCenter = std::vector<double>(2, 0.0);
m_spacecraftVelocity = std::vector<double>(3, 0.0);
m_sunPosition = std::vector<double>(3, 0.0);
- m_odtX = std::vector<double>(10, 0.0);
- m_odtY = std::vector<double>(10, 0.0);
+ m_distortionType = DistortionType::TRANSVERSE;
+ m_opticalDistCoeffs = std::vector<double>(20, 0.0);
m_transX = std::vector<double>(3, 0.0);
m_transY = std::vector<double>(3, 0.0);
m_iTransS = std::vector<double>(3, 0.0);
m_iTransL = std::vector<double>(3, 0.0);
m_boresight = std::vector<double>(3, 0.0);
m_parameterType = std::vector<csm::param::Type>(NUM_PARAMETERS, csm::param::REAL);
+ m_referencePointXyz.x = 0;
+ m_referencePointXyz.y = 0;
+ m_referencePointXyz.z = 0;
}
@@ -129,14 +134,15 @@ csm::ImageCoord UsgsAstroFrameSensorModel::groundToImage(
undistortedy = (f * (m[0][1] * xo + m[1][1] * yo + m[2][1] * zo)/denom) + m_linePp;
// Apply the distortion to the line/sample location and then convert back to line/sample
- double distortedx, distortedy;
- distortionFunction(undistortedx, undistortedy, distortedx, distortedy);
+ double distortedX, distortedY;
+ applyDistortion(undistortedx, undistortedy, distortedX, distortedY,
+ m_opticalDistCoeffs, m_distortionType);
// Convert distorted mm into line/sample
double sample, line;
- sample = m_iTransS[0] + m_iTransS[1] * distortedx + m_iTransS[2] * distortedy + m_ccdCenter[1];
- line = m_iTransL[0] + m_iTransL[1] * distortedx + m_iTransL[2] * distortedy + m_ccdCenter[0];
+ sample = m_iTransS[0] + m_iTransS[1] * distortedX + m_iTransS[2] * distortedY + m_ccdCenter[1];
+ line = m_iTransL[0] + m_iTransL[1] * distortedX + m_iTransL[2] * distortedY + m_ccdCenter[0];
return csm::ImageCoord(line, sample);
}
@@ -147,9 +153,17 @@ csm::ImageCoordCovar UsgsAstroFrameSensorModel::groundToImage(const csm::EcefCoo
double *achievedPrecision,
csm::WarningList *warnings) const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::groundToImage");
+ csm::EcefCoord gp;
+ gp.x = groundPt.x;
+ gp.y = groundPt.y;
+ gp.z = groundPt.z;
+
+ csm::ImageCoord ip = groundToImage(
+ gp, desiredPrecision, achievedPrecision, warnings);
+ csm::ImageCoordCovar result(ip.line, ip.samp);
+ // This is a partial, incorrect implementation to test if SocetGXP needs
+ // this method implemented in order to load the sensor.
+ return result;
}
@@ -179,18 +193,15 @@ csm::EcefCoord UsgsAstroFrameSensorModel::imageToGround(const csm::ImageCoord &i
x_camera = m_transX[0] + m_transX[1] * (lo - line_center) + m_transX[2] * (so - sample_center);
// Apply the distortion model (remove distortion)
- double undistorted_cameraX, undistorted_cameraY = 0.0;
-
- setFocalPlane(x_camera, y_camera, undistorted_cameraX, undistorted_cameraY);
+ double undistortedX, undistortedY;
+ removeDistortion(x_camera, y_camera, undistortedX, undistortedY,
+ m_opticalDistCoeffs,
+ m_distortionType);
// Now back from distorted mm to pixels
- double udx, udy; //distorted line and sample
- udx = undistorted_cameraX;
- udy = undistorted_cameraY;
-
- xl = m[0][0] * udx + m[0][1] * udy - m[0][2] * -m_focalLength;
- yl = m[1][0] * udx + m[1][1] * udy - m[1][2] * -m_focalLength;
- zl = m[2][0] * udx + m[2][1] * udy - m[2][2] * -m_focalLength;
+ xl = m[0][0] * undistortedX + m[0][1] * undistortedY - m[0][2] * - m_focalLength;
+ yl = m[1][0] * undistortedX + m[1][1] * undistortedY - m[1][2] * - m_focalLength;
+ zl = m[2][0] * undistortedX + m[2][1] * undistortedY - m[2][2] * - m_focalLength;
double x, y, z;
double xc, yc, zc;
@@ -210,9 +221,9 @@ csm::EcefCoordCovar UsgsAstroFrameSensorModel::imageToGround(const csm::ImageCoo
double heightVariance, double desiredPrecision,
double *achievedPrecision,
csm::WarningList *warnings) const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::imageToGround");
+ // This is an incomplete implementation to see if SocetGXP needs this method implemented.
+ csm::EcefCoordCovar result;
+ return result;
}
@@ -238,10 +249,11 @@ csm::EcefLocus UsgsAstroFrameSensorModel::imageToRemoteImagingLocus(const csm::I
double focalPlaneY = m_transY[0] + m_transY[1] * row + m_transY[2] * col;
// Distort
- double undistortedFocalPlaneX = focalPlaneX;
- double undistortedFocalPlaneY = focalPlaneY;
-
- setFocalPlane(focalPlaneX, focalPlaneY, undistortedFocalPlaneX, undistortedFocalPlaneY);
+ double undistortedFocalPlaneX, undistortedFocalPlaneY;
+ removeDistortion(focalPlaneX, focalPlaneY,
+ undistortedFocalPlaneX, undistortedFocalPlaneY,
+ m_opticalDistCoeffs,
+ m_distortionType);
// Get rotation matrix and transform to a body-fixed frame
double m[3][3];
@@ -533,9 +545,8 @@ const csm::CorrelationModel& UsgsAstroFrameSensorModel::getCorrelationModel() co
std::vector<double> UsgsAstroFrameSensorModel::getUnmodeledCrossCovariance(const csm::ImageCoord &pt1,
const csm::ImageCoord &pt2) const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getUnmodeledCrossCovariance");
+ // No unmodeled error
+ return std::vector<double>(4, 0.0);
}
@@ -550,52 +561,38 @@ std::string UsgsAstroFrameSensorModel::getModelName() const {
std::string UsgsAstroFrameSensorModel::getPedigree() const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getPedigree");
+ return "USGS_FRAMER";
}
std::string UsgsAstroFrameSensorModel::getImageIdentifier() const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getImageIdentifier");
+ return m_imageIdentifier;
}
void UsgsAstroFrameSensorModel::setImageIdentifier(const std::string& imageId,
csm::WarningList* warnings) {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::setImageIdentifier");
+ m_imageIdentifier = imageId;
}
std::string UsgsAstroFrameSensorModel::getSensorIdentifier() const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getSensorIdentifier");
+ return m_sensorName;
}
std::string UsgsAstroFrameSensorModel::getPlatformIdentifier() const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getPlatformIdentifier");
+ return m_platformName;
}
std::string UsgsAstroFrameSensorModel::getCollectionIdentifier() const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getCollectionIdentifier");
+ return m_collectionIdentifier;
}
std::string UsgsAstroFrameSensorModel::getTrajectoryIdentifier() const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getTrajectoryIdentifier");
+ return "";
}
@@ -610,15 +607,15 @@ std::string UsgsAstroFrameSensorModel::getSensorMode() const {
std::string UsgsAstroFrameSensorModel::getReferenceDateAndTime() const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getReferenceDateAndTime");
+ return "";
}
std::string UsgsAstroFrameSensorModel::getModelState() const {
json state = {
{"m_modelName", _SENSOR_MODEL_NAME},
+ {"m_sensorName", m_sensorName},
+ {"m_platformName", m_platformName},
{"m_focalLength" , m_focalLength},
{"m_iTransS", {m_iTransS[0], m_iTransS[1], m_iTransS[2]}},
{"m_iTransL", {m_iTransL[0], m_iTransL[1], m_iTransL[2]}},
@@ -642,10 +639,8 @@ std::string UsgsAstroFrameSensorModel::getModelState() const {
{"m_samplePp", m_samplePp},
{"m_minElevation", m_minElevation},
{"m_maxElevation", m_maxElevation},
- {"m_odtX", {m_odtX[0], m_odtX[1], m_odtX[2], m_odtX[3], m_odtX[4],
- m_odtX[5], m_odtX[6], m_odtX[7], m_odtX[8], m_odtX[9]}},
- {"m_odtY", {m_odtY[0], m_odtY[1], m_odtY[2], m_odtY[3], m_odtY[4],
- m_odtY[5], m_odtY[6], m_odtY[7], m_odtY[8], m_odtY[9]}},
+ {"m_distortionType", m_distortionType},
+ {"m_opticalDistCoeffs", m_opticalDistCoeffs},
{"m_originalHalfLines", m_originalHalfLines},
{"m_originalHalfSamples", m_originalHalfSamples},
{"m_spacecraftName", m_spacecraftName},
@@ -653,10 +648,12 @@ std::string UsgsAstroFrameSensorModel::getModelState() const {
{"m_ephemerisTime", m_ephemerisTime},
{"m_nLines", m_nLines},
{"m_nSamples", m_nSamples},
- {"m_currentParameterValue", {m_currentParameterValue[0], m_currentParameterValue[1],
- m_currentParameterValue[2], m_currentParameterValue[3],
- m_currentParameterValue[4], m_currentParameterValue[5],
- m_currentParameterValue[6]}},
+ {"m_currentParameterValue", m_currentParameterValue},
+ {"m_imageIdentifier", m_imageIdentifier},
+ {"m_collectionIdentifier", m_collectionIdentifier},
+ {"m_referencePointXyz", {m_referencePointXyz.x,
+ m_referencePointXyz.y,
+ m_referencePointXyz.z}},
{"m_currentParameterCovariance", m_currentParameterCovariance}
};
@@ -678,8 +675,8 @@ bool UsgsAstroFrameSensorModel::isValidModelState(const std::string& stringState
"m_ccdCenter",
"m_spacecraftVelocity",
"m_sunPosition",
- "m_odtX",
- "m_odtY",
+ "m_distortionType",
+ "m_opticalDistCoeffs",
"m_transX",
"m_transY",
"m_iTransS",
@@ -737,7 +734,6 @@ bool UsgsAstroFrameSensorModel::isValidIsd(const std::string& Isd, csm::WarningL
void UsgsAstroFrameSensorModel::replaceModelState(const std::string& stringState) {
json state = json::parse(stringState);
-
// The json library's .at() will except if key is missing
try {
m_modelName = state.at("m_modelName").get<std::string>();
@@ -753,12 +749,22 @@ void UsgsAstroFrameSensorModel::replaceModelState(const std::string& stringState
m_ccdCenter = state.at("m_ccdCenter").get<std::vector<double>>();
m_spacecraftVelocity = state.at("m_spacecraftVelocity").get<std::vector<double>>();
m_sunPosition = state.at("m_sunPosition").get<std::vector<double>>();
- m_odtX = state.at("m_odtX").get<std::vector<double>>();
- m_odtY = state.at("m_odtY").get<std::vector<double>>();
+ m_distortionType = (DistortionType)state.at("m_distortionType").get<int>();
+ m_opticalDistCoeffs = state.at("m_opticalDistCoeffs").get<std::vector<double>>();
m_transX = state.at("m_transX").get<std::vector<double>>();
m_transY = state.at("m_transY").get<std::vector<double>>();
m_iTransS = state.at("m_iTransS").get<std::vector<double>>();
m_iTransL = state.at("m_iTransL").get<std::vector<double>>();
+ m_imageIdentifier = state.at("m_imageIdentifier").get<std::string>();
+ m_platformName = state.at("m_platformName").get<std::string>();
+ m_sensorName = state.at("m_sensorName").get<std::string>();
+ m_collectionIdentifier = state.at("m_collectionIdentifier").get<std::string>();
+ std::vector<double> refpt = state.at("m_referencePointXyz").get<std::vector<double>>();
+ m_referencePointXyz.x = refpt[0];
+ m_referencePointXyz.y = refpt[1];
+ m_referencePointXyz.z = refpt[2];
+ m_currentParameterCovariance = state.at("m_currentParameterCovariance").get<std::vector<double>>();
+
// Leaving unused params commented out
// m_targetName = state.at("m_targetName").get<std::string>();
@@ -791,153 +797,105 @@ void UsgsAstroFrameSensorModel::replaceModelState(const std::string& stringState
std::string UsgsAstroFrameSensorModel::constructStateFromIsd(const std::string& jsonIsd, csm::WarningList* warnings) {
- auto metric_conversion = [](double val, std::string from, std::string to="m") {
- json typemap = {
- {"m", 0},
- {"km", 3}
- };
-
- // everything to lowercase
- std::transform(from.begin(), from.end(), from.begin(), ::tolower);
- std::transform(to.begin(), to.end(), to.begin(), ::tolower);
- return val*pow(10, typemap[from].get<int>() - typemap[to].get<int>());
- };
-
json isd = json::parse(jsonIsd);
json state = {};
+ csm::WarningList* parsingWarnings = new csm::WarningList;
- try {
- state["m_modelName"] = isd.at("name_model");
- std::cerr << "Model Name Parsed!" << std::endl;
-
- state["m_startingDetectorSample"] = isd.at("starting_detector_sample");
- state["m_startingDetectorLine"] = isd.at("starting_detector_line");
-
- std::cerr << "Detector Starting Pixel Parsed!" << std::endl;
-
- // get focal length
- {
- json jayson = isd.at("focal_length_model");
- json focal_length = jayson.at("focal_length");
- json epsilon = jayson.at("focal_epsilon");
-
- state["m_focalLength"] = focal_length;
- state["m_focalLengthEpsilon"] = epsilon;
+ state["m_modelName"] = getSensorModelName(isd, parsingWarnings);
+ state["m_imageIdentifier"] = getImageId(isd, parsingWarnings);
+ state["m_sensorName"] = getSensorName(isd, parsingWarnings);
+ state["m_platformName"] = getPlatformName(isd, parsingWarnings);
- std::cerr << "Focal Length Parsed!" << std::endl;
- }
-
- // get sensor_position
- {
- json jayson = isd.at("sensor_position");
- json positions = jayson.at("positions")[0];
- json velocities = jayson.at("velocities")[0];
- json unit = jayson.at("unit");
-
- unit = unit.get<std::string>();
- state["m_currentParameterValue"] = json();
- state["m_currentParameterValue"][0] = metric_conversion(positions[0].get<double>(), unit);
- state["m_currentParameterValue"][1] = metric_conversion(positions[1].get<double>(), unit);
- state["m_currentParameterValue"][2] = metric_conversion(positions[2].get<double>(), unit);
- state["m_spacecraftVelocity"] = velocities;
-
- std::cerr << "Sensor Location Parsed!" << std::endl;
- }
+ state["m_startingDetectorSample"] = getDetectorStartingSample(isd, parsingWarnings);
+ state["m_startingDetectorLine"] = getDetectorStartingLine(isd, parsingWarnings);
- // get sun_position
- // sun position is not strictly necessary, but is required for getIlluminationDirection.
- {
- json jayson = isd.at("sun_position");
- json positions = jayson.at("positions")[0];
- json unit = jayson.at("unit");
+ // get focal length
+ state["m_focalLength"] = getFocalLength(isd, parsingWarnings);
+ state["m_focalLengthEpsilon"] = getFocalLengthEpsilon(isd, parsingWarnings);
- unit = unit.get<std::string>();
- state["m_sunPosition"] = json();
- state["m_sunPosition"][0] = metric_conversion(positions[0].get<double>(), unit);
- state["m_sunPosition"][1] = metric_conversion(positions[1].get<double>(), unit);
- state["m_sunPosition"][2] = metric_conversion(positions[2].get<double>(), unit);
-
- std::cerr << "Sun Position Parsed!" << std::endl;
- }
- // get sensor_orientation quaternion
- {
- json jayson = isd.at("sensor_orientation");
- json quaternion = jayson.at("quaternions")[0];
+ state["m_currentParameterValue"] = json();
- state["m_currentParameterValue"][3] = quaternion[0];
- state["m_currentParameterValue"][4] = quaternion[1];
- state["m_currentParameterValue"][5] = quaternion[2];
- state["m_currentParameterValue"][6] = quaternion[3];
-
- std::cerr << "Sensor Orientation Parsed!" << std::endl;
- }
-
- // get optical_distortion
- {
- json jayson = isd.at("optical_distortion");
- std::vector<double> xDistortion = jayson.at("transverse").at("x");
- std::vector<double> yDistortion = jayson.at("transverse").at("y");
- xDistortion.resize(10, 0.0);
- yDistortion.resize(10, 0.0);
-
- state["m_odtX"] = xDistortion;
- state["m_odtY"] = yDistortion;
+ // get sensor_position
+ std::vector<double> position = getSensorPositions(isd, parsingWarnings);
+ if (!position.empty() && position.size() != 3) {
+ parsingWarnings->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Sensor position does not have 3 values.",
+ "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+ state["m_currentParameterValue"][0] = 0;
+ state["m_currentParameterValue"][1] = 0;
+ state["m_currentParameterValue"][2] = 0;
+ }
+ else {
+ state["m_currentParameterValue"] = position;
+ }
- std::cerr << "Distortion Parsed!" << std::endl;
- }
+ // get sensor_velocity
+ std::vector<double> velocity = getSensorVelocities(isd, parsingWarnings);
+ if (!velocity.empty() && velocity.size() != 3) {
+ parsingWarnings->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Sensor velocity does not have 3 values.",
+ "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+ }
+ else {
+ state["m_spacecraftVelocity"] = velocity;
+ }
- // get detector_center
- {
- json jayson = isd.at("detector_center");
- json sample = jayson.at("sample");
- json line = jayson.at("line");
- state["m_ccdCenter"][0] = line;
- state["m_ccdCenter"][1] = sample;
+ // get sun_position
+ // sun position is not strictly necessary, but is required for getIlluminationDirection.
+ state["m_sunPosition"] = getSunPositions(isd);
- std::cerr << "Detector Center Parsed!" << std::endl;
- }
+ // get sensor_orientation quaternion
+ std::vector<double> quaternion = getSensorOrientations(isd, parsingWarnings);
+ if (quaternion.size() != 4) {
+ parsingWarnings->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Sensor orientation quaternion does not have 4 values.",
+ "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+ }
+ else {
+ state["m_currentParameterValue"][3] = quaternion[0];
+ state["m_currentParameterValue"][4] = quaternion[1];
+ state["m_currentParameterValue"][5] = quaternion[2];
+ state["m_currentParameterValue"][6] = quaternion[3];
+ }
- // get radii
- {
- json jayson = isd.at("radii");
- json semiminor = jayson.at("semiminor");
- json semimajor = jayson.at("semimajor");
- json unit = jayson.at("unit");
+ // get optical_distortion
+ state["m_distortionType"] = getDistortionModel(isd);
+ state["m_opticalDistCoeffs"] = getDistortionCoeffs(isd);
- unit = unit.get<std::string>();
- state["m_minorAxis"] = metric_conversion(semiminor.get<double>(), unit);
- state["m_majorAxis"] = metric_conversion(semimajor.get<double>(), unit);
+ // get detector_center
+ state["m_ccdCenter"][0] = getDetectorCenterLine(isd, parsingWarnings);
+ state["m_ccdCenter"][1] = getDetectorCenterSample(isd, parsingWarnings);
- std::cerr << "Target Radii Parsed!" << std::endl;
- }
- // get reference_height
- {
- json reference_height = isd.at("reference_height");
- json maxheight = reference_height.at("maxheight");
- json minheight = reference_height.at("minheight");
- json unit = reference_height.at("unit");
+ // get radii
+ state["m_minorAxis"] = getSemiMinorRadius(isd, parsingWarnings);
+ state["m_majorAxis"] = getSemiMajorRadius(isd, parsingWarnings);
- unit = unit.get<std::string>();
- state["m_minElevation"] = metric_conversion(minheight.get<double>(), unit);
- state["m_maxElevation"] = metric_conversion(maxheight.get<double>(), unit);
- std::cerr << "Reference Height Parsed!" << std::endl;
- }
+ // get reference_height
+ state["m_minElevation"] = getMinHeight(isd, parsingWarnings);
+ state["m_maxElevation"] = getMaxHeight(isd, parsingWarnings);
- state["m_ephemerisTime"] = isd.at("center_ephemeris_time");
- state["m_nLines"] = isd.at("image_lines");
- state["m_nSamples"] = isd.at("image_samples");
- state["m_iTransL"] = isd.at("focal2pixel_lines");
+ state["m_ephemerisTime"] = getCenterTime(isd, parsingWarnings);
+ state["m_nLines"] = getTotalLines(isd, parsingWarnings);
+ state["m_nSamples"] = getTotalSamples(isd, parsingWarnings);
- state["m_iTransS"] = isd.at("focal2pixel_samples");
+ state["m_iTransL"] = getFocal2PixelLines(isd, parsingWarnings);
+ state["m_iTransS"] = getFocal2PixelSamples(isd, parsingWarnings);
- // We don't pass the pixel to focal plane transformation so invert the
- // focal plane to pixel transformation
+ // We don't pass the pixel to focal plane transformation so invert the
+ // focal plane to pixel transformation
+ try {
double determinant = state["m_iTransL"][1].get<double>() * state["m_iTransS"][2].get<double>() -
state["m_iTransL"][2].get<double>() * state["m_iTransS"][1].get<double>();
@@ -950,37 +908,47 @@ std::string UsgsAstroFrameSensorModel::constructStateFromIsd(const std::string&
state["m_transY"][2] = state["m_iTransS"][2].get<double>() / determinant;
state["m_transY"][0] = -(state["m_transY"][1].get<double>() * state["m_iTransL"][0].get<double>() +
state["m_transY"][2].get<double>() * state["m_iTransS"][0].get<double>());
-
- std::cerr << "Focal To Pixel Transformation Parsed!" << std::endl;
-
}
- catch(std::out_of_range& e) {
- throw csm::Error(csm::Error::SENSOR_MODEL_NOT_CONSTRUCTIBLE,
- "ISD missing necessary keywords to create sensor model: " + std::string(e.what()),
- "UsgsAstroFrameSensorModel::constructStateFromIsd");
+ catch (...) {
+ parsingWarnings->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not compute detector pixel to focal plane coordinate transformation.",
+ "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
}
- catch(...) {
+
+
+ state["m_referencePointXyz"] = std::vector<double>(3, 0.0);
+ state["m_currentParameterCovariance"] = std::vector<double>(NUM_PARAMETERS*NUM_PARAMETERS,0.0);
+ state["m_collectionIdentifier"] = "";
+
+
+ if (!parsingWarnings->empty()) {
+ if (warnings) {
+ warnings->insert(warnings->end(), parsingWarnings->begin(), parsingWarnings->end());
+ }
+ delete parsingWarnings;
+ parsingWarnings = nullptr;
throw csm::Error(csm::Error::SENSOR_MODEL_NOT_CONSTRUCTIBLE,
"ISD is invalid for creating the sensor model.",
"UsgsAstroFrameSensorModel::constructStateFromIsd");
}
- return state.dump();
+ delete parsingWarnings;
+ parsingWarnings = nullptr;
+ return state.dump();
}
+
csm::EcefCoord UsgsAstroFrameSensorModel::getReferencePoint() const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getReferencePoint");
+ return m_referencePointXyz;
}
void UsgsAstroFrameSensorModel::setReferencePoint(const csm::EcefCoord &groundPt) {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::setReferencePoint");
+ m_referencePointXyz = groundPt;
}
@@ -1006,28 +974,26 @@ std::string UsgsAstroFrameSensorModel::getParameterUnits(int index) const {
bool UsgsAstroFrameSensorModel::hasShareableParameters() const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::hasShareableParameters");
+ return false;
}
bool UsgsAstroFrameSensorModel::isParameterShareable(int index) const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::isParameterShareable");
+ return false;
}
csm::SharingCriteria UsgsAstroFrameSensorModel::getParameterSharingCriteria(int index) const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getParameterSharingCriteria");
+ // Parameter sharing is not supported for this sensor,
+ // all indices are out of range
+ throw csm::Error(
+ csm::Error::INDEX_OUT_OF_RANGE,
+ "Index out of range.",
+ "UsgsAstroLsSensorModel::getParameterSharingCriteria");
}
double UsgsAstroFrameSensorModel::getParameterValue(int index) const {
-
return m_currentParameterValue[index];
}
@@ -1049,46 +1015,48 @@ void UsgsAstroFrameSensorModel::setParameterType(int index, csm::param::Type pTy
double UsgsAstroFrameSensorModel::getParameterCovariance(int index1, int index2) const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getParameterCovariance");
+ int index = UsgsAstroFrameSensorModel::NUM_PARAMETERS * index1 + index2;
+ return m_currentParameterCovariance[index];
}
void UsgsAstroFrameSensorModel::setParameterCovariance(int index1, int index2, double covariance) {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::setParameterCovariance");
+ int index = UsgsAstroFrameSensorModel::NUM_PARAMETERS * index1 + index2;
+ m_currentParameterCovariance[index] = covariance;
}
int UsgsAstroFrameSensorModel::getNumGeometricCorrectionSwitches() const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getNumGeometricCorrectionSwitches");
+ return 0;
}
std::string UsgsAstroFrameSensorModel::getGeometricCorrectionName(int index) const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getGeometricCorrectionName");
+ // Since there are no geometric corrections, all indices are out of range
+ throw csm::Error(
+ csm::Error::INDEX_OUT_OF_RANGE,
+ "Index is out of range.",
+ "UsgsAstroLsSensorModel::getGeometricCorrectionName");
}
void UsgsAstroFrameSensorModel::setGeometricCorrectionSwitch(int index,
bool value,
csm::param::Type pType) {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::setGeometricCorrectionSwitch");
+ // Since there are no geometric corrections, all indices are out of range
+ throw csm::Error(
+ csm::Error::INDEX_OUT_OF_RANGE,
+ "Index is out of range.",
+ "UsgsAstroLsSensorModel::setGeometricCorrectionSwitch");
}
bool UsgsAstroFrameSensorModel::getGeometricCorrectionSwitch(int index) const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getGeometricCorrectionSwitch");
+ // Since there are no geometric corrections, all indices are out of range
+ throw csm::Error(
+ csm::Error::INDEX_OUT_OF_RANGE,
+ "Index is out of range.",
+ "UsgsAstroLsSensorModel::getGeometricCorrectionSwitch");
}
@@ -1096,9 +1064,14 @@ std::vector<double> UsgsAstroFrameSensorModel::getCrossCovarianceMatrix(
const GeometricModel &comparisonModel,
csm::param::Set pSet,
const GeometricModelList &otherModels) const {
- throw csm::Error(csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroFrameSensorModel::getCrossCovarianceMatrix");
+
+ // No correlation between models.
+ const std::vector<int>& indices = getParameterSetIndices(pSet);
+ size_t num_rows = indices.size();
+ const std::vector<int>& indices2 = comparisonModel.getParameterSetIndices(pSet);
+ size_t num_cols = indices.size();
+
+ return std::vector<double>(num_rows * num_cols, 0.0);
}
@@ -1191,180 +1164,6 @@ void UsgsAstroFrameSensorModel::losEllipsoidIntersect(
}
-/**
- * @brief Compute undistorted focal plane x/y.
- *
- * Computes undistorted focal plane (x,y) coordinates given a distorted focal plane (x,y)
- * coordinate. The undistorted coordinates are solved for using the Newton-Raphson
- * method for root-finding if the distortionFunction method is invoked.
- *
- * @param dx distorted focal plane x in millimeters
- * @param dy distorted focal plane y in millimeters
- * @param undistortedX The undistorted x coordinate, in millimeters.
- * @param undistortedY The undistorted y coordinate, in millimeters.
- *
- * @return if the conversion was successful
- * @todo Review the tolerance and maximum iterations of the root-
- * finding algorithm.
- * @todo Review the handling of non-convergence of the root-finding
- * algorithm.
- * @todo Add error handling for near-zero determinant.
-*/
-bool UsgsAstroFrameSensorModel::setFocalPlane(double dx,double dy,
- double &undistortedX,
- double &undistortedY ) const {
-
- // Solve the distortion equation using the Newton-Raphson method.
- // Set the error tolerance to about one millionth of a NAC pixel.
- const double tol = 1.4E-5;
-
- // The maximum number of iterations of the Newton-Raphson method.
- const int maxTries = 60;
-
- double x;
- double y;
- double fx;
- double fy;
- double Jxx;
- double Jxy;
- double Jyx;
- double Jyy;
-
- // Initial guess at the root
- x = dx;
- y = dy;
-
- distortionFunction(x, y, fx, fy);
-
-
- for (int count = 1; ((fabs(fx) +fabs(fy)) > tol) && (count < maxTries); count++) {
-
- this->distortionFunction(x, y, fx, fy);
-
- fx = dx - fx;
- fy = dy - fy;
-
- distortionJacobian(x, y, Jxx, Jxy, Jyx, Jyy);
-
- double determinant = Jxx * Jyy - Jxy * Jyx;
- if (fabs(determinant) < 1E-6) {
-
- undistortedX = x;
- undistortedY = y;
- //
- // Near-zero determinant. Add error handling here.
- //
- //-- Just break out and return with no convergence
- return false;
- }
-
- x = x + (Jyy * fx - Jxy * fy) / determinant;
- y = y + (Jxx * fy - Jyx * fx) / determinant;
- }
-
- if ( (fabs(fx) + fabs(fy)) <= tol) {
- // The method converged to a root.
- undistortedX = x;
- undistortedY = y;
- return true;
-
- }
- else {
- // The method did not converge to a root within the maximum
- // number of iterations. Return with no distortion.
- undistortedX = dx;
- undistortedY = dy;
- return false;
- }
- return true;
-}
-
-
-/**
- * @description Jacobian of the distortion function. The Jacobian was computed
- * algebraically from the function described in the distortionFunction
- * method.
- *
- * @param x
- * @param y
- * @param Jxx Partial_xx
- * @param Jxy Partial_xy
- * @param Jyx Partial_yx
- * @param Jyy Partial_yy
- */
-void UsgsAstroFrameSensorModel::distortionJacobian(double x, double y, double &Jxx, double &Jxy,
- double &Jyx, double &Jyy) const {
-
- double d_dx[10];
- d_dx[0] = 0;
- d_dx[1] = 1;
- d_dx[2] = 0;
- d_dx[3] = 2 * x;
- d_dx[4] = y;
- d_dx[5] = 0;
- d_dx[6] = 3 * x * x;
- d_dx[7] = 2 * x * y;
- d_dx[8] = y * y;
- d_dx[9] = 0;
- double d_dy[10];
- d_dy[0] = 0;
- d_dy[1] = 0;
- d_dy[2] = 1;
- d_dy[3] = 0;
- d_dy[4] = x;
- d_dy[5] = 2 * y;
- d_dy[6] = 0;
- d_dy[7] = x * x;
- d_dy[8] = 2 * x * y;
- d_dy[9] = 3 * y * y;
-
- Jxx = 0.0;
- Jxy = 0.0;
- Jyx = 0.0;
- Jyy = 0.0;
-
- for (int i = 0; i < 10; i++) {
- Jxx = Jxx + d_dx[i] * m_odtX[i];
- Jxy = Jxy + d_dy[i] * m_odtX[i];
- Jyx = Jyx + d_dx[i] * m_odtY[i];
- Jyy = Jyy + d_dy[i] * m_odtY[i];
- }
-}
-
-
-
-/**
- * @description Compute distorted focal plane (dx,dy) coordinate given an undistorted focal
- * plane (ux,uy) coordinate. This describes the third order Taylor approximation to the
- * distortion model.
- *
- * @param ux Undistored x
- * @param uy Undistored y
- * @param dx Result distorted x
- * @param dy Result distorted y
- */
-void UsgsAstroFrameSensorModel::distortionFunction(double ux, double uy, double &dx, double &dy) const {
-
- double f[10];
- f[0] = 1;
- f[1] = ux;
- f[2] = uy;
- f[3] = ux * ux;
- f[4] = ux * uy;
- f[5] = uy * uy;
- f[6] = ux * ux * ux;
- f[7] = ux * ux * uy;
- f[8] = ux * uy * uy;
- f[9] = uy * uy * uy;
-
- dx = 0.0;
- dy = 0.0;
- for (int i = 0; i < 10; i++) {
- dx = dx + f[i] * m_odtX[i];
- dy = dy + f[i] * m_odtY[i];
- }
-}
-
/***** Helper Functions *****/
double UsgsAstroFrameSensorModel::getValue(
diff --git a/src/UsgsAstroLsSensorModel.cpp b/src/UsgsAstroLsSensorModel.cpp
index c3b0be4..ec3d030 100644
--- a/src/UsgsAstroLsSensorModel.cpp
+++ b/src/UsgsAstroLsSensorModel.cpp
@@ -19,6 +19,8 @@
#define USGS_SENSOR_LIBRARY
#include "UsgsAstroLsSensorModel.h"
+#include "Utilities.h"
+#include "Distortion.h"
#include <algorithm>
#include <iostream>
@@ -58,53 +60,45 @@ const std::string UsgsAstroLsSensorModel::PARAMETER_NAME[] =
const std::string UsgsAstroLsSensorModel::_STATE_KEYWORD[] =
{
- "m_sensorModelName",
+ "m_modelName",
"m_imageIdentifier",
- "m_sensorType",
- "m_totalLines",
- "m_totalSamples",
- "m_offsetLines",
- "m_offsetSamples",
+ "m_sensorName",
+ "m_nLines",
+ "m_nSamples",
"m_platformFlag",
- "m_aberrFlag",
- "m_atmrefFlag",
"m_intTimeLines",
"m_intTimeStartTimes",
"m_intTimes",
"m_startingEphemerisTime",
"m_centerEphemerisTime",
"m_detectorSampleSumming",
- "m_startingSample",
+ "m_startingDetectorSample",
+ "m_startingDetectorLine",
"m_ikCode",
- "m_focal",
- "m_isisZDirection",
- "m_opticalDistCoef",
+ "m_focalLength",
+ "m_zDirection",
+ "m_distortionType",
+ "m_opticalDistCoeffs",
"m_iTransS",
"m_iTransL",
"m_detectorSampleOrigin",
"m_detectorLineOrigin",
- "m_detectorLineOffset",
- "m_mountingMatrix",
- "m_semiMajorAxis",
- "m_semiMinorAxis",
- "m_referenceDateAndTime",
+ "m_majorAxis",
+ "m_minorAxis",
"m_platformIdentifier",
"m_sensorIdentifier",
- "m_trajectoryIdentifier",
- "m_collectionIdentifier",
- "m_refElevation",
"m_minElevation",
"m_maxElevation",
"m_dtEphem",
"m_t0Ephem",
"m_dtQuat",
"m_t0Quat",
- "m_numEphem",
+ "m_numPositions",
"m_numQuaternions",
- "m_ephemPts",
- "m_ephemRates",
+ "m_positions",
+ "m_velocities",
"m_quaternions",
- "m_parameterVals",
+ "m_currentParameterValue",
"m_parameterType",
"m_referencePointXyz",
"m_gsd",
@@ -112,7 +106,6 @@ const std::string UsgsAstroLsSensorModel::_STATE_KEYWORD[] =
"m_halfSwath",
"m_halfTime",
"m_covariance",
- "m_imageFlipFlag"
};
const int UsgsAstroLsSensorModel::NUM_PARAM_TYPES = 4;
@@ -141,50 +134,44 @@ void UsgsAstroLsSensorModel::replaceModelState(const std::string &stateString )
int num_paramsSq = num_params * num_params;
m_imageIdentifier = j["m_imageIdentifier"].get<std::string>();
- m_sensorType = j["m_sensorType"];
- m_totalLines = j["m_totalLines"];
- m_totalSamples = j["m_totalSamples"];
- m_offsetLines = j["m_offsetLines"];
- m_offsetSamples = j["m_offsetSamples"];
+ m_sensorName = j["m_sensorName"];
+ m_nLines = j["m_nLines"];
+ m_nSamples = j["m_nSamples"];
m_platformFlag = j["m_platformFlag"];
- m_aberrFlag = j["m_aberrFlag"];
- m_atmRefFlag = j["m_atmRefFlag"];
+
m_intTimeLines = j["m_intTimeLines"].get<std::vector<double>>();
m_intTimeStartTimes = j["m_intTimeStartTimes"].get<std::vector<double>>();
m_intTimes = j["m_intTimes"].get<std::vector<double>>();
+
m_startingEphemerisTime = j["m_startingEphemerisTime"];
m_centerEphemerisTime = j["m_centerEphemerisTime"];
m_detectorSampleSumming = j["m_detectorSampleSumming"];
- m_startingSample = j["m_startingSample"];
+ m_startingSample = j["m_startingDetectorSample"];
m_ikCode = j["m_ikCode"];
- m_focal = j["m_focal"];
- m_isisZDirection = j["m_isisZDirection"];
+
+
+ m_focalLength = j["m_focalLength"];
+ m_zDirection = j["m_zDirection"];
+ m_distortionType = (DistortionType)j["m_distortionType"].get<int>();
+ m_opticalDistCoeffs = j["m_opticalDistCoeffs"].get<std::vector<double>>();
for (int i = 0; i < 3; i++) {
- m_opticalDistCoef[i] = j["m_opticalDistCoef"][i];
m_iTransS[i] = j["m_iTransS"][i];
m_iTransL[i] = j["m_iTransL"][i];
}
m_detectorSampleOrigin = j["m_detectorSampleOrigin"];
m_detectorLineOrigin = j["m_detectorLineOrigin"];
- m_detectorLineOffset = j["m_detectorLineOffset"];
- for (int i = 0; i < 9; i++) {
- m_mountingMatrix[i] = j["m_mountingMatrix"][i];
- }
- m_semiMajorAxis = j["m_semiMajorAxis"];
- m_semiMinorAxis = j["m_semiMinorAxis"];
- m_referenceDateAndTime = j["m_referenceDateAndTime"];
+ m_majorAxis = j["m_majorAxis"];
+ m_minorAxis = j["m_minorAxis"];
m_platformIdentifier = j["m_platformIdentifier"];
m_sensorIdentifier = j["m_sensorIdentifier"];
- m_trajectoryIdentifier = j["m_trajectoryIdentifier"];
- m_collectionIdentifier = j["m_collectionIdentifier"];
- m_refElevation = j["m_refElevation"];
m_minElevation = j["m_minElevation"];
m_maxElevation = j["m_maxElevation"];
m_dtEphem = j["m_dtEphem"];
m_t0Ephem = j["m_t0Ephem"];
m_dtQuat = j["m_dtQuat"];
m_t0Quat = j["m_t0Quat"];
- m_numEphem = j["m_numEphem"];
+ m_numPositions = j["m_numPositions"];
+
m_numQuaternions = j["m_numQuaternions"];
m_referencePointXyz.x = j["m_referencePointXyz"][0];
m_referencePointXyz.y = j["m_referencePointXyz"][1];
@@ -193,13 +180,14 @@ void UsgsAstroLsSensorModel::replaceModelState(const std::string &stateString )
m_flyingHeight = j["m_flyingHeight"];
m_halfSwath = j["m_halfSwath"];
m_halfTime = j["m_halfTime"];
- m_imageFlipFlag = j["m_imageFlipFlag"];
// Vector = is overloaded so explicit get with type required.
- m_ephemPts = j["m_ephemPts"].get<std::vector<double>>();
- m_ephemRates = j["m_ephemRates"].get<std::vector<double>>();
+
+ m_positions = j["m_positions"].get<std::vector<double>>();
+ m_velocities = j["m_velocities"].get<std::vector<double>>();
m_quaternions = j["m_quaternions"].get<std::vector<double>>();
- m_parameterVals = j["m_parameterVals"].get<std::vector<double>>();
+ m_currentParameterValue = j["m_currentParameterValue"].get<std::vector<double>>();
m_covariance = j["m_covariance"].get<std::vector<double>>();
+
for (int i = 0; i < num_params; i++) {
for (int k = 0; k < NUM_PARAM_TYPES; k++) {
if (j["m_parameterType"][i] == PARAM_STRING_ALL[k]) {
@@ -233,11 +221,11 @@ std::string UsgsAstroLsSensorModel::getModelNameFromModelState(
// If model name cannot be determined, return a blank string
std::string model_name;
- if (j.find("m_sensorModelName") != j.end()) {
- model_name = j["m_sensorModelName"];
+ 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_sensorModelName' key in the model state object.";
+ std::string aMessage = "No 'm_modelName' key in the model state object.";
std::string aFunction = "UsgsAstroLsPlugin::getModelNameFromModelState";
csm::Error csmErr(aErrorType, aMessage, aFunction);
throw(csmErr);
@@ -254,15 +242,13 @@ std::string UsgsAstroLsSensorModel::getModelNameFromModelState(
std::string UsgsAstroLsSensorModel::getModelState() const {
json state;
+ state["m_modelName"] = _SENSOR_MODEL_NAME;
+ state["m_startingDetectorSample"] = m_startingSample;
state["m_imageIdentifier"] = m_imageIdentifier;
- state["m_sensorType"] = m_sensorType;
- state["m_totalLines"] = m_totalLines;
- state["m_totalSamples"] = m_totalSamples;
- state["m_offsetLines"] = m_offsetLines;
- state["m_offsetSamples"] = m_offsetSamples;
+ state["m_sensorName"] = m_sensorName;
+ state["m_nLines"] = m_nLines;
+ state["m_nSamples"] = m_nSamples;
state["m_platformFlag"] = m_platformFlag;
- state["m_aberrFlag"] = m_aberrFlag;
- state["m_atmRefFlag"] = m_atmRefFlag;
state["m_intTimeLines"] = m_intTimeLines;
state["m_intTimeStartTimes"] = m_intTimeStartTimes;
state["m_intTimes"] = m_intTimes;
@@ -271,47 +257,41 @@ std::string UsgsAstroLsSensorModel::getModelState() const {
state["m_detectorSampleSumming"] = m_detectorSampleSumming;
state["m_startingSample"] = m_startingSample;
state["m_ikCode"] = m_ikCode;
- state["m_focal"] = m_focal;
- state["m_isisZDirection"] = m_isisZDirection;
- state["m_opticalDistCoef"] = std::vector<double>(m_opticalDistCoef, m_opticalDistCoef+3);
+ state["m_focalLength"] = m_focalLength;
+ state["m_zDirection"] = m_zDirection;
+ state["m_distortionType"] = m_distortionType;
+ state["m_opticalDistCoeffs"] = m_opticalDistCoeffs;
state["m_iTransS"] = std::vector<double>(m_iTransS, m_iTransS+3);
state["m_iTransL"] = std::vector<double>(m_iTransL, m_iTransL+3);
state["m_detectorSampleOrigin"] = m_detectorSampleOrigin;
state["m_detectorLineOrigin"] = m_detectorLineOrigin;
- state["m_detectorLineOffset"] = m_detectorLineOffset;
- state["m_mountingMatrix"] = std::vector<double>(m_mountingMatrix, m_mountingMatrix+9);
- state["m_semiMajorAxis"] = m_semiMajorAxis;
- state["m_semiMinorAxis"] = m_semiMinorAxis;
- state["m_referenceDateAndTime"] = m_referenceDateAndTime;
+ state["m_majorAxis"] = m_majorAxis;
+ state["m_minorAxis"] = m_minorAxis;
state["m_platformIdentifier"] = m_platformIdentifier;
state["m_sensorIdentifier"] = m_sensorIdentifier;
- state["m_trajectoryIdentifier"] = m_trajectoryIdentifier;
- state["m_collectionIdentifier"] = m_collectionIdentifier;
- state["m_refElevation"] = m_refElevation;
state["m_minElevation"] = m_minElevation;
state["m_maxElevation"] = m_maxElevation;
state["m_dtEphem"] = m_dtEphem;
state["m_t0Ephem"] = m_t0Ephem;
state["m_dtQuat"] = m_dtQuat;
state["m_t0Quat"] = m_t0Quat;
- state["m_numEphem"] = m_numEphem;
+ state["m_numPositions"] = m_numPositions;
state["m_numQuaternions"] = m_numQuaternions;
- state["m_ephemPts"] = m_ephemPts;
- state["m_ephemRates"] = m_ephemRates;
+ state["m_positions"] = m_positions;
+ state["m_velocities"] = m_velocities;
state["m_quaternions"] = m_quaternions;
- state["m_parameterVals"] = m_parameterVals;
+ state["m_currentParameterValue"] = m_currentParameterValue;
state["m_parameterType"] = m_parameterType;
state["m_gsd"] = m_gsd;
state["m_flyingHeight"] = m_flyingHeight;
state["m_halfSwath"] = m_halfSwath;
state["m_halfTime"] = m_halfTime;
state["m_covariance"] = m_covariance;
- state["m_imageFlipFlag"] = m_imageFlipFlag;
state["m_referencePointXyz"] = json();
- state["m_referencePointXyz"]["x"] = m_referencePointXyz.x;
- state["m_referencePointXyz"]["y"] = m_referencePointXyz.y;
- state["m_referencePointXyz"]["z"] = m_referencePointXyz.z;
+ state["m_referencePointXyz"][0] = m_referencePointXyz.x;
+ state["m_referencePointXyz"][1] = m_referencePointXyz.y;
+ state["m_referencePointXyz"][2] = m_referencePointXyz.z;
return state.dump();
}
@@ -331,14 +311,10 @@ void UsgsAstroLsSensorModel::reset()
_no_adjustment.assign(UsgsAstroLsSensorModel::NUM_PARAMETERS, 0.0);
m_imageIdentifier = ""; // 1
- m_sensorType = "USGSAstroLineScanner"; // 2
- m_totalLines = 0; // 3
- m_totalSamples = 0; // 4
- m_offsetLines = 0.0; // 7
- m_offsetSamples = 0.0; // 8
+ m_sensorName = "USGSAstroLineScanner"; // 2
+ m_nLines = 0; // 3
+ m_nSamples = 0; // 4
m_platformFlag = 1; // 9
- m_aberrFlag = 0; // 10
- m_atmRefFlag = 0; // 11
m_intTimeLines.clear();
m_intTimeStartTimes.clear();
m_intTimes.clear();
@@ -347,11 +323,10 @@ void UsgsAstroLsSensorModel::reset()
m_detectorSampleSumming = 1.0; // 15
m_startingSample = 1.0; // 16
m_ikCode = -85600; // 17
- m_focal = 350.0; // 18
- m_isisZDirection = 1.0; // 19
- m_opticalDistCoef[0] = 0.0; // 20
- m_opticalDistCoef[1] = 0.0; // 20
- m_opticalDistCoef[2] = 0.0; // 20
+ m_focalLength = 350.0; // 18
+ m_zDirection = 1.0; // 19
+ m_distortionType = DistortionType::RADIAL;
+ m_opticalDistCoeffs = std::vector<double>(3, 0.0);
m_iTransS[0] = 0.0; // 21
m_iTransS[1] = 0.0; // 21
m_iTransS[2] = 150.0; // 21
@@ -360,49 +335,34 @@ void UsgsAstroLsSensorModel::reset()
m_iTransL[2] = 0.0; // 22
m_detectorSampleOrigin = 2500.0; // 23
m_detectorLineOrigin = 0.0; // 24
- m_detectorLineOffset = 0.0; // 25
- m_mountingMatrix[0] = 1.0; // 26
- m_mountingMatrix[1] = 0.0; // 26
- m_mountingMatrix[2] = 0.0; // 26
- m_mountingMatrix[3] = 0.0; // 26
- m_mountingMatrix[4] = 1.0; // 26
- m_mountingMatrix[5] = 0.0; // 26
- m_mountingMatrix[6] = 0.0; // 26
- m_mountingMatrix[7] = 0.0; // 26
- m_mountingMatrix[8] = 1.0; // 26
- m_semiMajorAxis = 3400000.0; // 27
- m_semiMinorAxis = 3350000.0; // 28
- m_referenceDateAndTime = ""; // 30
+ m_majorAxis = 3400000.0; // 27
+ m_minorAxis = 3350000.0; // 28
m_platformIdentifier = ""; // 31
m_sensorIdentifier = ""; // 32
- m_trajectoryIdentifier = ""; // 33
- m_collectionIdentifier = ""; // 33
- m_refElevation = 30; // 34
m_minElevation = -8000.0; // 34
m_maxElevation = 8000.0; // 35
m_dtEphem = 2.0; // 36
m_t0Ephem = -70.0; // 37
m_dtQuat = 0.1; // 38
m_t0Quat = -40.0; // 39
- m_numEphem = 0; // 40
+ m_numPositions = 0; // 40
m_numQuaternions = 0; // 41
- m_ephemPts.clear(); // 42
- m_ephemRates.clear(); // 43
+ m_positions.clear(); // 42
+ m_velocities.clear(); // 43
m_quaternions.clear(); // 44
- m_parameterVals.assign(NUM_PARAMETERS,0.0);
+ m_currentParameterValue.assign(NUM_PARAMETERS,0.0);
m_parameterType.assign(NUM_PARAMETERS,csm::param::REAL);
- m_referencePointXyz.x = 0.0; // 47
- m_referencePointXyz.y = 0.0; // 47
- m_referencePointXyz.z = 0.0; // 47
- m_gsd = 1.0; // 48
- m_flyingHeight = 1000.0; // 49
- m_halfSwath = 1000.0; // 50
- m_halfTime = 10.0; // 51
+ m_referencePointXyz.x = 0.0;
+ m_referencePointXyz.y = 0.0;
+ m_referencePointXyz.z = 0.0;
+ m_gsd = 1.0;
+ m_flyingHeight = 1000.0;
+ m_halfSwath = 1000.0;
+ m_halfTime = 10.0;
m_covariance = std::vector<double>(NUM_PARAMETERS * NUM_PARAMETERS,0.0); // 52
- m_imageFlipFlag = 0; // 53
}
@@ -431,10 +391,10 @@ void UsgsAstroLsSensorModel::updateState()
// This routine will set some parameters not found in the ISD.
// Reference point (image center)
- double lineCtr = m_totalLines / 2.0;
- double sampCtr = m_totalSamples / 2.0;
+ double lineCtr = m_nLines / 2.0;
+ double sampCtr = m_nSamples / 2.0;
csm::ImageCoord ip(lineCtr, sampCtr);
- double refHeight = m_refElevation;
+ double refHeight = 0;
m_referencePointXyz = imageToGround(ip, refHeight);
// Compute ground sample distance
ip.line += 1;
@@ -453,11 +413,11 @@ void UsgsAstroLsSensorModel::updateState()
m_flyingHeight = sqrt(dx*dx + dy*dy + dz*dz);
// Compute half swath
- m_halfSwath = m_gsd * m_totalSamples / 2.0;
+ m_halfSwath = m_gsd * m_nSamples / 2.0;
// Compute half time duration
double fullImageTime = m_intTimeStartTimes.back() - m_intTimeStartTimes.front()
- + m_intTimes.back() * (m_totalLines - m_intTimeLines.back());
+ + m_intTimes.back() * (m_nLines - m_intTimeLines.back());
m_halfTime = fullImageTime / 2.0;
// Parameter covariance, hardcoded accuracy values
@@ -505,55 +465,22 @@ csm::ImageCoord UsgsAstroLsSensorModel::groundToImage(
csm::WarningList* warnings) const
{
// Search for the line, sample coordinate that viewed a given ground point.
- // This method uses an iterative bisection method to search for the image
+ // This method uses an iterative secant method to search for the image
// line.
- //
- // For a given search window, this routine involves projecting the
- // ground point onto the focal plane based on the instrument orientation
- // at the start and end of the search window. Then, it computes the focal
- // plane intersection at a mid-point of the search window. Then, it reduces
- // the search window based on the signs of the intersected line offsets from
- // the center of the ccd. For example, if the line offset is -145 at the
- // start of the window, 10 at the mid point, and 35 at the end of the search
- // window, the window will be reduced to the start of the old window to the
- // middle of the old window.
- //
- // In order to achieve faster convergence, the mid point is calculated
- // using the False Position Method instead of simple bisection. This method
- // uses the zero of the line between the two ends of the search window for
- // the mid point instead of a simple bisection. In most cases, this will
- // converge significantly faster, but it can be slower than simple bisection
- // in some situations.
-
- // Start bisection search on the image lines
- double sampCtr = m_totalSamples / 2.0;
- double firstTime = getImageTime(csm::ImageCoord(0.0, sampCtr));
- double lastTime = getImageTime(csm::ImageCoord(m_totalLines, sampCtr));
- double firstOffset = computeViewingPixel(firstTime, ground_pt, adj).line - 0.5;
- double lastOffset = computeViewingPixel(lastTime, ground_pt, adj).line - 0.5;
-
- // Check if both offsets have the same sign.
- // This means there is not guaranteed to be a zero.
- if ((firstOffset > 0) != (lastOffset < 0)) {
- throw csm::Warning(
- csm::Warning::IMAGE_COORD_OUT_OF_BOUNDS,
- "The image coordinate is out of bounds of the image size.",
- "UsgsAstroLsSensorModel::groundToImage");
- }
// Convert the ground precision to pixel precision so we can
// check for convergence without re-intersecting
csm::ImageCoord approxPoint;
computeLinearApproximation(ground_pt, approxPoint);
csm::ImageCoord approxNextPoint = approxPoint;
- if (approxNextPoint.line + 1 < m_totalLines) {
+ if (approxNextPoint.line + 1 < m_nLines) {
++approxNextPoint.line;
}
else {
--approxNextPoint.line;
}
- csm::EcefCoord approxIntersect = imageToGround(approxPoint, m_refElevation);
- csm::EcefCoord approxNextIntersect = imageToGround(approxNextPoint, m_refElevation);
+ csm::EcefCoord approxIntersect = imageToGround(approxPoint, 0);
+ csm::EcefCoord approxNextIntersect = imageToGround(approxNextPoint, 0);
double lineDX = approxNextIntersect.x - approxIntersect.x;
double lineDY = approxNextIntersect.y - approxIntersect.y;
double lineDZ = approxNextIntersect.z - approxIntersect.z;
@@ -561,32 +488,67 @@ csm::ImageCoord UsgsAstroLsSensorModel::groundToImage(
// Increase the precision by a small amount to ensure the desired precision is met
double pixelPrec = desired_precision / approxLineRes * 0.9;
- // Start bisection search for zero
+ // Start secant method search on the image lines
+ double sampCtr = m_nSamples / 2.0;
+ double firstTime = getImageTime(csm::ImageCoord(0.0, sampCtr));
+ double lastTime = getImageTime(csm::ImageCoord(m_nLines, sampCtr));
+ double firstOffset = computeViewingPixel(firstTime, ground_pt, adj, pixelPrec/2).line - 0.5;
+ double lastOffset = computeViewingPixel(lastTime, ground_pt, adj, pixelPrec/2).line - 0.5;
+
+ // Check if both offsets have the same sign.
+ // This means there is not guaranteed to be a zero.
+ if ((firstOffset > 0) != (lastOffset < 0)) {
+ throw csm::Warning(
+ csm::Warning::IMAGE_COORD_OUT_OF_BOUNDS,
+ "The image coordinate is out of bounds of the image size.",
+ "UsgsAstroLsSensorModel::groundToImage");
+ }
+
+ // Start secant method search
for (int it = 0; it < 30; it++) {
double nextTime = ((firstTime * lastOffset) - (lastTime * firstOffset))
/ (lastOffset - firstOffset);
- double nextOffset = computeViewingPixel(nextTime, ground_pt, adj).line - 0.5;
- // We're looking for a zero, so check that either firstLine and middleLine have
- // opposite signs, or middleLine and lastLine have opposite signs.
- if ((firstOffset > 0) == (nextOffset < 0)) {
- lastTime = nextTime;
- lastOffset = nextOffset;
+ // Because time across the image is not continuous, find the exposure closest
+ // to the computed nextTime and use that.
+
+ // I.E. if the computed nextTime is 0.3, and the middle exposure times for
+ // lines are 0.07, 0.17, 0.27, 0.37, and 0.47; then use 0.27 because it is
+ // the closest middle exposure time.
+ auto referenceTimeIt = std::upper_bound(m_intTimeStartTimes.begin(),
+ m_intTimeStartTimes.end(),
+ nextTime);
+ if (referenceTimeIt != m_intTimeStartTimes.begin()) {
+ --referenceTimeIt;
+ }
+ size_t referenceIndex = std::distance(m_intTimeStartTimes.begin(), referenceTimeIt);
+ double computedLine = (nextTime - m_intTimeStartTimes[referenceIndex]) / m_intTimes[referenceIndex]
+ + m_intTimeLines[referenceIndex] - 0.5; // subtract 0.5 for ISIS -> CSM pixel conversion
+ double closestLine = floor(computedLine + 0.5);
+ nextTime = getImageTime(csm::ImageCoord(closestLine, sampCtr));
+
+ double nextOffset = computeViewingPixel(nextTime, ground_pt, adj, pixelPrec/2).line - 0.5;
+
+ // remove the farthest away node
+ if (fabs(firstTime - nextTime) > fabs(lastTime - nextTime)) {
+ firstTime = nextTime;
+ firstOffset = nextOffset;
}
else {
- firstTime = nextTime;
- firstOffset = nextOffset;
+ lastTime = nextTime;
+ lastOffset = nextOffset;
}
if (fabs(lastOffset - firstOffset) < pixelPrec) {
break;
}
}
- // Check that the desired precision was met
+ // Avoid division by 0 if the first and last nodes are the same
+ double computedTime = firstTime;
+ if (fabs(lastOffset - firstOffset) > 10e-15) {
+ computedTime = ((firstTime * lastOffset) - (lastTime * firstOffset))
+ / (lastOffset - firstOffset);
+ }
- double computedTime = ((firstTime * lastOffset) - (lastTime * firstOffset))
- / (lastOffset - firstOffset);
- csm::ImageCoord calculatedPixel = computeViewingPixel(computedTime, ground_pt, adj);
- // The computed viewing line is the detector line, so we need to convert that to image lines
auto referenceTimeIt = std::upper_bound(m_intTimeStartTimes.begin(),
m_intTimeStartTimes.end(),
computedTime);
@@ -594,33 +556,21 @@ csm::ImageCoord UsgsAstroLsSensorModel::groundToImage(
--referenceTimeIt;
}
size_t referenceIndex = std::distance(m_intTimeStartTimes.begin(), referenceTimeIt);
- calculatedPixel.line += m_intTimeLines[referenceIndex] - 1
- + (computedTime - m_intTimeStartTimes[referenceIndex])
- / m_intTimes[referenceIndex];
- csm::EcefCoord calculatedPoint = imageToGround(calculatedPixel, m_refElevation);
+ double computedLine = (computedTime - m_intTimeStartTimes[referenceIndex]) / m_intTimes[referenceIndex]
+ + m_intTimeLines[referenceIndex] - 0.5; // subtract 0.5 for ISIS -> CSM pixel conversion
+ double closestLine = floor(computedLine + 0.5); // This assumes pixels are the interval [n, n+1)
+ computedTime = getImageTime(csm::ImageCoord(closestLine, sampCtr));
+ csm::ImageCoord calculatedPixel = computeViewingPixel(computedTime, ground_pt, adj, pixelPrec/2);
+ // The computed pioxel is the detector pixel, so we need to convert that to image lines
+ calculatedPixel.line += closestLine;
+
+ // Reintersect to ensure the image point actually views the ground point.
+ csm::EcefCoord calculatedPoint = imageToGround(calculatedPixel, 0);
double dx = ground_pt.x - calculatedPoint.x;
double dy = ground_pt.y - calculatedPoint.y;
double dz = ground_pt.z - calculatedPoint.z;
double len = dx * dx + dy * dy + dz * dz;
- // Check that the pixel is actually in the image
- if ((calculatedPixel.samp < 0) ||
- (calculatedPixel.samp > m_totalSamples)) {
- throw csm::Warning(
- csm::Warning::IMAGE_COORD_OUT_OF_BOUNDS,
- "The image coordinate is out of bounds of the image size.",
- "UsgsAstroLsSensorModel::groundToImage");
- }
-
- // If the final correction is greater than 10 meters,
- // the solution is not valid enough to report even with a warning
- if (len > 100.0) {
- throw csm::Error(
- csm::Error::ALGORITHM,
- "Did not converge.",
- "UsgsAstroLsSensorModel::groundToImage");
- }
-
if (achieved_precision) {
*achieved_precision = sqrt(len);
}
@@ -726,13 +676,6 @@ csm::EcefCoord UsgsAstroLsSensorModel::imageToGround(
losToEcf(
image_pt.line, image_pt.samp, _no_adjustment,
xc, yc, zc, vx, vy, vz, xl, yl, zl);
- if (m_aberrFlag == 1)
- {
- lightAberrationCorr(vx, vy, vz, xl, yl, zl, dxl, dyl, dzl);
- xl += dxl;
- yl += dyl;
- zl += dzl;
- }
double aPrec;
double x, y, z;
@@ -1125,7 +1068,7 @@ void UsgsAstroLsSensorModel::setParameterCovariance(
//***************************************************************************
std::string UsgsAstroLsSensorModel::getTrajectoryIdentifier() const
{
- return m_trajectoryIdentifier;
+ return "UNKNOWN";
}
//***************************************************************************
@@ -1133,14 +1076,10 @@ std::string UsgsAstroLsSensorModel::getTrajectoryIdentifier() const
//***************************************************************************
std::string UsgsAstroLsSensorModel::getReferenceDateAndTime() const
{
- if (m_referenceDateAndTime == "UNKNOWN")
- {
- throw csm::Error(
- csm::Error::UNSUPPORTED_FUNCTION,
- "Unsupported function",
- "UsgsAstroLsSensorModel::getReferenceDateAndTime");
- }
- return m_referenceDateAndTime;
+ throw csm::Error(
+ csm::Error::UNSUPPORTED_FUNCTION,
+ "Unsupported function",
+ "UsgsAstroLsSensorModel::getReferenceDateAndTime");
}
//***************************************************************************
@@ -1155,8 +1094,8 @@ double UsgsAstroLsSensorModel::getImageTime(
// CSM image convention: UL pixel center == (0.5, 0.5)
// USGS image convention: UL pixel center == (1.0, 1.0)
- double lineCSMFull = line1 + m_offsetLines;
- double lineUSGSFull = lineCSMFull + 0.5;
+ double lineCSMFull = line1;
+ double lineUSGSFull = floor(lineCSMFull) + 0.5;
// These calculation assumes that the values in the integration time
// vectors are in terms of ISIS' pixels
@@ -1225,7 +1164,7 @@ csm::EcefVector UsgsAstroLsSensorModel::getSensorVelocity(double time) const
//***************************************************************************
void UsgsAstroLsSensorModel::setParameterValue(int index, double value)
{
- m_parameterVals[index] = value;
+ m_currentParameterValue[index] = value;
}
//***************************************************************************
@@ -1233,7 +1172,7 @@ void UsgsAstroLsSensorModel::setParameterValue(int index, double value)
//***************************************************************************
double UsgsAstroLsSensorModel::getParameterValue(int index) const
{
- return m_parameterVals[index];
+ return m_currentParameterValue[index];
}
//***************************************************************************
@@ -1363,7 +1302,7 @@ csm::ImageCoord UsgsAstroLsSensorModel::getImageStart() const
//***************************************************************************
csm::ImageVector UsgsAstroLsSensorModel::getImageSize() const
{
- return csm::ImageVector(m_totalLines, m_totalSamples);
+ return csm::ImageVector(m_nLines, m_nSamples);
}
//---------------------------------------------------------------------------
@@ -1380,14 +1319,10 @@ csm::ImageVector UsgsAstroLsSensorModel::getImageSize() const
// int num_paramsSq = num_params * num_params;
//
// m_imageIdentifier = j["m_imageIdentifier"];
-// m_sensorType = j["m_sensorType"];
-// m_totalLines = j["m_totalLines"];
-// m_totalSamples = j["m_totalSamples"];
-// m_offsetLines = j["m_offsetLines"];
-// m_offsetSamples = j["m_offsetSamples"];
+// m_sensorName = j["m_sensorName"];
+// m_nLines = j["m_nLines"];
+// m_nSamples = j["m_nSamples"];
// m_platformFlag = j["m_platformFlag"];
-// m_aberrFlag = j["m_aberrFlag"];
-// m_atmRefFlag = j["m_atmrefFlag"];
// m_intTimeLines = j["m_intTimeLines"].get<std::vector<double>>();
// m_intTimeStartTimes = j["m_intTimeStartTimes"].get<std::vector<double>>();
// m_intTimes = j["m_intTimes"].get<std::vector<double>>();
@@ -1396,7 +1331,7 @@ csm::ImageVector UsgsAstroLsSensorModel::getImageSize() const
// m_detectorSampleSumming = j["m_detectorSampleSumming"];
// m_startingSample = j["m_startingSample"];
// m_ikCode = j["m_ikCode"];
-// m_focal = j["m_focal"];
+// m_focalLength = j["m_focalLength"];
// m_isisZDirection = j["m_isisZDirection"];
// for (int i = 0; i < 3; i++) {
// m_opticalDistCoef[i] = j["m_opticalDistCoef"][i];
@@ -1405,25 +1340,20 @@ csm::ImageVector UsgsAstroLsSensorModel::getImageSize() const
// }
// m_detectorSampleOrigin = j["m_detectorSampleOrigin"];
// m_detectorLineOrigin = j["m_detectorLineOrigin"];
-// m_detectorLineOffset = j["m_detectorLineOffset"];
// for (int i = 0; i < 9; i++) {
// m_mountingMatrix[i] = j["m_mountingMatrix"][i];
// }
-// m_semiMajorAxis = j["m_semiMajorAxis"];
-// m_semiMinorAxis = j["m_semiMinorAxis"];
-// m_referenceDateAndTime = j["m_referenceDateAndTime"];
+// m_majorAxis = j["m_majorAxis"];
+// m_minorAxis = j["m_minorAxis"];
// m_platformIdentifier = j["m_platformIdentifier"];
// m_sensorIdentifier = j["m_sensorIdentifier"];
-// m_trajectoryIdentifier = j["m_trajectoryIdentifier"];
-// m_collectionIdentifier = j["m_collectionIdentifier"];
-// m_refElevation = j["m_refElevation"];
// m_minElevation = j["m_minElevation"];
// m_maxElevation = j["m_maxElevation"];
// m_dtEphem = j["m_dtEphem"];
// m_t0Ephem = j["m_t0Ephem"];
// m_dtQuat = j["m_dtQuat"];
// m_t0Quat = j["m_t0Quat"];
-// m_numEphem = j["m_numEphem"];
+// m_numPositions = j["m_numPositions"];
// m_numQuaternions = j["m_numQuaternions"];
// m_referencePointXyz.x = j["m_referencePointXyz"][0];
// m_referencePointXyz.y = j["m_referencePointXyz"][1];
@@ -1432,12 +1362,11 @@ csm::ImageVector UsgsAstroLsSensorModel::getImageSize() const
// m_flyingHeight = j["m_flyingHeight"];
// m_halfSwath = j["m_halfSwath"];
// m_halfTime = j["m_halfTime"];
-// m_imageFlipFlag = j["m_imageFlipFlag"];
// // Vector = is overloaded so explicit get with type required.
-// m_ephemPts = j["m_ephemPts"].get<std::vector<double>>();
-// m_ephemRates = j["m_ephemRates"].get<std::vector<double>>();
+// m_positions = j["m_positions"].get<std::vector<double>>();
+// m_velocities = j["m_velocities"].get<std::vector<double>>();
// m_quaternions = j["m_quaternions"].get<std::vector<double>>();
-// m_parameterVals = j["m_parameterVals"].get<std::vector<double>>();
+// m_currentParameterValue = j["m_currentParameterValue"].get<std::vector<double>>();
// m_covariance = j["m_covariance"].get<std::vector<double>>();
//
// for (int i = 0; i < num_params; i++) {
@@ -1471,7 +1400,7 @@ UsgsAstroLsSensorModel::getValidImageRange() const
{
return std::pair<csm::ImageCoord, csm::ImageCoord>(
csm::ImageCoord(0.0, 0.0),
- csm::ImageCoord(m_totalLines, m_totalSamples));
+ csm::ImageCoord(m_nLines, m_nSamples));
}
//***************************************************************************
@@ -1580,7 +1509,7 @@ std::vector<double> UsgsAstroLsSensorModel::getUnmodeledCrossCovariance(
//***************************************************************************
std::string UsgsAstroLsSensorModel::getCollectionIdentifier() const
{
- return m_collectionIdentifier;
+ return "UNKNOWN";
}
//***************************************************************************
@@ -1642,14 +1571,14 @@ csm::Version UsgsAstroLsSensorModel::getVersion() const
csm::Ellipsoid UsgsAstroLsSensorModel::getEllipsoid() const
{
- return csm::Ellipsoid(m_semiMajorAxis, m_semiMinorAxis);
+ return csm::Ellipsoid(m_majorAxis, m_minorAxis);
}
void UsgsAstroLsSensorModel::setEllipsoid(
const csm::Ellipsoid &ellipsoid)
{
- m_semiMajorAxis = ellipsoid.getSemiMajorRadius();
- m_semiMinorAxis = ellipsoid.getSemiMinorRadius();
+ m_majorAxis = ellipsoid.getSemiMajorRadius();
+ m_minorAxis = ellipsoid.getSemiMinorRadius();
}
@@ -1658,9 +1587,46 @@ double UsgsAstroLsSensorModel::getValue(
int index,
const std::vector<double> &adjustments) const
{
- return m_parameterVals[index] + adjustments[index];
+ return m_currentParameterValue[index] + adjustments[index];
+}
+
+
+//***************************************************************************
+// Functions pulled out of losToEcf and computeViewingPixel
+// **************************************************************************
+
+void UsgsAstroLsSensorModel::getQuaternions(const double& time, double q[4]) const{
+ int nOrder = 8;
+ if (m_platformFlag == 0)
+ nOrder = 4;
+ int nOrderQuat = nOrder;
+ if (m_numQuaternions < 6 && nOrder == 8)
+ nOrderQuat = 4;
+ lagrangeInterp(
+ m_numQuaternions, &m_quaternions[0], m_t0Quat, m_dtQuat, time, 4, nOrderQuat, q);
+}
+
+
+void UsgsAstroLsSensorModel::calculateAttitudeCorrection(
+ const double& time,
+ const std::vector<double>& adj,
+ double attCorr[9]) const
+{
+ double aTime = time - m_t0Quat;
+ double euler[3];
+ double nTime = aTime / m_halfTime;
+ double nTime2 = nTime * nTime;
+ euler[0] =
+ (getValue(6, adj) + getValue(9, adj)* nTime + getValue(12, adj)* nTime2) / m_flyingHeight;
+ euler[1] =
+ (getValue(7, adj) + getValue(10, adj)* nTime + getValue(13, adj)* nTime2) / m_flyingHeight;
+ euler[2] =
+ (getValue(8, adj) + getValue(11, adj)* nTime + getValue(14, adj)* nTime2) / m_halfSwath;
+
+ calculateRotationMatrixFromEuler(euler, attCorr);
}
+
//***************************************************************************
// UsgsAstroLsSensorModel::losToEcf
//***************************************************************************
@@ -1674,161 +1640,67 @@ void UsgsAstroLsSensorModel::losToEcf(
double& vx, // output sensor x velocity
double& vy, // output sensor y velocity
double& vz, // output sensor z velocity
- double& xl, // output line-of-sight x coordinate
- double& yl, // output line-of-sight y coordinate
- double& zl) const // output line-of-sight z coordinate
+ double& bodyLookX, // output line-of-sight x coordinate
+ double& bodyLookY, // output line-of-sight y coordinate
+ double& bodyLookZ) const // output line-of-sight z coordinate
{
//# private_func_description
- // Computes image ray in ecf coordinate system.
-
+ // Computes image ray (look vector) in ecf coordinate system.
// Compute adjusted sensor position and velocity
double time = getImageTime(csm::ImageCoord(line, sample));
getAdjSensorPosVel(time, adj, xc, yc, zc, vx, vy, vz);
// CSM image image convention: UL pixel center == (0.5, 0.5)
// USGS image convention: UL pixel center == (1.0, 1.0)
- double sampleCSMFull = sample + m_offsetSamples;
- double sampleUSGSFull = sampleCSMFull + 0.5;
- double fractionalLine = line - floor(line) - 0.5;
-
- // Compute distorted image coordinates in mm
-
- double isisDetSample = (sampleUSGSFull - 1.0)
- * m_detectorSampleSumming + m_startingSample;
- double m11 = m_iTransL[1];
- double m12 = m_iTransL[2];
- double m21 = m_iTransS[1];
- double m22 = m_iTransS[2];
- double t1 = fractionalLine + m_detectorLineOffset
- - m_detectorLineOrigin - m_iTransL[0];
- double t2 = isisDetSample - m_detectorSampleOrigin - m_iTransS[0];
- double determinant = m11 * m22 - m12 * m21;
- double p11 = m11 / determinant;
- double p12 = -m12 / determinant;
- double p21 = -m21 / determinant;
- double p22 = m22 / determinant;
- double isisNatFocalPlaneX = p11 * t1 + p12 * t2;
- double isisNatFocalPlaneY = p21 * t1 + p22 * t2;
+ double sampleCSMFull = sample;
+ double sampleUSGSFull = sampleCSMFull;
+ double fractionalLine = line - floor(line);
+
+ // Compute distorted image coordinates in mm (sample, line on image (pixels) -> focal plane
+ std::tuple<double, double> natFocalPlane;
+ computeDistortedFocalPlaneCoordinates(fractionalLine, sampleUSGSFull, m_detectorSampleOrigin, m_detectorLineOrigin, m_detectorSampleSumming, m_startingSample, m_iTransS, m_iTransL, natFocalPlane);
// Remove lens distortion
- double isisFocalPlaneX = isisNatFocalPlaneX;
- double isisFocalPlaneY = isisNatFocalPlaneY;
- if (m_opticalDistCoef[0] != 0.0 ||
- m_opticalDistCoef[1] != 0.0 ||
- m_opticalDistCoef[2] != 0.0)
- {
- double rr = isisNatFocalPlaneX * isisNatFocalPlaneX
- + isisNatFocalPlaneY * isisNatFocalPlaneY;
- if (rr > 1.0E-6)
- {
- double dr = m_opticalDistCoef[0] + (rr * (m_opticalDistCoef[1]
- + rr * m_opticalDistCoef[2]));
- isisFocalPlaneX = isisNatFocalPlaneX * (1.0 - dr);
- isisFocalPlaneY = isisNatFocalPlaneY * (1.0 - dr);
- }
- }
+ double undistortedFocalPlaneX, undistortedFocalPlaneY;
+ removeDistortion(std::get<0>(natFocalPlane), std::get<1>(natFocalPlane),
+ undistortedFocalPlaneX, undistortedFocalPlaneY,
+ m_opticalDistCoeffs,
+ DistortionType::RADIAL);
- // Define imaging ray in image space
-
- double losIsis[3];
- losIsis[0] = -isisFocalPlaneX * m_isisZDirection;
- losIsis[1] = -isisFocalPlaneY * m_isisZDirection;
- losIsis[2] = -m_focal * (1.0 - getValue(15, adj) / m_halfSwath);
- double isisMag = sqrt(losIsis[0] * losIsis[0]
- + losIsis[1] * losIsis[1]
- + losIsis[2] * losIsis[2]);
- losIsis[0] /= isisMag;
- losIsis[1] /= isisMag;
- losIsis[2] /= isisMag;
-
- // Apply boresight correction
-
- double losApl[3];
- losApl[0] =
- m_mountingMatrix[0] * losIsis[0]
- + m_mountingMatrix[1] * losIsis[1]
- + m_mountingMatrix[2] * losIsis[2];
- losApl[1] =
- m_mountingMatrix[3] * losIsis[0]
- + m_mountingMatrix[4] * losIsis[1]
- + m_mountingMatrix[5] * losIsis[2];
- losApl[2] =
- m_mountingMatrix[6] * losIsis[0]
- + m_mountingMatrix[7] * losIsis[1]
- + m_mountingMatrix[8] * losIsis[2];
+ // Define imaging ray (look vector) in camera space
+ double cameraLook[3];
+ createCameraLookVector(undistortedFocalPlaneX, undistortedFocalPlaneY, m_zDirection, m_focalLength, cameraLook);
// Apply attitude correction
-
- double aTime = time - m_t0Quat;
- double euler[3];
- double nTime = aTime / m_halfTime;
- double nTime2 = nTime * nTime;
- euler[0] =
- (getValue(6, adj) + getValue(9, adj)* nTime + getValue(12, adj)* nTime2) / m_flyingHeight;
- euler[1] =
- (getValue(7, adj) + getValue(10, adj)* nTime + getValue(13, adj)* nTime2) / m_flyingHeight;
- euler[2] =
- (getValue(8, adj) + getValue(11, adj)* nTime + getValue(14, adj)* nTime2) / m_halfSwath;
- double cos_a = cos(euler[0]);
- double sin_a = sin(euler[0]);
- double cos_b = cos(euler[1]);
- double sin_b = sin(euler[1]);
- double cos_c = cos(euler[2]);
- double sin_c = sin(euler[2]);
- double plFromApl[9];
- plFromApl[0] = cos_b * cos_c;
- plFromApl[1] = -cos_a * sin_c + sin_a * sin_b * cos_c;
- plFromApl[2] = sin_a * sin_c + cos_a * sin_b * cos_c;
- plFromApl[3] = cos_b * sin_c;
- plFromApl[4] = cos_a * cos_c + sin_a * sin_b * sin_c;
- plFromApl[5] = -sin_a * cos_c + cos_a * sin_b * sin_c;
- plFromApl[6] = -sin_b;
- plFromApl[7] = sin_a * cos_b;
- plFromApl[8] = cos_a * cos_b;
-
- double losPl[3];
- losPl[0] = plFromApl[0] * losApl[0] + plFromApl[1] * losApl[1]
- + plFromApl[2] * losApl[2];
- losPl[1] = plFromApl[3] * losApl[0] + plFromApl[4] * losApl[1]
- + plFromApl[5] * losApl[2];
- losPl[2] = plFromApl[6] * losApl[0] + plFromApl[7] * losApl[1]
- + plFromApl[8] * losApl[2];
-
- // Apply rotation matrix from sensor quaternions
-
- int nOrder = 8;
- if (m_platformFlag == 0)
- nOrder = 4;
- int nOrderQuat = nOrder;
- if (m_numQuaternions < 6 && nOrder == 8)
- nOrderQuat = 4;
- double q[4];
- lagrangeInterp(
- m_numQuaternions, &m_quaternions[0], m_t0Quat, m_dtQuat,
- time, 4, nOrderQuat, q);
- double norm = sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]);
- q[0] /= norm;
- q[1] /= norm;
- q[2] /= norm;
- q[3] /= norm;
- double ecfFromPl[9];
- ecfFromPl[0] = q[0] * q[0] - q[1] * q[1] - q[2] * q[2] + q[3] * q[3];
- ecfFromPl[1] = 2 * (q[0] * q[1] - q[2] * q[3]);
- ecfFromPl[2] = 2 * (q[0] * q[2] + q[1] * q[3]);
- ecfFromPl[3] = 2 * (q[0] * q[1] + q[2] * q[3]);
- ecfFromPl[4] = -q[0] * q[0] + q[1] * q[1] - q[2] * q[2] + q[3] * q[3];
- ecfFromPl[5] = 2 * (q[1] * q[2] - q[0] * q[3]);
- ecfFromPl[6] = 2 * (q[0] * q[2] - q[1] * q[3]);
- ecfFromPl[7] = 2 * (q[1] * q[2] + q[0] * q[3]);
- ecfFromPl[8] = -q[0] * q[0] - q[1] * q[1] + q[2] * q[2] + q[3] * q[3];
-
-
- xl = ecfFromPl[0] * losPl[0] + ecfFromPl[1] * losPl[1]
- + ecfFromPl[2] * losPl[2];
- yl = ecfFromPl[3] * losPl[0] + ecfFromPl[4] * losPl[1]
- + ecfFromPl[5] * losPl[2];
- zl = ecfFromPl[6] * losPl[0] + ecfFromPl[7] * losPl[1]
- + ecfFromPl[8] * losPl[2];
+ double attCorr[9];
+ calculateAttitudeCorrection(time, adj, attCorr);
+
+ double correctedCameraLook[3];
+ correctedCameraLook[0] = attCorr[0] * cameraLook[0]
+ + attCorr[1] * cameraLook[1]
+ + attCorr[2] * cameraLook[2];
+ correctedCameraLook[1] = attCorr[3] * cameraLook[0]
+ + attCorr[4] * cameraLook[1]
+ + attCorr[5] * cameraLook[2];
+ correctedCameraLook[2] = attCorr[6] * cameraLook[0]
+ + attCorr[7] * cameraLook[1]
+ + attCorr[8] * cameraLook[2];
+
+// Rotate the look vector into the body fixed frame from the camera reference frame by applying the rotation matrix from the sensor quaternions
+ double quaternions[4];
+ getQuaternions(time, quaternions);
+ double cameraToBody[9];
+ calculateRotationMatrixFromQuaternions(quaternions, cameraToBody);
+
+ bodyLookX = cameraToBody[0] * correctedCameraLook[0]
+ + cameraToBody[1] * correctedCameraLook[1]
+ + cameraToBody[2] * correctedCameraLook[2];
+ bodyLookY = cameraToBody[3] * correctedCameraLook[0]
+ + cameraToBody[4] * correctedCameraLook[1]
+ + cameraToBody[5] * correctedCameraLook[2];
+ bodyLookZ = cameraToBody[6] * correctedCameraLook[0]
+ + cameraToBody[7] * correctedCameraLook[1]
+ + cameraToBody[8] * correctedCameraLook[2];
}
@@ -2028,7 +1900,7 @@ void UsgsAstroLsSensorModel::computeElevation(
// Compute elevation given xyz
// Requires semi-major-axis and eccentricity-square
const int MKTR = 10;
- double ecc_sqr = 1.0 - m_semiMinorAxis * m_semiMinorAxis / m_semiMajorAxis / m_semiMajorAxis;
+ double ecc_sqr = 1.0 - m_minorAxis * m_minorAxis / m_majorAxis / m_majorAxis;
double ep2 = 1.0 - ecc_sqr;
double d2 = x * x + y * y;
double d = sqrt(d2);
@@ -2045,7 +1917,7 @@ void UsgsAstroLsSensorModel::computeElevation(
{
hPrev = h;
tt = tanPhi * tanPhi;
- r = m_semiMajorAxis / sqrt(1.0 + ep2 * tt);
+ r = m_majorAxis / sqrt(1.0 + ep2 * tt);
zz = z + r * ecc_sqr * tanPhi;
n = r * sqrt(1.0 + tt);
h = sqrt(d2 + zz * zz) - n;
@@ -2063,7 +1935,7 @@ void UsgsAstroLsSensorModel::computeElevation(
{
hPrev = h;
cc = cotPhi * cotPhi;
- r = m_semiMajorAxis / sqrt(ep2 + cc);
+ r = m_majorAxis / sqrt(ep2 + cc);
dd = d - r * ecc_sqr * cotPhi;
nn = r * sqrt(1.0 + cc) * ep2;
h = sqrt(dd * dd + z * z) - nn;
@@ -2100,8 +1972,8 @@ void UsgsAstroLsSensorModel::losEllipsoidIntersect(
const int MKTR = 10;
double ap, bp, k;
- ap = m_semiMajorAxis + height;
- bp = m_semiMinorAxis + height;
+ ap = m_majorAxis + height;
+ bp = m_minorAxis + height;
k = ap * ap / (bp * bp);
// Solve quadratic equation for scale factor
@@ -2141,19 +2013,6 @@ void UsgsAstroLsSensorModel::losEllipsoidIntersect(
computeElevation(x, y, z, h, aPrec, desired_precision);
slope = -1;
- while (MKTR > ktr && fabs(height - h) > desired_precision)
- {
- sprev = scale;
- scale += slope * (height - h);
- x = xc + scale * xl;
- y = yc + scale * yl;
- z = zc + scale * zl;
- hprev = h;
- computeElevation(x, y, z, h, aPrec, desired_precision);
- slope = (sprev - scale) / (hprev - h);
- ktr++;
- }
-
achieved_precision = fabs(height - h);
}
@@ -2250,14 +2109,6 @@ void UsgsAstroLsSensorModel::imageToPlane(
losToEcf(line, sample, adj, xc, yc, zc, vx, vy, vz, xl, yl, zl);
- if (m_aberrFlag == 1)
- {
- lightAberrationCorr(vx, vy, vz, xl, yl, zl, dxl, dyl, dzl);
- xl += dxl;
- yl += dyl;
- zl += dzl;
- }
-
losPlaneIntersect(xc, yc, zc, xl, yl, zl, x, y, z, mode);
}
@@ -2279,10 +2130,10 @@ void UsgsAstroLsSensorModel::getAdjSensorPosVel(
if (m_platformFlag == 0)
nOrder = 4;
double sensPosNom[3];
- lagrangeInterp(m_numEphem, &m_ephemPts[0], m_t0Ephem, m_dtEphem,
+ lagrangeInterp(m_numPositions, &m_positions[0], m_t0Ephem, m_dtEphem,
time, 3, nOrder, sensPosNom);
double sensVelNom[3];
- lagrangeInterp(m_numEphem, &m_ephemRates[0], m_t0Ephem, m_dtEphem,
+ lagrangeInterp(m_numPositions, &m_velocities[0], m_t0Ephem, m_dtEphem,
time, 3, nOrder, sensVelNom);
// Compute rotation matrix from ICR to ECF
@@ -2351,8 +2202,12 @@ void UsgsAstroLsSensorModel::getAdjSensorPosVel(
csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel(
const double& time,
const csm::EcefCoord& groundPoint,
- const std::vector<double>& adj) const
+ const std::vector<double>& adj,
+ const double& desiredPrecision) const
{
+ // Helper function to compute the CCD pixel that views a ground point based
+ // on the exterior orientation at a given time.
+
// Get the exterior orientation
double xc, yc, zc, vx, vy, vz;
getAdjSensorPosVel(time, adj, xc, yc, zc, vx, vy, vz);
@@ -2363,69 +2218,27 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel(
double bodyLookZ = groundPoint.z - zc;
// Rotate the look vector into the camera reference frame
- int nOrder = 8;
- if (m_platformFlag == 0)
- nOrder = 4;
- int nOrderQuat = nOrder;
- if (m_numQuaternions < 6 && nOrder == 8)
- nOrderQuat = 4;
- double q[4];
- lagrangeInterp(
- m_numQuaternions, &m_quaternions[0], m_t0Quat, m_dtQuat,
- time, 4, nOrderQuat, q);
- double norm = sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]);
- // Divide by the negative norm for 0 through 2 to invert the quaternion
- q[0] /= -norm;
- q[1] /= -norm;
- q[2] /= -norm;
- q[3] /= norm;
+ double quaternions[4];
+ getQuaternions(time, quaternions);
double bodyToCamera[9];
- bodyToCamera[0] = q[0] * q[0] - q[1] * q[1] - q[2] * q[2] + q[3] * q[3];
- bodyToCamera[1] = 2 * (q[0] * q[1] - q[2] * q[3]);
- bodyToCamera[2] = 2 * (q[0] * q[2] + q[1] * q[3]);
- bodyToCamera[3] = 2 * (q[0] * q[1] + q[2] * q[3]);
- bodyToCamera[4] = -q[0] * q[0] + q[1] * q[1] - q[2] * q[2] + q[3] * q[3];
- bodyToCamera[5] = 2 * (q[1] * q[2] - q[0] * q[3]);
- bodyToCamera[6] = 2 * (q[0] * q[2] - q[1] * q[3]);
- bodyToCamera[7] = 2 * (q[1] * q[2] + q[0] * q[3]);
- bodyToCamera[8] = -q[0] * q[0] - q[1] * q[1] + q[2] * q[2] + q[3] * q[3];
+ calculateRotationMatrixFromQuaternions(quaternions, bodyToCamera);
+
+ // Apply transpose of matrix to rotate body->camera
double cameraLookX = bodyToCamera[0] * bodyLookX
- + bodyToCamera[1] * bodyLookY
- + bodyToCamera[2] * bodyLookZ;
- double cameraLookY = bodyToCamera[3] * bodyLookX
+ + bodyToCamera[3] * bodyLookY
+ + bodyToCamera[6] * bodyLookZ;
+ double cameraLookY = bodyToCamera[1] * bodyLookX
+ bodyToCamera[4] * bodyLookY
- + bodyToCamera[5] * bodyLookZ;
- double cameraLookZ = bodyToCamera[6] * bodyLookX
- + bodyToCamera[7] * bodyLookY
+ + bodyToCamera[7] * bodyLookZ;
+ double cameraLookZ = bodyToCamera[2] * bodyLookX
+ + bodyToCamera[5] * bodyLookY
+ bodyToCamera[8] * bodyLookZ;
// Invert the attitude correction
- double aTime = time - m_t0Quat;
- double euler[3];
- double nTime = aTime / m_halfTime;
- double nTime2 = nTime * nTime;
- euler[0] =
- (getValue(6, adj) + getValue(9, adj)* nTime + getValue(12, adj)* nTime2) / m_flyingHeight;
- euler[1] =
- (getValue(7, adj) + getValue(10, adj)* nTime + getValue(13, adj)* nTime2) / m_flyingHeight;
- euler[2] =
- (getValue(8, adj) + getValue(11, adj)* nTime + getValue(14, adj)* nTime2) / m_halfSwath;
- double cos_a = cos(euler[0]);
- double sin_a = sin(euler[0]);
- double cos_b = cos(euler[1]);
- double sin_b = sin(euler[1]);
- double cos_c = cos(euler[2]);
- double sin_c = sin(euler[2]);
double attCorr[9];
- attCorr[0] = cos_b * cos_c;
- attCorr[1] = -cos_a * sin_c + sin_a * sin_b * cos_c;
- attCorr[2] = sin_a * sin_c + cos_a * sin_b * cos_c;
- attCorr[3] = cos_b * sin_c;
- attCorr[4] = cos_a * cos_c + sin_a * sin_b * sin_c;
- attCorr[5] = -sin_a * cos_c + cos_a * sin_b * sin_c;
- attCorr[6] = -sin_b;
- attCorr[7] = sin_a * cos_b;
- attCorr[8] = cos_a * cos_b;
+ calculateAttitudeCorrection(time, adj, attCorr);
+
+ // Apply transpose of matrix to invert the attidue correction
double adjustedLookX = attCorr[0] * cameraLookX
+ attCorr[3] * cameraLookY
+ attCorr[6] * cameraLookZ;
@@ -2436,40 +2249,29 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel(
+ attCorr[5] * cameraLookY
+ attCorr[8] * cameraLookZ;
- // Invert the boresight correction
- double correctedLookX = m_mountingMatrix[0] * adjustedLookX
- + m_mountingMatrix[3] * adjustedLookY
- + m_mountingMatrix[6] * adjustedLookZ;
- double correctedLookY = m_mountingMatrix[1] * adjustedLookX
- + m_mountingMatrix[4] * adjustedLookY
- + m_mountingMatrix[7] * adjustedLookZ;
- double correctedLookZ = m_mountingMatrix[2] * adjustedLookX
- + m_mountingMatrix[5] * adjustedLookY
- + m_mountingMatrix[8] * adjustedLookZ;
-
// Convert to focal plane coordinate
- double lookScale = m_focal / correctedLookZ;
- double focalX = correctedLookX * lookScale;
- double focalY = correctedLookY * lookScale;
+ double lookScale = m_focalLength / adjustedLookZ;
+ double focalX = adjustedLookX * lookScale;
+ double focalY = adjustedLookY * lookScale;
+ double distortedFocalX, distortedFocalY;
- // TODO invert distortion here
- // We probably only want to try and invert the distortion if we are
- // reasonably close to the actual CCD because the distortion equations are
- // sometimes only well behaved close to the CCD.
+ // Invert distortion
+ applyDistortion(focalX, focalY, distortedFocalX, distortedFocalY,
+ m_opticalDistCoeffs, m_distortionType, desiredPrecision);
// Convert to detector line and sample
double detectorLine = m_iTransL[0]
- + m_iTransL[1] * focalX
- + m_iTransL[2] * focalY;
+ + m_iTransL[1] * distortedFocalX
+ + m_iTransL[2] * distortedFocalY;
double detectorSample = m_iTransS[0]
- + m_iTransS[1] * focalX
- + m_iTransS[2] * focalY;
+ + m_iTransS[1] * distortedFocalX
+ + m_iTransS[2] * distortedFocalY;
// Convert to image sample line
- double line = detectorLine + m_detectorLineOrigin - m_detectorLineOffset
- - m_offsetLines + 0.5;
+ double line = detectorLine + m_detectorLineOrigin;
double sample = (detectorSample + m_detectorSampleOrigin - m_startingSample)
- / m_detectorSampleSumming - m_offsetSamples + 0.5;
+ / m_detectorSampleSumming;
+
return csm::ImageCoord(line, sample);
}
@@ -2488,8 +2290,8 @@ void UsgsAstroLsSensorModel::computeLinearApproximation(
// Since this is valid only over image,
// don't let result go beyond the image border.
- double numRows = m_totalLines;
- double numCols = m_totalSamples;
+ double numRows = m_nLines;
+ double numCols = m_nSamples;
if (ip.line < 0.0) ip.line = 0.0;
if (ip.line > numRows) ip.line = numRows;
@@ -2498,8 +2300,8 @@ void UsgsAstroLsSensorModel::computeLinearApproximation(
}
else
{
- ip.line = m_totalLines / 2.0;
- ip.samp = m_totalSamples / 2.0;
+ ip.line = m_nLines / 2.0;
+ ip.samp = m_nSamples / 2.0;
}
}
@@ -2524,8 +2326,8 @@ void UsgsAstroLsSensorModel::setLinearApproximation()
}
- double numRows = m_totalLines;
- double numCols = m_totalSamples;
+ double numRows = m_nLines;
+ double numCols = m_nSamples;
csm::ImageCoord imagePt;
csm::EcefCoord gp[8];
@@ -2645,129 +2447,134 @@ double UsgsAstroLsSensorModel::determinant3x3(double mat[9]) const
std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imageSupportData, csm::WarningList *warnings) const
{
- // Instantiate UsgsAstroLineScanner sensor model
- json isd = json::parse(imageSupportData);
- json state;
-
- int num_params = NUM_PARAMETERS;
-
- state["m_imageIdentifier"] = "UNKNOWN";
- state["m_sensorType"] = "LINE_SCAN";
- state["m_totalLines"] = isd.at("image_lines");
- state["m_totalSamples"] = isd.at("image_samples");
- state["m_offsetLines"] = 0.0;
- state["m_offsetSamples"] = 0.0;
- state["m_platformFlag"] = 1;
- state["m_aberrFlag"] = 0;
- state["m_atmRefFlag"] = 0;
- state["m_centerEphemerisTime"] = isd.at("center_ephemeris_time");
- state["m_startingEphemerisTime"] = isd.at("starting_ephemeris_time");
-
- for (auto& scanRate : isd.at("line_scan_rate")) {
- state["m_intTimeLines"].push_back(scanRate[0]);
- state["m_intTimeStartTimes"].push_back(scanRate[1]);
- state["m_intTimes"].push_back(scanRate[2]);
- }
-
- state["m_detectorSampleSumming"] = isd.at("detector_sample_summing");
- state["m_startingSample"] = isd.at("detector_line_summing");
- state["m_ikCode"] = 0;
- state["m_focal"] = isd.at("focal_length_model").at("focal_length");
- state["m_isisZDirection"] = 1;
- state["m_opticalDistCoef"] = isd.at("optical_distortion").at("radial").at("coefficients");
- state["m_iTransS"] = isd.at("focal2pixel_samples");
- state["m_iTransL"] = isd.at("focal2pixel_lines");
-
- state["m_detectorSampleOrigin"] = isd.at("detector_center").at("sample");
- state["m_detectorLineOrigin"] = isd.at("detector_center").at("line");
- state["m_detectorLineOffset"] = 0;
-
- double cos_a = cos(0);
- double sin_a = sin(0);
- double cos_b = cos(0);
- double sin_b = sin(0);
- double cos_c = cos(0);
- double sin_c = sin(0);
-
- state["m_mountingMatrix"] = json::array();
- state["m_mountingMatrix"][0] = cos_b * cos_c;
- state["m_mountingMatrix"][1] = -cos_a * sin_c + sin_a * sin_b * cos_c;
- state["m_mountingMatrix"][2] = sin_a * sin_c + cos_a * sin_b * cos_c;
- state["m_mountingMatrix"][3] = cos_b * sin_c;
- state["m_mountingMatrix"][4] = cos_a * cos_c + sin_a * sin_b * sin_c;
- state["m_mountingMatrix"][5] = -sin_a * cos_c + cos_a * sin_b * sin_c;
- state["m_mountingMatrix"][6] = -sin_b;
- state["m_mountingMatrix"][7] = sin_a * cos_b;
- state["m_mountingMatrix"][8] = cos_a * cos_b;
-
- state["m_dtEphem"] = isd.at("dt_ephemeris");
- state["m_t0Ephem"] = isd.at("t0_ephemeris");
- state["m_dtQuat"] = isd.at("dt_quaternion");
- state["m_t0Quat"] = isd.at("t0_quaternion");
-
- state["m_numEphem"] = isd.at("sensor_position").at("positions").size();
- state["m_numQuaternions"] = isd.at("sensor_orientation").at("quaternions").size();
-
- for (auto& location : isd.at("sensor_position").at("positions")) {
- state["m_ephemPts"].push_back(location[0]);
- state["m_ephemPts"].push_back(location[1]);
- state["m_ephemPts"].push_back(location[2]);
- }
-
- for (auto& velocity : isd.at("sensor_position").at("velocities")) {
- state["m_ephemRates"].push_back(velocity[0]);
- state["m_ephemRates"].push_back(velocity[1]);
- state["m_ephemRates"].push_back(velocity[2]);
- }
-
- for (auto& quaternion : isd.at("sensor_orientation").at("quaternions")) {
- state["m_quaternions"].push_back(quaternion[0]);
- state["m_quaternions"].push_back(quaternion[1]);
- state["m_quaternions"].push_back(quaternion[2]);
- state["m_quaternions"].push_back(quaternion[3]);
- }
-
-
- state["m_parameterVals"] = std::vector<double>(NUM_PARAMETERS, 0.0);
- state["m_parameterVals"][15] = state["m_focal"].get<double>();
-
- // Set the ellipsoid
- state["m_semiMajorAxis"] = isd.at("radii").at("semimajor");
- state["m_semiMinorAxis"] = isd.at("radii").at("semiminor");
+ // Instantiate UsgsAstroLineScanner sensor model
+ json isd = json::parse(imageSupportData);
+ json state = {};
+
+ csm::WarningList* parsingWarnings = new csm::WarningList;
+
+ int num_params = NUM_PARAMETERS;
+
+ state["m_modelName"] = getSensorModelName(isd, parsingWarnings);
+ state["m_imageIdentifier"] = getImageId(isd, parsingWarnings);
+ state["m_sensorName"] = getSensorName(isd, parsingWarnings);
+ state["m_platformName"] = getPlatformName(isd, parsingWarnings);
+
+ state["m_focalLength"] = getFocalLength(isd, parsingWarnings);
+
+ state["m_nLines"] = getTotalLines(isd, parsingWarnings);
+ state["m_nSamples"] = getTotalSamples(isd, parsingWarnings);
+
+ state["m_iTransS"] = getFocal2PixelSamples(isd, parsingWarnings);
+ state["m_iTransL"] = getFocal2PixelLines(isd, parsingWarnings);
+
+ state["m_platformFlag"] = 1;
+ state["m_ikCode"] = 0;
+ state["m_zDirection"] = 1;
+
+ state["m_distortionType"] = getDistortionModel(isd, parsingWarnings);
+ state["m_opticalDistCoeffs"] = getDistortionCoeffs(isd, parsingWarnings);
+
+ // Zero computed state values
+ state["m_referencePointXyz"] = std::vector<double>(3, 0.0);
+
+ // leave these be for now.
+ state["m_gsd"] = 1.0;
+ state["m_flyingHeight"] = 1000.0;
+ state["m_halfSwath"] = 1000.0;
+ state["m_halfTime"] = 10.0;
+
+ state["m_centerEphemerisTime"] = getCenterTime(isd, parsingWarnings);
+ state["m_startingEphemerisTime"] = getStartingTime(isd, parsingWarnings);
+
+ state["m_intTimeLines"] = getIntegrationStartLines(isd, parsingWarnings);
+ state["m_intTimeStartTimes"] = getIntegrationStartTimes(isd, parsingWarnings);
+ state["m_intTimes"] = getIntegrationTimes(isd, parsingWarnings);
+
+ state["m_detectorSampleSumming"] = getSampleSumming(isd, parsingWarnings);
+ state["m_startingDetectorSample"] = getDetectorStartingSample(isd, parsingWarnings);
+ state["m_startingDetectorLine"] = getDetectorStartingLine(isd, parsingWarnings);
+ state["m_detectorSampleOrigin"] = getDetectorCenterSample(isd, parsingWarnings);
+ state["m_detectorLineOrigin"] = getDetectorCenterLine(isd, parsingWarnings);
+
+ // These are exlusive to LineScanners, leave them here for now.
+ try {
+ state["m_dtEphem"] = isd.at("dt_ephemeris");
+ }
+ catch(...) {
+ parsingWarnings->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "dt_ephemeris not in ISD",
+ "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+ }
+
+ try {
+ state["m_t0Ephem"] = isd.at("t0_ephemeris");
+ }
+ catch(...) {
+ parsingWarnings->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "t0_ephemeris not in ISD",
+ "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+ }
+
+ try{
+ state["m_dtQuat"] = isd.at("dt_quaternion");
+ }
+ catch(...) {
+ parsingWarnings->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "dt_quaternion not in ISD",
+ "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+ }
+
+ try{
+ state["m_t0Quat"] = isd.at("t0_quaternion");
+ }
+ catch(...) {
+ parsingWarnings->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "t0_quaternion not in ISD",
+ "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+ }
+
+ std::vector<double> positions = getSensorPositions(isd, parsingWarnings);
+ state["m_positions"] = positions;
+ state["m_numPositions"] = positions.size();
+
+ state["m_velocities"] = getSensorVelocities(isd, parsingWarnings);
+
+ std::vector<double> quaternions = getSensorOrientations(isd, parsingWarnings);
+ state["m_quaternions"] = quaternions;
+ state["m_numQuaternions"] = quaternions.size();
+
+ state["m_currentParameterValue"] = std::vector<double>(NUM_PARAMETERS, 0.0);
+
+ // get radii
+ state["m_minorAxis"] = getSemiMinorRadius(isd, parsingWarnings);
+ state["m_majorAxis"] = getSemiMajorRadius(isd, parsingWarnings);
+
+ // set identifiers
+ state["m_platformIdentifier"] = getPlatformName(isd, parsingWarnings);
+ state["m_sensorIdentifier"] = getSensorName(isd, parsingWarnings);
+
+ // get reference_height
+ state["m_minElevation"] = getMinHeight(isd, parsingWarnings);
+ state["m_maxElevation"] = getMaxHeight(isd, parsingWarnings);
+
+
+ // Default to identity covariance
+ state["m_covariance"] =
+ std::vector<double>(NUM_PARAMETERS * NUM_PARAMETERS, 0.0);
+ for (int i = 0; i < NUM_PARAMETERS; i++) {
+ state["m_covariance"][i * NUM_PARAMETERS + i] = 1.0;
+ }
- // Now finish setting the state data from the ISD read in
-
- // set identifiers
- state["m_referenceDateAndTime"] = "UNKNOWN";
- state["m_platformIdentifier"] = isd.at("name_platform");
- state["m_sensorIdentifier"] = isd.at("name_sensor");
- state["m_trajectoryIdentifier"] = "UNKNOWN";
- state["m_collectionIdentifier"] = "UNKNOWN";
-
- // Ground elevations
- state["m_refElevation"] = 0.0;
- state["m_minElevation"] = isd.at("reference_height").at("minheight");
- state["m_maxElevation"] = isd.at("reference_height").at("maxheight");
-
- // Zero ateter values
- for (int i = 0; i < NUM_PARAMETERS; i++) {
- state["m_ateterType"][i] = csm::param::REAL;
- }
-
- // Default to identity covariance
- state["m_covariance"] =
- std::vector<double>(NUM_PARAMETERS * NUM_PARAMETERS, 0.0);
- for (int i = 0; i < NUM_PARAMETERS; i++) {
- state["m_covariance"][i * NUM_PARAMETERS + i] = 1.0;
- }
- // Zero computed state values
- state["m_referencePointXyz"] = std::vector<double>(3, 0.0);
- state["m_gsd"] = 1.0;
- state["m_flyingHeight"] = 1000.0;
- state["m_halfSwath"] = 1000.0;
- state["m_halfTime"] = 10.0;
- state["m_imageFlipFlag"] = 0;
// The state data will still be updated when a sensor model is created since
// some state data is notin the ISD and requires a SM to compute them.
return state.dump();
diff --git a/src/UsgsAstroPlugin.cpp b/src/UsgsAstroPlugin.cpp
index fc67553..563f24d 100644
--- a/src/UsgsAstroPlugin.cpp
+++ b/src/UsgsAstroPlugin.cpp
@@ -28,9 +28,8 @@ using json = nlohmann::json;
// Declaration of static variables
const std::string UsgsAstroPlugin::_PLUGIN_NAME = "UsgsAstroPluginCSM";
const std::string UsgsAstroPlugin::_MANUFACTURER_NAME = "UsgsAstrogeology";
-const std::string UsgsAstroPlugin::_RELEASE_DATE = "20170425";
+const std::string UsgsAstroPlugin::_RELEASE_DATE = "20190222";
const int UsgsAstroPlugin::_N_SENSOR_MODELS = 2;
-const int UsgsAstroPlugin::_NUM_ISD_KEYWORDS = 21;
const std::string UsgsAstroPlugin::_ISD_KEYWORD[] =
{
@@ -58,12 +57,6 @@ const std::string UsgsAstroPlugin::_ISD_KEYWORD[] =
"sun_position"
};
-// const json UsgsAstroPlugin::MODEL_KEYWORDS = {
-// {UsgsAstroFrameSensorModel::_SENSOR_MODEL_NAME, UsgsAstroFrameSensorModel::_STATE_KEYWORD},
-// {UsgsAstroLsSensorModel::_SENSOR_MODEL_NAME, UsgsAstroLsSensorModel::_STATE_KEYWORD}
-// };
-
-
// Static Instance of itself
const UsgsAstroPlugin UsgsAstroPlugin::m_registeredPlugin;
@@ -154,12 +147,15 @@ std::string UsgsAstroPlugin::loadImageSupportData(const csm::Isd &imageSupportDa
std::string imageFilename = imageSupportDataOriginal.filename();
size_t lastIndex = imageFilename.find_last_of(".");
std::string baseName = imageFilename.substr(0, lastIndex);
+ lastIndex = baseName.find_last_of(DIR_DELIMITER_STR);
+ std::string filename = baseName.substr(lastIndex + 1);
std::string isdFilename = baseName.append(".json");
try {
std::ifstream isd_sidecar(isdFilename);
json jsonisd;
isd_sidecar >> jsonisd;
+ jsonisd["image_identifier"] = filename;
return jsonisd.dump();
} catch (...) {
@@ -205,7 +201,7 @@ bool UsgsAstroPlugin::canISDBeConvertedToModelState(const csm::Isd &imageSupport
std::string UsgsAstroPlugin::getStateFromISD(csm::Isd imageSupportData) const {
std::string stringIsd = loadImageSupportData(imageSupportData);
json jsonIsd = json::parse(stringIsd);
- return convertISDToModelState(imageSupportData, jsonIsd.at("modelName"));
+ return convertISDToModelState(imageSupportData, jsonIsd.at("name_model"));
}
@@ -221,7 +217,6 @@ std::string UsgsAstroPlugin::convertISDToModelState(const csm::Isd &imageSupport
csm::Model *UsgsAstroPlugin::constructModelFromISD(const csm::Isd &imageSupportDataOriginal,
const std::string &modelName,
csm::WarningList *warnings) const {
-
std::string stringIsd = loadImageSupportData(imageSupportDataOriginal);
if (modelName == UsgsAstroFrameSensorModel::_SENSOR_MODEL_NAME) {
@@ -263,7 +258,7 @@ csm::Model *UsgsAstroPlugin::constructModelFromState(const std::string& modelSta
csm::WarningList *warnings) const {
json state = json::parse(modelState);
- std::string modelName = state["modelName"];
+ std::string modelName = state["m_modelName"];
if (modelName == UsgsAstroFrameSensorModel::_SENSOR_MODEL_NAME) {
UsgsAstroFrameSensorModel* model = new UsgsAstroFrameSensorModel();
diff --git a/src/Utilities.cpp b/src/Utilities.cpp
new file mode 100644
index 0000000..1812d66
--- /dev/null
+++ b/src/Utilities.cpp
@@ -0,0 +1,729 @@
+#include "Utilities.h"
+
+using json = nlohmann::json;
+
+// Calculates a rotation matrix from Euler angles
+void calculateRotationMatrixFromEuler(
+ double euler[],
+ double rotationMatrix[])
+{
+ double cos_a = cos(euler[0]);
+ double sin_a = sin(euler[0]);
+ double cos_b = cos(euler[1]);
+ double sin_b = sin(euler[1]);
+ double cos_c = cos(euler[2]);
+ double sin_c = sin(euler[2]);
+
+ rotationMatrix[0] = cos_b * cos_c;
+ rotationMatrix[1] = -cos_a * sin_c + sin_a * sin_b * cos_c;
+ rotationMatrix[2] = sin_a * sin_c + cos_a * sin_b * cos_c;
+ rotationMatrix[3] = cos_b * sin_c;
+ rotationMatrix[4] = cos_a * cos_c + sin_a * sin_b * sin_c;
+ rotationMatrix[5] = -sin_a * cos_c + cos_a * sin_b * sin_c;
+ rotationMatrix[6] = -sin_b;
+ rotationMatrix[7] = sin_a * cos_b;
+ rotationMatrix[8] = cos_a * cos_b;
+}
+
+
+// uses a quaternion to calclate a rotation matrix.
+void calculateRotationMatrixFromQuaternions(
+ double q[4],
+ double rotationMatrix[9])
+{
+ double norm = sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]);
+ q[0] /= norm;
+ q[1] /= norm;
+ q[2] /= norm;
+ q[3] /= norm;
+
+ rotationMatrix[0] = q[0] * q[0] - q[1] * q[1] - q[2] * q[2] + q[3] * q[3];
+ rotationMatrix[1] = 2 * (q[0] * q[1] - q[2] * q[3]);
+ rotationMatrix[2] = 2 * (q[0] * q[2] + q[1] * q[3]);
+ rotationMatrix[3] = 2 * (q[0] * q[1] + q[2] * q[3]);
+ rotationMatrix[4] = -q[0] * q[0] + q[1] * q[1] - q[2] * q[2] + q[3] * q[3];
+ rotationMatrix[5] = 2 * (q[1] * q[2] - q[0] * q[3]);
+ rotationMatrix[6] = 2 * (q[0] * q[2] - q[1] * q[3]);
+ rotationMatrix[7] = 2 * (q[1] * q[2] + q[0] * q[3]);
+ rotationMatrix[8] = -q[0] * q[0] - q[1] * q[1] + q[2] * q[2] + q[3] * q[3];
+}
+
+void computeDistortedFocalPlaneCoordinates(
+ const double& line,
+ const double& sample,
+ const double& sampleOrigin,
+ const double& lineOrigin,
+ const double& sampleSumming,
+ const double& startingSample,
+ const double iTransS[],
+ const double iTransL[],
+ std::tuple<double, double>& natFocalPlane)
+{
+ double detSample = sample * sampleSumming + startingSample;
+ double m11 = iTransL[1];
+ double m12 = iTransL[2];
+ double m21 = iTransS[1];
+ double m22 = iTransS[2];
+ double t1 = line - lineOrigin - iTransL[0];
+ double t2 = detSample - sampleOrigin - iTransS[0];
+ double determinant = m11 * m22 - m12 * m21;
+ double p11 = m11 / determinant;
+ double p12 = -m12 / determinant;
+ double p21 = -m21 / determinant;
+ double p22 = m22 / determinant;
+
+ std::get<0>(natFocalPlane) = p11 * t1 + p12 * t2;
+ std::get<1>(natFocalPlane) = p21 * t1 + p22 * t2;
+};
+
+// Define imaging ray in image space (In other words, create a look vector in camera space)
+void createCameraLookVector(
+ const double& undistortedFocalPlaneX,
+ const double& undistortedFocalPlaneY,
+ const double& zDirection,
+ const double& focalLength,
+ double cameraLook[])
+{
+ cameraLook[0] = -undistortedFocalPlaneX * zDirection;
+ cameraLook[1] = -undistortedFocalPlaneY * zDirection;
+ cameraLook[2] = -focalLength;
+ double magnitude = sqrt(cameraLook[0] * cameraLook[0]
+ + cameraLook[1] * cameraLook[1]
+ + cameraLook[2] * cameraLook[2]);
+ cameraLook[0] /= magnitude;
+ cameraLook[1] /= magnitude;
+ cameraLook[2] /= magnitude;
+}
+
+double metric_conversion(double val, std::string from, std::string to) {
+ json typemap = {
+ {"m", 0},
+ {"km", 3}
+ };
+
+ // everything to lowercase
+ std::transform(from.begin(), from.end(), from.begin(), ::tolower);
+ std::transform(to.begin(), to.end(), to.begin(), ::tolower);
+ return val*pow(10, typemap[from].get<int>() - typemap[to].get<int>());
+}
+
+std::string getSensorModelName(json isd, csm::WarningList *list) {
+ std::string name = "";
+ try {
+ name = isd.at("name_model");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the sensor model name.",
+ "Utilities::getSensorModelName()"));
+ }
+ }
+ return name;
+}
+
+std::string getImageId(json isd, csm::WarningList *list) {
+ std::string id = "";
+ try {
+ id = isd.at("image_identifier");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the image identifier.",
+ "Utilities::getImageId()"));
+ }
+ }
+ return id;
+}
+
+std::string getSensorName(json isd, csm::WarningList *list) {
+ std::string name = "";
+ try {
+ name = isd.at("name_sensor");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the sensor name.",
+ "Utilities::getSensorName()"));
+ }
+ }
+ return name;
+}
+
+std::string getPlatformName(json isd, csm::WarningList *list) {
+ std::string name = "";
+ try {
+ name = isd.at("name_platform");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the platform name.",
+ "Utilities::getPlatformName()"));
+ }
+ }
+ return name;
+}
+
+int getTotalLines(json isd, csm::WarningList *list) {
+ int lines = 0;
+ try {
+ lines = isd.at("image_lines");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the number of lines in the image.",
+ "Utilities::getTotalLines()"));
+ }
+ }
+ return lines;
+}
+
+int getTotalSamples(json isd, csm::WarningList *list) {
+ int samples = 0;
+ try {
+ samples = isd.at("image_samples");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the number of samples in the image.",
+ "Utilities::getTotalSamples()"));
+ }
+ }
+ return samples;
+}
+
+double getStartingTime(json isd, csm::WarningList *list) {
+ double time = 0.0;
+ try {
+ time = isd.at("starting_ephemeris_time");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the image start time.",
+ "Utilities::getStartingTime()"));
+ }
+ }
+ return time;
+}
+
+double getCenterTime(json isd, csm::WarningList *list) {
+ double time = 0.0;
+ try {
+ time = isd.at("center_ephemeris_time");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the center image time.",
+ "Utilities::getCenterTime()"));
+ }
+ }
+ return time;
+}
+
+std::vector<double> getIntegrationStartLines(json isd, csm::WarningList *list) {
+ std::vector<double> lines;
+ try {
+ for (auto& scanRate : isd.at("line_scan_rate")) {
+ lines.push_back(scanRate[0]);
+ }
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the integration start lines in the integration time table.",
+ "Utilities::getIntegrationStartLines()"));
+ }
+ }
+ return lines;
+}
+
+std::vector<double> getIntegrationStartTimes(json isd, csm::WarningList *list) {
+ std::vector<double> times;
+ try {
+ for (auto& scanRate : isd.at("line_scan_rate")) {
+ times.push_back(scanRate[1]);
+ }
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the integration start times in the integration time table.",
+ "Utilities::getIntegrationStartTimes()"));
+ }
+ }
+ return times;
+}
+
+std::vector<double> getIntegrationTimes(json isd, csm::WarningList *list) {
+ std::vector<double> times;
+ try {
+ for (auto& scanRate : isd.at("line_scan_rate")) {
+ times.push_back(scanRate[2]);
+ }
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the integration times in the integration time table.",
+ "Utilities::getIntegrationTimes()"));
+ }
+ }
+ return times;
+}
+
+int getSampleSumming(json isd, csm::WarningList *list) {
+ int summing = 0;
+ try {
+ summing = isd.at("detector_sample_summing");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the sample direction detector pixel summing.",
+ "Utilities::getSampleSumming()"));
+ }
+ }
+ return summing;
+}
+
+int getLineSumming(json isd, csm::WarningList *list) {
+ int summing = 0;
+ try {
+ summing = isd.at("detector_line_summing");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the line direction detector pixel summing.",
+ "Utilities::getLineSumming()"));
+ }
+ }
+ return summing;
+}
+
+double getFocalLength(json isd, csm::WarningList *list) {
+ double length = 0.0;
+ try {
+ length = isd.at("focal_length_model").at("focal_length");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the focal length.",
+ "Utilities::getFocalLength()"));
+ }
+ }
+ return length;
+}
+
+double getFocalLengthEpsilon(json isd, csm::WarningList *list) {
+ double epsilon = 0.0;
+ try {
+ epsilon = isd.at("focal_length_model").at("focal_epsilon");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the focal length uncertainty.",
+ "Utilities::getFocalLengthEpsilon()"));
+ }
+ }
+ return epsilon;
+}
+
+std::vector<double> getFocal2PixelLines(json isd, csm::WarningList *list) {
+ std::vector<double> transformation;
+ try {
+ transformation = isd.at("focal2pixel_lines").get<std::vector<double>>();
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the focal plane coordinate to detector lines transformation.",
+ "Utilities::getFocal2PixelLines()"));
+ }
+ }
+ return transformation;
+}
+
+std::vector<double> getFocal2PixelSamples(json isd, csm::WarningList *list) {
+ std::vector<double> transformation;
+ try {
+ transformation = isd.at("focal2pixel_samples").get<std::vector<double>>();
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the focal plane coordinate to detector samples transformation.",
+ "Utilities::getFocal2PixelSamples()"));
+ }
+ }
+ return transformation;
+}
+
+double getDetectorCenterLine(json isd, csm::WarningList *list) {
+ double line;
+ try {
+ line = isd.at("detector_center").at("line");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the detector center line.",
+ "Utilities::getDetectorCenterLine()"));
+ }
+ }
+ return line;
+}
+
+double getDetectorCenterSample(json isd, csm::WarningList *list) {
+ double sample;
+ try {
+ sample = isd.at("detector_center").at("sample");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the detector center sample.",
+ "Utilities::getDetectorCenterSample()"));
+ }
+ }
+ return sample;
+}
+
+
+double getDetectorStartingLine(json isd, csm::WarningList *list) {
+ double line;
+ try {
+ line = isd.at("starting_detector_line");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the detector starting line.",
+ "Utilities::getDetectorStartingLine()"));
+ }
+ }
+ return line;
+}
+
+double getDetectorStartingSample(json isd, csm::WarningList *list) {
+ double sample;
+ try {
+ sample = isd.at("starting_detector_sample");
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the detector starting sample.",
+ "Utilities::getDetectorStartingSample()"));
+ }
+ }
+ return sample;
+}
+
+double getMinHeight(json isd, csm::WarningList *list) {
+ double height = 0.0;
+ try {
+ json referenceHeight = isd.at("reference_height");
+ json minHeight = referenceHeight.at("minheight");
+ json unit = referenceHeight.at("unit");
+ height = metric_conversion(minHeight.get<double>(), unit.get<std::string>());
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the minimum height above the reference ellipsoid.",
+ "Utilities::getMinHeight()"));
+ }
+ }
+ return height;
+}
+
+double getMaxHeight(json isd, csm::WarningList *list) {
+ double height = 0.0;
+ try {
+ json referenceHeight = isd.at("reference_height");
+ json maxHeight = referenceHeight.at("maxheight");
+ json unit = referenceHeight.at("unit");
+ height = metric_conversion(maxHeight.get<double>(), unit.get<std::string>());
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the maximum height above the reference ellipsoid.",
+ "Utilities::getMaxHeight()"));
+ }
+ }
+ return height;
+}
+
+double getSemiMajorRadius(json isd, csm::WarningList *list) {
+ double radius = 0.0;
+ try {
+ json radii = isd.at("radii");
+ json semiMajor = radii.at("semimajor");
+ json unit = radii.at("unit");
+ radius = metric_conversion(semiMajor.get<double>(), unit.get<std::string>());
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the reference ellipsoid semimajor radius.",
+ "Utilities::getSemiMajorRadius()"));
+ }
+ }
+ return radius;
+}
+
+double getSemiMinorRadius(json isd, csm::WarningList *list) {
+ double radius = 0.0;
+ try {
+ json radii = isd.at("radii");
+ json semiMinor = radii.at("semiminor");
+ json unit = radii.at("unit");
+ radius = metric_conversion(semiMinor.get<double>(), unit.get<std::string>());
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the reference ellipsoid semiminor radius.",
+ "Utilities::getSemiMinorRadius()"));
+ }
+ }
+ return radius;
+}
+
+DistortionType getDistortionModel(json isd, csm::WarningList *list) {
+ json distoriton_subset = isd.at("optical_distortion");
+
+ json::iterator it = distoriton_subset.begin();
+
+ std::string distortion = (std::string)it.key();
+
+ if (distortion.compare("transverse") == 0) {
+ return DistortionType::TRANSVERSE;
+ }
+ else if (distortion.compare("radial") == 0) {
+ return DistortionType::RADIAL;
+ }
+
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the distortion model.",
+ "Utilities::getDistortionModel()"));
+ }
+ return DistortionType::TRANSVERSE;
+}
+
+std::vector<double> getDistortionCoeffs(json isd, csm::WarningList *list) {
+ std::vector<double> coefficients;
+
+ DistortionType distortion = getDistortionModel(isd);
+
+ switch (distortion) {
+ case DistortionType::TRANSVERSE: {
+ try {
+ std::vector<double> coefficientsX, coefficientsY;
+
+ coefficientsX = isd.at("optical_distortion").at("transverse").at("x").get<std::vector<double>>();
+ coefficientsX.resize(10, 0.0);
+
+ coefficientsY = isd.at("optical_distortion").at("transverse").at("y").get<std::vector<double>>();
+ coefficientsY.resize(10, 0.0);
+
+ coefficients = coefficientsX;
+
+ coefficients.insert(coefficients.end(), coefficientsY.begin(), coefficientsY.end());
+ return coefficients;
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse a set of transverse distortion model coefficients.",
+ "Utilities::getDistortion()"));
+ }
+ coefficients = std::vector<double>(20, 0.0);
+ }
+ }
+ break;
+ case DistortionType::RADIAL: {
+ try {
+ coefficients = isd.at("optical_distortion").at("radial").at("coefficients").get<std::vector<double>>();
+
+ return coefficients;
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the radial distortion model coefficients.",
+ "Utilities::getDistortion()"));
+ }
+ coefficients = std::vector<double>(3, 0.0);
+ }
+ }
+ }
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the distortion model coefficients.",
+ "Utilities::getDistortion()"));
+ }
+
+ return coefficients;
+}
+
+std::vector<double> getSunPositions(json isd, csm::WarningList *list) {
+ std::vector<double> positions;
+ try {
+ json jayson = isd.at("sun_position");
+ json unit = jayson.at("unit");
+ for (auto& location : jayson.at("positions")) {
+ positions.push_back(metric_conversion(location[0].get<double>(), unit.get<std::string>()));
+ positions.push_back(metric_conversion(location[1].get<double>(), unit.get<std::string>()));
+ positions.push_back(metric_conversion(location[2].get<double>(), unit.get<std::string>()));
+ }
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the sun positions.",
+ "Utilities::getSunPositions()"));
+ }
+ }
+ return positions;
+}
+
+std::vector<double> getSensorPositions(json isd, csm::WarningList *list) {
+ std::vector<double> positions;
+ try {
+ json jayson = isd.at("sensor_position");
+ json unit = jayson.at("unit");
+ for (auto& location : jayson.at("positions")) {
+ positions.push_back(metric_conversion(location[0].get<double>(), unit.get<std::string>()));
+ positions.push_back(metric_conversion(location[1].get<double>(), unit.get<std::string>()));
+ positions.push_back(metric_conversion(location[2].get<double>(), unit.get<std::string>()));
+ }
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the sensor positions.",
+ "Utilities::getSensorPositions()"));
+ }
+ }
+ return positions;
+}
+
+std::vector<double> getSensorVelocities(json isd, csm::WarningList *list) {
+ std::vector<double> velocities;
+ try {
+ json jayson = isd.at("sensor_position");
+ json unit = jayson.at("unit");
+ for (auto& velocity : jayson.at("velocities")) {
+ velocities.push_back(metric_conversion(velocity[0].get<double>(), unit.get<std::string>()));
+ velocities.push_back(metric_conversion(velocity[1].get<double>(), unit.get<std::string>()));
+ velocities.push_back(metric_conversion(velocity[2].get<double>(), unit.get<std::string>()));
+ }
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the sensor velocities.",
+ "Utilities::getSensorVelocities()"));
+ }
+ }
+ return velocities;
+}
+
+std::vector<double> getSensorOrientations(json isd, csm::WarningList *list) {
+ std::vector<double> quaternions;
+ try {
+ for (auto& quaternion : isd.at("sensor_orientation").at("quaternions")) {
+ quaternions.push_back(quaternion[0]);
+ quaternions.push_back(quaternion[1]);
+ quaternions.push_back(quaternion[2]);
+ quaternions.push_back(quaternion[3]);
+ }
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the sensor orientations.",
+ "Utilities::getSensorOrientations()"));
+ }
+ }
+ return quaternions;
+}
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 3388d81..3948c05 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -3,7 +3,11 @@ cmake_minimum_required(VERSION 3.10)
# Link runCSMCameraModelTests with what we want to test and the GTest and pthread library
add_executable(runCSMCameraModelTests
PluginTests.cpp
- FrameCameraTests.cpp)
+ FrameCameraTests.cpp
+ LineScanCameraTests.cpp
+ DistortionTests.cpp
+ ISDParsingTests.cpp
+ UtilitiesTests.cpp)
if(WIN32)
option(CMAKE_USE_WIN32_THREADS_INIT "using WIN32 threads" ON)
option(gtest_disable_pthreads "Disable uses of pthreads in gtest." ON)
diff --git a/tests/DistortionTests.cpp b/tests/DistortionTests.cpp
new file mode 100644
index 0000000..26715ef
--- /dev/null
+++ b/tests/DistortionTests.cpp
@@ -0,0 +1,151 @@
+#include "Distortion.h"
+
+#include <gtest/gtest.h>
+
+#include "Fixtures.h"
+
+// NOTE: The imagePt format is (Lines,Samples)
+
+INSTANTIATE_TEST_CASE_P(JacobianTest,ImageCoordParameterizedTest,
+ ::testing::Values(csm::ImageCoord(2.5,2.5),csm::ImageCoord(7.5,7.5)));
+
+TEST_P(ImageCoordParameterizedTest, JacobianTest) {
+ std::vector<double> transverseDistortionCoeffs = {1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0,
+ 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0};
+
+ double Jxx,Jxy,Jyx,Jyy;
+
+ csm::ImageCoord imagePt = GetParam();
+ double jacobian[4];
+ distortionJacobian(imagePt.samp, imagePt.line, jacobian, transverseDistortionCoeffs);
+
+ // Jxx * Jyy - Jxy * Jyx
+ double determinant = fabs(jacobian[0] * jacobian[3] - jacobian[1] * jacobian[2]);
+ EXPECT_GT(determinant,1e-3);
+}
+
+TEST(Transverse, Jacobian1) {
+ csm::ImageCoord imagePt(7.5, 7.5);
+
+ std::vector<double> transverseDistortionCoeffs = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0};
+
+ double jacobian[4];
+ distortionJacobian(imagePt.samp, imagePt.line, jacobian, transverseDistortionCoeffs);
+
+ EXPECT_NEAR(jacobian[0],56.25,1e-8 );
+ EXPECT_NEAR(jacobian[1],112.5,1e-8);
+ EXPECT_NEAR(jacobian[2],56.25,1e-8);
+ EXPECT_NEAR(jacobian[3],281.25,1e-8);
+}
+
+TEST(Transverse, distortMe_AlternatingOnes) {
+ csm::ImageCoord imagePt(7.5, 7.5);
+
+ std::vector<double> transverseDistortionCoeffs = {1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,
+ 0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0};
+
+ double dx, dy;
+ computeTransverseDistortion(imagePt.samp, imagePt.line, dx, dy, transverseDistortionCoeffs);
+
+ EXPECT_NEAR(dx,908.5,1e-8 );
+ EXPECT_NEAR(dy,963.75,1e-8);
+}
+
+TEST(Transverse, distortMe_AllCoefficientsOne) {
+ csm::ImageCoord imagePt(7.5, 7.5);
+
+ std::vector<double> transverseDistortionCoeffs = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
+ 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
+
+ double dx, dy;
+ computeTransverseDistortion(imagePt.samp, imagePt.line, dx, dy, transverseDistortionCoeffs);
+
+ EXPECT_NEAR(dx,1872.25,1e-8 );
+ EXPECT_NEAR(dy,1872.25,1e-8);
+}
+
+TEST(transverse, removeDistortion1) {
+ csm::ImageCoord imagePt(7.5, 7.5);
+ double ux, uy;
+
+ std::vector<double> transverseDistortionCoeffs = {0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
+ 0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
+
+ removeDistortion(imagePt.samp, imagePt.line, ux, uy, transverseDistortionCoeffs, DistortionType::TRANSVERSE);
+
+ EXPECT_NEAR(imagePt.samp, 7.5, 1e-8 );
+ EXPECT_NEAR(imagePt.line, 7.5, 1e-8);
+ EXPECT_NEAR(imagePt.line, ux, 1e-8);
+ EXPECT_NEAR(imagePt.samp, uy, 1e-8);
+}
+
+TEST(transverse, removeDistortion_AllCoefficientsOne) {
+ csm::ImageCoord imagePt(1872.25, 1872.25);
+ double ux, uy;
+
+ std::vector<double> transverseDistortionCoeffs = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
+ 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
+
+ removeDistortion(imagePt.samp, imagePt.line, ux, uy, transverseDistortionCoeffs, DistortionType::TRANSVERSE);
+
+ // The Jacobian is singular, so the setFocalPlane should break out of it's iteration and
+ // returns the same distorted coordinates that were passed in.
+ EXPECT_NEAR(ux, imagePt.samp,1e-8 );
+ EXPECT_NEAR(uy, imagePt.line,1e-8);
+}
+
+TEST(transverse, removeDistortion_AlternatingOnes) {
+ csm::ImageCoord imagePt(963.75, 908.5);
+ double ux, uy;
+
+ std::vector<double> transverseDistortionCoeffs = {1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,
+ 0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0};
+
+ removeDistortion(imagePt.samp, imagePt.line, ux, uy, transverseDistortionCoeffs, DistortionType::TRANSVERSE);
+
+ EXPECT_NEAR(ux, 7.5, 1e-8 );
+ EXPECT_NEAR(uy, 7.5, 1e-8);
+}
+
+TEST(Radial, testRemoveDistortion) {
+ csm::ImageCoord imagePt(0.0, 4.0);
+
+ double ux, uy;
+ std::vector<double> coeffs = {0, 0, 0};
+
+ removeDistortion(imagePt.samp, imagePt.line, ux, uy, coeffs, DistortionType::RADIAL);
+
+ EXPECT_NEAR(ux, 4, 1e-8);
+ EXPECT_NEAR(uy, 0, 1e-8);
+}
+
+// If coeffs are 0 then this will have the same result as removeDistortion
+// with 0 distortion coefficients
+TEST(Radial, testInverseDistortion){
+ csm::ImageCoord imagePt(0.0, 4.0);
+
+ double dx, dy;
+ double desiredPrecision = 0.01;
+ std::vector<double> coeffs = {0, 0, 0};
+
+ applyDistortion(imagePt.samp, imagePt.line, dx, dy, coeffs,
+ DistortionType::RADIAL, desiredPrecision);
+
+ EXPECT_NEAR(dx,4,1e-8);
+ EXPECT_NEAR(dy,0,1e-8);
+}
+
+TEST(Radial, testInverseOnesCoeffs){
+ csm::ImageCoord imagePt(0.0, 4.0);
+
+ double dx, dy;
+ double desiredPrecision = 0.01;
+ std::vector<double> coeffs = {1, 1, 1};
+
+ applyDistortion(imagePt.samp, imagePt.line, dx, dy, coeffs,
+ DistortionType::RADIAL, desiredPrecision);
+
+ EXPECT_NEAR(dx,4,1e-8);
+ EXPECT_NEAR(dy,0,1e-8);
+}
diff --git a/tests/Fixtures.h b/tests/Fixtures.h
index be28df3..e372b12 100644
--- a/tests/Fixtures.h
+++ b/tests/Fixtures.h
@@ -3,6 +3,7 @@
#include "UsgsAstroPlugin.h"
#include "UsgsAstroFrameSensorModel.h"
+#include "UsgsAstroLsSensorModel.h"
#include <json.hpp>
@@ -32,6 +33,8 @@ inline void jsonToIsd(json &object, csm::Isd &isd, std::string prefix="") {
}
}
+//////////Framing Camera Sensor Model Fixtures//////////
+
class FrameSensorModel : public ::testing::Test {
protected:
csm::Isd isd;
@@ -77,6 +80,8 @@ class ConstVelLineScanIsdTest : public ::testing::Test {
}
};
+class ImageCoordParameterizedTest : public ::testing::TestWithParam<csm::ImageCoord> {};
+
class FramerParameterizedTest : public ::testing::TestWithParam<csm::ImageCoord> {
protected:
@@ -139,6 +144,62 @@ class FrameStateTest : public ::testing::Test {
}
};
+//////////Line Scan Camera Sensor Model Fixtures//////////
+
+class ConstVelocityLineScanSensorModel : public ::testing::Test {
+ protected:
+ csm::Isd isd;
+ UsgsAstroLsSensorModel *sensorModel;
+
+ void SetUp() override {
+ sensorModel = NULL;
+
+ isd.setFilename("data/constVelocityLineScan.img");
+ UsgsAstroPlugin cameraPlugin;
+
+ csm::Model *model = cameraPlugin.constructModelFromISD(
+ isd,
+ "USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL");
+ sensorModel = dynamic_cast<UsgsAstroLsSensorModel *>(model);
+
+ ASSERT_NE(sensorModel, nullptr);
+ }
+
+ void TearDown() override {
+ if (sensorModel) {
+ delete sensorModel;
+ sensorModel = NULL;
+ }
+ }
+};
+
+class ConstAngularVelocityLineScanSensorModel : public ::testing::Test {
+ protected:
+ csm::Isd isd;
+ UsgsAstroLsSensorModel *sensorModel;
+
+ void SetUp() override {
+ sensorModel = NULL;
+
+ isd.setFilename("data/constAngularVelocityLineScan.img");
+ UsgsAstroPlugin cameraPlugin;
+
+ csm::Model *model = cameraPlugin.constructModelFromISD(
+ isd,
+ "USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL");
+ sensorModel = dynamic_cast<UsgsAstroLsSensorModel *>(model);
+
+ ASSERT_NE(sensorModel, nullptr);
+ }
+
+ void TearDown() override {
+ if (sensorModel) {
+ delete sensorModel;
+ sensorModel = NULL;
+ }
+ }
+};
+
#endif
diff --git a/tests/FrameCameraTests.cpp b/tests/FrameCameraTests.cpp
index cab179e..d41b461 100644
--- a/tests/FrameCameraTests.cpp
+++ b/tests/FrameCameraTests.cpp
@@ -7,34 +7,6 @@
#include "Fixtures.h"
using json = nlohmann::json;
-
-INSTANTIATE_TEST_CASE_P(JacobianTest,FramerParameterizedTest,
- ::testing::Values(csm::ImageCoord(2.5,2.5),csm::ImageCoord(7.5,7.5)));
-
-TEST_P(FramerParameterizedTest, JacobianTest) {
-
- UsgsAstroFrameSensorModel* sensorModel = createModel(isd);
- std::string modelState = sensorModel->getModelState();
- auto state = json::parse(modelState);
-
- state["m_odtX"] = {1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0};
- state["m_odtY"] = {0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0};
- sensorModel->replaceModelState(state.dump());
-
- double Jxx,Jxy,Jyx,Jyy;
- ASSERT_NE(sensorModel, nullptr);
-
- csm::ImageCoord imagePt1 = GetParam();
-// cout << "[" << imagePt1.samp << "," << imagePt1.line << "]"<< endl;
- sensorModel->distortionJacobian(imagePt1.samp, imagePt1.line, Jxx, Jxy,Jyx,Jyy);
-
- double determinant = fabs(Jxx*Jyy - Jxy*Jyx);
- EXPECT_GT(determinant,1e-3);
-
- delete sensorModel;
- sensorModel=NULL;
-}
-
// NOTE: The imagePt format is (Lines,Samples)
TEST_F(FrameSensorModel, State) {
@@ -86,6 +58,13 @@ TEST_F(FrameSensorModel, OffBody3) {
EXPECT_NEAR(groundPt.z, 14.99325304, 1e-8);
}
+TEST_F(FrameSensorModel, getReferencePoint) {
+ csm::EcefCoord groundPt = sensorModel->getReferencePoint();
+ EXPECT_EQ(groundPt.x, 0.0);
+ EXPECT_EQ(groundPt.y, 0.0);
+ EXPECT_EQ(groundPt.z, 0.0);
+}
+
TEST_F(FrameSensorModel, OffBody4) {
csm::ImageCoord imagePt(15.0, 15.0);
csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0);
@@ -94,133 +73,8 @@ TEST_F(FrameSensorModel, OffBody4) {
EXPECT_NEAR(groundPt.z, -14.99325304, 1e-8);
}
-TEST_F(FrameSensorModel, setFocalPlane1) {
- csm::ImageCoord imagePt(7.5, 7.5);
- double ux,uy;
-
- std::string modelState = sensorModel->getModelState();
- auto state = json::parse(modelState);
-
- state["m_odtX"] = {0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
- state["m_odtY"] = {0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
- sensorModel->replaceModelState(state.dump());
-
- ASSERT_NE(sensorModel, nullptr);
- sensorModel->setFocalPlane(imagePt.samp, imagePt.line, ux, uy);
- EXPECT_NEAR(imagePt.samp,7.5,1e-8 );
- EXPECT_NEAR(imagePt.line,7.5,1e-8);
- delete sensorModel;
- sensorModel = NULL;
-}
-
-
-TEST_F(FrameSensorModel, Jacobian1) {
- csm::ImageCoord imagePt(7.5, 7.5);
-
- std::string modelState = sensorModel->getModelState();
- auto state = json::parse(modelState);
- state["m_odtX"] = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0};
- state["m_odtY"] = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0};
- sensorModel->replaceModelState(state.dump());
-
- double Jxx,Jxy,Jyx,Jyy;
-
- ASSERT_NE(sensorModel, nullptr);
- sensorModel->distortionJacobian(imagePt.samp, imagePt.line, Jxx, Jxy,Jyx,Jyy);
-
- EXPECT_NEAR(Jxx,56.25,1e-8 );
- EXPECT_NEAR(Jxy,112.5,1e-8);
- EXPECT_NEAR(Jyx,56.25,1e-8);
- EXPECT_NEAR(Jyy,281.25,1e-8);
-
- delete sensorModel;
- sensorModel = NULL;
-}
-
-
-TEST_F(FrameSensorModel, distortMe_AllCoefficientsOne) {
- csm::ImageCoord imagePt(7.5, 7.5);
-
- std::string modelState = sensorModel->getModelState();
- auto state = json::parse(modelState);
- state["m_odtX"] = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
- state["m_odtY"] = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
- sensorModel->replaceModelState(state.dump());
-
- double dx,dy;
- ASSERT_NE(sensorModel, nullptr);
- sensorModel->distortionFunction(imagePt.samp, imagePt.line,dx,dy );
-
- EXPECT_NEAR(dx,1872.25,1e-8 );
- EXPECT_NEAR(dy,1872.25,1e-8);
-
- delete sensorModel;
- sensorModel = NULL;
-}
-
-TEST_F(FrameSensorModel, setFocalPlane_AllCoefficientsOne) {
- csm::ImageCoord imagePt(1872.25, 1872.25);
-
- std::string modelState = sensorModel->getModelState();
- auto state = json::parse(modelState);
- state["m_odtX"] = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
- state["m_odtY"] = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
- sensorModel->replaceModelState(state.dump());
-
- double ux,uy;
- ASSERT_NE(sensorModel, nullptr);
- sensorModel->setFocalPlane(imagePt.samp, imagePt.line,ux,uy );
-
- // The Jacobian is singular, so the setFocalPlane should break out of it's iteration and
- // returns the same distorted coordinates that were passed in.
- EXPECT_NEAR(ux,imagePt.samp,1e-8 );
- EXPECT_NEAR(uy,imagePt.line,1e-8);
-
- delete sensorModel;
- sensorModel = NULL;
-}
-
-
-TEST_F(FrameSensorModel, distortMe_AlternatingOnes) {
- csm::ImageCoord imagePt(7.5, 7.5);
-
- std::string modelState = sensorModel->getModelState();
- auto state = json::parse(modelState);
- state["m_odtX"] = {1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0};
- state["m_odtY"] = {0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0};
- sensorModel->replaceModelState(state.dump());
-
- double dx,dy;
- ASSERT_NE(sensorModel, nullptr);
- sensorModel->distortionFunction(imagePt.samp, imagePt.line,dx,dy );
-
- EXPECT_NEAR(dx,908.5,1e-8 );
- EXPECT_NEAR(dy,963.75,1e-8);
-
- delete sensorModel;
- sensorModel = NULL;
-}
-
-
-TEST_F(FrameSensorModel, setFocalPlane_AlternatingOnes) {
- csm::ImageCoord imagePt(963.75, 908.5);
-
- std::string modelState = sensorModel->getModelState();
- auto state = json::parse(modelState);
- state["m_odtX"] = {1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0};
- state["m_odtY"] = {0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0};
- sensorModel->replaceModelState(state.dump());
-
- double ux,uy;
- ASSERT_NE(sensorModel, nullptr);
-
- sensorModel->setFocalPlane(imagePt.samp, imagePt.line,ux,uy );
-
- EXPECT_NEAR(ux,7.5,1e-8 );
- EXPECT_NEAR(uy,7.5,1e-8);
-
- delete sensorModel;
- sensorModel = NULL;
+TEST_F(FrameSensorModel, getImageIdentifier) {
+ EXPECT_EQ("simpleFramerISD", sensorModel->getImageIdentifier());
}
diff --git a/tests/ISDParsingTests.cpp b/tests/ISDParsingTests.cpp
new file mode 100644
index 0000000..3c2c062
--- /dev/null
+++ b/tests/ISDParsingTests.cpp
@@ -0,0 +1,313 @@
+#include "Utilities.h"
+
+#include <gtest/gtest.h>
+
+#include "Fixtures.h"
+
+TEST(MetricConversion, DistanceConversion) {
+ EXPECT_EQ(1, metric_conversion(1000, "m", "km"));
+ EXPECT_EQ(1000, metric_conversion(1000, "m", "m"));
+ EXPECT_EQ(1000, metric_conversion(1, "km", "m"));
+ EXPECT_EQ(1, metric_conversion(1, "km", "km"));
+ EXPECT_EQ(0, metric_conversion(0, "m", "km"));
+ EXPECT_EQ(0, metric_conversion(0, "m", "m"));
+ EXPECT_EQ(0, metric_conversion(0, "km", "m"));
+ EXPECT_EQ(0, metric_conversion(0, "km", "km"));
+ EXPECT_EQ(-1, metric_conversion(-1000, "m", "km"));
+ EXPECT_EQ(-1000, metric_conversion(-1000, "m", "m"));
+ EXPECT_EQ(-1000, metric_conversion(-1, "km", "m"));
+ EXPECT_EQ(-1, metric_conversion(-1, "km", "km"));
+}
+
+TEST(ISDParsing, ModelName) {
+ json isd = {
+ {"name_model", "Test"}
+ };
+ EXPECT_EQ("Test", getSensorModelName(isd));
+}
+
+TEST(ISDParsing, ImageIdentifier) {
+ json isd = {
+ {"image_identifier", "Test"}
+ };
+ EXPECT_EQ("Test", getImageId(isd));
+}
+
+TEST(ISDParsing, SensorName) {
+ json isd = {
+ {"name_sensor", "Test"}
+ };
+ EXPECT_EQ("Test", getSensorName(isd));
+}
+
+TEST(ISDParsing, PlatformName) {
+ json isd = {
+ {"name_platform", "Test"}
+ };
+ EXPECT_EQ("Test", getPlatformName(isd));
+}
+
+TEST(ISDParsing, TotalLines) {
+ json isd = {
+ {"image_lines", 16}
+ };
+ EXPECT_EQ(16, getTotalLines(isd));
+}
+
+TEST(ISDParsing, TotalSamples) {
+ json isd = {
+ {"image_samples", 16}
+ };
+ EXPECT_EQ(16, getTotalSamples(isd));
+}
+
+TEST(ISDParsing, StartTime) {
+ json isd = {
+ {"starting_ephemeris_time", -10}
+ };
+ EXPECT_EQ(-10, getStartingTime(isd));
+}
+
+TEST(ISDParsing, CenterTime) {
+ json isd = {
+ {"center_ephemeris_time", -10}
+ };
+ EXPECT_EQ(-10, getCenterTime(isd));
+}
+
+TEST(ISDParsing, IntegrationStartLines) {
+ json isd = {
+ {"line_scan_rate", {
+ {0, 1, 2},
+ {3, 4, 5},
+ {6, 7, 8}}
+ }
+ };
+ std::vector<double> startLines = {0, 3, 6};
+ EXPECT_EQ(startLines, getIntegrationStartLines(isd));
+}
+
+TEST(ISDParsing, IntegrationStartTimes) {
+ json isd = {
+ {"line_scan_rate", {
+ {0, 1, 2},
+ {3, 4, 5},
+ {6, 7, 8}}
+ }
+ };
+ std::vector<double> startTimes = {1, 4, 7};
+ EXPECT_EQ(startTimes, getIntegrationStartTimes(isd));
+}
+
+TEST(ISDParsing, IntegrationTimes) {
+ json isd = {
+ {"line_scan_rate", {
+ {0, 1, 2},
+ {3, 4, 5},
+ {6, 7, 8}}
+ }
+ };
+ std::vector<double> times = {2, 5, 8};
+ EXPECT_EQ(times, getIntegrationTimes(isd));
+}
+
+TEST(ISDParsing, SampleSumming) {
+ json isd = {
+ {"detector_sample_summing", 4}
+ };
+ EXPECT_EQ(4, getSampleSumming(isd));
+}
+
+TEST(ISDParsing, LineSumming) {
+ json isd = {
+ {"detector_line_summing", 4}
+ };
+ EXPECT_EQ(4, getLineSumming(isd));
+}
+
+TEST(ISDParsing, FocalLength) {
+ json isd = {
+ {"focal_length_model", {
+ {"focal_length", 2}}
+ }
+ };
+ EXPECT_EQ(2, getFocalLength(isd));
+}
+
+TEST(ISDParsing, FocalLengthEpsilon) {
+ json isd = {
+ {"focal_length_model", {
+ {"focal_epsilon", 0.1}}
+ }
+ };
+ EXPECT_EQ(0.1, getFocalLengthEpsilon(isd));
+}
+
+TEST(ISDParsing, Focal2PixelLines) {
+ json isd = {
+ {"focal2pixel_lines", {0, 1, 2}}
+ };
+ std::vector<double> coefficients = {0, 1, 2};
+ EXPECT_EQ(coefficients, getFocal2PixelLines(isd));
+}
+
+TEST(ISDParsing, Focal2PixelSamples) {
+ json isd = {
+ {"focal2pixel_samples", {0, 1, 2}}
+ };
+ std::vector<double> coefficients = {0, 1, 2};
+ EXPECT_EQ(coefficients, getFocal2PixelSamples(isd));
+}
+
+TEST(ISDParsing, DetectorCenterLine) {
+ json isd = {
+ {"detector_center", {
+ {"line", 2}}
+ }
+ };
+ EXPECT_EQ(2, getDetectorCenterLine(isd));
+}
+
+TEST(ISDParsing, DetectorCenterSample) {
+ json isd = {
+ {"detector_center", {
+ {"sample", 3}}
+ }
+ };
+ EXPECT_EQ(3, getDetectorCenterSample(isd));
+}
+
+TEST(ISDParsing, DetectorStartingLine) {
+ json isd = {
+ {"starting_detector_line", 1}
+ };
+ EXPECT_EQ(1, getDetectorStartingLine(isd));
+}
+
+TEST(ISDParsing, DetectorStartingSample) {
+ json isd = {
+ {"starting_detector_sample", 2}
+ };
+ EXPECT_EQ(2, getDetectorStartingSample(isd));
+}
+
+TEST(ISDParsing, MinHeight) {
+ json isd = {
+ {"reference_height", {
+ {"minheight", -1},
+ {"unit", "km"}}
+ }
+ };
+ EXPECT_EQ(-1000, getMinHeight(isd));
+}
+
+TEST(ISDParsing, MaxHeight) {
+ json isd = {
+ {"reference_height", {
+ {"maxheight", 10},
+ {"unit", "km"}}
+ }
+ };
+ EXPECT_EQ(10000, getMaxHeight(isd));
+}
+
+TEST(ISDParsing, SemiMajor) {
+ json isd = {
+ {"radii", {
+ {"semimajor", 2},
+ {"unit", "km"}}
+ }
+ };
+ EXPECT_EQ(2000, getSemiMajorRadius(isd));
+}
+
+TEST(ISDParsing, SemiMinor) {
+ json isd = {
+ {"radii", {
+ {"semiminor", 1},
+ {"unit", "km"}}
+ }
+ };
+ EXPECT_EQ(1000, getSemiMinorRadius(isd));
+}
+
+TEST(ISDParsing, TransverseDistortion) {
+ json isd = {
+ {"optical_distortion", {
+ {"transverse", {
+ {"y", {-11, 21, 24}},
+ {"x", {-1, 2, 4}}}}
+ }
+ }
+ };
+ std::vector<double> coefficients = {-1, 2, 4, 0, 0, 0, 0, 0, 0, 0,
+ -11, 21, 24, 0, 0, 0, 0, 0, 0, 0};
+ EXPECT_EQ(coefficients, getDistortionCoeffs(isd));
+}
+
+TEST(ISDParsing, Radial) {
+ json isd = {
+ {"optical_distortion", {
+ {"radial", {
+ {"coefficients", {0, 1, 2}}}}
+ }
+ }
+ };
+ std::vector<double> coefficients = {0, 1, 2};
+ EXPECT_EQ(coefficients, getDistortionCoeffs(isd));
+}
+
+TEST(ISDParsing, SunPosition) {
+ json isd = {
+ {"sun_position", {
+ {"positions", {
+ {0, 1, 2},
+ {3, 4, 5},
+ {6, 7, 8}}},
+ {"unit", "km"}}
+ }
+ };
+ std::vector<double> positions = {0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000};
+ EXPECT_EQ(positions, getSunPositions(isd));
+}
+
+TEST(ISDParsing, SensorPosition) {
+ json isd = {
+ {"sensor_position", {
+ {"positions", {
+ {0, 1, 2},
+ {3, 4, 5},
+ {6, 7, 8}}},
+ {"unit", "km"}}
+ }
+ };
+ std::vector<double> positions = {0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000};
+ EXPECT_EQ(positions, getSensorPositions(isd));
+}
+
+TEST(ISDParsing, SensorVelocities) {
+ json isd = {
+ {"sensor_position", {
+ {"velocities", {
+ {0, 1, 2},
+ {3, 4, 5},
+ {6, 7, 8}}},
+ {"unit", "km"}}
+ }
+ };
+ std::vector<double> velocity = {0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000};
+ EXPECT_EQ(velocity, getSensorVelocities(isd));
+}
+
+TEST(ISDParsing, SensorOrientations) {
+ json isd = {
+ {"sensor_orientation", {
+ {"quaternions", {
+ {0, 1, 2, 3},
+ {4, 5, 6, 7},
+ {8, 9, 10, 11}}}}
+ }
+ };
+ std::vector<double> rotations = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11};
+ EXPECT_EQ(rotations, getSensorOrientations(isd));
+}
diff --git a/tests/LineScanCameraTests.cpp b/tests/LineScanCameraTests.cpp
new file mode 100644
index 0000000..627b570
--- /dev/null
+++ b/tests/LineScanCameraTests.cpp
@@ -0,0 +1,127 @@
+#define _USE_MATH_DEFINES
+
+#include "Fixtures.h"
+#include "UsgsAstroPlugin.h"
+#include "UsgsAstroLsSensorModel.h"
+
+#include <json.hpp>
+#include <gtest/gtest.h>
+
+#include <math.h>
+#include <iostream>
+
+using json = nlohmann::json;
+
+// TODO all commented out expects are failing and need to either have updated numbers
+// for the line scanner test cases, or we need to figure out why the line scanner
+// is not honoring this functionality.
+
+
+TEST_F(ConstVelocityLineScanSensorModel, State) {
+ std::string modelState = sensorModel->getModelState();
+ sensorModel->replaceModelState(modelState);
+
+ // When this is different, the output is very hard to parse
+ // TODO implement JSON diff for gtest
+
+ EXPECT_EQ(sensorModel->getModelState(), modelState);
+}
+
+// Fly by tests
+
+TEST_F(ConstVelocityLineScanSensorModel, Center) {
+ csm::ImageCoord imagePt(8.5, 8.0);
+ csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0);
+ EXPECT_DOUBLE_EQ(groundPt.x, 10.0);
+ EXPECT_DOUBLE_EQ(groundPt.y, 0.0);
+ EXPECT_DOUBLE_EQ(groundPt.z, 0.0);
+}
+
+TEST_F(ConstVelocityLineScanSensorModel, Inversion) {
+ csm::ImageCoord imagePt(8.5, 8);
+ csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0);
+ csm::ImageCoord imageReprojPt = sensorModel->groundToImage(groundPt);
+
+ EXPECT_DOUBLE_EQ(imagePt.line, imageReprojPt.line);
+ EXPECT_DOUBLE_EQ(imagePt.samp, imageReprojPt.samp);
+}
+
+TEST_F(ConstVelocityLineScanSensorModel, OffBody) {
+ csm::ImageCoord imagePt(4.5, 4.0);
+ csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0);
+ EXPECT_NEAR(groundPt.x, 0.063995905, 1e-8);
+ EXPECT_NEAR(groundPt.y, -7.999488033, 1e-8);
+ EXPECT_NEAR(groundPt.z, 8, 1e-8);
+}
+
+TEST_F(ConstVelocityLineScanSensorModel, ProximateImageLocus) {
+ csm::ImageCoord imagePt(8.5, 8.0);
+ csm::EcefCoord groundPt(5, 5, 5);
+ csm::EcefLocus locus = sensorModel->imageToProximateImagingLocus(imagePt, groundPt);
+ EXPECT_DOUBLE_EQ(locus.direction.x, -1.0);
+ EXPECT_DOUBLE_EQ(locus.direction.y, 0.0);
+ EXPECT_DOUBLE_EQ(locus.direction.z, 0.0);
+ EXPECT_DOUBLE_EQ(locus.point.x, 5.0);
+ EXPECT_DOUBLE_EQ(locus.point.y, 0.0);
+ EXPECT_DOUBLE_EQ(locus.point.z, 0.0);
+}
+
+TEST_F(ConstVelocityLineScanSensorModel, RemoteImageLocus) {
+ csm::ImageCoord imagePt(8.5, 8.0);
+ csm::EcefLocus locus = sensorModel->imageToRemoteImagingLocus(imagePt);
+ EXPECT_DOUBLE_EQ(locus.direction.x, -1.0);
+ EXPECT_DOUBLE_EQ(locus.direction.y, 0.0);
+ EXPECT_DOUBLE_EQ(locus.direction.z, 0.0);
+ EXPECT_DOUBLE_EQ(locus.point.x, 10.0);
+ EXPECT_DOUBLE_EQ(locus.point.y, 0.0);
+ EXPECT_DOUBLE_EQ(locus.point.z, 0.0);
+}
+
+// Pan tests
+
+TEST_F(ConstAngularVelocityLineScanSensorModel, Center) {
+ csm::ImageCoord imagePt(8.5, 8.0);
+ csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0);
+ EXPECT_DOUBLE_EQ(groundPt.x, 10.0);
+ EXPECT_DOUBLE_EQ(groundPt.y, 0.0);
+ EXPECT_DOUBLE_EQ(groundPt.z, 0.0);
+}
+
+TEST_F(ConstAngularVelocityLineScanSensorModel, Inversion) {
+ csm::ImageCoord imagePt(8.5, 8);
+ csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0);
+ csm::ImageCoord imageReprojPt = sensorModel->groundToImage(groundPt);
+
+ EXPECT_DOUBLE_EQ(imagePt.line, imageReprojPt.line);
+ EXPECT_DOUBLE_EQ(imagePt.samp, imageReprojPt.samp);
+}
+
+TEST_F(ConstAngularVelocityLineScanSensorModel, OffBody) {
+ csm::ImageCoord imagePt(4.5, 4.0);
+ csm::EcefCoord groundPt = sensorModel->imageToGround(imagePt, 0.0);
+ EXPECT_NEAR(groundPt.x, 4.15414478, 1e-8);
+ EXPECT_NEAR(groundPt.y, -7.98311067, 1e-8);
+ EXPECT_NEAR(groundPt.z, 63.82129588, 1e-8);
+}
+
+TEST_F(ConstVelocityLineScanSensorModel, calculateAttitudeCorrection) {
+ std::vector<double> adj;
+ double attCorr[9];
+ adj.resize(15, 0);
+ // Pi/2 with simply compensating for member variable m_flyingHeight in UsgsAstroLsSensorModel
+ adj[7] = (M_PI / 2) * 990.0496255790623081338708;
+ sensorModel->calculateAttitudeCorrection(999.5, adj, attCorr);
+
+ // EXPECT_NEARs are used here instead of EXPECT_DOUBLE_EQs because index 0 and 8 of the matrix
+ // are evaluating to 6.12...e-17. This is too small to be worried about here, but
+ // EXPECT_DOUBLE_EQ is too sensitive.
+ EXPECT_NEAR(attCorr[0], 0, 1e-8);
+ EXPECT_NEAR(attCorr[1], 0, 1e-8);
+ EXPECT_NEAR(attCorr[2], 1, 1e-8);
+ EXPECT_NEAR(attCorr[3], 0, 1e-8);
+ EXPECT_NEAR(attCorr[4], 1, 1e-8);
+ EXPECT_NEAR(attCorr[5], 0, 1e-8);
+ EXPECT_NEAR(attCorr[6], -1, 1e-8);
+ EXPECT_NEAR(attCorr[7], 0, 1e-8);
+ EXPECT_NEAR(attCorr[8], 0, 1e-8);
+}
diff --git a/tests/PluginTests.cpp b/tests/PluginTests.cpp
index 2252e88..62da7c5 100644
--- a/tests/PluginTests.cpp
+++ b/tests/PluginTests.cpp
@@ -21,7 +21,7 @@ TEST(PluginTests, ManufacturerName) {
TEST(PluginTests, ReleaseDate) {
UsgsAstroPlugin testPlugin;
- EXPECT_EQ("20170425", testPlugin.getReleaseDate());
+ EXPECT_EQ("20190222", testPlugin.getReleaseDate());
}
TEST(PluginTests, NumModels) {
@@ -82,6 +82,14 @@ TEST_F(FrameIsdTest, Constructible) {
"USGS_ASTRO_FRAME_SENSOR_MODEL"));
}
+TEST_F(FrameIsdTest, ConstructibleFromState) {
+ UsgsAstroPlugin testPlugin;
+ std::string modelState = testPlugin.getStateFromISD(isd);
+ EXPECT_TRUE(testPlugin.canModelBeConstructedFromState(
+ "USGS_ASTRO_FRAME_SENSOR_MODEL",
+ modelState));
+}
+
TEST_F(FrameIsdTest, NotConstructible) {
UsgsAstroPlugin testPlugin;
isd.setFilename("data/constVelocityLineScan.img");
@@ -110,13 +118,13 @@ TEST_F(FrameIsdTest, ConstructValidCamera) {
TEST_F(FrameIsdTest, ConstructInValidCamera) {
UsgsAstroPlugin testPlugin;
- isd.setFilename("data/constVelocityLineScan.img");
+ isd.setFilename("data/empty.img");
csm::Model *cameraModel = NULL;
try {
testPlugin.constructModelFromISD(
isd,
"USGS_ASTRO_FRAME_SENSOR_MODEL",
- NULL);
+ nullptr);
FAIL() << "Expected an error";
}
catch(csm::Error &e) {
@@ -137,6 +145,14 @@ TEST_F(ConstVelLineScanIsdTest, Constructible) {
"USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL"));
}
+TEST_F(ConstVelLineScanIsdTest, ConstructibleFromState) {
+ UsgsAstroPlugin testPlugin;
+ std::string modelState = testPlugin.getStateFromISD(isd);
+ EXPECT_TRUE(testPlugin.canModelBeConstructedFromState(
+ "USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL",
+ modelState));
+}
+
TEST_F(ConstVelLineScanIsdTest, NotConstructible) {
UsgsAstroPlugin testPlugin;
isd.setFilename("data/simpleFramerISD.img");
@@ -163,13 +179,13 @@ TEST_F(ConstVelLineScanIsdTest, ConstructValidCamera) {
TEST_F(ConstVelLineScanIsdTest, ConstructInValidCamera) {
UsgsAstroPlugin testPlugin;
- isd.setFilename("data/simpleFramerISD.img");
+ isd.setFilename("data/empty.img");
csm::Model *cameraModel = NULL;
try {
testPlugin.constructModelFromISD(
isd,
"USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL",
- NULL);
+ nullptr);
FAIL() << "Expected an error";
}
diff --git a/tests/UtilitiesTests.cpp b/tests/UtilitiesTests.cpp
new file mode 100644
index 0000000..5473dc5
--- /dev/null
+++ b/tests/UtilitiesTests.cpp
@@ -0,0 +1,69 @@
+#define _USE_MATH_DEFINES
+
+#include "Fixtures.h"
+#include "UsgsAstroPlugin.h"
+#include "Utilities.h"
+
+#include <json.hpp>
+#include <gtest/gtest.h>
+
+#include <math.h>
+
+using json = nlohmann::json;
+
+
+TEST(UtilitiesTests, calculateRotationMatrixFromEuler) {
+ double euler[3], rotationMatrix[9];
+ euler[0] = 0;
+ euler[1] = M_PI/2;
+ euler[2] = 0;
+ calculateRotationMatrixFromEuler(euler, rotationMatrix);
+
+ // EXPECT_NEARs are used here instead of EXPECT_DOUBLE_EQs because index 0 and 8 of the matrix
+ // are evaluating to 6.12...e-17. This is too small to be worried about here, but
+ // EXPECT_DOUBLE_EQ is too sensitive.
+ EXPECT_NEAR(rotationMatrix[0], 0, 1e-8);
+ EXPECT_NEAR(rotationMatrix[1], 0, 1e-8);
+ EXPECT_NEAR(rotationMatrix[2], 1, 1e-8);
+ EXPECT_NEAR(rotationMatrix[3], 0, 1e-8);
+ EXPECT_NEAR(rotationMatrix[4], 1, 1e-8);
+ EXPECT_NEAR(rotationMatrix[5], 0, 1e-8);
+ EXPECT_NEAR(rotationMatrix[6], -1, 1e-8);
+ EXPECT_NEAR(rotationMatrix[7], 0, 1e-8);
+ EXPECT_NEAR(rotationMatrix[8], 0, 1e-8);
+}
+
+TEST(UtilitiesTests, calculateRotationMatrixFromQuaternions) {
+ double q[4], rotationMatrix[9];
+ q[0] = 0;
+ q[1] = -1/sqrt(2);
+ q[2] = 0;
+ q[3] = 1/sqrt(2);
+ calculateRotationMatrixFromQuaternions(q, rotationMatrix);
+ EXPECT_DOUBLE_EQ(rotationMatrix[0], 0);
+ EXPECT_DOUBLE_EQ(rotationMatrix[1], 0);
+ EXPECT_DOUBLE_EQ(rotationMatrix[2], -1);
+ EXPECT_DOUBLE_EQ(rotationMatrix[3], 0);
+ EXPECT_DOUBLE_EQ(rotationMatrix[4], 1);
+ EXPECT_DOUBLE_EQ(rotationMatrix[5], 0);
+ EXPECT_DOUBLE_EQ(rotationMatrix[6], 1);
+ EXPECT_DOUBLE_EQ(rotationMatrix[7], 0);
+ EXPECT_DOUBLE_EQ(rotationMatrix[8], 0);
+}
+
+TEST(UtilitiesTests, computeDistortedFocalPlaneCoordinates) {
+ double iTransS[] = {0.0, 0.0, 10.0};
+ double iTransL[] = {0.0, 10.0, 0.0};
+ std::tuple<double, double> natFocalPlane;
+ computeDistortedFocalPlaneCoordinates(0.5, 4, 8, 0.5, 1, 0, iTransS, iTransL, natFocalPlane);
+ EXPECT_DOUBLE_EQ(std::get<0> (natFocalPlane), 0);
+ EXPECT_DOUBLE_EQ(std::get<1> (natFocalPlane), -0.4);
+}
+
+TEST(UtilitiesTests, createCameraLookVector) {
+ double cameraLook[3];
+ createCameraLookVector(0, -0.4, 1, 50, cameraLook);
+ EXPECT_NEAR(cameraLook[0], 0, 1e-8);
+ EXPECT_NEAR(cameraLook[1], 0.007999744, 1e-8);
+ EXPECT_NEAR(cameraLook[2], -0.999968001, 1e-8);
+}
diff --git a/tests/data/constAngularVelocityLineScan.json b/tests/data/constAngularVelocityLineScan.json
new file mode 100644
index 0000000..c5c991c
--- /dev/null
+++ b/tests/data/constAngularVelocityLineScan.json
@@ -0,0 +1,92 @@
+{
+ "name_model" : "USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL",
+ "name_platform" : "TEST_PLATFORM",
+ "name_sensor" : "TEST_SENSOR",
+ "center_ephemeris_time": 1000.0,
+ "starting_ephemeris_time": 999.2,
+ "line_scan_rate": [
+ [0.5, -0.8, 0.1]
+ ],
+ "detector_sample_summing": 1,
+ "detector_line_summing": 1,
+ "t0_ephemeris": -0.8,
+ "dt_ephemeris": 0.2,
+ "t0_quaternion": -0.8,
+ "dt_quaternion": 0.2,
+ "focal2pixel_lines": [0.0, 10.0, 0.0],
+ "focal2pixel_samples": [0.0, 0.0, 10.0],
+ "focal_length_model": {
+ "focal_length": 50
+ },
+ "image_lines": 16,
+ "image_samples": 16,
+ "detector_center" : {
+ "line" : 0.5,
+ "sample" : 8.0
+ },
+ "interpolation_method": "lagrange",
+ "optical_distortion": {
+ "radial": {
+ "coefficients": [0.0, 0.0, 0.0]
+ }
+ },
+ "radii": {
+ "semimajor": 10,
+ "semiminor": 10,
+ "unit":"m"
+ },
+ "reference_height": {
+ "maxheight": 1,
+ "minheight": -1,
+ "unit": "m"
+ },
+ "sensor_position": {
+ "unit": "m",
+ "positions": [
+ [1000.0, 0.0, 0.0],
+ [1000.0, 0.0, 0.0],
+ [1000.0, 0.0, 0.0],
+ [1000.0, 0.0, 0.0],
+ [1000.0, 0.0, 0.0],
+ [1000.0, 0.0, 0.0],
+ [1000.0, 0.0, 0.0],
+ [1000.0, 0.0, 0.0],
+ [1000.0, 0.0, 0.0]
+ ],
+ "velocities": [
+ [0.0, 10.0, 0.0],
+ [0.0, 10.0, 0.0],
+ [0.0, 10.0, 0.0],
+ [0.0, 10.0, 0.0],
+ [0.0, 10.0, 0.0],
+ [0.0, 10.0, 0.0],
+ [0.0, 10.0, 0.0],
+ [0.0, 10.0, 0.0],
+ [0.0, 10.0, 0.0]
+ ]
+ },
+ "sensor_orientation": {
+ "quaternions": [
+ [0.0, -0.660435174378928, 0, 0.750883067090392],
+ [0.0, -0.672364256957188, 0, 0.740220444169444],
+ [0.0, -0.68412121795764, 0, 0.729368328857344],
+ [0.0, -0.695703047662478, 0, 0.718329499236346],
+ [0.0, -0.707106781186547, 0.0, 0.707106781186547],
+ [0.0, -0.718329499236346, 0, 0.695703047662478],
+ [0.0, -0.729368328857344, 0, 0.68412121795764],
+ [0.0, -0.740220444169444, 0, 0.672364256957188],
+ [0.0, -0.750883067090392, 0, 0.660435174378928]
+ ]
+ },
+ "starting_detector_line": 0,
+ "starting_detector_sample": 0,
+ "sun_position": {
+ "positions": [
+ [100.0, 100.0, 0.0]
+ ],
+ "velocities": [
+ [0.0, 0.0, 0.0]
+ ],
+ "unit": "m"
+ }
+}
diff --git a/tests/data/constVelocityLineScan.json b/tests/data/constVelocityLineScan.json
index 7f910ef..142a4b6 100644
--- a/tests/data/constVelocityLineScan.json
+++ b/tests/data/constVelocityLineScan.json
@@ -3,27 +3,26 @@
"name_platform" : "TEST_PLATFORM",
"name_sensor" : "TEST_SENSOR",
"center_ephemeris_time": 1000.0,
- "starting_ephemeris_time": 992.0,
+ "starting_ephemeris_time": 999.2,
"line_scan_rate": [
- [0.5, 992.0, 0.1]
+ [0.5, -0.8, 0.1]
],
"detector_sample_summing": 1,
"detector_line_summing": 1,
- "t0_ephemeris": -8.0,
+ "t0_ephemeris": -0.8,
"dt_ephemeris": 0.2,
- "t0_quaternion": -8.0,
+ "t0_quaternion": -0.8,
"dt_quaternion": 0.2,
"focal2pixel_lines": [0.0, 10.0, 0.0],
"focal2pixel_samples": [0.0, 0.0, 10.0],
"focal_length_model": {
- "focal_length": 50,
- "focal_epsilon": 1.0
+ "focal_length": 50
},
"image_lines": 16,
"image_samples": 16,
"detector_center" : {
"line" : 0.5,
- "sample" : 7.5
+ "sample" : 8.0
},
"interpolation_method": "lagrange",
"optical_distortion": {
diff --git a/tests/data/empty.json b/tests/data/empty.json
new file mode 100644
index 0000000..0967ef4
--- /dev/null
+++ b/tests/data/empty.json
@@ -0,0 +1 @@
+{}
diff --git a/tests/data/simpleFramerISD.json b/tests/data/simpleFramerISD.json
index 63fb799..b22c3e8 100644
--- a/tests/data/simpleFramerISD.json
+++ b/tests/data/simpleFramerISD.json
@@ -1,5 +1,8 @@
{
"name_model" : "USGS_ASTRO_FRAME_SENSOR_MODEL",
+ "name_sensor": "SENSOR",
+ "name_platform": "CarryMe",
+ "image_identifier": "TEST_IMAGE",
"center_ephemeris_time": 50,
"dt_ephemeris": 100,
"focal2pixel_lines": [0.0, 0.0, 10.0],
--
GitLab