SensorUtil Additions (#5406) (24817381) · Commits · aflab / astrogeology / ISIS3

CHANGELOG.md

+3 −0

Original line number	Diff line number	Diff line
		@@ -45,6 +45,9 @@ release.
		- Fixed <i>noproj</i> bug where missing shapemodel-related keywords (RayTraceEngine, BulletParts, Tolerance) are dropped when the output label is created. This resulted in the Bullet collision detection engine not being used. Issue: [#5377](https://github.com/USGS-Astrogeology/ISIS3/issues/5377)
		- Fixed ISIS failing to expand env variables with an "_" in them. [#5402](https://github.com/DOI-USGS/ISIS3/pull/5402)

		### Added
		- Added 8 new functions to the Sensor Utility Library: Slant Distance, Target Center Distance, Right Ascension Declination, Local Solar Time, Line Resolution, Sample Resolution, Pixel Resolution, and Solar Longitude.

		## [8.1.0] - 2024-01-08

		### Changed

SensorUtilities/include/MathUtils.h

+81 −2

Original line number	Diff line number	Diff line
		@@ -5,6 +5,8 @@

		namespace SensorUtilities {

		const double RAD2DEG = 57.29577951308232087679815481;

		/**
		* Structure for a 2-dimensional spherical ground point.
		* Latitude and longitude are in radians.
		@@ -29,6 +31,17 @@ namespace SensorUtilities {
		};


		/**
		* Structure for right ascension and declination angles.
		* Right Ascenion angle in degrees.
		* Declination angle in degrees.
		*/
		struct RaDec {
		double ra;
		double dec;
		};


		/**
		* Structure for a point in an image.
		* The line and sample origin is the upper-left corner at (0, 0).
		@@ -50,6 +63,11 @@ namespace SensorUtilities {
		double y;
		double z;

		/**
		* Default Vec constructor needed for python wrapper.
		*/
		Vec() {};

		/**
		* Construct a vector from 3 values
		*/
		@@ -68,6 +86,15 @@ namespace SensorUtilities {
		operator std::vector<double>() const;
		};

		/**
		* Structure for a 3 by 3 Matrix.
		*/
		struct Matrix {
		Vec a;
		Vec b;
		Vec c;
		};


		bool operator==(const GroundPt2D& lhs, const GroundPt2D& rhs);

		@@ -129,6 +156,16 @@ namespace SensorUtilities {
		double distance(Vec start, Vec stop);


		/**
		* Converts coordinate from radians to degrees
		*
		* @param radianLatLon coordinate in radians
		*
		* @return Coordinate in degrees
		*/
		GroundPt2D radiansToDegrees(GroundPt2D radianLatLon);


		/**
		* Convert spherical coordinates to rectangular coordinates
		*
		@@ -149,7 +186,49 @@ namespace SensorUtilities {
		*/
		GroundPt3D rectToSpherical(Vec rectangular);



		/**
		* Computes and returns the ground azimuth between the ground point and
		* another point of interest, such as the subspacecraft point or the
		* subsolar point.
		*
		* @param groundPt Ground point latitude and longitude
		* @param subPt latitude and longitude of a point of interest (e.g. subspacecraft or subsolar)
		*
		* @return double The azimuth in degrees
		*/
		double groundAzimuth(GroundPt2D groundPt, GroundPt2D subPt);


		/**
		* Find the component of a vector that is perpendicular to a second vector.
		*/
		Vec perpendicularVec(Vec aVec, Vec bVec);


		/**
		* Compute the cross product of two Vecs.
		*/
		Vec crossProduct(Vec aVec, Vec bVec);


		/**
		* Multiply a scalar and a Vec.
		*/
		Vec scaleVector(Vec vec, double scalar);


		/**
		* Find the unit vector along a Vec.
		*/
		Vec unitVector(Vec vec);


		/**
		* Multiply a Matrix with a Vec.
		*/
		Vec matrixVecProduct(Matrix mat, Vec vec);
		}

		#endif
		No newline at end of file

SensorUtilities/include/SensorUtilities.h

+96 −3

Original line number	Diff line number	Diff line
		@@ -9,7 +9,6 @@

		namespace SensorUtilities {


		/**
		* Structure for storing the state of an observer at a specific image coordinate
		* and time.
		@@ -85,8 +84,30 @@ namespace SensorUtilities {
		* Get the position for the illumination source at a given time.
		*/
		virtual Vec position(double time) = 0;
		/**
		* Get the velocity for the illumination source at a given time.
		*/
		virtual Vec velocity(double time) = 0;
		};

		/**
		* Interface for the location of the body source.
		* Implementations of this interface operate in object space.
		*/
		class Body {
		public:
		/**
		* Get the rotation matrix for the body source at a given time.
		*/
		virtual std::vector<double> rotation(double time) = 0;

		/**
		* J2000 position vector to BodyFixed position vector.
		*
		* @param pos J2000 position to get as a bodyFixed vector
		*/
		virtual Vec fixedVector(Vec pos) = 0;
		};

		/**
		* Compute the phase angle at an image point.
		@@ -127,6 +148,22 @@ namespace SensorUtilities {
		double illuminationDistance(const ImagePt &imagePoint, Sensor sensor, Shape shape, Illuminator *illuminator);


		/**
		* Return the distance between the spacecraft and surface point.
		*
		* @return the distance in meters
		*/
		double slantDistance(ImagePt imagePoint, Sensor sensor, Shape shape);


		/**
		* Return the distance between the spacecraft and center point of a body.
		*
		* @return the distance in meters
		*/
		double targetCenterDistance(ImagePt &imagePoint, Sensor sensor, Body body);


		/**
		* Compute the latitude and longitude on the body below the sensor
		* when an image coordinate was observed.
		@@ -173,7 +210,63 @@ namespace SensorUtilities {
		* @return The local radius in meters
		*/
		double localRadius(const GroundPt2D &groundPt, Shape *shape, double maxRadius=1e9);
		}

		#endif

		/**
		* Computes the right ascension angle (sky longitude) and declination angle (sky latitude).
		*
		* @returns right ascension angle in degrees and declination angle in degrees.
		*/
		RaDec rightAscensionDeclination(ImagePt &imagePoint, Sensor *sensor);


		/**
		* Computes the local solar time in hours
		*
		* @param imagePoint The line and sample
		* @param sensor The sensor model used to calculate ground coordinates
		* @param shape The shape to compute the local radius of
		* @param illuminator The illumination model
		*
		* @returns local solar time in hours
		*/
		double localSolarTime(ImagePt &imagePoint, Sensor sensor, Shape shape, Illuminator *illuminator);

		/**
		* @brief Returns the line resolution at the current position in meters.
		*
		* @return @b double The line resolution
		*/
		double lineResolution(const ImagePt &imagePoint, Sensor sensor, Shape shape, double focalLength, double pixelPitch, double lineScaleFactor=1);


		/**
		* @brief Returns the sample resolution at the current position in meters.
		*
		* @return @b double The sample resolution
		*/
		double sampleResolution(const ImagePt &imagePoint, Sensor sensor, Shape shape, double focalLength, double pixelPitch, double sampleScaleFactor=1);


		/**
		* @brief Returns the pixel resolution at the current position in meters/pixel.
		*
		* @return @b double The pixel resolution
		*/
		double pixelResolution(const ImagePt &imagePoint, Sensor sensor, Shape shape, double focalLength, double pixelPitch, double lineScaleFactor=1, double sampleScaleFactor=1);


		/**
		*
		* Computes the solar longitude for the given ephemeris time.
		*
		* @param imagePoint The line and sample
		* @param sensor The sensor model used to calculate ground coordinates
		* @param illuminator The illumination model
		* @param body The target body
		*
		* @returns solar longitude in degrees.
		*/
		double solarLongitude(ImagePt &imagePoint, Sensor sensor, Illuminator illuminator, Body *body);
		}
		#endif

SensorUtilities/src/MathUtils.cpp

+154 −1

Original line number	Diff line number	Diff line
		@@ -79,6 +79,18 @@ namespace SensorUtilities {
		}


		GroundPt2D radiansToDegrees(GroundPt2D radianLatLon) {
		double degreeLon = radianLatLon.lon;
		if (degreeLon < 0) {
		degreeLon += 2 * M_PI;
		}

		degreeLon *= RAD2DEG;
		double degreeLat = radianLatLon.lat * RAD2DEG;
		return {degreeLat, degreeLon};
		}


		Vec sphericalToRect(GroundPt3D spherical) {
		return {
		spherical.radius * cos(spherical.lat) * cos(spherical.lon),
		@@ -99,5 +111,146 @@ namespace SensorUtilities {
		};
		}

		double groundAzimuth(GroundPt2D groundPt, GroundPt2D subPt) {
		double a;
		double b;

		if (groundPt.lat >= 0.0) {
		a = (90.0 - subPt.lat) * M_PI / 180.0;
		b = (90.0 - groundPt.lat) * M_PI / 180.0;
		}
		else {
		a = (90.0 + subPt.lat) * M_PI / 180.0;
		b = (90.0 + groundPt.lat) * M_PI / 180.0;
		}

		GroundPt2D cs;
		GroundPt2D cg;

		cs.lon = subPt.lon;
		cg.lon = groundPt.lon;

		if (cs.lon > cg.lon) {
		if ((cs.lon - cg.lon) > 180.0) {
		while ((cs.lon - cg.lon) > 180.0) cs.lon = cs.lon - 360.0;
		}
		}
		if (cg.lon > cs.lon) {
		if ((cg.lon-cs.lon) > 180.0) {
		while ((cg.lon-cs.lon) > 180.0) cg.lon = cg.lon - 360.0;
		}
		}

		int quad;
		if (subPt.lat > groundPt.lat) {
		if (cs.lon < cg.lon) {
		quad = 2;
		}
		else {
		quad = 1;
		}
		}
		else if (subPt.lat < groundPt.lat) {
		if (cs.lon < cg.lon) {
		quad = 3;
		}
		else {
		quad = 4;
		}
		}
		else {
		if (cs.lon > cg.lon) {
		quad = 1;
		}
		else if (cs.lon < cg.lon) {
		quad = 2;
		}
		else {
		return 0.0;
		}
		}

		double C = (cg.lon - cs.lon) * M_PI / 180.0;
		if (C < 0) C = -C;

		double c = acos( cos(a) * cos(b) + sin(a) * sin(b) * cos(C));

		double azimuth = 0.0;

		if (sin(b) == 0.0 \|\| sin(c) == 0.0) {
		return azimuth;
		}

		double intermediate = (cos(a) - cos(b) * cos(c)) / (sin(b) * sin(c));
		if (intermediate < -1.0) {
		intermediate = -1.0;
		}
		else if (intermediate > 1.0) {
		intermediate = 1.0;
		}

		double A = acos(intermediate) * 180.0 / M_PI;

		if (groundPt.lat >= 0.0) {
		if (quad == 1 \|\| quad == 4) {
		azimuth = A;
		}
		else if (quad == 2 \|\| quad == 3) {
		azimuth = 360.0 - A;
		}
		}
		else {
		if (quad == 1 \|\| quad == 4) {
		azimuth = 180.0 - A;
		}
		else if (quad == 2 \|\| quad == 3) {
		azimuth = 180.0 + A;
		}
		}
		return azimuth;
		}

		Vec crossProduct(Vec aVec, Vec bVec) {
		double x = aVec.y * bVec.z - aVec.z * bVec.y;
		double y = aVec.z * bVec.x - aVec.x * bVec.z;
		double z = aVec.x * bVec.y - aVec.y * bVec.x;
		return {x, y, z};
		}

		Vec unitVector(Vec vec) {
		double mag = magnitude(vec);
		return {vec.x / mag, vec.y / mag, vec.z / mag};
		}

		Vec scaleVector(Vec vec, double scalar) {
		return {vec.x * scalar, vec.y * scalar, vec.z * scalar};
		}

		Vec perpendicularVec(Vec aVec, Vec bVec) {
		if (magnitude(aVec) == 0){
		return bVec;
		}

		Vec aNorm = unitVector(aVec);
		Vec bNorm = unitVector(bVec);

		double angle = bNorm.x * aNorm.x + bNorm.y * aNorm.y + bNorm.z * aNorm.z;
		double aMag = magnitude(aVec);
		double magP = angle * aMag;

		Vec p = {bNorm.x * magP, bNorm.y * magP, bNorm.z * magP};
		Vec q = aVec - p;

		return q;
		}

		Vec matrixVecProduct(Matrix mat, Vec vec) {
		Vec result;

		result.x = mat.a.x * vec.x + mat.a.y * vec.y + mat.a.z * vec.z;
		result.y = mat.b.x * vec.x + mat.b.y * vec.y + mat.b.z * vec.z;
		result.z = mat.c.x * vec.x + mat.c.y * vec.y + mat.c.z * vec.z;

		return result;
		}
		}
		No newline at end of file

SensorUtilities/src/SensorUtilities.cpp

+119 −1

Original line number	Diff line number	Diff line
		#include "SensorUtilities.h"
		#include <cmath>
		#include <iostream>


		namespace SensorUtilities {
		double phaseAngle(const ImagePt &imagePoint, Sensor sensor, Shape shape, Illuminator *illuminator) {
		@@ -27,6 +30,20 @@ namespace SensorUtilities {
		}


		double slantDistance(ImagePt imagePoint, Sensor sensor, Shape shape) {
		ObserverState sensorState = sensor->getState(imagePoint);
		Intersection intersect = shape->intersect(sensorState.sensorPos, sensorState.lookVec);

		return distance(sensorState.sensorPos, intersect.groundPt);
		}


		double targetCenterDistance(ImagePt &imagePoint, Sensor sensor, Body body) {
		ObserverState sensorState = sensor->getState(imagePoint);
		Vec sB = body->fixedVector(sensorState.sensorPos);
		return distance(sB, {0,0,0});
		}

		double illuminationDistance(const ImagePt &imagePoint, Sensor sensor, Shape shape, Illuminator *illuminator) {
		ObserverState sensorState = sensor->getState(imagePoint);
		Intersection intersect = shape->intersect(sensorState.sensorPos, sensorState.lookVec);
		@@ -79,4 +96,105 @@ namespace SensorUtilities {
		Intersection intersect = shape->intersect(position, lookVec);
		return magnitude(intersect.groundPt);
		}


		RaDec rightAscensionDeclination(ImagePt &imagePoint, Sensor *sensor) {
		ObserverState sensorState = sensor->getState(imagePoint);
		GroundPt3D sphericalPt = rectToSpherical(sensorState.j2000LookVec);

		GroundPt2D sphericalPtRadians = {sphericalPt.lat, sphericalPt.lon};
		GroundPt2D sphericalPtDegrees = radiansToDegrees(sphericalPtRadians);

		return {sphericalPtDegrees.lon, sphericalPtDegrees.lat};
		}


		double localSolarTime(ImagePt &imagePoint, Sensor sensor, Shape shape, Illuminator *illuminator) {
		ObserverState sensorState = sensor->getState(imagePoint);

		GroundPt2D subSolarPt = subSolarPoint(imagePoint, sensor, illuminator);
		GroundPt2D subSolarPtDegrees = radiansToDegrees(subSolarPt);

		Intersection intersect = shape->intersect(sensorState.sensorPos, sensorState.lookVec);
		GroundPt3D sphericalPt = rectToSpherical(intersect.groundPt);
		GroundPt2D sphericalPtRadians = {sphericalPt.lat, sphericalPt.lon};
		GroundPt2D sphericalPtDegrees = radiansToDegrees(sphericalPtRadians);

		double lst = sphericalPtDegrees.lon - subSolarPtDegrees.lon + 180.0;
		lst = lst / 15.0; // 15 degrees per hour
		if (lst < 0.0) lst += 24.0;
		if (lst > 24.0) lst -= 24.0;
		return lst;
		}

		double lineResolution(const ImagePt &imagePoint, Sensor sensor, Shape shape, double focalLength, double pixelPitch, double lineScaleFactor) {
		ObserverState sensorState = sensor->getState(imagePoint);
		Intersection intersect = shape->intersect(sensorState.sensorPos, sensorState.lookVec);

		double dist = distance(sensorState.sensorPos, intersect.groundPt) * 1000.0;

		return (dist / (focalLength / pixelPitch)) * lineScaleFactor;
		}


		double sampleResolution(const ImagePt &imagePoint, Sensor sensor, Shape shape, double focalLength, double pixelPitch, double sampleScaleFactor) {
		ObserverState sensorState = sensor->getState(imagePoint);
		Intersection intersect = shape->intersect(sensorState.sensorPos, sensorState.lookVec);

		double dist = distance(sensorState.sensorPos, intersect.groundPt) * 1000.0;

		return (dist / (focalLength / pixelPitch)) * sampleScaleFactor;
		}


		double pixelResolution(const ImagePt &imagePoint, Sensor sensor, Shape shape, double focalLength, double pixelPitch, double lineScaleFactor, double sampleScaleFactor) {
		double lineRes = lineResolution(imagePoint, sensor, shape, focalLength, pixelPitch, lineScaleFactor);
		double sampRes = sampleResolution(imagePoint, sensor, shape, focalLength, pixelPitch, sampleScaleFactor);
		if (lineRes < 0.0) return 0.0;
		if (sampRes < 0.0) return 0.0;
		return (lineRes + sampRes) / 2.0;
		}


		double solarLongitude(ImagePt &imagePoint, Sensor sensor, Illuminator illuminator, Body *body) {
		ObserverState sensorState = sensor->getState(imagePoint);

		Vec illumPos = illuminator->position(sensorState.time);
		Vec illumVel = illuminator->velocity(sensorState.time);

		std::vector<double> bodyRot = body->rotation(sensorState.time);

		Vec sunAv = unitVector(crossProduct(illumPos, illumVel));
		double npole[3];
		for (int i = 0; i < 3; i++) {
		npole[i] = bodyRot[6+i];
		}

		Vec z = sunAv;
		Vec x = unitVector(crossProduct(npole, z));
		Vec y = unitVector(crossProduct(z, x));

		Matrix trans;
		trans.a.x = x.x;
		trans.b.x = y.x;
		trans.c.x = z.x;

		trans.a.y = x.y;
		trans.b.y = y.y;
		trans.c.y = z.y;

		trans.a.z = x.z;
		trans.b.z = y.z;
		trans.c.z = z.z;

		Vec pos = matrixVecProduct(trans, illumPos);
		GroundPt3D spherical = rectToSpherical(pos);

		double longitude360 = spherical.lon * RAD2DEG;

		if (longitude360 != 360.0) {
		longitude360 -= 360 * floor(longitude360 / 360);
		}
		return longitude360;
		}
		}