Preliminary camera model changes and faked distortion for Hayabusa 2 (#688) (fe79ae4f) · Commits · aflab / astrogeology / ISIS3

isis/src/hayabusa2/objs/Hyb2OncCamera/Hyb2OncCamera.cpp

+43 −22

Original line number	Diff line number	Diff line
		/** This is free and unencumbered software released into the public domain.

		The authors of ISIS do not claim copyright on the contents of this file.
		For more details about the LICENSE terms and the AUTHORS, you will
		find files of those names at the top level of this repository. **/

		/* SPDX-License-Identifier: CC0-1.0 */
		/**
		* @file
		*
		* 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 "Hyb2OncCamera.h"

		@@ -75,7 +87,7 @@ namespace Isis {
		// Setup focal plane map
		CameraFocalPlaneMap *focalMap = new CameraFocalPlaneMap(this, naifIkCode());

		// BORESIGHT SAMPLE AND LINE still need to be added to the IAK
		// Set the boresight sample and line
		double bLines = Spice::getDouble("INS" + toString(naifIkCode()) + "_BORESIGHT_LINE");
		double bSamples = Spice::getDouble("INS" + toString(naifIkCode()) + "_BORESIGHT_SAMPLE");

		@@ -90,8 +102,17 @@ namespace Isis {
		detMap->SetDetectorSampleSumming(binning);

		// Setup distortion map (use default for now)
		// Level L2a and L2b images have not been corrected for optical distortion (false)
		// Level L2d and up images have been corrected for optical distortion (true)
		// Create an appropriate DistortionMap based on the processing level IDed in the ingestion app
		if ((!inst.hasKeyword("DistortionCorrection")) \|\|
		(inst.hasKeyword("DistortionCorrection") && toBool(inst["DistortionCorrection"]) == false)) {
		CameraDistortionMap *distortionMap = new Hyb2OncDistortionMap(this);
		distortionMap->SetDistortion(naifIkCode());
		}
		else {
		new CameraDistortionMap(this);
		}

		// Setup the ground and sky map
		new CameraGroundMap(this);

isis/src/hayabusa2/objs/Hyb2OncCamera/Hyb2OncCamera.h

+25 −8

Original line number	Diff line number	Diff line
		#ifndef Hyb2OncCamera_h
		#define Hyb2OncCamera_h

		/** This is free and unencumbered software released into the public domain.

		The authors of ISIS do not claim copyright on the contents of this file.
		For more details about the LICENSE terms and the AUTHORS, you will
		find files of those names at the top level of this repository. **/

		/* SPDX-License-Identifier: CC0-1.0 */
		/**
		* @file
		*
		* 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 "FramingCamera.h"

		@@ -22,6 +34,11 @@ namespace Isis {
		*
		* @internal
		* @history 2017-07-07 Kristin Berry - Original version
		* @history 2018-07-05 Lucille Le Corre - Added support for sub-frame images
		* @history 2018-11-03 Stuart Sides - Changed sub-frame/cropped image support to use the
		* AlphaCube class. The ingestion program creates the
		* AlphaCube lable group when necessary. Modified to work
		* with L2c and L2d image (i.e., undistorted)
		*
		*/
		class Hyb2OncCamera : public FramingCamera {

isis/src/hayabusa2/objs/Hyb2OncCamera/Hyb2OncDistortionMap.cpp

+98 −65

Original line number	Diff line number	Diff line
		/** This is free and unencumbered software released into the public domain.

		The authors of ISIS do not claim copyright on the contents of this file.
		For more details about the LICENSE terms and the AUTHORS, you will
		find files of those names at the top level of this repository. **/

		/* SPDX-License-Identifier: CC0-1.0 */
		/**
		* @file
		* $Revision: 1.4 $
		* $Date: 2008/02/21 16:04:33 $
		*
		* 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 <QDebug>
		#include <QtGlobal>
		@@ -60,28 +75,21 @@ namespace Isis {
		p_focalPlaneX = dx;
		p_focalPlaneY = dy;

		// reducing to principal point offset (xp,yp)
		double x = dx;// - p_xp;
		double y = dy;// - p_yp;
		//
		// Get the distance from the focal plane center and if we are close
		// then skip the distortion (this prevents division by zero)
		double r = (x * x) + (y * y);
		double r2 = r*r;
		double x = dx;
		double y = dy;

		double r2 = (x * x) + (y * y);
		double r4 = r2*r2;
		if (r <= 1.0E-6) {
		p_undistortedFocalPlaneX = dx;
		p_undistortedFocalPlaneY = dy;
		return true;
		}

		// apply distortion correction
		// r = x^2 + y^2
		// rprime = L0r + L1r^3 + L2*r^5, where Li are distortion coeffs
		// Apply distortion correction
		// https://www.darts.isas.jaxa.jp/pub/hayabusa2/onc_bundle/browse/
		// r2 = x^2 + y^2
		// rprime = L1r + L2r^3 + L3*r^5, where Li are distortion coeffs
		// "dr" is rprime divided by r, used to reduce operations
		double dr = p_odk[0] + p_odk[1] * r2 + p_odk[2] * r4;
		p_undistortedFocalPlaneX = dr * x;
		p_undistortedFocalPlaneY = dr * y;

		return true;
		}

		@@ -101,52 +109,77 @@ namespace Isis {
		* @param uy undistorted focal plane y in millimeters
		*
		* @return if the conversion was successful
		* @see SetDistortion
		* @todo Generalize polynomial equation
		* @todo Figure out a better solution for divergence condition
		*/
		bool Hyb2OncDistortionMap::SetUndistortedFocalPlane(const double ux,
		const double uy) {

		// Image coordinates prior to introducing distortion
		p_undistortedFocalPlaneX = ux;
		p_undistortedFocalPlaneY = uy;

		// Compute the distance from the focal plane center and if we are
		// close to the center then no distortion is required

		bool converged = false;
		int iteration = 0;
		double tolMilliMeters = p_camera->PixelPitch() / 100.0;
		double x = ux;
		double y = uy;
		double r = (x * x) + (y * y);
		if (r <= 1.0E-6) {
		// TEMP GET IT COMPILING REMOVE BEFORE REAL WORK
		if (1 == 1) {
		p_focalPlaneX = ux;
		p_focalPlaneY = uy;
		return true;
		}
		double rPrevious = r;

		while (!converged && qAbs(r - rPrevious) > tolMilliMeters) {
		double r2 = r*r;
		double r4 = r2*r2;
		double dr = p_odk[0] + p_odk[1] * r2 + p_odk[2] * r4;
		rPrevious = r;
		x = dr * x;
		y = dr * y;
		r = xx + yy;

		iteration++;
		if (iteration > 50) {
		converged = false;
		double xt = ux;
		double yt = uy;

		double xx, yy, rr, rrrr, dr;
		double xdistortion, ydistortion;
		double xdistorted, ydistorted;
		double xprevious, yprevious;

		xprevious = 1000000.0;
		yprevious = 1000000.0;

		double tolerance = 0.000001;

		bool bConverged = false;

		// Iterating to introduce distortion...
		// We stop when the difference between distorted coordinates
		// in successive iterations is below the given tolerance
		for (int i = 0; i < 50; i++) {
		xx = xt * xt;
		yy = yt * yt;
		rr = xx + yy;
		rrrr = rr * rr;

		// Radial distortion
		// dr is the radial distortion contribution
		dr = p_odk[0] + p_odk[1] * rr + p_odk[2] * rrrr;

		// Distortion at the current point location
		xdistortion = xt * dr;
		ydistortion = yt * dr;

		// updated image coordinates
		xt = ux - xdistortion;
		yt = uy - ydistortion;

		// distorted point corrected for principal point
		xdistorted = xt;
		ydistorted = yt;

		// check for convergence
		if ((fabs(xt - xprevious) < tolerance) && (fabs(yt - yprevious) < tolerance)) {
		bConverged = true;
		break;
		}
		}

		converged = true;
		p_focalPlaneX = x;
		p_focalPlaneY = y;
		xprevious = xt;
		yprevious = yt;
		}

		return converged;
		if (bConverged) {
		p_focalPlaneX = xdistorted;
		p_focalPlaneY = ydistorted;
		}

		return bConverged;
		}
		}

isis/src/hayabusa2/objs/Hyb2OncCamera/Hyb2OncDistortionMap.h

+4 −0

Original line number	Diff line number	Diff line
		@@ -26,6 +26,10 @@ namespace Isis {
		* @internal
		* @history 2017-07-11 Jeannie Backer and Ian Humphrey - Original version.
		*
		* @history 2018-12-11 Stuart Sides - Loop in
		* SetUndistortedFocalPlane was not being exicuted.
		* Modified to use both itteration and convergence.
		*
		*/
		class Hyb2OncDistortionMap : public CameraDistortionMap {
		public:

isis/src/hayabusa2/objs/Hyb2OncCamera/unitTest.cpp

+78 −39

Original line number	Diff line number	Diff line
		/** This is free and unencumbered software released into the public domain.

		The authors of ISIS do not claim copyright on the contents of this file.
		For more details about the LICENSE terms and the AUTHORS, you will
		find files of those names at the top level of this repository. **/

		/* SPDX-License-Identifier: CC0-1.0 */

		/**
		* @file
		*
		* 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 <QDebug>

		#include <iomanip>
		@@ -23,9 +34,14 @@ find files of those names at the top level of this repository. **/
		using namespace std;
		using namespace Isis;

		void testCamera(Cube &c, double knownLat, double knownLon);
		void testCamera(Cube &c, double knownLat, double knownLon,
		double s1, double l1,
		double s2, double l2,
		double s3, double l3,
		double s4, double l4);
		void testLineSamp(Camera *cam, double sample, double line);


		/**
		* Unit Test for the Hayabusa2 ONC camera.
		*
		@@ -45,26 +61,33 @@ int main(void) {
		qDebug() << "";
		qDebug() << "----------------------------------------------";
		qDebug() << "Test for Telecopic Camera...";
		double knownLat = -54.63487131147738;
		double knownLon = 40.43436155430055;
		Cube c("$ISISTESTDATA/isis/src/hayabusa2/unitTestData/hyb2_onc_20151204_041012_tbf_l2a.fit.cub", "r");
		testCamera(c, knownLat, knownLon);
		double knownLat = -47.23506367622795;
		double knownLon = 45.06713880290044;
		Cube c("$hayabusa2/testData/hyb2_onc_20151204_041012_tbf_l2a.fit.cub", "r");
		testCamera(c, knownLat, knownLon,
		512.5, 512.5, 602.0, 334.0, 378.0, 557.0, 594.0, 580.0);
		// 362.0, 352.0, 602.0, 334.0, 378.0, 557.0, 594.0, 580.0);

		qDebug() << "";
		qDebug() << "----------------------------------------------";
		qDebug() << "Test for W1 Camera...";
		knownLat = -50.11857108654684;
		knownLon = 91.03535388676204;
		Cube w1("$ISISTESTDATA/isis/src/hayabusa2/unitTestData/hyb2_onc_20151204_045429_w1f_l2a.fit_crop.cub", "r");
		testCamera(w1, knownLat, knownLon);
		knownLat = -50.66777190122887;
		knownLon = 97.68522302461859;
		Cube w1("$hayabusa2/testData/hyb2_onc_20151204_045429_w1f_l2a.fit_crop.cub", "r");
		testCamera(w1, knownLat, knownLon,
		21.0, 20.0, 31.0, 11.0, 16.0, 29.0, 32.0, 28.0);
		// 16.0, 14.0, 31.0, 11.0, 16.0, 29.0, 32.0, 28.0);

		qDebug() << "";
		qDebug() << "----------------------------------------------";
		qDebug() << "Test for W2 Camera...";
		knownLat = 25.38911363842043;
		knownLon = 90.86547761107917;
		Cube w2("$ISISTESTDATA/isis/src/hayabusa2/unitTestData/hyb2_onc_20151203_072958_w2f_l2a.fit_crop.cub", "r");
		testCamera(w2, knownLat, knownLon);
		// knownLat = 4.676892803564044;
		// knownLon = -12.46121470106279;
		knownLat = 30.06610722049293;
		knownLon = 78.40416492096558;
		Cube w2("$hayabusa2/testData/hyb2_onc_20151203_072958_w2f_l2a.fit_crop.cub", "r");
		testCamera(w2, knownLat, knownLon,
		51.0, 42.0, 173.0, 21.0, 54.0, 149.0, 174.0, 155.0);

		}
		catch (IException &e) {
		@@ -74,7 +97,11 @@ int main(void) {
		}


		void testCamera(Cube &c, double knownLat, double knownLon) {
		void testCamera(Cube &c, double knownLat, double knownLon,
		double s1, double l1,
		double s2, double l2,
		double s3, double l3,
		double s4, double l4) {
		Hyb2OncCamera cam = (Hyb2OncCamera ) CameraFactory::Create(c);
		qDebug() << "FileName: " << FileName(c.fileName()).name();
		qDebug() << "CK Frame: " << cam->instrumentRotation()->Frame();
		@@ -108,16 +135,16 @@ void testCamera(Cube &c, double knownLat, double knownLon) {

		// Test all four corners to make sure the conversions are right
		qDebug() << "For upper left corner ...";
		testLineSamp(cam, 1.0, 1.0);
		testLineSamp(cam, s1, l1);

		qDebug() << "For upper right corner ...";
		testLineSamp(cam, cam->Samples(), 1.0);
		testLineSamp(cam, s2, l2);

		qDebug() << "For lower left corner ...";
		testLineSamp(cam, 1.0, cam->Lines());
		testLineSamp(cam, s3, l3);

		qDebug() << "For lower right corner ...";
		testLineSamp(cam, cam->Samples(), cam->Lines());
		testLineSamp(cam, s4, l4);

		qDebug() << "For center pixel position ...";

		@@ -147,16 +174,27 @@ void testLineSamp(Camera *cam, double sample, double line) {
		bool success = cam->SetImage(sample, line);

		if (success) {
		success = cam->SetUniversalGround(cam->UniversalLatitude(), cam->UniversalLongitude());
		double lat = cam->UniversalLatitude();
		double lon = cam->UniversalLongitude();
		success = cam->SetUniversalGround(lat, lon);
		// success = cam->SetUniversalGround(cam->UniversalLatitude(), cam->UniversalLongitude());
		}

		if (success) {
		double deltaSamp = sample - cam->Sample();
		double deltaLine = line - cam->Line();
		if(fabs(deltaSamp) < 0.001) deltaSamp = 0;
		if(fabs(deltaLine) < 0.001) deltaLine = 0;
		qDebug() << "DeltaSample = " << deltaSamp;
		qDebug() << "DeltaLine = " << deltaLine;
		if(fabs(deltaSamp) < 0.001) {
		qDebug() << "Delta Sample less than 0.001";
		}
		else {
		qDebug() << "Delta Sample larger than expected " << deltaSamp;
		}
		if(fabs(deltaLine) < 0.001) {
		qDebug() << "Delta Line less than 0.001";
		}
		else {
		qDebug() << "Delta Line larger than expected " << deltaLine;
		}
		qDebug() << "";
		}
		else {
		@@ -165,3 +203,4 @@ void testLineSamp(Camera *cam, double sample, double line) {
		qDebug() << "";
		}
		}