Implemented Rosetta OSIRIS distortion maps. Fixes #4496. (6132911b) · Commits · aflab / astrogeology / ISIS3

isis/src/rosetta/objs/RosettaOsirisCamera/Camera.plugin

+2 −2

Original line number	Diff line number	Diff line
		Group = ROSETTA-ORBITER/OSINAC
		Version = 1
		Version = 2
		Library = RosettaOsirisCamera
		Routine = RosettaOsirisCameraPlugin
		EndGroup

		Group = ROSETTA-ORBITER/OSIWAC
		Version = 1
		Version = 2
		Library = RosettaOsirisCamera
		Routine = RosettaOsirisCameraPlugin
		EndGroup

isis/src/rosetta/objs/RosettaOsirisCamera/RosettaOsirisCamera.cpp

+53 −11

Original line number	Diff line number	Diff line
		@@ -20,14 +20,19 @@

		#include "RosettaOsirisCamera.h"

		#include <QDebug>
		#include <QFile>
		#include <QString>

		#include "CameraDetectorMap.h"
		#include "CameraDistortionMap.h"
		#include "CameraFocalPlaneMap.h"
		#include "CameraGroundMap.h"
		#include "CameraSkyMap.h"
		#include "IString.h"
		#include "iTime.h"
		#include "FileName.h"
		#include "NaifStatus.h"
		#include "Preference.h"

		using namespace std;

		@@ -70,18 +75,10 @@ namespace Isis {

		// Setup focal plane map. The class will read data from the instrument addendum kernel to pull
		// out the affine transforms from detector samp,line to focal plane x,y.
		// cout << "Setting up FocalPlaneMap...\n";
		CameraFocalPlaneMap *focalMap = new CameraFocalPlaneMap(this, naifIkCode());

		// The boresight position recorded in the IK is zero-based and therefore needs to be adjusted
		// for ISIS
		// Don't know if this is true for OSIRIS. For now, we'll keep as is and see if things look off -Sasha
		double boresightSample = Spice::getDouble("INS" + ikCode + "_BORESIGHT",0) + 1.0;
		double boresightLine = Spice::getDouble("INS" + ikCode + "_BORESIGHT",1) + 1.0;
		focalMap->SetDetectorOrigin(boresightSample,boresightLine); //Presumably, don't need to worry about z (?)

		new CameraDetectorMap(this);
		new CameraDistortionMap(this);
		RosettaOsirisCameraDistortionMap* distortionMap = new RosettaOsirisCameraDistortionMap(this);

		// Setup the ground and sky map
		new CameraGroundMap(this);
		@@ -96,6 +93,23 @@ namespace Isis {
		// double stop = getClockTime(clockStopCount).Et();
		double exposureTime = (double) inst["ExposureDuration"];

		// Setup the distortion map
		PvlGroup &BandBin = lab.findGroup("BandBin", Pvl::Traverse);
		QString filterNumber = BandBin["FilterNumber"];
		initDistortion(ikCode, distortionMap);
		distortionMap->setPixelPitch(pixelPitch);

		// The boresight position depends on the filter. They are all defined as
		// offsets from the middle of the ccd.
		double referenceSample = Spice::getDouble("INS" + ikCode + "_BORESIGHT",0) + 1.0;
		double referenceLine = Spice::getDouble("INS" + ikCode + "_BORESIGHT",1) + 1.0;
		// The offsets in the IAK are based on the S/C frame, not the camera frame
		// For now, do not adjust based on filter. -JAM
		// referenceSample += Spice::getDouble("INS" + ikCode + "_FILTER_" + filterNumber + "_DX");
		// referenceLine += Spice::getDouble("INS" + ikCode + "_FILTER_" + filterNumber + "_DY");
		focalMap->SetDetectorOrigin(referenceSample, referenceLine);
		distortionMap->setBoresight(referenceSample, referenceLine);

		iTime centerTime = start + (exposureTime / 2.0);
		setTime( centerTime );

		@@ -106,6 +120,7 @@ namespace Isis {
		return;
		}


		/**
		* Returns the shutter open and close times. The LORRI camera doesn't use a shutter to start and
		* end an observation, but this function is being used to get the observation start and end times,
		@@ -128,6 +143,33 @@ namespace Isis {
		double exposureDuration) {
		return FramingCamera::ShutterOpenCloseTimes(time, exposureDuration);
		}


		/**
		* Initialize the distortion map using the paramters from the NAIF SPICE kernels.
		*
		* @param ikCode The NAIF IK code of the instrument
		* @param[out] distortionMap The distortion map that will be initialized
		*/
		void RosettaOsirisCamera::initDistortion(QString ikCode,
		RosettaOsirisCameraDistortionMap *distortionMap) {

		// Initialize matrices
		LinearAlgebra::Matrix toUnDistX = LinearAlgebra::zeroMatrix(4, 4);
		LinearAlgebra::Matrix toUnDistY = LinearAlgebra::zeroMatrix(4, 4);

		// Fill matrices from the kernels
		for (int i = 0; i < 4; i++) {
		for (int j = 0; j < 4; j++) {
		toUnDistX(i, j) = Spice::getDouble("INS" + ikCode + "_TO_UNDISTORTED_X", 4 * i + j);
		toUnDistY(i, j) = Spice::getDouble("INS" + ikCode + "_TO_UNDISTORTED_Y", 4 * i + j);
		}
		}

		// Save the matrices
		distortionMap->setUnDistortedXMatrix(toUnDistX);
		distortionMap->setUnDistortedYMatrix(toUnDistY);
		}
		}

		/**

isis/src/rosetta/objs/RosettaOsirisCamera/RosettaOsirisCamera.h

+10 −2

Original line number	Diff line number	Diff line
		@@ -22,6 +22,11 @@

		#include "FramingCamera.h"

		#include <QXmlStreamReader>

		#include "LinearAlgebra.h"
		#include "RosettaOsirisCameraDistortionMap.h"

		namespace Isis {
		/**
		* This is the camera model for the Osiris NAC Framing Camera
		@@ -32,7 +37,8 @@ namespace Isis {
		* @author 2015-05-21 Sasha Brownsberger
		*
		* @internal
		*
		* @history 2015-05-21 Sasha Brownsberger - Original Version.
		* @history 2017-06-02 Jesse Mapel - Added a distortion map Fixes #4496.
		*/
		class RosettaOsirisCamera : public FramingCamera {
		public:
		@@ -81,6 +87,8 @@ namespace Isis {
		* Kernel Reference ID
		*/
		virtual int SpkReferenceId() const { return (1); }

		void initDistortion(QString ikCode, RosettaOsirisCameraDistortionMap *distortionMap);
		};
		};
		#endif

isis/src/rosetta/objs/RosettaOsirisCamera/RosettaOsirisCamera.truth

+36 −0

Original line number	Diff line number	Diff line
		Unit Test for RosettaOsirisCameraDistortionMap...
		Create default distortion map

		Set distorted coordinates
		Distorted coordinates: (0.5, 0.75)
		Undistorted coordinates: (0.5, 0.75)

		Set undistorted coordinates
		Distorted coordinates: (0.25, 0.1)
		Undistorted coordinates: (0.25, 0.1)

		Modify the coefficient matrices
		Set distorted coordinates
		Distorted coordinates: (1.0, 1.0)
		Undistorted coordinates: (10.0, -10.0)

		Set distorted coordinates
		Distorted coordinates: (0.0, 1.0)
		Undistorted coordinates: (1.0, -10.0)

		Set distorted coordinates
		Distorted coordinates: (1.0, 0.0)
		Undistorted coordinates: (1.0, -1.0)

		Set undistorted coordinates
		Distorted coordinates: (-1.04607821983755e-08, -0.85196984405198)
		Undistorted coordinates: (1.0, 1.0)

		Set undistorted coordinates
		Distorted coordinates: (0.71109275805717, -0.85196984405198)
		Undistorted coordinates: (0.0, 1.0)

		Set undistorted coordinates
		Distorted coordinates: (0.71109275805717, -0.85196984405198)
		Undistorted coordinates: (1.0, 0.0)

		Unit Test for RosettaOsirisCamera...
		FileName: n20100710t154539230id20f22.cub
		CK Frame: -226111

isis/src/rosetta/objs/RosettaOsirisCamera/RosettaOsirisCameraDistortionMap.cpp

0 → 100644

+374 −0

Original line number	Diff line number	Diff line
		/**
		* @file
		* $Revision: 1.5 $
		* $Date: 2010/05/12 23:28:12 $
		*
		* Unless noted otherwise, the portions of Isis written by the USGS are
		* public domain. See individual third-party library and package descriptions
		* for intellectual property information, user agreements, and related
		* information.
		*
		* Although Isis has been used by the USGS, no warranty, expressed or
		* implied, is made by the USGS as to the accuracy and functioning of such
		* software and related material nor shall the fact of distribution
		* constitute any such warranty, and no responsibility is assumed by the
		* USGS in connection therewith.
		*
		* For additional information, launch
		* $ISISROOT/doc//documents/Disclaimers/Disclaimers.html
		* in a browser or see the Privacy & Disclaimers page on the Isis website,
		* http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on
		* http://www.usgs.gov/privacy.html.
		*/

		#include "RosettaOsirisCameraDistortionMap.h"

		using namespace std;

		namespace Isis {
		/**
		* Create a camera distortion map. This class maps between distorted
		* and undistorted focal plane x/y's. The default mapping is the
		* identity, that is, the focal plane x/y and undistorted focal plane
		* x/y will be identical.
		*
		* @param parent the parent camera that will use this distortion map
		*
		*/
		RosettaOsirisCameraDistortionMap::RosettaOsirisCameraDistortionMap(Camera *parent) :
		CameraDistortionMap(parent) {
		// Initialize to the identity transform
		m_toUnDistortedX = LinearAlgebra::zeroMatrix(4, 4);
		m_toUnDistortedY = LinearAlgebra::zeroMatrix(4, 4);
		m_toUnDistortedX(0, 1) = 1.0;
		m_toUnDistortedY(1, 0) = 1.0;

		m_boresightSample = 0.0;
		m_boresightLine = 0.0;
		m_pixelPitch = 1.0;
		}

		/**
		* Compute undistorted focal plane x/y given a distorted focal plane x/y.
		*
		* The distortion is modeled by pixelspace polynomials. The polynomials use
		* zero-based pixel space with the origin at the top left corner of the
		* image, so the input focal plane coordinates are converted to pixel
		* coordinates using the boresight location and pixel pitch. After
		* computation, they are converted back into focal plane coordinates by the
		* inverse process.
		*
		* Given a set of distorted pixel coordinates (dx, dy), the undistorted pixel
		* coordinates (ux, uy) are computed as:
		* @f[ (ux, uy) = F(dx, dy) = ( \sum_{i=0}^3 \sum_{j=0}^3 C_{i,j}^x dx^i dy^j,
		* \sum_{i=0}^3 \sum_{j=0}^3 C_{i,j}^y dx^i dy^j) @f] where @f$ C_{i,j}^y @f$
		* and @f$ C_{i,j}^y @f$ are the @f$ (i,j)^{\text{th}} @f$ coefficients of
		* the @f$ x @f$ and @f$ y @f$ polynomials respectively.
		*
		* This calculation is actually performed as follows:
		* @f[ (ux, uy) =
		* \left( C^x \begin{bmatrix} 1 \\ dx \\ dx^2 \\ dx^3 \end{bmatrix}
		* \cdot \begin{bmatrix} 1 \\ dy \\ dy^2 \\ dy^3 \end{bmatrix},
		* C^y \begin{bmatrix} 1 \\ dx \\ dx^2 \\ dx^3 \end{bmatrix}
		* \cdot \begin{bmatrix} 1 \\ dy \\ dy^2 \\ dy^3 \end{bmatrix} \right)
		* @f]
		* where @f$ C^x @f$ and @f$ C^y @f$ are the @f$ x @f$ and @f$ y @f$
		* coefficient matrices respectively.
		*
		* @param dx distorted focal plane x in millimeters
		* @param dy distorted focal plane y in millimeters
		*
		* @return @b bool if the conversion was successful
		*/
		bool RosettaOsirisCameraDistortionMap::SetFocalPlane(const double dx, const double dy) {
		p_focalPlaneX = dx;
		p_focalPlaneY = dy;

		// The equations are in pixel coordinates so convert
		double dxPixel = focalXToLine(dx);
		double dyPixel = focalYToSample(dy);

		LinearAlgebra::Vector xTerms = LinearAlgebra::zeroVector(4);
		xTerms(0) = 1.0;
		xTerms(1) = dxPixel;
		xTerms(2) = dxPixel*dxPixel;
		xTerms(3) = dxPixeldxPixeldxPixel;

		LinearAlgebra::Vector yTerms = LinearAlgebra::zeroVector(4);
		yTerms(0) = 1.0;
		yTerms(1) = dyPixel;
		yTerms(2) = dyPixel*dyPixel;
		yTerms(3) = dyPixeldyPixeldyPixel;

		double ux = LinearAlgebra::dotProduct(
		LinearAlgebra::multiply(m_toUnDistortedX, xTerms),
		yTerms);
		double uy = LinearAlgebra::dotProduct(
		LinearAlgebra::multiply(m_toUnDistortedY, xTerms),
		yTerms);

		p_undistortedFocalPlaneX = lineToFocalX(ux);
		p_undistortedFocalPlaneY = sampleToFocalY(uy);

		return true;
		}

		/**
		* Compute distorted focal plane x/y given an undistorted focal plane x/y.
		*
		* The conversion is performed using Newton's method to find distorted
		* coordinates whose undistorted coordinates are within @f$ 10^{-7} @f$
		* pixels of the input undistorted coordinates. The input undistorted
		* coordinates are used as an initial guess for the distorted coordinates.
		*
		* Given a set of undistorted pixel coordinates (ux, uy), the object function
		* is:
		* @f[
		* G(dx, dy) = (ux, uy) - F(dx, dy)
		* @f]
		* where @f$ F @f$ is the transformation from distorted to undistorted pixel
		* coordinates.
		*
		* Then, the negative jacobian is:
		* @f[
		* -J_G(dx, dy) = \begin{bmatrix}
		* C^x \begin{bmatrix} 0 \\ 1 \\ 2dx \\ 3dx^2 \end{bmatrix}
		* \cdot \begin{bmatrix} 1 \\ dy \\ dy^2 \\ dy^3 \end{bmatrix} &
		* C^x \begin{bmatrix} 1 \\ dx \\ dx^2 \\ dx^3 \end{bmatrix}
		* \cdot \begin{bmatrix} 0 \\ 1 \\ 2dy \\ 3dy^2 \end{bmatrix} \\
		* C^y \begin{bmatrix} 0 \\ 1 \\ 2dx \\ 3dx^2 \end{bmatrix}
		* \cdot \begin{bmatrix} 1 \\ dy \\ dy^2 \\ dy^3 \end{bmatrix} &
		* C^y \begin{bmatrix} 1 \\ dx \\ dx^2 \\ dx^3 \end{bmatrix}
		* \cdot \begin{bmatrix} 0 \\ 1 \\ 2dy \\ 3dy^2 \end{bmatrix} \end{bmatrix}
		* @f]
		*
		* @param ux undistorted focal plane x in millimeters
		* @param uy undistorted focal plane y in millimeters
		*
		* @return @b bool if the conversion was successful
		*/
		bool RosettaOsirisCameraDistortionMap::SetUndistortedFocalPlane(const double ux, const double uy) {
		// image coordinates prior to introducing distortion
		p_undistortedFocalPlaneX = ux;
		p_undistortedFocalPlaneY = uy;

		// The equations are in pixel coordinates so convert
		double uxPixel = focalXToLine(ux);
		double uyPixel = focalYToSample(uy);

		// Loop setup
		double tolerance = 1e-7;
		int iteration = 1;
		int maxIterations = 20;
		bool done = false;

		LinearAlgebra::Vector undistortedCoordinate = LinearAlgebra::zeroVector(2);
		undistortedCoordinate(0) = uxPixel;
		undistortedCoordinate(1) = uyPixel;
		// Use the undistorted coordinate as the initial point
		LinearAlgebra::Vector distortedCoordinate = undistortedCoordinate;

		LinearAlgebra::Vector xTerms = LinearAlgebra::zeroVector(4);
		LinearAlgebra::Vector yTerms = LinearAlgebra::zeroVector(4);
		LinearAlgebra::Vector delXTerms = LinearAlgebra::zeroVector(4);
		LinearAlgebra::Vector delYTerms = LinearAlgebra::zeroVector(4);

		LinearAlgebra::Vector objectFuncValue = LinearAlgebra::zeroVector(2);
		LinearAlgebra::Vector updateStep = LinearAlgebra::zeroVector(2);
		LinearAlgebra::Matrix negJacobian = LinearAlgebra::zeroMatrix(2, 2);
		LinearAlgebra::Matrix negInvJacobian = LinearAlgebra::zeroMatrix(2, 2);

		// Compute the term vectors
		// These are recomputed at the END of each iteration, so compute them
		// before the first iteration
		xTerms(0) = 1.0;
		xTerms(1) = distortedCoordinate(0);
		xTerms(2) = distortedCoordinate(0)*distortedCoordinate(0);
		xTerms(3) = distortedCoordinate(0)distortedCoordinate(0)distortedCoordinate(0);

		yTerms(0) = 1.0;
		yTerms(1) = distortedCoordinate(1);
		yTerms(2) = distortedCoordinate(1)*distortedCoordinate(1);
		yTerms(3) = distortedCoordinate(1)distortedCoordinate(1)distortedCoordinate(1);

		delXTerms(1) = 1.0;
		delXTerms(2) = 2.0*distortedCoordinate(0);
		delXTerms(3) = 3.0distortedCoordinate(0)distortedCoordinate(0);

		delYTerms(1) = 1.0;
		delYTerms(2) = 2.0*distortedCoordinate(1);
		delYTerms(3) = 3.0distortedCoordinate(1)distortedCoordinate(1);

		// Compute the object function, the distance between the redistorted and undistorted
		objectFuncValue(0) = undistortedCoordinate(0)
		- LinearAlgebra::dotProduct( LinearAlgebra::multiply(m_toUnDistortedX, xTerms),
		yTerms );
		objectFuncValue(1) = undistortedCoordinate(1)
		- LinearAlgebra::dotProduct( LinearAlgebra::multiply(m_toUnDistortedY, xTerms),
		yTerms );

		while(!done) {

		// Compute the negative jacobian
		// The variable terms and object function value have already been
		// computed priort to checking for convergence at the end of the
		// previous iteration. These are not needed to check for convergence
		// so, compute them now.
		negJacobian(0, 0) = LinearAlgebra::dotProduct( LinearAlgebra::multiply(m_toUnDistortedX, delXTerms),
		yTerms );
		negJacobian(0, 1) = LinearAlgebra::dotProduct( LinearAlgebra::multiply(m_toUnDistortedX, xTerms),
		delYTerms );
		negJacobian(1, 0) = LinearAlgebra::dotProduct( LinearAlgebra::multiply(m_toUnDistortedY, delXTerms),
		yTerms );
		negJacobian(1, 1) = LinearAlgebra::dotProduct( LinearAlgebra::multiply(m_toUnDistortedY, xTerms),
		delYTerms );

		// Invert the negative jacobian
		// If it is not invertible, then fail
		double det = negJacobian(0, 0) * negJacobian(1, 1) - negJacobian(1, 0) * negJacobian(0, 1);
		if ( fabs(det) < 1e-15 ) {
		return false;
		}
		negInvJacobian(0, 0) = negJacobian(1, 1) / det;
		negInvJacobian(0, 1) = - negJacobian(0, 1) / det;
		negInvJacobian(1, 0) = - negJacobian(1, 0) / det;
		negInvJacobian(1, 1) = negJacobian(0, 0) / det;

		// Compute the update step
		updateStep = LinearAlgebra::multiply(negInvJacobian, objectFuncValue);

		// Update and check for convergence
		distortedCoordinate += updateStep;

		// Compute the term vectors
		xTerms(0) = 1.0;
		xTerms(1) = distortedCoordinate(0);
		xTerms(2) = distortedCoordinate(0)*distortedCoordinate(0);
		xTerms(3) = distortedCoordinate(0)distortedCoordinate(0)distortedCoordinate(0);

		yTerms(0) = 1.0;
		yTerms(1) = distortedCoordinate(1);
		yTerms(2) = distortedCoordinate(1)*distortedCoordinate(1);
		yTerms(3) = distortedCoordinate(1)distortedCoordinate(1)distortedCoordinate(1);

		delXTerms(1) = 1.0;
		delXTerms(2) = 2.0*distortedCoordinate(0);
		delXTerms(3) = 3.0distortedCoordinate(0)distortedCoordinate(0);

		delYTerms(1) = 1.0;
		delYTerms(2) = 2.0*distortedCoordinate(1);
		delYTerms(3) = 3.0distortedCoordinate(1)distortedCoordinate(1);

		// Compute the object function, the distance between the redistorted and undistorted
		objectFuncValue(0) = undistortedCoordinate(0)
		- LinearAlgebra::dotProduct( LinearAlgebra::multiply(m_toUnDistortedX, xTerms),
		yTerms );
		objectFuncValue(1) = undistortedCoordinate(1)
		- LinearAlgebra::dotProduct( LinearAlgebra::multiply(m_toUnDistortedY, xTerms),
		yTerms );

		if ( LinearAlgebra::magnitude(objectFuncValue) < tolerance \|\|
		iteration > maxIterations ) {
		done = true;
		}
		iteration++;
		}

		p_focalPlaneX = lineToFocalX( distortedCoordinate(0) );
		p_focalPlaneY = sampleToFocalY( distortedCoordinate(1) );

		return true;
		}


		/**
		* Set the matrix for converting from distorted to undistorted samples.
		*
		* @param xMat The conversion matrix
		*/
		void RosettaOsirisCameraDistortionMap::setUnDistortedXMatrix(LinearAlgebra::Matrix xMat) {
		m_toUnDistortedX = xMat;
		}


		/**
		* Set the matrix for converting from distorted to undistorted lines.
		*
		* @param yMat The conversion matrix
		*/
		void RosettaOsirisCameraDistortionMap::setUnDistortedYMatrix(LinearAlgebra::Matrix yMat) {
		m_toUnDistortedY = yMat;
		}


		/**
		* Set the boresight location for converting from focal plane coordinates to
		* pixel coordinates.
		*
		* @param sample The sample location of the boresight
		* @param line The line location of the boresight
		*/
		void RosettaOsirisCameraDistortionMap::setBoresight(double sample, double line) {
		m_boresightSample = sample;
		m_boresightLine = line;
		}


		/**
		* Set the pixel pitch for converting from focal plane coordinates to
		* pixel coordinates.
		*
		* @param pitch The pixel pitch in mm per pixel.
		*/
		void RosettaOsirisCameraDistortionMap::setPixelPitch(double pitch) {
		m_pixelPitch = pitch;
		}


		/**
		* Convert a focal plane x coordinate to a pixel space line coordinate.
		*
		* @param x The focal plane x coordinate in mm
		*
		* @return @b double The pixel space line coordinate
		*/
		double RosettaOsirisCameraDistortionMap::focalXToLine(double x) {
		return ( x / m_pixelPitch + m_boresightLine );
		}


		/**
		* Convert a focal plane y coordinate to a pixel space sample coordinate.
		*
		* @param y The focal plane y coordinate in mm
		*
		* @return @b double The pixel space sample coordinate
		*/
		double RosettaOsirisCameraDistortionMap::focalYToSample(double y) {
		return ( y / m_pixelPitch + m_boresightSample);
		}


		/**
		* Convert pixel space line coordinate to a focal plane x coordinate.
		*
		* @param line The pixel space line coordinate
		*
		* @return @b double The focal plane x coordinate in mm
		*/
		double RosettaOsirisCameraDistortionMap::lineToFocalX(double line) {
		return ( ( line - m_boresightLine ) * m_pixelPitch );
		}


		/**
		* Convert pixel space sample coordinate to a focal plane y coordinate.
		*
		* @param sample The pixel space sample coordinate
		*
		* @return @b double The focal plane y coordinate in mm
		*/
		double RosettaOsirisCameraDistortionMap::sampleToFocalY(double sample) {
		return ( ( sample - m_boresightSample) * m_pixelPitch );
		}
		}