Fixed a bug in the ProjectionFactory class which causes maps to have an... (3c4e1655) · Commits · aflab / astrogeology / ISIS3

isis/src/base/objs/ProjectionFactory/ProjectionFactory.cpp

+78 −24

Original line number	Diff line number	Diff line
		@@ -22,8 +22,10 @@

		#include "ProjectionFactory.h"

		#include <cmath>

		#include <cfloat>
		#include <cmath>
		#include <iomanip>

		#include "Camera.h"
		#include "Cube.h"
		@@ -109,9 +111,6 @@ namespace Isis {
		}





		/**
		* This method returns a pointer to a RingPlaneProjection object. The projection is
		* intialized using information contained in a Label object. The information
		@@ -178,6 +177,7 @@ namespace Isis {
		}
		}


		/**
		* This method creates a map projection for a cube given a label. The label
		* must contain all the proper mapping information (radii, projection name,
		@@ -269,46 +269,99 @@ namespace Isis {
		throw IException(IException::Unknown, msg, _FILEINFO_);
		}


		// Convert upperleft coordinate to units of pixel
		// Truncate it to the nearest whole pixel (floor/ceil)
		// Convert it back to meters. But don't do this if
		// the X/Y position is already close to a whole pixel because
		// the floor/ceil function could cause an extra pixel to be added
		// just due to machine precision issues
		if (fabs(fmod(minX, pixelResolution)) > 1.0e-6) {
		if (pixelResolution - fabs(fmod(minX, pixelResolution)) > 1.0e-6) {
		double sampleOffset = floor(minX / pixelResolution);
		minX = sampleOffset * pixelResolution;

		bool flipX = false;
		bool flipY = false;

		double minXFlipped = minX;
		double maxXFlipped = maxX;
		double minYFlipped = minY;
		double maxYFlipped = maxY;

		//New range is (-1)*[0,maxX] = [-maxX, 0]

		if (minX == 0) {
		minXFlipped = -maxX;
		maxXFlipped = 0;
		flipX = true;

		}

		if (fabs(fmod(minXFlipped, pixelResolution)) > 1.0e-6) {
		if (pixelResolution - fabs(fmod(minXFlipped, pixelResolution)) > 1.0e-6) {
		double sampleOffset = floor(minXFlipped / pixelResolution);

		minXFlipped = sampleOffset * pixelResolution;

		}
		}


		// make sure that the distance from minX to maxX is at least one pixel wide
		// so we have at least one sample in the created cube
		if (maxX < minX + pixelResolution) {
		maxX = minX + pixelResolution;


		if (maxXFlipped < minXFlipped + pixelResolution) {
		maxXFlipped = minXFlipped + pixelResolution;
		}

		if (fabs(fmod(maxY, pixelResolution)) > 1.0e-6) {
		if (pixelResolution - fabs(fmod(maxY, pixelResolution)) > 1.0e-6) {
		double lineOffset = -1.0 * ceil(maxY / pixelResolution);
		maxY = -1.0 * lineOffset * pixelResolution;

		//New range is (-1)*[minY,0] = [0,-minY]
		if (maxY == 0) {
		maxYFlipped = -minY;
		minYFlipped = 0;
		flipY = true;

		}

		if (fabs(fmod(maxYFlipped, pixelResolution)) > 1.0e-6) {
		if (abs(pixelResolution - fabs(fmod(maxYFlipped, pixelResolution))) > 1.0e-6) {
		double lineOffset = ceil(maxYFlipped / pixelResolution);
		maxYFlipped = lineOffset * pixelResolution;

		}
		}


		// make sure that the distance from minY to maxY is at least one pixel wide
		// so we have at least one line in the created cube
		if (minY > maxY - pixelResolution) {
		minY = maxY - pixelResolution;
		if (minYFlipped > maxYFlipped - pixelResolution) {
		minYFlipped = maxYFlipped - pixelResolution;
		}

		// Determine the number of samples and lines
		samples = (int)((maxX - minX) / pixelResolution + 0.5);
		lines = (int)((maxY - minY) / pixelResolution + 0.5);

		samples = (int)((maxXFlipped - minXFlipped) / pixelResolution + 0.5);
		lines = (int)((maxYFlipped - minYFlipped) / pixelResolution + 0.5);


		// Set the upper left corner and add to the labels
		upperLeftX = minX;
		if (flipX) {
		upperLeftX = 0;

		}
		else {
		upperLeftX = minXFlipped;
		}
		mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerX", toString(upperLeftX)),
		Isis::Pvl::Replace);

		upperLeftY = maxY;
		if (flipY) {
		upperLeftY = 0;

		}
		else {
		upperLeftY = maxYFlipped;

		}

		mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerY", toString(upperLeftY)),
		Isis::Pvl::Replace);

		@@ -367,6 +420,7 @@ namespace Isis {
		return (Isis::Projection *) proj;
		}


		/**
		* This method creates a projection for a cube to a ring plane given a label.
		* Currently this is utilized only for projecting images of rings to the ring
		@@ -388,7 +442,7 @@ namespace Isis {
		* the size of the output cube to the size of the input cube.
		* Defaults to true.
		*
		* @return A pointer to a Projection object.
		* @return @b Projection* A pointer to a Projection object.
		*
		*/
		Isis::Projection *ProjectionFactory::RingsCreateForCube(Isis::Pvl &label,
		@@ -577,7 +631,7 @@ namespace Isis {
		* and returned.
		* @param cam An initialized camera model
		*
		* @return A pointer to a Projection object.
		* @return @b Projection* A pointer to a Projection object.
		*
		*/
		Isis::Projection *ProjectionFactory::CreateForCube(Isis::Pvl &label,
		@@ -828,7 +882,7 @@ namespace Isis {
		* and returned.
		* @param cam An initialized camera model
		*
		* @return A pointer to a Projection object.
		* @return Projection* A pointer to a Projection object.
		*
		*/
		Isis::Projection *ProjectionFactory::RingsCreateForCube(Isis::Pvl &label,

isis/src/base/objs/ProjectionFactory/ProjectionFactory.h

+3 −0

Original line number	Diff line number	Diff line
		@@ -84,6 +84,9 @@ namespace Isis {
		* at least one line. References #775.
		* @history 2014-01-16 Kimberly Oyama - Updated the error messages for unsupported projections
		* to include the source of the error. Fixes #988.
		* @history 2016-09-02 Tyler Wilson - Fixed a bug in the CreateForCube function because it was
		* producing an unequal number of lines/samples for maps from projections that are symmetric about
		* the prime meridian or the equator. Fixes #2245.
		*/
		class ProjectionFactory {
		public:

isis/src/base/objs/ProjectionFactory/ProjectionFactory.truth

+12 −0

Original line number	Diff line number	Diff line
		@@ -73,3 +73,15 @@ Test for ProjectionFactory's Create with unsupported projection
		ERROR Unsupported projection, unable to find plugin for [UnsupportedProjection].
		ERROR Unable to find PVL group [UnsupportedProjection] in file [isis/lib/Projection.plugin].


		Test for ProjectionFactory's Create method when the maximum latitude = 0.0
		South Map Samples = 8420
		North Map Samples = 8420
		South Map Lines = 1404
		North Map Lines = 1404

		Test for ProjectionFactory's Create method when the maximum longitude = 0.0
		West Map Samples = 2105
		East Map Samples = 2105
		West Map Lines = 1404
		East Map Lines = 1404

isis/src/base/objs/ProjectionFactory/unitTest.cpp

+143 −0

Original line number	Diff line number	Diff line
		@@ -134,6 +134,149 @@ int main(int argc, char *argv[]) {
		catch(IException &e) {
		ReportError( e.toString() );
		}


		try {

		cout << endl << "Test for ProjectionFactory's Create method when the maximum latitude = 0.0"
		<< endl;

		Pvl southMap;
		int southLines,southSamples;
		southMap.addGroup(PvlGroup("Mapping"));
		PvlGroup &mapGroupSouth = southMap.findGroup("Mapping");
		mapGroupSouth += PvlKeyword("ProjectionName", "Equirectangular");
		mapGroupSouth += PvlKeyword("CenterLongitude", toString(0.0));
		mapGroupSouth += PvlKeyword("CenterLatitude", toString(0.0));
		mapGroupSouth += PvlKeyword("EquatorialRadius", toString(13400.0));
		mapGroupSouth += PvlKeyword("PolarRadius", toString(9200.0));
		mapGroupSouth += PvlKeyword("LatitudeType", "Planetocentric");
		mapGroupSouth += PvlKeyword("LongitudeDirection", "PositiveEast");
		mapGroupSouth += PvlKeyword("LongitudeDomain", toString(360));
		mapGroupSouth += PvlKeyword("MinimumLatitude", toString(-60.0));
		mapGroupSouth += PvlKeyword("MaximumLatitude", toString(0.0));
		mapGroupSouth += PvlKeyword("MinimumLongitude", toString(0.0));
		mapGroupSouth += PvlKeyword("MaximumLongitude", toString(360.0));
		mapGroupSouth += PvlKeyword("PixelResolution", toString(10.0));


		Pvl northMap;
		int northLines,northSamples;
		northMap.addGroup(PvlGroup("Mapping"));
		PvlGroup &mapGroupNorth = northMap.findGroup("Mapping");
		mapGroupNorth += PvlKeyword("ProjectionName", "Equirectangular");
		mapGroupNorth += PvlKeyword("CenterLongitude", toString(0.0));
		mapGroupNorth += PvlKeyword("CenterLatitude", toString(0.0));
		mapGroupNorth += PvlKeyword("EquatorialRadius", toString(13400.0));
		mapGroupNorth += PvlKeyword("PolarRadius", toString(9200.0));
		mapGroupNorth += PvlKeyword("LatitudeType", "Planetocentric");
		mapGroupNorth += PvlKeyword("LongitudeDirection", "PositiveEast");
		mapGroupNorth += PvlKeyword("LongitudeDomain", toString(360));
		mapGroupNorth += PvlKeyword("MinimumLatitude", toString(0.0));
		mapGroupNorth += PvlKeyword("MaximumLatitude", toString(60.0));
		mapGroupNorth += PvlKeyword("MinimumLongitude", toString(0.0));
		mapGroupNorth += PvlKeyword("MaximumLongitude", toString(360.0));
		mapGroupNorth += PvlKeyword("PixelResolution", toString(10.0));


		TProjection projNorth = (TProjection ) ProjectionFactory::CreateForCube(southMap,
		southSamples, southLines);
		TProjection projSouth = (TProjection ) ProjectionFactory::CreateForCube(northMap,
		northSamples, northLines);

		projNorth->SetWorld(245.0, 355.0);
		projSouth->SetWorld(245.0, 355.0);



		cout << "South Map Samples = " << southSamples << endl;
		cout << "North Map Samples = " << northSamples << endl;
		cout << "South Map Lines = " << southLines << endl;
		cout << "North Map Lines = " << northLines << endl;




		}

		catch (IException &e) {

		ReportError (e.toString() );

		}


		try {

		cout << endl << "Test for ProjectionFactory's Create method when the maximum longitude = 0.0"
		<< endl;

		Pvl eastMap;
		int eastLines,eastSamples;
		eastMap.addGroup(PvlGroup("Mapping"));
		PvlGroup &mapGroupEast = eastMap.findGroup("Mapping");
		mapGroupEast += PvlKeyword("ProjectionName", "Equirectangular");
		mapGroupEast += PvlKeyword("CenterLongitude", toString(0.0));
		mapGroupEast += PvlKeyword("CenterLatitude", toString(0.0));
		mapGroupEast += PvlKeyword("EquatorialRadius", toString(13400.0));
		mapGroupEast += PvlKeyword("PolarRadius", toString(9200.0));
		mapGroupEast += PvlKeyword("LatitudeType", "Planetocentric");
		mapGroupEast += PvlKeyword("LongitudeDirection", "PositiveEast");
		mapGroupEast += PvlKeyword("LongitudeDomain", toString(360));
		mapGroupEast += PvlKeyword("MinimumLatitude", toString(0.0));
		mapGroupEast += PvlKeyword("MaximumLatitude", toString(60.0));
		mapGroupEast += PvlKeyword("MinimumLongitude", toString(0.0));
		mapGroupEast += PvlKeyword("MaximumLongitude", toString(90.0));
		mapGroupEast += PvlKeyword("PixelResolution", toString(10.0));


		Pvl westMap;
		int westLines,westSamples;
		westMap.addGroup(PvlGroup("Mapping"));
		PvlGroup &mapGroupWest = westMap.findGroup("Mapping");
		mapGroupWest += PvlKeyword("ProjectionName", "Equirectangular");
		mapGroupWest += PvlKeyword("CenterLongitude", toString(0.0));
		mapGroupWest += PvlKeyword("CenterLatitude", toString(0.0));
		mapGroupWest += PvlKeyword("EquatorialRadius", toString(13400.0));
		mapGroupWest += PvlKeyword("PolarRadius", toString(9200.0));
		mapGroupWest += PvlKeyword("LatitudeType", "Planetocentric");
		mapGroupWest += PvlKeyword("LongitudeDirection", "PositiveEast");
		mapGroupWest += PvlKeyword("LongitudeDomain", toString(360));
		mapGroupWest += PvlKeyword("MinimumLatitude", toString(0.0));
		mapGroupWest += PvlKeyword("MaximumLatitude", toString(60.0));
		mapGroupWest += PvlKeyword("MinimumLongitude", toString(-90.0));
		mapGroupWest += PvlKeyword("MaximumLongitude", toString(0.0));
		mapGroupWest += PvlKeyword("PixelResolution", toString(10.0));


		TProjection projEast = (TProjection ) ProjectionFactory::CreateForCube(eastMap,
		eastSamples, eastLines);
		TProjection projWest = (TProjection ) ProjectionFactory::CreateForCube(westMap,
		westSamples, westLines);

		projEast->SetWorld(245.0, 355.0);
		projWest->SetWorld(245.0, 355.0);



		cout << "West Map Samples = " << westSamples << endl;
		cout << "East Map Samples = " << eastSamples << endl;
		cout << "West Map Lines = " << westLines << endl;
		cout << "East Map Lines = " << eastLines << endl;




		}

		catch (IException &e) {

		ReportError (e.toString() );

		}



		}

		void doit(Pvl &lab) {