C++ ColorMap::writeFits方法代码示例

ColorMap* getAggregatedARMap(const ColorMap* map, Real cleaningFactor, Real aggregationFactor, const string& projection)
{
	// We convert the factors from arcsec to pixels
	cleaningFactor /= sqrt(map->PixelArea());
	aggregationFactor /= sqrt(map->PixelArea());
	
	ColorMap* aggregated = new ColorMap(map);

	#if defined DEBUG
	aggregated->writeFits(filenamePrefix + "pure.fits");
	#endif
	
	if(projection == "exact")
	{
		/*! Clean the color map to remove very small components (like protons)*/
		aggregated->erodeCircularProjected(cleaningFactor, 0);
	
		#if defined DEBUG
		aggregated->writeFits(filenamePrefix + "eroded.fits");
		#endif
	
		/*! Aggregate the blobs together */
		aggregated->dilateCircularProjected(cleaningFactor + aggregationFactor, 0);
	
		#if defined DEBUG
		aggregated->writeFits(filenamePrefix + "dilated.fits");
		#endif
	
		/*! Give back the original size */
		aggregated->erodeCircularProjected(aggregationFactor, 0);
	
		#if defined DEBUG
		aggregated->writeFits(filenamePrefix + "closed.fits");
		#endif
	}
	else
	{
		/*! Apply the projection */
		if(projection == "equirectangular")
		{
			aggregated->equirectangular_projection(map, false);
			#if defined DEBUG
			aggregated->writeFits(filenamePrefix + "equirectangular_projection.fits");
			#endif
			// We adjust the factors because in the projection the pixel size changes
			cleaningFactor *= (2./3.);
			aggregationFactor *= (2./3.);
		}
		else if(projection == "lambert")
		{
			aggregated->Lambert_cylindrical_projection(map, false);
			#if defined DEBUG
			aggregated->writeFits(filenamePrefix + "Lambert_cylindrical_projection.fits");
			#endif
			// We adjust the factors because in the projection the pixel size changes
			cleaningFactor *= (2./3.);
			aggregationFactor *= (2./3.);
		}
		else if(projection == "sinusoidal")
		{
			aggregated->sinusoidal_projection(map, false);
			#if defined DEBUG
			aggregated->writeFits(filenamePrefix + "sinusoidal_projection.fits");
			#endif
			// We adjust the factors because in the projection the pixel size changes
			cleaningFactor *= (2./3.);
			aggregationFactor *= (2./3.);
		}
		else if(projection != "none")
		{
			cerr<<"Unknown projection type "<<projection<<endl;
			exit(EXIT_FAILURE);
		}
	
		/*! Clean the color map to remove very small components (like protons)*/
		aggregated->erodeCircular(cleaningFactor, 0);
	
		#if defined DEBUG
		aggregated->writeFits(filenamePrefix + "eroded.fits");
		#endif
	
		/*! Aggregate the blobs together */
		aggregated->dilateCircular(cleaningFactor + aggregationFactor, 0);
	
		#if defined DEBUG
		aggregated->writeFits(filenamePrefix + "dilated.fits");
		#endif
	
		/*! Give back the original size */
		aggregated->erodeCircular(aggregationFactor, 0);
	
		#if defined DEBUG
		aggregated->writeFits(filenamePrefix + "closed.fits");
		#endif
	
		/* Apply the deprojection */
		if(projection == "equirectangular")
		{
			ColorMap* projeted = new ColorMap(aggregated);
			aggregated->equirectangular_deprojection(projeted, false);
//.........这里部分代码省略.........

int main(int argc, const char **argv)
{
	// We declare our program description
	string programDescription = "This Program does attribution and segmentation.";
	programDescription+="\nVersion: 3.0";
	programDescription+="\nAuthor: Benjamin Mampaey, [email protected]";
	
	programDescription+="\nCompiled on "  __DATE__  " with options :";
	programDescription+="\nNUMBERCHANNELS: " + toString(NUMBERCHANNELS);
	#if defined DEBUG
	programDescription+="\nDEBUG: ON";
	#endif
	#if defined EXTRA_SAFE
	programDescription+="\nEXTRA_SAFE: ON";
	#endif
	#if defined VERBOSE
	programDescription+="\nVERBOSE: ON";
	#endif
	programDescription+="\nEUVPixelType: " + string(typeid(EUVPixelType).name());
	programDescription+="\nReal: " + string(typeid(Real).name());
	
	// We define our program parameters
	ArgParser args(programDescription);
	
	args("segmentation") = Classifier::segmentationParameters();
	args("classification") = Classifier::classificationParameters();
	
	args["config"] = ArgParser::ConfigurationFile('C');
	args["help"] = ArgParser::Help('h');
	
	args["type"] = ArgParser::Parameter("SPoCA2", 'T', "The type of classifier to use for the attribution.\nPossible values: FCM, PFCM, PCM, PCM2, SPoCA, SPoCA2");
	args["imageType"] = ArgParser::Parameter("Unknown", 'I', "The type of the images.\nPossible values: EIT, EUVI, AIA, SWAP");
	args["imagePreprocessing"] = ArgParser::Parameter("ALC", 'P', "The steps of preprocessing to apply to the sun images.\nCan be any combination of the following:\n NAR=zz.z (Nullify pixels above zz.z*radius)\n ALC (Annulus Limb Correction)\n DivMedian (Division by the median)\n TakeSqrt (Take the square root)\n TakeLog (Take the log)\n TakeAbs (Take the absolute value)\n DivMode (Division by the mode)\n DivExpTime (Division by the Exposure Time)\n ThrMin=zz.z (Threshold intensities to minimum zz.z)\n ThrMax=zz.z (Threshold intensities to maximum zz.z)\n ThrMinPer=zz.z (Threshold intensities to minimum the zz.z percentile)\n ThrMaxPer=zz.z (Threshold intensities to maximum the zz.z percentile)\n Smooth=zz.z (Binomial smoothing of zz.z arcsec)");
	args["registerImages"] = ArgParser::Parameter(false, 'r', "Set to register/align the images when running multi channel attribution.");
	args["centersFile"] = ArgParser::Parameter("centers.txt", 'c', "The name of the file containing the centers. If it it not provided the centers will be initialized randomly.");
	args["computeEta"] = ArgParser::Parameter(false, 'e', "If the enters file do not contain the values for Eta or if you want to force Eta to be recomputed (slow!).");
	args["stats"] = ArgParser::Parameter(false, 's', "Set to compute stats about the generated maps.");
	args["statsPreprocessing"] = ArgParser::Parameter("NAR=0.95", 'P', "The steps of preprocessing to apply to the sun images.\nCan be any combination of the following:\n NAR=zz.z (Nullify pixels above zz.z*radius)\n ALC (Annulus Limb Correction)\n DivMedian (Division by the median)\n TakeSqrt (Take the square root)\n TakeLog (Take the log)\n TakeAbs (Take the absolute value)\n DivMode (Division by the mode)\n DivExpTime (Division by the Exposure Time)\n ThrMin=zz.z (Threshold intensities to minimum zz.z)\n ThrMax=zz.z (Threshold intensities to maximum zz.z)\n ThrMinPer=zz.z (Threshold intensities to minimum the zz.z percentile)\n ThrMaxPer=zz.z (Threshold intensities to maximum the zz.z percentile)\n Smooth=zz.z (Binomial smoothing of zz.z arcsec)");
	args["output"] = ArgParser::Parameter(".", 'O', "The name for the output file or of a directory.");
	args["uncompressed"] = ArgParser::Parameter(false, 'u', "Set this flag if you want results maps to be uncompressed.");
	args["fitsFile"] = ArgParser::RemainingPositionalParameters("Path to a fits file", NUMBERCHANNELS, NUMBERCHANNELS);
	
	// We parse the arguments
	try
	{
		args.parse(argc, argv);
	}
	catch ( const invalid_argument& error)
	{
		cerr<<"Error : "<<error.what()<<endl;
		cerr<<args.help_message(argv[0])<<endl;
		return EXIT_FAILURE;
	}
	
	// We setup the output directory
	string outputDirectory;
	string outputFile = args["output"];
	if (isDir(outputFile))
	{
		outputDirectory = outputFile;
	}
	else
	{
		outputDirectory = getPath(outputFile);
		if (! isDir(outputDirectory))
		{
			cerr<<"Error : "<<outputDirectory<<" is not a directory!"<<endl;
			return EXIT_FAILURE;
		}
	}
	
	// We read and preprocess the sun images
	deque<string> imagesFilenames = args.RemainingPositionalArguments();
	vector<EUVImage*> images;
	for (unsigned p = 0; p < imagesFilenames.size(); ++p)
	{
		EUVImage* image = getImageFromFile(args["imageType"], imagesFilenames[p]);
		image->preprocessing(args["imagePreprocessing"]);
		
		#if defined DEBUG
			image->getHeader().set("IPREPROC", args["imagePreprocessing"], "Image Preprocessing");
			image->writeFits(outputDirectory + "/" + stripPath(stripSuffix(imagesFilenames[p])) + ".preprocessed.fits");
		#endif
		images.push_back(image);
	}
	
	// We verify the images are aligned and we register them
	for(unsigned p = 1; p < images.size(); ++p)
	{
		string dissimilarity = checkSimilar(images[0], images[p]);
		if(! dissimilarity.empty())
		{
			if(args["registerImages"])
			{
				#if defined VERBOSE
				cout<<"Image "<<imagesFilenames[p]<<" will be registered to image "<<imagesFilenames[0]<<endl;
				#endif
				images[p]->align(images[0]);
				#if defined DEBUG
				images[p]->writeFits(outputDirectory + "/" + stripPath(stripSuffix(imagesFilenames[p])) + ".registered.fits");
//.........这里部分代码省略.........