[casacore] r21231 committed - Removed files that are now part of CASA's imageanalysis package

[casa...@googlecode.com]

Revision: 21231
Author: gervandiepen
Date: Mon Apr 2 05:08:24 2012
Log: Removed files that are now part of CASA's imageanalysis package

http://code.google.com/p/casacore/source/detail?r=21231

Deleted:
/trunk/images/IO
/trunk/images/Images/ImageAnalysis.cc
/trunk/images/Images/ImageAnalysis.h
/trunk/images/Images/ImageCollapser.cc
/trunk/images/Images/ImageCollapser.h
/trunk/images/Images/ImageFitter.cc
/trunk/images/Images/ImageFitter.h
/trunk/images/Images/ImageInputProcessor.cc
/trunk/images/Images/ImageInputProcessor.h
/trunk/images/Images/ImageProfileFit.cc
/trunk/images/Images/ImageProfileFit.h
/trunk/images/Images/ImageProfileFitter.cc
/trunk/images/Images/ImageProfileFitter.h
/trunk/images/Images/test/dImageProfileFit.cc
/trunk/images/Images/test/dImageProfileFit.run
/trunk/images/Images/test/tImageCollapser.cc
/trunk/images/Images/test/tImageFitter.cc
/trunk/images/Images/test/tImageInputProcessor.cc
/trunk/images/Images/test/tImageProfileFitter.cc
Modified:
/trunk/images/CMakeLists.txt
/trunk/images/Images/ImageUtilities2.tcc
/trunk/images/Images/test/CMakeLists.txt
/trunk/images/apps/CMakeLists.txt

=======================================
--- /trunk/images/Images/ImageAnalysis.cc Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,7169 +0,0 @@
-//# ImageAnalysis.cc: Image analysis and handling tool
-//# Copyright (C) 1995-2007
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This library is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU Library General Public License as published
by
-//# the Free Software Foundation; either version 2 of the License, or (at
your
-//# option) any later version.
-//#
-//# This library is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
-//# License for more details.
-//#
-//# You should have received a copy of the GNU Library General Public
License
-//# along with this library; if not, write to the Free Software Foundation,
-//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-//
-#include <casa/aips.h>
-#include <casa/iostream.h>
-#include <casa/sstream.h>
-#include <casa/Arrays/ArrayIO.h>
-#include <casa/Arrays/ArrayMath.h>
-#include <casa/Arrays/ArrayUtil.h>
-#include <casa/Arrays/MaskedArray.h>
-#include <casa/Arrays/MaskArrMath.h>
-#include <casa/BasicMath/Random.h>
-#include <casa/BasicSL/String.h>
-#include <casa/Containers/Record.h>
-#include <casa/Exceptions/Error.h>
-#include <casa/fstream.h>
-#include <casa/Logging/LogFilter.h>
-#include <casa/Logging/LogIO.h>
-#include <casa/Logging/LogOrigin.h>
-#include <casa/OS/Directory.h>
-#include <casa/OS/EnvVar.h>
-#include <casa/OS/File.h>
-#include <casa/OS/HostInfo.h>
-#include <casa/OS/RegularFile.h>
-#include <casa/OS/SymLink.h>
-#include <casa/Quanta/QuantumHolder.h>
-#include <casa/Quanta/Quantum.h>
-#include <casa/Utilities/Assert.h>
-#include <casa/Utilities/Regex.h>
-#include <measures/Measures/MeasureHolder.h>
-#include <measures/Measures/MDirection.h>
-#include <measures/Measures/MEpoch.h>
-#include <measures/Measures/MDoppler.h>
-#include <measures/Measures/MFrequency.h>
-#include <measures/Measures/MPosition.h>
-#include <components/ComponentModels/ComponentList.h>
-#include <components/ComponentModels/ComponentShape.h>
-#include <components/ComponentModels/GaussianShape.h>
-#include <components/ComponentModels/SkyComponent.h>
-#include <components/ComponentModels/SkyCompRep.h>
-#include <components/ComponentModels/TwoSidedShape.h>
-#include <components/SpectralComponents/SpectralElement.h>
-#include <coordinates/Coordinates/CoordinateSystem.h>
-#include <coordinates/Coordinates/StokesCoordinate.h>
-#include <coordinates/Coordinates/CoordinateUtil.h>
-#include <coordinates/Coordinates/DirectionCoordinate.h>
-#include <coordinates/Coordinates/SpectralCoordinate.h>
-#include <coordinates/Coordinates/GaussianConvert.h>
-#include <coordinates/Coordinates/LinearCoordinate.h>
-#include <images/Images/ComponentImager.h>
-#include <images/Images/Image2DConvolver.h>
-#include <images/Images/ImageCollapser.h>
-#include <images/Images/ImageConcat.h>
-#include <images/Images/ImageConvolver.h>
-#include <images/Images/ImageDecomposer.h>
-#include <images/Images/ImageExpr.h>
-#include <images/Images/ImageExprParse.h>
-#include <images/Images/ImageFFT.h>
-#include <images/Images/ImageFITSConverter.h>
-#include <images/Images/ImageHistograms.h>
-#include <images/Images/ImageMoments.h>
-#include <images/Images/ImageMetaData.h>
-#include <images/Images/ImageOpener.h>
-#include <images/Regions/ImageRegion.h>
-#include <images/Images/ImageRegrid.h>
-#include <images/Images/ImageSourceFinder.h>
-#include <images/Images/ImageStatistics.h>
-#include <images/Images/ImageSummary.h>
-#include <images/Images/ImageTwoPtCorr.h>
-#include <images/Images/ImageUtilities.h>
-#include <images/Images/LELImageCoord.h>
-#include <images/Images/PagedImage.h>
-#include <images/Images/RebinImage.h>
-#include <images/Regions/RegionManager.h>
-#include <images/Images/SepImageConvolver.h>
-#include <images/Images/SubImage.h>
-#include <images/Images/TempImage.h>
-#include <images/Regions/WCLELMask.h>
-#include <images/Images/ImageAnalysis.h>
-#include <lattices/LatticeMath/Fit2D.h>
-#include <lattices/LatticeMath/LatticeFit.h>
-#include <lattices/Lattices/LatticeAddNoise.h>
-#include <lattices/Lattices/LatticeExprNode.h>
-#include <lattices/Lattices/LatticeExprNode.h>
-#include <lattices/Lattices/LatticeIterator.h>
-#include <lattices/Lattices/LatticeRegion.h>
-#include <lattices/Lattices/LatticeSlice1D.h>
-#include <lattices/Lattices/LatticeUtilities.h>
-#include <lattices/Lattices/LCBox.h>
-#include <lattices/Lattices/LCSlicer.h>
-#include <lattices/Lattices/MaskedLatticeIterator.h>
-#include <lattices/Lattices/PixelCurve1D.h>
-#include <lattices/Lattices/RegionType.h>
-#include <lattices/Lattices/TiledLineStepper.h>
-#include <scimath/Fitting/LinearFitSVD.h>
-#include <scimath/Functionals/Polynomial.h>
-#include <scimath/Mathematics/VectorKernel.h>
-#include <tables/LogTables/NewFile.h>
-#include <images/Images/FITSImage.h>
-#include <images/Images/MIRIADImage.h>
-
-#include <casa/namespace.h>
-
-namespace casa { //# name space casa begins
-
-ImageAnalysis::ImageAnalysis() :
- pImage_p(0), pStatistics_p(0), pHistograms_p(0),
- pOldStatsRegionRegion_p(0), pOldStatsMaskRegion_p(0),
- pOldHistRegionRegion_p(0), pOldHistMaskRegion_p(0) {
-
- // Register the functions to create a FITSImage or MIRIADImage object.
- FITSImage::registerOpenFunction();
- MIRIADImage::registerOpenFunction();
-
- itsLog = new LogIO();
-
-}
-
-ImageAnalysis::ImageAnalysis(const ImageInterface<Float>* inImage) :
- pImage_p(inImage->cloneII()), itsLog(new LogIO()), pStatistics_p(0),
- pHistograms_p(0), pOldStatsRegionRegion_p(0),
- pOldStatsMaskRegion_p(0), pOldHistRegionRegion_p(0),
- pOldHistMaskRegion_p(0) {}
-
-ImageAnalysis::ImageAnalysis(ImageInterface<Float>* inImage, const Bool
cloneInputPointer) :
-itsLog(new LogIO()), pStatistics_p(0), pHistograms_p(0),
- pOldStatsRegionRegion_p(0),
- pOldStatsMaskRegion_p(0), pOldHistRegionRegion_p(0),
- pOldHistMaskRegion_p(0) {
- pImage_p = cloneInputPointer ? inImage->cloneII() : inImage;
-}
-
-ImageAnalysis::~ImageAnalysis() {
-
- if (pImage_p != 0) {
- if((pImage_p->isPersistent()) && ((pImage_p->imageType())
== "PagedImage")){
- ImageOpener::ImageTypes type =
ImageOpener::imageType(pImage_p->name());
- if (type == ImageOpener::AIPSPP) {
- (static_cast<PagedImage<Float>
*>(pImage_p))->table().relinquishAutoLocks(True);
- (static_cast<PagedImage<Float> *>(pImage_p))->table().unlock();
- }
- }
-
-
- delete pImage_p;
- pImage_p = 0;
- }
- deleteHistAndStats();
-}
-
-Bool ImageAnalysis::toRecord(RecordInterface& rec) {
-
- if (pImage_p != 0) {
- String err;
- return pImage_p->toRecord(err, rec);
-
- }
-
- return False;
-}
-
-Bool ImageAnalysis::fromRecord(const RecordInterface& rec, const String&
name) {
-
- Bool retval = False;
- String err;
- if (name != "") {
- TempImage<Float> tempim;
- retval = tempim.fromRecord(err, rec);
- if (retval) {
- {
- PagedImage<Float> diskim(tempim.shape(), tempim.coordinates(),
- name);
- diskim.copyData(tempim);
- // go out of context hence flush to disk
- }
- retval = open(name);
- }
- } else {
- if (itsLog == 0)
- itsLog = new LogIO();
- if (pImage_p != 0) {
- delete pImage_p;
- *itsLog << LogOrigin("ImageAnalysis", "fromRecord");
- *itsLog << LogIO::WARN
- << "Image is already open, disconnecting first"
- << LogIO::POST;
- pImage_p = 0;
- }
- pImage_p = new TempImage<Float> ();
- retval = pImage_p->fromRecord(err, rec);
-
- }
-
- return retval;
-
-}
-
-Bool ImageAnalysis::open(const String& infile) {
- Bool rstat = True;
-
- if (itsLog == 0)
- itsLog = new LogIO();
-
- *itsLog << LogOrigin("ImageAnalysis", "open");
-
- // Check whether infile exists
- if (infile.empty()) {
- *itsLog << LogIO::WARN << "File string is empty" << LogIO::POST;
- return false;
- }
- File thefile(infile);
- if (!thefile.exists()) {
- *itsLog << LogIO::WARN << "File [" << infile << "] does not exist."
- << LogIO::POST;
- return false;
- }
-
- // Generally used if the image is already closed !b
- if (pImage_p != 0) {
- delete pImage_p;
- *itsLog << LogIO::WARN << "Image is already open, closing first"
- << LogIO::POST;
- pImage_p = 0;
- }
-
- // Open input image. We don't handle an Image tool because
- // we would get a bit confused as to who owns the pointer
- ImageUtilities::openImage(pImage_p, infile, *itsLog);
-
- // Ensure that we reconstruct the statistics and histograms objects
- deleteHistAndStats();
- return rstat;
-}
-
-Bool ImageAnalysis::detached() {
- if (pImage_p == 0)
- return True;
- return False;
-
-}
-
-Bool ImageAnalysis::addnoise(const String& type, const Vector<Double>&
pars,
- Record& region, const Bool zeroIt) {
- bool rstat(False);
- *itsLog << LogOrigin("ImageAnalysis", "addnoise");
-
- // Make SubImage
- String mask;
- ImageRegion* pRegionRegion = 0;
- ImageRegion* pMaskRegion = 0;
- SubImage<Float> subImage = SubImage<Float>::createSubImage(
- pRegionRegion, pMaskRegion, *pImage_p,
- region,
- mask, itsLog, True
- );
- delete pRegionRegion;
- delete pMaskRegion;
-
- // Zero subimage if requested
- if (zeroIt)
- subImage.set(0.0);
-
- // Do it
- Random::Types typeNoise = Random::asType(type);
- LatticeAddNoise lan(typeNoise, pars);
- lan.add(subImage);
- //
- deleteHistAndStats();
- rstat = true;
-
- return rstat;
-}
-
-ImageInterface<Float> *
-ImageAnalysis::imagecalc(const String& outfile, const String& expr,
- const Bool overwrite) {
-
- *itsLog << LogOrigin("ImageAnalysis", "imagecalc");
-
- Record regions;
-
- // Check output file name
- if (!outfile.empty() && !overwrite) {
- NewFile validfile;
- String errmsg;
- if (!validfile.valueOK(outfile, errmsg)) {
- *itsLog << errmsg << LogIO::EXCEPTION;
- }
- }
-
- // Get LatticeExprNode (tree) from parser. Substitute possible
- // object-id's by a sequence number, while creating a
- // LatticeExprNode for it. Convert the GlishRecord containing
- // regions to a PtrBlock<const ImageRegion*>.
- if (expr.empty()) {
- *itsLog << "You must specify an expression" << LogIO::EXCEPTION;
- }
-
- Block<LatticeExprNode> temps;
- String exprName;
- // String newexpr = substituteOID (temps, exprName, expr);
- String newexpr = expr;
- PtrBlock<const ImageRegion*> tempRegs;
- makeRegionBlock(tempRegs, regions, *itsLog);
- LatticeExprNode node = ImageExprParse::command(newexpr, temps, tempRegs);
-
- // Get the shape of the expression
- const IPosition shapeOut = node.shape();
-
- // Get the CoordinateSystem of the expression
- const LELAttribute attr = node.getAttribute();
- const LELLattCoordBase* lattCoord = &(attr.coordinates().coordinates());
- if (!lattCoord->hasCoordinates() || lattCoord->classname()
- != "LELImageCoord") {
- *itsLog << "Images in expression have no coordinates"
- << LogIO::EXCEPTION;
- }
- const LELImageCoord* imCoord =
- dynamic_cast<const LELImageCoord*> (lattCoord);
-
- AlwaysAssert (imCoord != 0, AipsError);
- CoordinateSystem cSysOut = imCoord->coordinates();
-
- // Create LatticeExpr create mask if needed
- LatticeExpr<Float> latEx(node);
-
- // Construct output image - an ImageExpr or a PagedImage
- if (outfile.empty()) {
- pImage_p = new ImageExpr<Float> (latEx, exprName);
- if (pImage_p == 0) {
- *itsLog << "Failed to create ImageExpr" << LogIO::EXCEPTION;
- }
- } else {
- *itsLog << LogIO::NORMAL << "Creating image `" << outfile
- << "' of shape " << shapeOut << LogIO::POST;
- pImage_p = new PagedImage<Float> (shapeOut, cSysOut, outfile);
- if (pImage_p == 0) {
- *itsLog << "Failed to create PagedImage" << LogIO::EXCEPTION;
- }
-
- // Make mask if needed, and copy data and mask
- if (latEx.isMasked()) {
- String maskName("");
- makeMask(*pImage_p, maskName, False, True, *itsLog, True);
- }
- LatticeUtilities::copyDataAndMask(*itsLog, *pImage_p, latEx);
- }
-
- // Copy miscellaneous stuff over
- pImage_p->setMiscInfo(imCoord->miscInfo());
- pImage_p->setImageInfo(imCoord->imageInfo());
- //
- if (expr.contains("spectralindex")) {
- pImage_p->setUnits("");
- } else if (expr.contains(Regex("pa\\(*"))) {
- pImage_p->setUnits("deg");
- Vector<Int> newstokes(1);
- newstokes = Stokes::Pangle;
- StokesCoordinate scOut(newstokes);
- CoordinateSystem cSys = pImage_p->coordinates();
- Int iStokes = cSys.findCoordinate(Coordinate::STOKES, -1);
- cSys.replaceCoordinate(scOut, iStokes);
- pImage_p->setCoordinateInfo(cSys);
- } else {
- pImage_p->setUnits(imCoord->unit());
- }
-
- // Logger not yet available
- // pImage_p->appendLog(imCoord->logger());
-
- // Delete the ImageRegions (by using an empty GlishRecord).
- makeRegionBlock(tempRegs, Record(), *itsLog);
-
- return pImage_p;
-
-}
-
-ImageInterface<Float> *
-ImageAnalysis::imageconcat(const String& outfile,
- const Vector<String>& inFiles, const Int axis, const Bool relax,
- const Bool tempclose, const Bool overwrite) {
- *itsLog << LogOrigin("ImageAnalysis", "imageconcat");
-
- // There could be wild cards embedded in our list so expand them out
- Vector<String> expInNames = Directory::shellExpand(inFiles, False);
- if (expInNames.nelements() <= 1) {
- *itsLog << "You must give at least two valid input images"
- << LogIO::EXCEPTION;
- }
- *itsLog << LogIO::NORMAL << "Number of expanded file names = "
- << expInNames.nelements() << LogIO::POST;
-
- // Verify output file
- String outFile(outfile);
- if (!outFile.empty() && !overwrite) {
- NewFile validfile;
- String errmsg;
- if (!validfile.valueOK(outFile, errmsg)) {
- *itsLog << errmsg << LogIO::EXCEPTION;
- }
- }
-
- // Find spectral axis of first image
- PtrHolder<ImageInterface<Float> > im;
- ImageUtilities::openImage(im, expInNames(0), *itsLog);
-
- // PagedImage<Float> im(expInNames(0), MaskSpecifier(True));
- CoordinateSystem cSys = im.ptr()->coordinates();
- Int iAxis = axis;
- if (iAxis < 0) {
- iAxis = CoordinateUtil::findSpectralAxis(cSys);
- if (iAxis < 0) {
- *itsLog << "Could not find a spectral axis in first input image"
- << LogIO::EXCEPTION;
- }
- }
-
- // Create concatenator. Use holder so if exceptions, the ImageConcat
- // object gets cleaned up
- uInt axis2 = uInt(iAxis);
- PtrHolder<ImageConcat<Float> > pConcat(new ImageConcat<Float> (axis2,
- tempclose));
-
- // Set first image
- pConcat.ptr()->setImage(*(im.ptr()), relax);
-
- // Set the other images. We may run into the open file limit.
- for (uInt i = 1; i < expInNames.nelements(); i++) {
- Bool doneOpen = False;
- try {
- PtrHolder<ImageInterface<Float> > im2;
- ImageUtilities::openImage(im2, expInNames(i), *itsLog);
- // PagedImage<Float> im2(expInNames(i), MaskSpecifier(True));
- doneOpen = True;
- pConcat.ptr()->setImage(*(im2.ptr()), relax);
- } catch (AipsError x) {
- if (!doneOpen) {
- *itsLog << "Failed to open file " << expInNames(i) << endl;
- *itsLog
- << "This may mean you have too many files open simultaneously"
- << endl;
- *itsLog
- << "Try using tempclose=T in the imageconcat constructor"
- << LogIO::EXCEPTION;
- } else {
- *itsLog << x.getMesg() << LogIO::EXCEPTION;
- }
- }
- }
- //
- if (!outFile.empty()) {
- // Construct output image and give it a mask if needed
- pImage_p = new PagedImage<Float> (pConcat.ptr()->shape(),
- pConcat.ptr()->coordinates(), outFile);
- if (!pImage_p) {
- *itsLog << "Failed to create PagedImage" << LogIO::EXCEPTION;
- }
- *itsLog << LogIO::NORMAL << "Creating image '" << outfile
- << "' of shape " << pImage_p->shape() << LogIO::POST;
- //
- if (pConcat.ptr()->isMasked()) {
- String maskName("");
- makeMask(*pImage_p, maskName, False, True, *itsLog, True);
- }
-
- // Copy to output
- LatticeUtilities::copyDataAndMask(*itsLog, *pImage_p, *(pConcat.ptr()));
- ImageUtilities::copyMiscellaneous(*pImage_p, *(pConcat.ptr()));
- } else {
- pImage_p = pConcat.ptr()->cloneII();
- }
- return pImage_p;
-}
-
-Bool ImageAnalysis::imagefromarray(const String& outfile,
- Array<Float> & pixelsArray, const Record& csys, const Bool linear,
- const Bool overwrite, const Bool log) {
-
- Bool rstat = False;
- try {
- *itsLog << LogOrigin("ImageAnalysis", "imagefromarray");
-
- String error;
- if (csys.nfields() != 0) {
- PtrHolder<CoordinateSystem> cSys(makeCoordinateSystem(csys,
- pixelsArray.shape()));
- CoordinateSystem* pCS = cSys.ptr();
- if (!make_image(error, outfile, *pCS, pixelsArray.shape(), *itsLog,
- log, overwrite)) {
- *itsLog << error << LogIO::EXCEPTION;
- }
- } else {
- // Make default CoordinateSystem
- CoordinateSystem cSys = CoordinateUtil::makeCoordinateSystem(
- pixelsArray.shape(), linear);
- centreRefPix(cSys, pixelsArray.shape());
- if (!make_image(error, outfile, cSys, pixelsArray.shape(), *itsLog,
- log, overwrite)) {
- *itsLog << error << LogIO::EXCEPTION;
- }
- }
-
- // Fill image
- pImage_p->putSlice(pixelsArray, IPosition(pixelsArray.ndim(), 0),
- IPosition(pixelsArray.ndim(), 1));
- rstat = True;
- } catch (AipsError x) {
- *itsLog << LogIO::SEVERE << "Exception Reported: " << x.getMesg()
- << LogIO::POST;
- }
- return rstat;
-}
-
-Bool ImageAnalysis::imagefromascii(const String& outfile, const String&
infile,
- const Vector<Int>& shape, const String& sep, const Record& csys,
- const Bool linear, const Bool overwrite) {
- // The glish code ignored sep (assumed to be ' ') so will we to
- Bool rstat = False;
-
- try {
- *itsLog << LogOrigin("ImageAnalysis", "imagefromascii");
-
- Path filePath(infile);
- String fileName = filePath.expandedName();
-
- ifstream inFile(fileName.c_str());
- if (!inFile) {
- *itsLog << LogIO::SEVERE << "Cannot open " << infile << LogIO::POST;
- return false;
- }
-
- IPosition shp(shape);
- uInt n = shp.product();
- uInt nx = shp(0);
- Vector<Float> a(n, 0.0);
- int idx = 0;
- string line;
- vector<string> line2(2 * nx);
- uInt iline = 0;
- uInt nl = 1;
- while (nl > 0) {
- getline(inFile, line, '\n');
- nl = split(line, &(line2[0]), 2 * nx, sep);
- if (nl > 0) {
- if (nl != nx) {
- *itsLog << LogIO::SEVERE << "Length of line " << iline
- << " is " << nl << " but should be " << nx
- << LogIO::POST;
- return false;
- }
- for (uInt i = 0; i < nx; i++) {
- a[idx + i] = atof(line2[i].c_str());
- }
- idx += nx;
- iline += 1;
- }
- }
- Vector<Float> vec(n);
- for (uInt i = 0; i < n; i++)
- vec[i] = a[i];
- Array<Float> pixels(vec.reform(IPosition(shape)));
- rstat = imagefromarray(outfile, pixels, csys, linear, overwrite);
-
- } catch (AipsError x) {
- *itsLog << LogIO::SEVERE << "Exception Reported: " << x.getMesg()
- << LogIO::POST;
- }
- return rstat;
-}
-
-Bool ImageAnalysis::imagefromfits(const String& outfile,
- const String& fitsfile, const Int whichrep, const Int whichhdu,
- const Bool zeroBlanks, const Bool overwrite) {
- Bool rstat = False;
- try {
- *itsLog << LogOrigin("ImageAnalysis", "imagefromfits");
-
- // Check output file
- if (!overwrite && !outfile.empty()) {
- NewFile validfile;
- String errmsg;
- if (!validfile.valueOK(outfile, errmsg)) {
- *itsLog << errmsg << LogIO::EXCEPTION;
- }
- }
- //
- if (whichrep < 0) {
- *itsLog << "The Coordinate Representation index must be non-negative"
- << LogIO::EXCEPTION;
- }
- //
- ImageInterface<Float>* pOut = 0;
- String error;
- ImageFITSConverter::FITSToImage(pOut, error, outfile, fitsfile,
- whichrep, whichhdu, HostInfo::memoryFree() / 1024, overwrite,
- zeroBlanks);
- //
- if (pOut == 0) {
- *itsLog << error << LogIO::EXCEPTION;
- }
- pImage_p = pOut;
- rstat = True;
- } catch (AipsError x) {
- *itsLog << LogIO::SEVERE << "Exception Reported: " << x.getMesg()
- << LogIO::POST;
- }
- return rstat;
-}
-
-Bool ImageAnalysis::imagefromimage(const String& outfile, const String&
infile,
- Record& region, const String& Mask, const bool dropdeg,
- const bool overwrite) {
- Bool rstat = False;
- try {
- *itsLog << LogOrigin("ImageAnalysis", "imagefromimage");
-
- // Open
- PtrHolder<ImageInterface<Float> > inImage;
- ImageUtilities::openImage(inImage, infile, *itsLog);
- ImageInterface<Float>* pInImage = inImage.ptr();
- //
- // Convert region from Glish record to ImageRegion.
- // Convert mask to ImageRegion and make SubImage.
- //
- ImageRegion* pRegionRegion = 0;
- ImageRegion* pMaskRegion = 0;
- AxesSpecifier axesSpecifier;
- if (dropdeg)
- axesSpecifier = AxesSpecifier(False);
- SubImage<Float> subImage = SubImage<Float>::createSubImage(
- pRegionRegion, pMaskRegion,
- *pInImage, region,
- Mask, itsLog, True, axesSpecifier
- );
- delete pRegionRegion;
- delete pMaskRegion;
-
- // Create output image
-
- if (outfile.empty()) {
- pImage_p = new SubImage<Float> (subImage);
- rstat = True;
- } else {
- if (!overwrite) {
- NewFile validfile;
- String errmsg;
- if (!validfile.valueOK(outfile, errmsg)) {
- *itsLog << errmsg << LogIO::EXCEPTION;
- }
- }
- //
- *itsLog << LogIO::NORMAL << "Creating image '" << outfile
- << "' of shape " << subImage.shape() << LogIO::POST;
- pImage_p = new PagedImage<Float> (subImage.shape(),
- subImage.coordinates(), outfile);
- if (pImage_p == 0) {
- *itsLog << "Failed to create PagedImage" << LogIO::EXCEPTION;
- }
- ImageUtilities::copyMiscellaneous(*pImage_p, *pInImage);
-
- // Make output mask if required
-
- if (subImage.isMasked()) {
- String maskName("");
- makeMask(*pImage_p, maskName, False, True, *itsLog, True);
- }
-
- // Copy data and mask
-
- LatticeUtilities::copyDataAndMask(*itsLog, *pImage_p, subImage);
- rstat = True;
- }
- } catch (AipsError x) {
- *itsLog << LogIO::SEVERE << "Exception Reported: " << x.getMesg()
- << LogIO::POST;
- }
- return rstat;
-}
-
-Bool ImageAnalysis::imagefromshape(const String& outfile,
- const Vector<Int>& shapeV, const Record& coordinates,
- const Bool linear, const Bool overwrite, const Bool log) {
- Bool rstat = False;
- try {
- *itsLog << LogOrigin("ImageAnalysis", "imagefromshape");
-
- // Some protection
- if (shapeV.nelements() == 0) {
- *itsLog << "The shape is invalid" << LogIO::EXCEPTION;
- }
- for (uInt i = 0; i < shapeV.nelements(); i++) {
- if (shapeV(i) <= 0) {
- *itsLog << "The shape is invalid" << LogIO::EXCEPTION;
- }
- }
-
- // Make with supplied CoordinateSystem if record not empty
- String error;
- if (coordinates.nfields() > 0) { //
- PtrHolder<CoordinateSystem> pCS(makeCoordinateSystem(coordinates,
- shapeV));
- if (!make_image(error, outfile, *(pCS.ptr()), shapeV, *itsLog, log,
- overwrite)) {
- *itsLog << error << LogIO::EXCEPTION;
- }
- } else {
- // Make default CoordinateSystem
- CoordinateSystem cSys = CoordinateUtil::makeCoordinateSystem(
- shapeV, linear);
- centreRefPix(cSys, shapeV);
- if (!make_image(error, outfile, cSys, shapeV, *itsLog, log,
- overwrite)) {
- *itsLog << error << LogIO::EXCEPTION;
- }
- }
- pImage_p->set(0.0);
- rstat = True;
- } catch (AipsError x) {
- *itsLog << LogIO::SEVERE << "Exception Reported: " << x.getMesg()
- << LogIO::POST;
- }
- return rstat;
-}
-
-Bool ImageAnalysis::adddegaxes(const String& outfile,
PtrHolder<ImageInterface<
- Float> >& outImage, const Bool direction, const Bool spectral,
- const String& stokes, const Bool linear, const Bool tabular,
- const Bool overwrite) {
-
- *itsLog << LogOrigin("ImageAnalysis", "adddegaxes");
-
- String outFile(outfile);
- String sStokes(stokes);
-
- ImageUtilities::addDegenerateAxes(*itsLog, outImage, *pImage_p, outFile,
- direction, spectral, sStokes, linear, tabular, overwrite);
-
- return True;
-
-}
-
-ImageInterface<Float> *
-ImageAnalysis::convolve(const String& outFile, Array<Float>& kernelArray,
- const String& kernelFileName, const Double in_scale, Record& region,
- String& mask, const Bool overwrite, const Bool) {
-
- *itsLog << LogOrigin("ImageAnalysis", "convolve");
-
- //Need to deal with the string part
- // String kernelFileName(kernel.toString());
- if (mask == "[]")
- mask = "";
- Bool autoScale;
- Double scale(in_scale);
-
- if (scale > 0) {
- autoScale = False;
- } else {
- autoScale = True;
- scale = 1.0;
- }
-
- // Check output file name
- if (!outFile.empty() && !overwrite) {
- NewFile validfile;
- String errmsg;
- if (!validfile.valueOK(outFile, errmsg)) {
- *itsLog << errmsg << LogIO::EXCEPTION;
- }
- }
-
- // Convert region from Glish record to ImageRegion. Convert mask
- // to ImageRegion and make SubImage.
- ImageRegion* pRegionRegion = 0;
- ImageRegion* pMaskRegion = 0;
- SubImage<Float> subImage = SubImage<Float>::createSubImage(
- pRegionRegion, pMaskRegion, *pImage_p,
- region,
- mask, itsLog, False
- );
- delete pRegionRegion;
- delete pMaskRegion;
-
- // Create output image
- IPosition outShape = subImage.shape();
- PtrHolder<ImageInterface<Float> > imOut;
- if (outFile.empty()) {
- *itsLog << LogIO::NORMAL << "Creating (temp)image of shape "
- << outShape << LogIO::POST;
- imOut.set(new TempImage<Float> (outShape, subImage.coordinates()));
- } else {
- *itsLog << LogIO::NORMAL << "Creating image '" << outFile
- << "' of shape " << outShape << LogIO::POST;
- imOut.set(new PagedImage<Float> (outShape, subImage.coordinates(),
- outFile));
- }
- ImageInterface<Float>* pImOut = imOut.ptr()->cloneII();
-
- // Make the convolver
- ImageConvolver<Float> aic;
- Bool copyMisc = True;
- Bool warnOnly = True;
- ImageConvolver<Float>::ScaleTypes scaleType(ImageConvolver<Float>::NONE);
- if (autoScale) {
- scaleType = ImageConvolver<Float>::AUTOSCALE;
- } else {
- scaleType = ImageConvolver<Float>::SCALE;
- }
- if (kernelFileName.empty()) {
- if (kernelArray.nelements() > 1) {
- aic.convolve(*itsLog, *pImOut, subImage, kernelArray, scaleType,
- scale, copyMisc);
- } else {
- *itsLog << "Kernel array dimensions are invalid"
- << LogIO::EXCEPTION;
- }
- } else {
- if (!Table::isReadable(kernelFileName)) {
- *itsLog << LogIO::SEVERE << "kernel image " << kernelFileName
- << " is not available " << LogIO::POST;
- return 0;
- }
- PagedImage<Float> kernelImage(kernelFileName);
- aic.convolve(*itsLog, *pImOut, subImage, kernelImage, scaleType, scale,
- copyMisc, warnOnly);
- }
-
- return pImOut;
-}
-
-Record*
-ImageAnalysis::boundingbox(const Record& Region) {
- *itsLog << LogOrigin("ImageAnalysis", "boundingbox");
- // Find the bounding box of this region
- Record tmpR(Region);
- const ImageRegion* pRegion = ImageRegion::fromRecord(
- 0, pImage_p->coordinates(), pImage_p->shape(),
- tmpR
- );
- LatticeRegion latRegion =
pRegion->toLatticeRegion(pImage_p->coordinates(),
- pImage_p->shape());
- //
- Slicer sl = latRegion.slicer();
- IPosition blc(sl.start()); // 1-rel for Glish
- IPosition trc(sl.end());
- IPosition inc(sl.stride());
- IPosition length(sl.length());
- RecordDesc outRecDesc;
- outRecDesc.addField("blc", TpArrayInt);
- outRecDesc.addField("trc", TpArrayInt);
- outRecDesc.addField("inc", TpArrayInt);
- outRecDesc.addField("bbShape", TpArrayInt);
- outRecDesc.addField("regionShape", TpArrayInt);
- outRecDesc.addField("imageShape", TpArrayInt);
- outRecDesc.addField("blcf", TpString);
- outRecDesc.addField("trcf", TpString);
- Record *outRec = new Record(outRecDesc);
- outRec->define("blc", blc.asVector());
- outRec->define("trc", trc.asVector());
- outRec->define("inc", inc.asVector());
- outRec->define("bbShape", (trc - blc + 1).asVector());
- outRec->define("regionShape", length.asVector());
- outRec->define("imageShape", pImage_p->shape().asVector());
- //
- CoordinateSystem cSys(pImage_p->coordinates());
- outRec->define("blcf", CoordinateUtil::formatCoordinate(blc, cSys)); //
0-rel for use in C++
- outRec->define("trcf", CoordinateUtil::formatCoordinate(trc, cSys));
- return outRec;
-}
-
-String ImageAnalysis::brightnessunit() {
- String rstat;
- *itsLog << LogOrigin("ImageAnalysis", "brightnessunit");
- rstat = pImage_p->units().getName();
- return rstat;
-}
-
-Bool ImageAnalysis::calc(const String& expr) {
-
- *itsLog << LogOrigin("ImageAnalysis", "calc");
- Record regions;
-
- // Get LatticeExprNode (tree) from parser
- // Convert the GlishRecord containing regions to a
- // PtrBlock<const ImageRegion*>.
- if (expr.empty()) {
- *itsLog << "You must specify an expression" << LogIO::EXCEPTION;
- return False;
- }
-
- Block<LatticeExprNode> temps;
- // String exprName;
- // String newexpr = substituteOID (temps, exprName, expr);
- String newexpr = expr;
- PtrBlock<const ImageRegion*> tempRegs;
- makeRegionBlock(tempRegs, regions, *itsLog);
- LatticeExprNode node = ImageExprParse::command(newexpr, temps, tempRegs);
-
- // Delete the ImageRegions (by using an empty GlishRecord)
- makeRegionBlock(tempRegs, Record(), *itsLog);
-
- // Get the shape of the expression and check it matches that
- // of the output image
- if (!node.isScalar()) {
- const IPosition shapeOut = node.shape();
- if (!pImage_p->shape().isEqual(shapeOut)) {
- *itsLog << LogIO::SEVERE
- << "The shape of the expression does not conform " << endl;
- *itsLog << "with the shape of the output image" << LogIO::POST;
- *itsLog << "Expression shape = " << shapeOut << endl;
- *itsLog << "Image shape = " << pImage_p->shape()
- << LogIO::EXCEPTION;
- }
- }
-
- // Get the CoordinateSystem of the expression and check it matches
- // that of the output image
- if (!node.isScalar()) {
- const LELAttribute attr = node.getAttribute();
- const LELLattCoordBase* lattCoord = &(attr.coordinates().coordinates());
- if (!lattCoord->hasCoordinates() || lattCoord->classname()
- != "LELImageCoord") {
- // We assume here that the output coordinates are ok
- *itsLog << LogIO::WARN
- << "Images in expression have no coordinates"
- << LogIO::POST;
- } else {
- const LELImageCoord* imCoord =
- dynamic_cast<const LELImageCoord*> (lattCoord);
- AlwaysAssert (imCoord != 0, AipsError);
- const CoordinateSystem& cSysOut = imCoord->coordinates();
- if (!pImage_p->coordinates().near(cSysOut)) {
- // Since the output image has coordinates, and the shapes
- // have conformed, just issue a warning
- *itsLog << LogIO::WARN
- << "The coordinates of the expression do not conform "
- << endl;
- *itsLog << "with the coordinates of the output image" << endl;
- *itsLog << "Proceeding with output image coordinates"
- << LogIO::POST;
- }
- }
- }
- // Make a LatticeExpr and see if it is masked
- Bool exprIsMasked = node.isMasked();
- if (exprIsMasked) {
- if (!pImage_p->isMasked()) {
- // The image does not have a default mask set. So try and make it one.
- String maskName("");
- makeMask(*pImage_p, maskName, True, True, *itsLog, True);
- }
- }
- // Evaluate the expression and fill the output image and mask
- if (node.isScalar()) {
- LatticeExprNode node2 = toFloat(node);
- // If the scalar value is masked, there is nothing
- // to do.
- if (!exprIsMasked) {
- Float value = node2.getFloat();
- if (pImage_p->isMasked()) {
- // We implement with a LEL expression of the form
- // iif(mask(image)", value, image)
- LatticeExprNode node3 = iif(mask(*pImage_p), node2, *pImage_p);
- pImage_p->copyData(LatticeExpr<Float> (node3));
- } else {
- // Just set all values to the scalar. There is no mask to
- // worry about.
- pImage_p->set(value);
- }
- }
- } else {
- if (pImage_p->isMasked()) {
- // We implement with a LEL expression of the form
- // iif(mask(image)", expr, image)
***The diff for this file has been truncated for email.***
=======================================
--- /trunk/images/Images/ImageAnalysis.h Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,760 +0,0 @@
-//# ImageAnalysis.h: Image analysis and handling tool
-//# Copyright (C) 2007
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or (at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-
-#ifndef _IMAGEANALYSIS__H__
-#define _IMAGEANALYSIS__H__
-
-
-//# put includes here
-#include <coordinates/Coordinates/CoordinateSystem.h>
-#include <lattices/LatticeMath/Fit2D.h>
-#include <casa/Quanta.h>
-#include <casa/Utilities/PtrHolder.h>
-#include <measures/Measures/Stokes.h>
-#include <images/Images/ImageFit1D.h>
-#include <images/Images/ImageInfo.h>
-#include <images/Images/ImageInterface.h>
-#include <components/ComponentModels/ComponentType.h>
-#include <casa/Arrays/AxesSpecifier.h>
-#include <measures/Measures/Stokes.h>
-
-namespace casa {
-
-class DirectionCoordinate;
-class LogIO;
-class SkyComponent;
-class Record;
-class Fit2D;
-class ImageRegion;
-class ComponentList;
-template<class T> class Array;
-template<class T> class Block;
-template<class T> class PtrBlock;
-template<class T> class Flux;
-template<class T> class ImageStatistics;
-template<class T> class ImageHistograms;
-template<class T> class MaskedArray;
-template<class T> class Quantum;
-template<class T> class SubLattice;
-template<class T> class SubImage;
-template<class T> class Vector;
-
-// <summary>
-// Image analysis and handling tool
-// </summary>
-
-// <synopsis>
-// This the casapy image tool.
-// One time it should be merged with pyrap's image tool ImageProxy.
-// </synopsis>
-
-class ImageAnalysis
-{
- public:
-
- ImageAnalysis();
-
- //ImageInterface constructor
- ImageAnalysis(const ImageInterface<Float>* inImage);
-
- //Use this constructor with cloneInputPointer=False if you want this
object
- // to take over management of the input pointer. The input pointer
will be deleted
- // when this object is destroyed.
- ImageAnalysis(ImageInterface<Float>* inImage, const Bool
cloneInputPointer);
-
- virtual ~ImageAnalysis();
-
- Bool addnoise(const String& type, const Vector<Double>& pars,
- Record& region, const Bool zero = False);
-
- ImageInterface<Float> * imagecalc(const String& outfile,
- const String& pixels,
- const Bool overwrite = False);
-
- ImageInterface<Float> * imageconcat(const String& outfile,
- const Vector<String>& infiles,
- const Int axis,
- const Bool relax = False,
- const Bool tempclose = True,
- const Bool overwrite = False);
-
- Bool imagefromarray(const String& outfile, Array<Float>& pixels,
- const Record& csys, const Bool linear = False,
- const Bool overwrite = False, const Bool log =
True);
-
- Bool imagefromascii(const String& outfile, const String& infile,
- const Vector<Int>& shape, const String& sep,
- const Record& csys, const Bool linear = False,
- const Bool overwrite = False);
-
- Bool imagefromfits(const String& outfile, const String& infile,
- const Int whichrep = 0, const Int whichhdu = 0,
- const Bool zeroblanks = False,
- const Bool overwrite = False);
-
- Bool imagefromimage(const String& outfile, const String& infile,
- Record& region, const String& mask,
- const Bool dropdeg = False,
- const Bool overwrite = False);
-
- Bool imagefromshape(const String& outfile, const Vector<Int>& shape,
- const Record& csys, const Bool linear = True,
- const Bool overwrite = False, const Bool log =
True);
-
- Bool adddegaxes(const String& outfile,
- PtrHolder<ImageInterface<Float> >& ph,
- const Bool direction,
- const Bool spectral,
- const String& stokes,
- const Bool linear = False,
- const Bool tabular = False,
- const Bool overwrite = False);
-
- ImageInterface<Float> * convolve(const String& outfile,
- Array<Float>& kernel,
- const String& kernImage,
- const Double scale,
- Record& region, String& mask,
- const Bool overwrite = False,
- const Bool async = False);
-
- Record* boundingbox(const Record& region);
-
- String brightnessunit();
-
- Bool calc(const String& pixels);
-
- // regions should be a Record of Records having different regions
-
- Bool calcmask(const String& mask, Record& regions, const String& name,
- const Bool asdefault = True);
-
-
- ImageInterface<Float> * continuumsub(const String& outline,
- const String& outcont, Record&
region,
- const Vector<int>& channels,
- const String& pol = "",
- const Int fitorder = 0,
- const Bool overwrite = false);
-
- // the output <src>fakeBeam</src> indicates if there was no beam in
the header and a fake one
- // was assumed to do the conversion.
- Quantity convertflux(Bool& fakeBeam, const Quantity& value, const
Quantity& major,
- const Quantity& minor,
- const String& type = "Gaussian",
- const Bool topeak = True,
- Bool supressNoBeamWarnings = False);
-
- ImageInterface<Float>* convolve2d(
- const String& outfile, const Vector<Int>& axes,
- const String& type, const Quantity& major,
- const Quantity& minor, const Quantity& pa,
- Double scale, Record& region, const String& mask,
- const Bool overwrite = False
- );
-
- CoordinateSystem coordsys(const Vector<int>& axes);
-
- CoordinateSystem csys(const Vector<int>& axes);
-
- Record* coordmeasures(Quantity& intensity, Record& direction,
- Record& frequency, Record& velocity,
- const Vector<double>& pixel);
-
- Matrix<Float> decompose(Record& region, const String& mask,
- const Bool simple = false,
- const Double threshold = -1,
- const Int ncontour = 11,
- const Int minrange = 1,
- const Int naxis = 2,
- const Bool fit = True,
- const Double maxrms = -1,
- const Int maxretry = -1,
- const Int maxiter = 256,
- const Double convcriteria = 0.0001);
-
- Matrix<Float> decompose(Matrix<Int>& blcs, Matrix<Int>& trcs, Record&
region, const String& mask,
- const Bool simple = false,
- const Double threshold = -1,
- const Int ncontour = 11,
- const Int minrange = 1,
- const Int naxis = 2,
- const Bool fit = True,
- const Double maxrms = -1,
- const Int maxretry = -1,
- const Int maxiter = 256,
- const Double convcriteria = 0.0001);
-
- Record deconvolvecomponentlist(Record& complist);
-
- Bool remove(Bool verbose=true);
-
- Bool fft(const String& real, const String& imag, const String& amp,
- const String& phase, const Vector<Int>& axes, Record& region,
- const String& mask);
-
- Record findsources(const Int nmax, const Double cutoff, Record& region,
- const String& mask, const Bool point = True,
- const Int width = 5, const Bool negfind = False);
-
- // <src>subImage</src> will contain the subimage of the original image
- // on which the fit is performed.
- Vector<ImageFit1D<Float> > fitallprofiles(
- Record& region, SubImage<Float>& subImage, String& xUnit,
- const Int axis, const String& mask,
- const Int ngauss, const Int poly,
- const String& weightsImageName = "", const String& fit = "",
- const String& resid = ""
- ) const;
-
- // <src>subImage</src> will contain the subimage of the original image
- // on which the fit is performed.
- ImageFit1D<Float> fitprofile(
- Record& regionRecord, SubImage<Float>& subImage,
- String& xUnit,
- const uInt axis, const String& mask,
- const Record& estimate, const uInt ngauss = 0,
- const Int poly = -1, const String& modelName = "",
- const String& residName = "", const Bool fit = True,
- const String weightsImageName = ""
- );
-
- ImageInterface<Float>* fitpolynomial(const String& residfile,
- const String& fitfile,
- const String& sigmafile,
- const Int axis, const Int order,
- Record& region, const String&
mask,
- const bool overwrite = false);
-
- // adding residImage parameter. If set to blank, do not write residual
iamge,
- // else write to the file given.
- ComponentList fitsky(
- Array<Float>& pixels, Array<Bool>& pixelmask,
- Bool& converged, Record& inputStats, Record& residStats,
- Double& chiSquared, Record& region, const uInt& chan,
- const String& stokesString, const String& mask,
- const Vector<String>& models, Record& inputEstimate,
- const Vector<String>& fixedparams,
- const Vector<Float>& includepix,
- const Vector<Float>& excludepix,
- const Bool fit=True,
- const Bool deconvolve=False, const Bool list=True,
- const String& residImage="", const String& modelImageName=""
- );
-
- Bool getchunk(Array<Float>& pixel, Array<Bool>& pixmask,
- const Vector<Int>& blc, const Vector<Int>& trc,
- const Vector<Int>& inc, const Vector<Int>& axes,
- const Bool list = False, const Bool dropdeg = False,
- const bool getmask = False);
-
- Bool getregion(Array<Float>& pixels, Array<Bool>& pixmask, Record&
region,
- const Vector<Int>& axes, const String& mask,
- const Bool list = False, const Bool dropdeg = False,
- const bool getmask = False);
-
- Record* getslice(const Vector<Double>& x, const Vector<Double>& y,
- const Vector<Int>& axes, const Vector<Int>& coord,
- const Int npts = 0, const String& method = "linear");
-
- ImageInterface<Float> * hanning(const String& outfile, Record& region,
- const String& mask, const Int axis =
-10,
- const Bool drop = True,
- const bool overwrite = False);
-
- Vector<Bool> haslock();
-
- Bool histograms(Record& histout, const Vector<Int>& axes, Record&
region,
- const String& mask, const Int nbins,
- const Vector<Double>& includepix, const Bool gauss,
- const Bool cumu, const Bool log, const Bool list,
- const String& plotter, const Int nx, const Int ny,
- const Vector<Int>& size, const Bool force = False,
- const Bool disk = False);
-
- Vector<String> history(const Bool list = False, const Bool browse =
True);
-
- ImageInterface<Float> * insert(const String& infile, Record& region,
- const Vector<double>& locate);
-
- // Bool isopen();
-
- Bool ispersistent();
-
- Bool lock(const Bool writelock = False, const Int nattempts = 0);
-
- Bool makecomplex(const String& outfile, const String& imag, Record&
region,
- const Bool overwrite = False);
-
- Vector<String> maskhandler(const String& op,const Vector<String>& nam);
-
- Record miscinfo();
-
- Bool modify(Record& model, Record& region , const String& mask,
- const Bool subtract = True, const Bool list = True);
-
- Record maxfit(Record& region, const Bool point, const Int width = 5,
- const Bool negfind = False, const Bool list = True);
-
- ImageInterface<Float> * moments(const Vector<Int>& moments, const Int
axis,
- Record& region, const String& mask,
- const Vector<String>& method,
- const Vector<Int>& smoothaxes,
- const Vector<String>& smoothtypes,
- const Vector<Quantity>& smoothwidths,
- const Vector<Float>& includepix,
- const Vector<Float>& excludepix,
- const Double peaksnr, const Double
stddev,
- const String& doppler = "RADIO",
- const String& outfile = "",
- const String& smoothout = "",
- const String& plotter = "/NULL",
- const Int nx = 1, const Int ny = 1,
- const Bool yind = False,
- const Bool overwrite = False,
- const Bool drop = True);
-
- String name(const Bool strippath = False);
-
- Bool open(const String& infile);
-
- Record* pixelvalue(const Vector<Int>& pixel);
- void pixelValue (Bool& offImage, Quantum<Double>& value, Bool& mask,
- Vector<Int>& pos) const;
-
- Bool putchunk(const Array<Float>& pixels, const Vector<Int>& blc,
- const Vector<Int>& inc, const Bool list = False,
- const Bool locking = True, const Bool replicate = False);
-
- Bool putregion(const Array<Float>& pixels, const Array<Bool>&
pixelmask,
- Record& region, const Bool list = False,
- const Bool usemask = True,
- const Bool locking = True, const Bool replicate =
False);
-
- ImageInterface<Float> * rebin(const String& outfile,
- const Vector<Int>& bin, Record& region,
- const String& mask, const Bool dropdeg,
- const Bool overwrite = False);
-
- //regrids to a given coordinate system...one uses a record that is
- //converted to a CoordinateSytem
-
- ImageInterface<Float> * regrid(const String& outfile,
- const Vector<Int>& shape,
- const Record& csys, const Vector<Int>&
axes,
- Record& region, const String& mask,
- const String& method = "linear",
- const Int decimate = 10,
- const Bool replicate = False,
- const Bool doref = True,
- const Bool dropdeg = False,
- const Bool overwrite = False,
- const Bool force = False);
-
- ImageInterface<Float> * regrid(const String& outfile,
- const Vector<Int>& shape,
- const CoordinateSystem& csys,
- const Vector<Int>& axes,
- Record& region, const String& mask,
- const String& method = "linear",
- const Int decimate = 10,
- const Bool replicate = False,
- const Bool doref = True,
- const Bool dropdeg = False,
- const Bool overwrite = False,
- const Bool force = False);
-
- ImageInterface<Float> * rotate(const String& outfile,
- const Vector<int>& shape,
- const Quantity& pa, Record& region,
- const String& mask,
- const String& method = "cubic",
- const Int decimate = 0,
- const Bool replicate = False,
- const Bool dropdeg = False,
- const Bool overwrite = False);
-
- Bool rename(const String& name, const Bool overwrite = False);
-
- Bool replacemaskedpixels(const String& pixels, Record& region,
- const String& mask, const Bool update = False,
- const Bool list = False);
-
- Record restoringbeam();
-
- ImageInterface<Float> * sepconvolve(const String& outfile,
- const Vector<Int>& axes,
- const Vector<String>& types,
- const Vector<Quantity>& widths,
- Double scale,
- Record& region,
- const String& mask,
- const bool overwrite = False);
-
- Bool set(const String& pixels, const Int pixelmask,
- Record& region, const Bool list = false);
-
- Bool setbrightnessunit(const String& unit);
-
- bool setcoordsys(const Record& csys);
-
- bool sethistory(const String& origin, const Vector<String>& history);
-
- bool setmiscinfo(const Record& info);
-
- Vector<Int> shape();
-
- Bool setrestoringbeam(const Quantity& major, const Quantity& minor,
- const Quantity& pa, const Record& beam,
- const Bool remove = False, const Bool log =
True);
-
- Bool statistics(
- Record& statsout, const Vector<Int>& axes, Record& region,
- const String& mask, const Vector<String>& plotstats,
- const Vector<Float>& includepix, const Vector<Float>& excludepix,
- const String& plotter = "/NULL", const Int nx = 1,
- const Int ny = 1, const Bool list = True,
- const Bool force = False, const Bool disk = False,
- const Bool robust = False, const Bool verbose = True
- );
-
- bool twopointcorrelation(const String& outfile, Record& region,
- const String& mask, const Vector<Int>& axes,
- const String& method = "structurefunction",
- const Bool overwrite = False);
-
- ImageInterface<Float> * subimage(const String& outfile, Record& region,
- const String& mask,
- const Bool dropdeg = False,
- const Bool overwrite = False,
- const Bool list = True);
-
- Vector<String> summary(Record& header, const String& doppler = "RADIO",
- const Bool list = True,
- const Bool pixelorder = True);
-
- Bool tofits(const String& outfile, const Bool velocity, const Bool
optical,
- const Int bitpix, const Double minpix, const Double maxpix,
- Record& region, const String& mask,
- const Bool overwrite = False,
- const Bool dropdeg = False, const Bool deglast = False,
- const Bool dropstokes = False, const Bool stokeslast =
False,
- const Bool wavelength = False);
-
- Bool toASCII(const String& outfile, Record& region, const String& mask,
- const String& sep = " ", const String& format = "%e",
- const Double maskvalue = -999, const Bool overwrite =
False);
-
-
- Vector<Double> topixel(Record& value);
-
- Record toworld(const Vector<double>& value, const String& format
= "n");
-
- Bool unlock();
-
- Bool detached();
-
- Record setregion(const Vector<Int>& blc, const Vector<Int>& trc,
- const String& infile = "");
-
- Record setboxregion(const Vector<Double>& blc, const Vector<Double>&
trc,
- const Bool frac = False, const String& infile
= "");
-
- //make test image...cube or 2d (default)
- bool maketestimage(const String& outfile="", const Bool
overwrite=False,
- const String& imagetype="2d");
-
- ImageInterface<Float> * newimage(const String& infile,
- const String& outfile,
- Record& region,
- const String& Mask,
- const bool dropdeg = False,
- const bool overwrite = False);
-
- ImageInterface<Float> * newimagefromfile(const String& fileName);
-
- ImageInterface<Float> * newimagefromarray(const String& outfile,
- Array<Float> & pixelsArray,
- const Record& csys,
- const Bool linear = False,
- const Bool overwrite = False,
- const Bool log = True);
-
- ImageInterface<Float> * newimagefromshape(const String& outfile,
- const Vector<Int>& shape,
- const Record& csys,
- const Bool linear = True,
- const Bool overwrite = False,
- const Bool log = True);
-
- ImageInterface<Float> * newimagefromfits(const String& outfile,
- const String& infile,
- const Int whichrep = 0,
- const Int whichhdu = 0,
- const Bool zeroblanks = False,
- const Bool overwrite = False);
-
- Record* echo(Record& v, const Bool godeep = False);
-
- //Functions to get you back a spectral profile at direction position
x, y.
- //x, y are to be in the world coord value or pixel value...user
specifies
- //by parameter xytype ("world" or "pixel").
- //On success returns true
- //return value of profile is in zyaxisval, zxaxisval contains the
spectral
- //values at which zyaxisval is evaluated its in the spectral type
- //specified by specaxis...possibilities
are "pixel", "freq", "radiovel", or "opticalvel"
- //(the code looks for the keywords "pixel", "freq", "vel", "optical",
and "radio"
- // in the string)
- // if "vel" is found but no "radio" or "optical", the full
relativistic velocity
- // is generated (MFrequency::RELATIVISTIC)
- // xunits determines the units of the x-axis values...default is "GHz"
for
- // freq and "km/s" for vel
- //PLEASE note that the returned value of zyaxisval are the units of
the image
- //specframe can be a valid frame from MFrequency...i.e LSRK, LSRD
etc...
- Bool getFreqProfile(const Vector<Double>& xy,
- Vector<Float>& zxaxisval, Vector<Float>& zyaxisval,
- const String& xytype="world",
- const String& specaxis="freq",
- const Int& whichStokes=0,
- const Int& whichTabular=0,
- const Int& whichLinear=0,
- const String& xunits="",
- const String& specframe="");
-
- //how about using this ?
- //for x.shape(xn) & y shape(yn)
- //if xn == yn == 1, single point
- //if xn == yn == 2, rectangle
- //if (xn == yn) > 2, polygon
- Bool getFreqProfile(const Vector<Double>& x,
- const Vector<Double>& y,
- Vector<Float>& zxaxisval, Vector<Float>& zyaxisval,
- const String& xytype="world",
- const String& specaxis="freq",
- const Int& whichStokes=0,
- const Int& whichTabular=0,
- const Int& whichLinear=0,
- const String& xunits="",
- const String& specframe="");
-
- // Return a record of the associates ImageInterface
- Bool toRecord(RecordInterface& rec);
- // Create a pagedimage if imagename is not "" else create a tempimage
- Bool fromRecord(const RecordInterface& rec, const String&
imagename="");
-
- // Deconvolve from beam
- // return True if the deconvolved source is a point source, False
otherwise
- // The returned value of successFit is True if the deconvolution could
be performed, False otherwise.
- casa::Bool
- deconvolveFromBeam(Quantity& majorFit,
- Quantity& minorFit,
- Quantity& paFit,
- Bool& successFit,
- const Vector<Quantity>& beam);
-
- // get the associated ImageInterface object
- const ImageInterface<Float>* getImage() const;
-
- private:
-
- ImageInterface<Float>* pImage_p;
- LogIO * itsLog;
-
-
-
- // Having private version of IS and IH means that they will
- // only recreate storage images if they have to
-
- ImageStatistics<casa::Float>* pStatistics_p;
- ImageHistograms<casa::Float>* pHistograms_p;
- //
- IPosition last_chunk_shape_p;
- ImageRegion* pOldStatsRegionRegion_p;
- casa::ImageRegion* pOldStatsMaskRegion_p;
- casa::ImageRegion* pOldHistRegionRegion_p;
- casa::ImageRegion* pOldHistMaskRegion_p;
- casa::Bool oldStatsStorageForce_p, oldHistStorageForce_p;
-
-
-
- // Center refpix apart from STokes
- void centreRefPix (casa::CoordinateSystem& cSys,
- const casa::IPosition& shape) const;
-
- // Convert types
- casa::ComponentType::Shape convertModelType (casa::Fit2D::Types
typeIn) const;
-
- // Delete private ImageStatistics and ImageHistograms objects
- bool deleteHistAndStats();
-
- // Convert error parameters from pixel to world and insert in
SkyComponent
- void encodeSkyComponentError (casa::LogIO& os,
- casa::SkyComponent& sky,
- casa::Double fluxRatio,
- const casa::ImageInterface<casa::Float>&
subIm,
- const casa::Vector<casa::Double>&
parameters,
- const casa::Vector<casa::Double>& errors,
- casa::Stokes::StokesTypes stokes,
- casa::Bool xIsLong) const;
-
- // Hanning smooth a vector
- void hanning_smooth (casa::Array<casa::Float>& out,
- casa::Array<casa::Bool>& maskOut,
- const casa::Vector<casa::Float>& in,
- const casa::Array<casa::Bool>& maskIn,
- casa::Bool isMasked) const;
-
-
-// Make a new image with given CS
- casa::Bool make_image(casa::String &error, const casa::String &image,
- const casa::CoordinateSystem& cSys,
- const casa::IPosition& shape,
- casa::LogIO& os, casa::Bool log=casa::True,
- casa::Bool overwrite=casa::False);
-
- // If file name empty make TempImage (allowTemp=T) or do nothing.
- // Otherwise, make a PagedImage from file name and copy mask and
- // misc from inimage. Returns T if image made, F if not
- casa::Bool
- makeExternalImage (casa::PtrHolder<casa::ImageInterface<casa::Float>
>& image,
- const casa::String& fileName,
- const casa::CoordinateSystem& cSys,
- const casa::IPosition& shape,
- const casa::ImageInterface<casa::Float>& inImage,
- casa::LogIO& os, casa::Bool overwrite=casa::False,
- casa::Bool allowTemp=casa::False,
- casa::Bool copyMask=casa::True) const;
-
- // Make a mask and define it in the image.
- casa::Bool makeMask(casa::ImageInterface<casa::Float>& out,
- casa::String& maskName,
- casa::Bool init, casa::Bool makeDefault,
- casa::LogIO& os, casa::Bool list) const;
- /*
- // Make a SubImage from a region and a WCLELMask string
- SubImage<Float> makeSubImage(
- ImageRegion*& pRegionRegion, ImageRegion*& pMaskRegion,
- ImageInterface<Float>& inImage, const Record& theRegion,
- const String& mask, LogIO *os, Bool writableIfPossible,
- const AxesSpecifier& axesSpecifier=casa::AxesSpecifier()
- );
-*/
-// See if the combination of the 'region' and 'mask' ImageRegions have
changed
- casa::Bool haveRegionsChanged (casa::ImageRegion* pNewRegionRegion,
- casa::ImageRegion* pNewMaskRegion,
- casa::ImageRegion* pOldRegionRegion,
- casa::ImageRegion* pOldMaskRegion)
const;
-
-// Convert a Record to a CoordinateSystem
- casa::CoordinateSystem*
- makeCoordinateSystem(const casa::Record& cSys,
- const casa::IPosition& shape) const;
-
- // Make a block of regions from a Record
- void makeRegionBlock(casa::PtrBlock<const casa::ImageRegion*>& regions,
- const casa::Record& Regions,
- casa::LogIO& logger);
-
- // Set the cache
- void set_cache(const casa::IPosition& chunk_shape) const;
-
- // Make an estimate for single parameter fit models
- casa::Vector<casa::Double>
- singleParameterEstimate (casa::Fit2D& fitter,
- casa::Fit2D::Types modelType,
- const casa::MaskedArray<casa::Float>&
pixels,
- casa::Float minVal, casa::Float maxVal,
- const casa::IPosition& minPos,
- const casa::IPosition& maxPos,
- casa::Stokes::StokesTypes stokes,
- const casa::ImageInterface<casa::Float>& im,
- casa::Bool xIsLong,
- casa::LogIO& os) const;
-
- // Prints an error message if the image DO is detached and returns
True.
- //bool detached() const;
-
- // Convert object-id's in the expression to LatticeExprNode objects.
- // It returns a string where the object-id's are placed by $n.
- // That string can be parsed by ImageExprParse.
- // Furthermore it fills the string exprName with the expression
- // where the object-id's are replaced by the image names.
- // Note that an image name can be an expression in itself, so
- // this string is not suitable for the ImageExprParse.
- //casa::String substituteOID (casa::Block<casa::LatticeExprNode>&
nodes,
- // casa::String& exprName,
- // const casa::String& expr) const;
-
-
- // Some helper functions that needs to be in casa namespace coordsys
-
- Record toWorldRecord (const Vector<Double>& pixel,
- const String& format) const;
-
- Record worldVectorToRecord (const Vector<Double>& world,
- Int c, const String& format,
- Bool isAbsolute, Bool showAsAbsolute)
const;
-
- Record worldVectorToMeasures(const Vector<Double>& world,
- Int c, Bool abs) const;
-
- void trim (Vector<Double>& inout,
- const Vector<Double>& replace) const;
-
- //return a vector of the spectral axis values in units requested
- //e.g "vel", "fre" or "pix"..specVal has to be sized already
- Bool getSpectralAxisVal(const String& specaxis, Vector<Float>& specVal,
- const CoordinateSystem& cSys, const String&
xunits,
- const String& freqFrame="");
-
-
- // Set the include and/or exclude pixel range for fitsky(). The
algorithm is:
- // if both includepix and excludepix are set, throw and exception
- // if neither is set and stokes=="I", all pixels >= 0 are included and
all < 0 are excluded
- // if neither is set and stokes != "I", all pixels <= 0 are included
and all pixels > 0 are excluded
- // if include/exclude has one element then -abs(element(0)) to
abs(element(0)) are included/excluded
- // if include/exclude has two elements, then the range betwenn those
elements is included/excluded
-
- void _setFitSkyIncludeExclude(
- const Vector<Float>& includepix, const Vector<Float>& excludepix,
- Fit2D& fitter
- ) const;
-
- void _fitskyExtractBeam(
- Vector<Double>& parameters, const ImageInfo& imageInfo,
- const Bool xIsLong, const CoordinateSystem& cSys
- ) const;
-
- // Writes a residual image and gets the stats from that image, and then
- // deletes the image if the user did not request it be saved.
- Record _fitskyWriteResidualAndGetStats(
- const SubImage<Float>& subImage, const Array<Float>& residPixels,
String residImageName
- ) const;
-
-};
-
-} // casac namespace
-#endif
-
=======================================
--- /trunk/images/Images/ImageCollapser.cc Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,409 +0,0 @@
-//# tSubImage.cc: Test program for class SubImage
-//# Copyright (C) 1998,1999,2000,2001,2003
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or (at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#include <images/Images/ImageCollapser.h>
-
-#include <casa/Arrays/ArrayMath.h>
-#include <casa/Containers/HashMap.h>
-#include <casa/Containers/HashMapIter.h>
-#include <casa/OS/Directory.h>
-#include <casa/OS/RegularFile.h>
-#include <casa/OS/SymLink.h>
-#include <images/Images/PagedImage.h>
-#include <images/Images/SubImage.h>
-#include <images/Images/TempImage.h>
-
-namespace casa {
-
- HashMap<uInt, String> *ImageCollapser::_funcNameMap = 0;
- HashMap<uInt, String> *ImageCollapser::_minMatchMap = 0;
- HashMap<uInt, Float (*)(const Array<Float>&)> *ImageCollapser::_funcMap =
0;
-
- ImageCollapser::ImageCollapser(
- String aggString, const String& imagename,
- const String& region, const String& box,
- const String& chanInp, const String& stokes,
- const String& maskInp, const uInt axis,
- const String& outname, const Bool overwrite
- ) : _log(new LogIO()), _image(0),
- _chan(chanInp), _stokesString(stokes), _mask(maskInp),
- _outname(outname),
- _overwrite(overwrite), _destructImage(True),
- _invertAxesSelection(False),
- _axes(Vector<uInt>(1, axis)), _aggType(UNKNOWN) {
- _construct(aggString, imagename, box, region);
- }
-
- ImageCollapser::ImageCollapser(
- String aggString, const String& imagename,
- const String& region, const String& box,
- const String& chanInp, const String& stokes,
- const String& maskInp, const Vector<uInt>& axes,
- const String& outname, const Bool overwrite
- ) : _log(new LogIO()), _image(0),
- _chan(chanInp), _stokesString(stokes), _mask(maskInp),
- _outname(outname),
- _overwrite(overwrite), _destructImage(True),
- _invertAxesSelection(False),
- _axes(axes), _aggType(UNKNOWN) {
- _construct(aggString, imagename, box, region);
- }
-
- ImageCollapser::ImageCollapser(
- String aggString, const ImageInterface<Float> * const image,
- const String& region, const String& box,
- const String& chanInp, const String& stokes,
- const String& maskInp, const uInt axis,
- const String& outname, const Bool overwrite
- ) : _log(new LogIO()), _image(image->cloneII()),
- _chan(chanInp), _stokesString(stokes), _mask(maskInp),
- _outname(outname),
- _overwrite(overwrite), _destructImage(True),
- _invertAxesSelection(False),
- _axes(Vector<uInt>(1, axis)), _aggType(UNKNOWN) {
- _construct(aggString, _image, box, region);
- }
-
- ImageCollapser::ImageCollapser(
- String aggString, const ImageInterface<Float> * const image,
- const String& region, const String& box,
- const String& chanInp, const String& stokes,
- const String& maskInp, const Vector<uInt> axes,
- const String& outname, const Bool overwrite
- ) : _log(new LogIO()), _image(image->cloneII()),
- _chan(chanInp), _stokesString(stokes), _mask(maskInp),
- _outname(outname),
- _overwrite(overwrite), _destructImage(True),
- _invertAxesSelection(False),
- _axes(axes), _aggType(UNKNOWN) {
- _construct(aggString, _image, box, region);
- }
-
- ImageCollapser::ImageCollapser(
- const ImageInterface<Float> * const image,
- const Vector<uInt>& axes, const Bool invertAxesSelection,
- const AggregateType aggregateType,
- const String& outname, const Bool overwrite
- ) : _log(new LogIO()), _image(image->cloneII()), _regionRecord(Record()),
- _mask(""), _outname(outname),
- _overwrite(overwrite), _destructImage(False),
- _invertAxesSelection(invertAxesSelection),
- _axes(axes), _aggType(aggregateType) {
- LogOrigin logOrigin("ImageCollapser", __FUNCTION__);
- *_log << logOrigin;
- if (_aggType == UNKNOWN) {
- *_log << "UNKNOWN aggregateType not allowed" << LogIO::EXCEPTION;
- }
- if (! _image) {
- *_log << "Cannot use a null image pointer with this constructor"
- << LogIO::EXCEPTION;
- }
- _invert();
- *_log << logOrigin;
- }
-
- ImageCollapser::~ImageCollapser() {
- delete _log;
- if (_destructImage) {
- delete _image;
- }
- }
-
- ImageInterface<Float>* ImageCollapser::collapse(const Bool wantReturn)
const {
- *_log << LogOrigin("ImageCollapser", __FUNCTION__);
- ImageRegion *imageRegion = 0;
- ImageRegion *maskRegion = 0;
- SubImage<Float> subImage = SubImage<Float>::createSubImage(
- imageRegion, maskRegion, *_image,
- _regionRecord, _mask, _log, False
- );
- delete imageRegion;
- delete maskRegion;
- IPosition inShape = subImage.shape();
- // Set the compressed axis reference pixel and reference value
- CoordinateSystem outCoords(subImage.coordinates());
- Vector<Double> blc, trc;
- IPosition pixblc(inShape.nelements(), 0);
- IPosition pixtrc = inShape - 1;
- if(
- ! outCoords.toWorld(blc, pixblc)
- || ! outCoords.toWorld(trc, pixtrc)
- ) {
- *_log << "Could not set new coordinate values" << LogIO::EXCEPTION;
- }
- Vector<Double> refValues = outCoords.referenceValue();
- Vector<Double> refPixels = outCoords.referencePixel();
- IPosition outShape = inShape;
- IPosition shape(outShape.nelements(), 1);
-
- for (Vector<uInt>::const_iterator iter=_axes.begin(); iter !=
_axes.end(); iter++) {
- uInt i = *iter;
- refValues[i] = (blc[i] + trc[i])/2;
- refPixels[i] = 0;
- outShape[i] = 1;
- shape[i] = inShape[i];
- }
-
- if (! outCoords.setReferenceValue(refValues)) {
- *_log << "Unable to set reference value" << LogIO::EXCEPTION;
- }
- if (! outCoords.setReferencePixel(refPixels)) {
- *_log << "Unable to set reference pixel" << LogIO::EXCEPTION;
- }
-
- ImageInterface<Float> *outImage = 0;
- if (_outname.empty()) {
- outImage = new TempImage<Float>(outShape, outCoords);
- }
- else {
- File out(_outname);
- if (out.exists()) {
- // remove file if it exists which prevents emission of
- // file is already open in table cache exceptions
- if (_overwrite) {
- if (out.isDirectory()) {
- Directory dir(_outname);
- dir.removeRecursive();
- }
- else if (out.isRegular()) {
- RegularFile reg(_outname);
- reg.remove();
- }
- else if (out.isSymLink()) {
- SymLink link(_outname);
- link.remove();
- }
- }
- else {
- // The only way this block can be entered is if a file by this name
- // has been written between the checking of inputs in the
constructor
- // call and the call of this method.
- *_log << "File " << _outname << " exists but overwrite is false so
it cannot be overwritten"
- << LogIO::EXCEPTION;
- }
- }
- outImage = new PagedImage<Float>(outShape, outCoords, _outname);
- }
- if (_aggType == ZERO) {
- Array<Float> zeros(outShape, 0.0);
- outImage->put(zeros);
- }
- else {
- Float (*function)(const Array<Float>&) = (*funcMap())(_aggType);
- Array<Float> data = subImage.get(False);
- for (uInt i=0; i<outShape.product(); i++) {
- IPosition start = toIPositionInArray(i, outShape);
- IPosition end = start + shape - 1;
- Slicer s(start, end, Slicer::endIsLast);
- outImage->putAt(function(data(s)), start);
- }
- }
- if (! _outname.empty()) {
- outImage->flush();
- }
- if (! wantReturn) {
- delete outImage;
- outImage = 0;
- }
- return outImage;
- }
-
- const HashMap<uInt, Float (*)(const Array<Float>&)>*
ImageCollapser::funcMap() {
- if (! _funcMap) {
- _funcMap = new HashMap<uInt, Float (*)(const
Array<Float>&)>(casa::mean<Float>);
- _funcMap->define((uInt)AVDEV, casa::avdev);
- _funcMap->define((uInt)MAX, casa::max);
- _funcMap->define((uInt)MEAN, casa::mean);
- _funcMap->define((uInt)MEDIAN, casa::median);
- _funcMap->define((uInt)MIN, casa::min);
- _funcMap->define((uInt)RMS, casa::rms);
- _funcMap->define((uInt)STDDEV, casa::stddev);
- _funcMap->define((uInt)SUM, casa::sum);
- _funcMap->define((uInt)VARIANCE, casa::variance);
- }
- return _funcMap;
- }
-
- const HashMap<uInt, String>* ImageCollapser::funcNameMap() {
- if (! _funcNameMap) {
- _funcNameMap = new HashMap<uInt, String>;
- _funcNameMap->define((uInt)AVDEV, "avdev");
- _funcNameMap->define((uInt)MAX, "max");
- _funcNameMap->define((uInt)MEAN, "mean");
- _funcNameMap->define((uInt)MEDIAN, "median");
- _funcNameMap->define((uInt)MIN, "min");
- _funcNameMap->define((uInt)RMS, "rms");
- _funcNameMap->define((uInt)STDDEV, "stddev");
- _funcNameMap->define((uInt)SUM, "sum");
- _funcNameMap->define((uInt)VARIANCE, "variance");
- _funcNameMap->define((uInt)ZERO, "zero");
- }
- return _funcNameMap;
- }
-
- const HashMap<uInt, String>* ImageCollapser::minMatchMap() {
- if (! _minMatchMap) {
- _minMatchMap = new HashMap<uInt, String>;
- _minMatchMap->define((uInt)AVDEV, "a");
- _minMatchMap->define((uInt)MAX, "ma");
- _minMatchMap->define((uInt)MEAN, "mea");
- _minMatchMap->define((uInt)MEDIAN, "med");
- _minMatchMap->define((uInt)MIN, "mi");
- _minMatchMap->define((uInt)RMS, "r");
- _minMatchMap->define((uInt)STDDEV, "st");
- _minMatchMap->define((uInt)SUM, "su");
- _minMatchMap->define((uInt)VARIANCE, "v");
- _minMatchMap->define((uInt)ZERO, "z");
- }
- return _minMatchMap;
- }
-
- Vector<ImageInputProcessor::OutputStruct>
ImageCollapser::_getOutputStruct() {
- Vector<ImageInputProcessor::OutputStruct> outputs(0);
- _outname.trim();
- if (! _outname.empty()) {
- ImageInputProcessor::OutputStruct outputImage;
- outputImage.label = "output image";
- outputImage.outputFile = &_outname;
- outputImage.required = True;
- outputImage.replaceable = _overwrite;
- outputs.resize(1);
- outputs[0] = outputImage;
- }
- return outputs;
- }
-
- void ImageCollapser::_finishConstruction() {
- for (
- Vector<uInt>::const_iterator iter=_axes.begin();
- iter != _axes.end(); iter++
- ) {
- if (*iter >= _image->ndim()) {
- *_log << "Specified zero-based axis (" << *iter
- << ") must be less than the number of axes in " <<
_image->name()
- << "(" << _image->ndim() << LogIO::EXCEPTION;
- }
- }
- _invert();
- }
-
- void ImageCollapser::_construct(
- String& aggString, const String& imagename,
- const String& box, const String& regionName
- ) {
- LogOrigin logOrigin("ImageCollapser", __FUNCTION__);
- _setAggregateType(aggString);
- *_log << logOrigin;
-
- String diagnostics;
- Vector<ImageInputProcessor::OutputStruct> outputs =
_getOutputStruct();
- Vector<ImageInputProcessor::OutputStruct> *outputPtr =
outputs.size() > 0
- ? &outputs
- : 0;
- ImageInputProcessor inputProcessor;
- inputProcessor.process(
- _image, _regionRecord, diagnostics,
- outputPtr, _stokesString, imagename,
- 0, regionName, box, _chan,
- RegionManager::USE_ALL_STOKES,
- False, 0
- );
- _finishConstruction();
- }
-
- void ImageCollapser::_construct(
- String& aggString, const ImageInterface<Float> *image,
- const String& box, const String& regionName
- ) {
- LogOrigin logOrigin("ImageCollapser", __FUNCTION__);
- _setAggregateType(aggString);
- *_log << logOrigin;
-
- String diagnostics;
- Vector<ImageInputProcessor::OutputStruct> outputs =
_getOutputStruct();
- Vector<ImageInputProcessor::OutputStruct> *outputPtr =
outputs.size() > 0
- ? &outputs
- : 0;
- ImageInputProcessor inputProcessor;
- inputProcessor.process(
- _regionRecord, diagnostics,
- outputPtr, _stokesString, image,
- 0, regionName, box, _chan,
- RegionManager::USE_ALL_STOKES,
- False, 0
- );
- _finishConstruction();
- }
-
- void ImageCollapser::_setAggregateType(String& aggString) {
- LogOrigin logOrigin("ImageCollapser", __FUNCTION__);
- *_log << logOrigin;
- if (aggString.empty()) {
- *_log << "Aggregate function name is not specified and it must
be."
- << LogIO::EXCEPTION;
- }
- aggString.downcase();
- const HashMap<uInt, String> *funcNamePtr = funcNameMap();
- ConstHashMapIter<uInt, String> iter = *minMatchMap();
- for (iter.toStart(); ! iter.atEnd(); iter++) {
- uInt key = iter.getKey();
- String minMatch = iter.getVal();
- String funcName = (*funcNamePtr)(key);
- if (
- aggString.startsWith(minMatch)
- && funcName.startsWith(aggString)
- ) {
- _aggType = (AggregateType)key;
- return;
- }
- }
- *_log << "Unknown aggregate function specified by " << aggString <<
LogIO::EXCEPTION;
- }
-
- void ImageCollapser::_invert() {
- if (_invertAxesSelection) {
- Vector<uInt> newAxes(_image->ndim() - _axes.size(), 0);
- uInt index=0;
- for (uInt i=0; i<_image->ndim(); i++) {
- Bool found = False;
- for (uInt j=0; j<_axes.size(); j++) {
- if (i == _axes[j]) {
- found = True;
- break;
- }
- }
- if (! found) {
- newAxes[index] = i;
- index++;
- }
- }
- _axes.resize(newAxes.size());
- _axes = newAxes;
- }
- }
-}
-
=======================================
--- /trunk/images/Images/ImageCollapser.h Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,180 +0,0 @@
-//# tSubImage.cc: Test program for class SubImage
-//# Copyright (C) 1998,1999,2000,2001,2003
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or (at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#ifndef IMAGES_IMAGECOLLAPSER_H
-#define IMAGES_IMAGECOLLAPSER_H
-
-#include <casa/Containers/HashMap.h>
-#include <casa/Logging/LogIO.h>
-#include <images/Images/ImageInputProcessor.h>
-#include <images/Images/ImageInterface.h>
-
-#include <casa/namespace.h>
-
-namespace casa {
-
- class ImageCollapser {
- // <summary>
- // Top level interface which allows collapsing of images along a
single axis. An aggregate method
- // (average, sum, etc) is applied to the collapsed pixels.
- // </summary>
-
- // <reviewed reviewer="" date="" tests="" demos="">
- // </reviewed>
-
- // <prerequisite>
- // </prerequisite>
-
- // <etymology>
- // Collapses image.
- // </etymology>
-
- // <synopsis>
- // High level interface for collapsing an image along a single
axis.
- // </synopsis>
-
- // <example>
- // <srcblock>
- // ImageCollapser collapser();
- // collapser.collapse();
- // </srcblock>
- // </example>
-
- public:
-
- enum AggregateType {
- AVDEV,
- MAX,
- MEAN,
- MEDIAN,
- MIN,
- RMS,
- STDDEV,
- SUM,
- VARIANCE,
- // set all pixels in output image to 0
- ZERO,
- UNKNOWN
- };
-
- // if <src>outname</src> is empty, no image will be written
- // if <src>overwrite</src> is True, if image already exists it will
be removed
- // if <src>overwrite</src> is False, if image already exists
exception will be thrown
- //
- // <group>
- ImageCollapser(
- String aggString,const String& imagename,
- const String& region, const String& box,
- const String& chanInp, const String& stokes,
- const String& maskInp, const uInt axis,
- const String& outname, const Bool overwrite
- );
-
- ImageCollapser(
- String aggString,const String& imagename,
- const String& region, const String& box,
- const String& chanInp, const String& stokes,
- const String& maskInp, const Vector<uInt>& axes,
- const String& outname, const Bool overwrite
- );
-
- ImageCollapser(
- String aggString, const ImageInterface<Float> * const image,
- const String& region, const String& box,
- const String& chanInp, const String& stokes,
- const String& maskInp, const uInt axis,
- const String& outname, const Bool overwrite
- );
-
- ImageCollapser(
- String aggString, const ImageInterface<Float> * const image,
- const String& region, const String& box,
- const String& chanInp, const String& stokes,
- const String& maskInp, const Vector<uInt> axes,
- const String& outname, const Bool overwrite
- );
-
- ImageCollapser(
- const ImageInterface<Float> * const image,
- const Vector<uInt>& axes, const Bool invertAxesSelection,
- const AggregateType aggregateType,
- const String& outname, const Bool overwrite
- );
- // </group>
-
- // destructor
- ~ImageCollapser();
-
- // perform the collapse. If <src>wantReturn</src> is True, return
a pointer to the
- // collapsed image. The returned pointer is created via new(); it
is the caller's
- // responsibility to delete the returned pointer. If
<src>wantReturn</src> is False,
- // a NULL pointer is returned and pointer deletion is performed
internally.
- ImageInterface<Float>* collapse(const Bool wantReturn) const;
-
- static const HashMap<uInt, Float (*)(const Array<Float>&)>*
funcMap();
- static const HashMap<uInt, String>* funcNameMap();
- static const HashMap<uInt, String>* minMatchMap();
-
- private:
- LogIO *_log;
- ImageInterface<Float> * _image;
- Record _regionRecord;
- String _chan, _stokesString, _mask, _outname;
- Bool _overwrite, _destructImage, _invertAxesSelection;
- Vector<uInt> _axes;
- AggregateType _aggType;
-
- static HashMap<uInt, Float (*)(const Array<Float>&)> *_funcMap;
- static HashMap<uInt, String> *_funcNameMap;
- static HashMap<uInt, String> *_minMatchMap;
-
- // disallow default constructor
- ImageCollapser();
-
- // does the lion's share of constructing the object, ie checks
validity of
- // inputs, etc.
- void _construct(
- String& aggString, const String& imagename,
- const String& box, const String& regionName
- );
-
- void _construct(
- String& aggString, const ImageInterface<Float> *image,
- const String& box, const String& regionName
- );
-
- void _setAggregateType(String& aggString);
-
- void _invert();
-
- Vector<ImageInputProcessor::OutputStruct> _getOutputStruct();
-
- void _finishConstruction();
- };
-}
-
-#endif
=======================================
--- /trunk/images/Images/ImageFitter.cc Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,966 +0,0 @@
-//# tSubImage.cc: Test program for class SubImage
-//# Copyright (C) 1998,1999,2000,2001,2003
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or (at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#include <casa/Containers/HashMap.h>
-#include <casa/Containers/HashMapIter.h>
-
-#include <casa/IO/FilebufIO.h>
-#include <casa/IO/FiledesIO.h>
-
-#include <casa/Quanta/Quantum.h>
-#include <casa/Quanta/Unit.h>
-#include <casa/Quanta/UnitMap.h>
-#include <casa/Quanta/MVAngle.h>
-#include <casa/Quanta/MVTime.h>
-#include <casa/OS/Time.h>
-#include <casa/Utilities/Precision.h>
-
-#include <components/ComponentModels/Flux.h>
-#include <components/ComponentModels/SpectralModel.h>
-
-#include <images/IO/FitterEstimatesFileParser.h>
-#include <images/Images/ImageAnalysis.h>
-#include <images/Images/ImageFitter.h>
-#include <images/Images/ImageMetaData.h>
-#include <images/Images/ImageStatistics.h>
-#include <images/Images/FITSImage.h>
-#include <images/Images/MIRIADImage.h>
-#include <images/Regions/RegionManager.h>
-
-#include <images/Regions/WCUnion.h>
-#include <images/Regions/WCBox.h>
-#include <lattices/Lattices/LCRegion.h>
-
-#include <components/ComponentModels/SkyComponent.h>
-#include <components/ComponentModels/ComponentShape.h>
-
-#include <components/ComponentModels/GaussianShape.h>
-
-namespace casa {
-
- ImageFitter::ImageFitter(
- const String& imagename, const String& region, const String& box,
- const uInt chanInp, const String& stokes,
- const String& maskInp, const Vector<Float>& includepix,
- const Vector<Float>& excludepix, const String& residualInp,
- const String& modelInp, const String& estimatesFilename,
- const String& logfile, const Bool& append,
- const String& newEstimatesInp, const String& compListName,
- const CompListWriteControl writeControl
- ) : _log(new LogIO()), _image(0), _chan(chanInp),
_stokesString(stokes),
- _mask(maskInp), _residual(residualInp),_model(modelInp),
_logfileName(logfile),
- regionString(""), estimatesString(""),
_newEstimatesFileName(newEstimatesInp),
- _compListName(compListName),
- includePixelRange(includepix), excludePixelRange(excludepix),
- estimates(), fixed(0), logfileAppend(append), _fitConverged(False),
- fitDone(False), _noBeam(False), _deleteImageOnDestruct(True),
_peakIntensities(),
- _writeControl(writeControl) {
- *_log << LogOrigin("ImageFitter", __FUNCTION__);
- *_log << LogIO::WARN << "This constructor has been deprecated and
will be removed "
- << "in an upcoming version. Please modify existing code not to
use another constructor"
- << LogIO::POST;
- _construct(imagename, box, region, 0, estimatesFilename);
- }
-
- ImageFitter::ImageFitter(
- const String& imagename, const Record* regionPtr, const String&
box,
- const uInt chanInp, const String& stokes,
- const String& maskInp, const Vector<Float>& includepix,
- const Vector<Float>& excludepix, const String& residualInp,
- const String& modelInp, const String& estimatesFilename,
- const String& logfile, const Bool& append,
- const String& newEstimatesInp, const String& compListName,
- const CompListWriteControl writeControl
- ) : _log(new LogIO()), _image(0), _chan(chanInp),
_stokesString(stokes),
- _mask(maskInp), _residual(residualInp),_model(modelInp),
_logfileName(logfile),
- regionString(""), estimatesString(""),
_newEstimatesFileName(newEstimatesInp),
- _compListName(compListName),
- includePixelRange(includepix), excludePixelRange(excludepix),
- estimates(), fixed(0), logfileAppend(append), _fitConverged(False),
- fitDone(False), _noBeam(False), _deleteImageOnDestruct(True),
_peakIntensities(),
- _writeControl(writeControl) {
- *_log << LogOrigin("ImageFitter", __FUNCTION__);
- *_log << LogIO::WARN << "This constructor has been deprecated and
will be removed "
- << "in an upcoming version. Please modify existing code not to
use another constructor"
- << LogIO::POST;
- _construct(imagename, box, "", regionPtr, estimatesFilename);
- }
-
- ImageFitter::ImageFitter(
- const ImageInterface<Float>*& image, const String& region, const
String& box,
- const uInt chanInp, const String& stokes,
- const String& maskInp, const Vector<Float>& includepix,
- const Vector<Float>& excludepix, const String& residualInp,
- const String& modelInp, const String& estimatesFilename,
- const String& logfile, const Bool& append,
- const String& newEstimatesInp, const String& compListName,
- const CompListWriteControl writeControl
- ) : _log(new LogIO()), _image(image->cloneII()), _chan(chanInp),
_stokesString(stokes),
- _mask(maskInp), _residual(residualInp),_model(modelInp),
_logfileName(logfile),
- regionString(""), estimatesString(""),
_newEstimatesFileName(newEstimatesInp),
- _compListName(compListName),
- includePixelRange(includepix), excludePixelRange(excludepix),
- estimates(), fixed(0), logfileAppend(append), _fitConverged(False),
- fitDone(False), _noBeam(False), _deleteImageOnDestruct(False),
_peakIntensities(),
- _writeControl(writeControl) {
- _construct(_image, box, region, 0, estimatesFilename);
- }
-
- ImageFitter::ImageFitter(
- const ImageInterface<Float>*& image, const Record* regionPtr,
const String& box,
- const uInt chanInp, const String& stokes,
- const String& maskInp, const Vector<Float>& includepix,
- const Vector<Float>& excludepix, const String& residualInp,
- const String& modelInp, const String& estimatesFilename,
- const String& logfile, const Bool& append,
- const String& newEstimatesInp, const String& compListName,
- const CompListWriteControl writeControl
- ) : _log(new LogIO()), _image(image->cloneII()), _chan(chanInp),
_stokesString(stokes),
- _mask(maskInp), _residual(residualInp),_model(modelInp),
_logfileName(logfile),
- regionString(""), estimatesString(""),
_newEstimatesFileName(newEstimatesInp),
- _compListName(compListName),
- includePixelRange(includepix), excludePixelRange(excludepix),
- estimates(), fixed(0), logfileAppend(append), _fitConverged(False),
- fitDone(False), _noBeam(False), _deleteImageOnDestruct(False),
_peakIntensities(),
- _writeControl(writeControl) {
- _construct(_image, box, "", regionPtr, estimatesFilename);
- }
-
-
- ImageFitter::~ImageFitter() {
- delete _log;
- if (_deleteImageOnDestruct) {
- delete _image;
- }
- }
-
- ComponentList ImageFitter::fit() {
- *_log << LogOrigin("ImageFitter", __FUNCTION__);
- Array<Float> residPixels;
- Array<Bool> residMask;
- Bool converged;
- uInt ngauss = estimates.nelements() > 0 ? estimates.nelements() :
1;
- Vector<String> models(ngauss);
- models.set("gaussian");
- Bool fit = True;
- Bool deconvolve = False;
- Bool list = True;
- String errmsg;
- Record estimatesRecord;
-
- estimates.toRecord(errmsg, estimatesRecord);
- ImageAnalysis myImage(_image);
- Vector<String> allowFluxUnits(1);
- allowFluxUnits[0] = "Jy.km/s";
- FluxRep<Double>::setAllowedUnits(allowFluxUnits);
-
- try {
- _results = myImage.fitsky(
- residPixels, residMask, converged,
- inputStats, _residStats, chiSquared,
- _regionRecord, _chan, _stokesString,
- _mask, models, estimatesRecord, fixed,
- includePixelRange, excludePixelRange,
- fit, deconvolve, list, _residual, _model
- );
- }
- catch (AipsError err) {
- *_log << LogIO::WARN << "Fit failed to converge because of
exception: "
- << err.getMesg() << LogIO::POST;
- converged = false;
- }
- FluxRep<Double>::clearAllowedUnits();
-
- fitDone = True;
- _fitConverged = converged;
- if(converged) {
- _setFluxes();
- _setSizes();
- _writeCompList();
- }
- String resultsString = _resultsToString();
- if (converged && ! _newEstimatesFileName.empty()) {
- _writeNewEstimatesFile();
- }
- *_log << LogIO::NORMAL << resultsString << LogIO::POST;
- if (! _logfileName.empty()) {
- _writeLogfile(resultsString);
- }
- return _results;
- }
-
- Bool ImageFitter::converged() const {
- if (!fitDone) {
- throw AipsError("fit has not yet been performed");
- }
- return _fitConverged;
- }
-
- void ImageFitter::_getStandardDeviations(Double& inputStdDev, Double&
residStdDev) const {
- inputStdDev = _getStatistic("sigma", inputStats);
- residStdDev = _getStatistic("sigma", _residStats);
- }
-
- void ImageFitter::_getRMSs(Double& inputRMS, Double& residRMS) const {
- inputRMS = _getStatistic("rms", inputStats);
- residRMS = _getStatistic("rms", _residStats);
- }
-
- Double ImageFitter::_getStatistic(const String& type, const Record&
stats) const {
- Vector<Double> statVec;
- stats.get(stats.fieldNumber(type), statVec);
- return statVec[0];
- }
-
- Vector<ImageInputProcessor::OutputStruct> ImageFitter::_getOutputs() {
- LogOrigin logOrigin("ImageFitter", __FUNCTION__);
- *_log << logOrigin;
-
- ImageInputProcessor::OutputStruct residualIm;
- residualIm.label = "residual image";
- residualIm.outputFile = &_residual;
- residualIm.required = False;
- residualIm.replaceable = True;
- ImageInputProcessor::OutputStruct modelIm;
- modelIm.label = "model image";
- modelIm.outputFile = &_model;
- modelIm.required = False;
- modelIm.replaceable = True;
- ImageInputProcessor::OutputStruct newEstFile;
- newEstFile.label = "new estiamtes file";
- newEstFile.outputFile = &_newEstimatesFileName;
- newEstFile.required = False;
- newEstFile.replaceable = True;
- ImageInputProcessor::OutputStruct logFile;
- logFile.label = "log file";
- logFile.outputFile = &_logfileName;
- logFile.required = False;
- logFile.replaceable = True;
-
- Vector<ImageInputProcessor::OutputStruct> outputs(4);
- outputs[0] = residualIm;
- outputs[1] = modelIm;
- outputs[2] = newEstFile;
- outputs[3] = logFile;
-
- return outputs;
- // Vector<Coordinate::Type> reqCoordTypes(1,
Coordinate::DIRECTION);
- //reqCoordTypes[0] = Coordinate::DIRECTION;
- // return inputProcessor;
- }
-
- void ImageFitter::_construct(
- const String& imagename, const String& box, const String&
regionName,
- const Record* regionPtr, const String& estimatesFilename
- ) {
- LogOrigin logOrigin("ImageFitter", __FUNCTION__);
- *_log << logOrigin;
- *_log << LogIO::WARN
- << "THIS CONSTRUCTOR IS DEPRECATED. PLEASE USE A CONSTRUCTOR
WHICH "
- << "TAKES A VALID IMAGE INTERFACE POINTER" << LogIO::POST;
- Vector<ImageInputProcessor::OutputStruct> outputs = _getOutputs();
- //ImageInputProcessor inputProcessor = _beginConstruction();
-
- String diagnostics;
- Vector<Coordinate::Type> reqCoordTypes(1, Coordinate::DIRECTION);
- ImageInputProcessor inputProcessor;
-
- inputProcessor.process(
- _image, _regionRecord, diagnostics, &outputs,
- _stokesString, imagename, regionPtr,
- regionName, box, String::toString(_chan),
- RegionManager::USE_FIRST_STOKES, False,
- &reqCoordTypes
- );
- _finishConstruction(estimatesFilename);
- }
-
- void ImageFitter::_construct(
- const ImageInterface<Float> *image, const String& box, const
String& regionName,
- const Record* regionPtr, const String& estimatesFilename
- ) {
- LogOrigin logOrigin("ImageFitter", __FUNCTION__);
- *_log << logOrigin;
- Vector<ImageInputProcessor::OutputStruct> outputs = _getOutputs();
- //ImageInputProcessor inputProcessor = _beginConstruction();
-
- String diagnostics;
- Vector<Coordinate::Type> reqCoordTypes(1, Coordinate::DIRECTION);
- ImageInputProcessor inputProcessor;
-
- inputProcessor.process(
- _regionRecord, diagnostics, &outputs,
- _stokesString, image, regionPtr,
- regionName, box, String::toString(_chan),
- RegionManager::USE_FIRST_STOKES, False,
- &reqCoordTypes
- );
- _finishConstruction(estimatesFilename);
- }
-
- void ImageFitter::_finishConstruction(const String& estimatesFilename)
{
- LogOrigin logOrigin("ImageFitter", __FUNCTION__);
- *_log << logOrigin;
- // <todo> kludge because Flux class is really only made for I, Q,
U, and V stokes
-
- String iquv = "IQUV";
- _kludgedStokes = (iquv.index(_stokesString) == String::npos) ||
_stokesString.empty()
- ? "I" : _stokesString;
- // </todo>
- //_doRegion(box, regionName, regionPtr);
- if(estimatesFilename.empty()) {
- *_log << LogIO::NORMAL << "No estimates file specified, so will
attempt to find and fit one gaussian."
- << LogIO::POST;
- }
- else {
- FitterEstimatesFileParser parser(estimatesFilename, *_image);
- estimates = parser.getEstimates();
- estimatesString = parser.getContents();
- fixed = parser.getFixed();
- Record rec;
- String errmsg;
- estimates.toRecord(errmsg, rec);
- *_log << LogIO::NORMAL << "File " << estimatesFilename << " has "
<< estimates.nelements()
- << " specified, so will attempt to fit that many gaussians " <<
LogIO::POST;
- }
- }
-
- String ImageFitter::_resultsToString() {
- ostringstream summary;
- summary << "****** Fit performed at " << Time().toString()
<< "******" << endl << endl;
- summary << "Input parameters ---" << endl;
- summary << " --- imagename: " << _image->name() <<
endl;
- summary << " --- region: " << regionString << endl;
- summary << " --- channel: " << _chan << endl;
- summary << " --- stokes: " << _stokesString <<
endl;
- summary << " --- mask: " << _mask << endl;
- summary << " --- include pixel ragne: " << includePixelRange <<
endl;
- summary << " --- exclude pixel ragne: " << excludePixelRange <<
endl;
- summary << " --- initial estimates: " << estimatesString <<
endl;
-
- if (converged()) {
- if (_noBeam) {
- *_log << LogIO::WARN << "Flux density not reported because "
- << "there is no clean beam in image header so these quantities
cannot "
- << "be calculated" << LogIO::POST;
- }
- summary << _statisticsToString() << endl;
- for (uInt i = 0; i < _results.nelements(); i++) {
- summary << "Fit on " << _image->name(True) << " component " << i <<
endl;
- summary <<
_results.component(i).positionToString(&(_image->coordinates())) << endl;
- summary << _sizeToString(i) << endl;
- summary << _fluxToString(i) << endl;
- summary << _spectrumToString(i) << endl;
- }
- }
- else {
- summary << "*** FIT FAILED ***" << endl;
- }
- return summary.str();
- }
-
- String ImageFitter::_statisticsToString() const {
- ostringstream stats;
- // TODO It is not clear how this chi squared value is calculated and
atm it does not
- // appear to be useful, so don't report it. In the future,
investigate more deeply
- // how it is calculated and see if a useful value for reporting can
be derived from
- // it.
- // stats << " --- Chi-squared of fit " << chiSquared << endl;
- stats << "Input and residual image statistics (to be used as a rough
guide only as to goodness of fit)" << endl;
- Double inputStdDev, residStdDev, inputRMS, residRMS;
- _getStandardDeviations(inputStdDev, residStdDev);
- _getRMSs(inputRMS, residRMS);
- String unit = _fluxDensities[0].getUnit();
- stats << " --- Standard deviation of input image " <<
inputStdDev << " " << unit << endl;
- stats << " --- Standard deviation of residual image " <<
residStdDev << " " << unit << endl;
- stats << " --- RMS of input image " << inputRMS << " " << unit
<< endl;
- stats << " --- RMS of residual image " << residRMS << " " <<
unit << endl;
- return stats.str();
- }
-
- void ImageFitter::_setFluxes() {
- uInt ncomps = _results.nelements();
-
- _fluxDensities.resize(ncomps);
- _fluxDensityErrors.resize(ncomps);
-
- _peakIntensities.resize(ncomps);
- _peakIntensityErrors.resize(ncomps);
-
- _majorAxes.resize(ncomps);
- _minorAxes.resize(ncomps);
- Vector<Quantity> fluxQuant;
-
- ImageAnalysis ia;
- ia.open(_image->name());
-
- Double rmsPeak = Vector<Double>(_residStats.asArrayDouble("rms"))[0];
- Quantity rmsPeakError(
- rmsPeak,
- _image->units()
- );
-
- ImageMetaData md(*_image);
- Quantity resArea;
- Bool found = False;
- Quantity intensityToFluxConversion =
_image->units().getName().contains("/beam")
- ? Quantity(1.0, "beam")
- : Quantity(1.0, "pixel");
-
- if (intensityToFluxConversion.getUnit() == "beam") {
- if(md.getBeamArea(resArea)) {
- found = True;
- }
- else {
- *_log << LogIO::WARN
- << "Image units are per beam but beam area could not "
- << "be determined. Assume beam area is pixel area."
- << LogIO::POST;
- }
- }
- if (! found) {
- if (! md.getDirectionPixelArea(resArea)) {
- *_log << LogIO::EXCEPTION
- << "Pixel area could not be determined";
- }
- }
-
- uInt polNum = 0;
- for(uInt i=0; i<ncomps; i++) {
- _results.getFlux(fluxQuant, i);
- // TODO there is probably a better way to get the flux component we
want...
- Vector<String> polarization = _results.getStokes(i);
- for (uInt j=0; j<polarization.size(); j++) {
- if (polarization[j] == _kludgedStokes) {
- _fluxDensities[i] = fluxQuant[j];
- std::complex<double> error = _results.component(i).flux().errors()[j];
- _fluxDensityErrors[i].setValue(sqrt(error.real()*error.real() +
error.imag()*error.imag()));
- _fluxDensityErrors[i].setUnit(_fluxDensities[i].getUnit());
- polNum = j;
- break;
- }
- }
- const ComponentShape* compShape = _results.getShape(i);
-
- AlwaysAssert(compShape->type() == ComponentType::GAUSSIAN,
AipsError);
- _majorAxes[i] = (static_cast<const GaussianShape
*>(compShape))->majorAxis();
- _minorAxes[i] = (static_cast<const GaussianShape
*>(compShape))->minorAxis();
- _peakIntensities[i] = ia.convertflux(
- _noBeam, _fluxDensities[i], _majorAxes[i],
_minorAxes[i], "Gaussian", True, True
- );
-
- rmsPeakError.convert(_peakIntensities[i].getUnit());
- Double rmsPeakErrorValue = rmsPeakError.getValue();
- Double peakErrorFromFluxErrorValue = (
- _peakIntensities[i]*_fluxDensityErrors[i]/_fluxDensities[i]
- ).getValue();
-
- _peakIntensityErrors[i].setValue(
- max(
- rmsPeakErrorValue,
- peakErrorFromFluxErrorValue
- )
- );
- _peakIntensityErrors[i].setUnit(_image->units());
- if (rmsPeakErrorValue > peakErrorFromFluxErrorValue) {
- const GaussianShape *gaussShape = static_cast<const
GaussianShape *>(compShape);
- Quantity compArea = gaussShape->getArea();
- Quantity rmsFluxError = rmsPeakError*compArea/resArea;
- rmsFluxError.convert(_fluxDensityErrors[i].getUnit());
- _fluxDensityErrors[i].setValue(
- max(
- _fluxDensityErrors[i].getValue(),
- rmsFluxError.getValue()
- )
- );
- Vector<std::complex<double> > errors(4,
std::complex<double>(0, 0));
- errors[polNum] =
std::complex<double>(_fluxDensityErrors[i].getValue(), 0);
- _results.component(i).flux().setErrors(errors);
-
- }
- }
- }
-
- void ImageFitter::_setSizes() {
- uInt ncomps = _results.nelements();
-
- _positionAngles.resize(ncomps);
- _majorAxisErrors.resize(ncomps);
- _minorAxisErrors.resize(ncomps);
- _positionAngleErrors.resize(ncomps);
- Double rmsPeak = Vector<Double>(_residStats.asArrayDouble("rms"))[0];
-
- Quantity rmsPeakError(
- rmsPeak,
- _image->units()
- );
-
- Vector<Quantity> beam = _image->imageInfo().restoringBeam();
- Bool hasBeam = beam.nelements() == 3;
-
- Quantity xBeam;
- Quantity yBeam;
- Quantity paBeam;
-
- if (hasBeam) {
- xBeam = beam[0];
- yBeam = beam[1];
- paBeam = beam[2];
- }
- else {
- ImageMetaData md(*_image);
- Int dirCoordNumber = md.directionCoordinateNumber();
- Vector<Double> pixInc =
_image->coordinates().directionCoordinate(dirCoordNumber).increment();
- xBeam = Quantity(pixInc[0], "rad");
- yBeam = Quantity(pixInc[1], "rad");
- paBeam = Quantity(0, "rad");
- }
-
- for(uInt i=0; i<_results.nelements(); i++) {
- const ComponentShape* compShape = _results.getShape(i);
- AlwaysAssert(compShape->type() == ComponentType::GAUSSIAN,
AipsError);
- _positionAngles[i] = (static_cast<const GaussianShape
*>(compShape))->positionAngle();
- _majorAxisErrors[i] = (static_cast<const GaussianShape
*>(compShape))->majorAxisError();
- _minorAxisErrors[i] = (static_cast<const GaussianShape
*>(compShape))->minorAxisError();
- _positionAngleErrors[i] = (static_cast<const GaussianShape
*>(compShape))->positionAngleError();
- Double signalToNoise =
fabs((_peakIntensities[i]/rmsPeakError).getValue());
-
- Quantity paRelToBeam = _positionAngles[i] - paBeam;
- paRelToBeam.convert("rad");
-
- xBeam.convert(_majorAxisErrors[i].getUnit());
- yBeam.convert(_majorAxisErrors[i].getUnit());
- Double xBeamVal = xBeam.getValue();
- Double yBeamVal = yBeam.getValue();
-
- Double cosPA = cos(paRelToBeam.getValue());
- Double sinPA = sin(paRelToBeam.getValue());
-
- // angles are measured from north (y direction).
- _majorAxisErrors[i].setValue(
- max(
- _majorAxisErrors[i].getValue(),
- sqrt(
- (
- xBeamVal*sinPA * xBeamVal*sinPA
- + yBeamVal*cosPA * yBeamVal*cosPA
- )/signalToNoise
- )
- )
- );
- _minorAxisErrors[i].setValue(
- max(
- _minorAxisErrors[i].getValue(),
- sqrt(
- (
- xBeamVal*cosPA * xBeamVal*cosPA
- + yBeamVal*sinPA * yBeamVal*sinPA
- )/signalToNoise
- )
- )
- );
-
- Double posAngleRad = _positionAngles[i].getValue(Unit("rad"));
- Quantity posAngErrorFromSN = _positionAngles[i] * sqrt(
- _majorAxisErrors[i]/_majorAxes[i] *
_majorAxisErrors[i]/_majorAxes[i]
- + _minorAxisErrors[i]/_minorAxes[i] *
_minorAxisErrors[i]/_minorAxes[i]
- );
- posAngErrorFromSN *= 1/(1 + posAngleRad*posAngleRad);
- posAngErrorFromSN.convert(_positionAngleErrors[i].getUnit());
- _positionAngleErrors[i].setValue(
- max(_positionAngleErrors[i].getValue(),
posAngErrorFromSN.getValue())
- );
-
- _majorAxisErrors[i].convert(_majorAxes[i].getUnit());
- _minorAxisErrors[i].convert(_minorAxes[i].getUnit());
- _positionAngleErrors[i].convert(_positionAngles[i].getUnit());
-
- GaussianShape* newShape = dynamic_cast<GaussianShape
*>(compShape->clone());
-
- newShape->setErrors(
- _majorAxisErrors[i], _minorAxisErrors[i],
- _positionAngleErrors[i]
- );
-
- // newShape->setErrors(newErrors);
-
- // set the position uncertainties
-
- Quantity latError = compShape->refDirectionErrorLat();
- Quantity longError = compShape->refDirectionErrorLong();
-
- paBeam.convert("rad");
- Double cosPaBeam = cos(paBeam.getValue());
- Double sinPaBeam = sin(paBeam.getValue());
-
- Quantity longErrorFromSN = sqrt(
- xBeam*sinPaBeam*xBeam*sinPaBeam +
yBeam*cosPaBeam*yBeam*cosPaBeam
- )/(2*signalToNoise);
- Quantity latErrorFromSN = sqrt(
- xBeam*cosPaBeam*xBeam*cosPaBeam +
yBeam*sinPaBeam*yBeam*sinPaBeam
- )/(2*signalToNoise);
-
- longErrorFromSN.convert(longError.getUnit());
- latErrorFromSN.convert(latError.getUnit());
- longError.setValue(
- max(longError.getValue(), longErrorFromSN.getValue())
- );
- latError.setValue(
- max(latError.getValue(), latErrorFromSN.getValue())
- );
- newShape->setRefDirectionError(latError, longError);
-
- Vector<Int> index(1, i);
- _results.setShape(index, *newShape);
- }
- }
-
- String ImageFitter::_sizeToString(const uInt compNumber) const {
- ostringstream size;
- const ComponentShape* compShape = _results.getShape(compNumber);
- AlwaysAssert(compShape->type() == ComponentType::GAUSSIAN, AipsError);
-
- Vector<Quantum<Double> > beam = _image->imageInfo().restoringBeam();
- Bool hasBeam = beam.nelements() == 3;
- size << "Image component size";
- if (hasBeam) {
- size << " (convolved with beam)";
- }
- size << " ---" << endl;
- size << compShape->sizeToString() << endl;
- if (hasBeam) {
- Quantity maj = _majorAxes[compNumber];
- Quantity min = _minorAxes[compNumber];
- Quantity pa = _positionAngles[compNumber];
- const GaussianShape *gaussShape = static_cast<const GaussianShape
*>(compShape);
- Quantity emaj = gaussShape->majorAxisError();
- Quantity emin = gaussShape->minorAxisError();
- Quantity epa = gaussShape->positionAngleError();
-
- size << "Clean beam size ---" << endl;
- size << TwoSidedShape::sizeToString(beam[0], beam[1], beam[2],
False) << endl;
- // Quantity femaj = emaj/maj;
- // Quantity femin = emin/min;
- Bool fitSuccess = False;
- Vector<Quantity> best(3);
- best[0] = maj;
- best[1] = min;
- best[2] = pa;
- Quantity majFit;
- Quantity minFit;
- Quantity paFit;
-
- Vector<Quantity> bestFit(3);
- Bool isPointSource = ImageUtilities::deconvolveFromBeam(
- bestFit[0], bestFit[1], bestFit[2], fitSuccess, *_log, best, beam
- );
- size << "Image component size (deconvolved from beam) ---" << endl;
-
- if(fitSuccess) {
- if (isPointSource) {
- size << " Component is a point source" << endl;
- }
- else {
- Vector<Quantity> majRange(2, maj - emaj);
- majRange[1] = maj + emaj;
- Vector<Quantity> minRange(2, min - emin);
- minRange[1] = min + emin;
- Vector<Quantity> paRange(2, pa - epa);
- paRange[1] = pa + epa;
- Vector<Quantity> sourceIn(3);
- for (uInt i=0; i<2; i++) {
- sourceIn[0] = majRange[i];
- for (uInt j=0; j<2; j++) {
- sourceIn[1] = minRange[j];
- for (uInt k=0; k<2; k++) {
- sourceIn[2] = paRange[k];
- minFit = Quantity();
- majFit = Quantity();
- paFit = Quantity();
- isPointSource = ImageUtilities::deconvolveFromBeam(
- majFit, minFit, paFit, fitSuccess,
- *_log, sourceIn, beam, False
- );
- if (fitSuccess) {
- Quantity errMaj = bestFit[0] - majFit;
- errMaj.convert(emaj.getUnit());
- Quantity errMin = bestFit[1] - minFit;
- errMin.convert(emin.getUnit());
- Quantity errPA = bestFit[2] - paFit;
- errPA.convert("deg");
- errPA.setValue(fmod(errPA.getValue(), 180.0));
- errPA.convert(epa.getUnit());
- emaj.setValue(
- max(emaj.getValue(), abs(errMaj.getValue()))
- );
- emin.setValue(
- max(emin.getValue(), abs(errMin.getValue()))
- );
- epa.setValue(
- max(epa.getValue(), abs(errPA.getValue()))
- );
- }
- }
- }
- }
- size << TwoSidedShape::sizeToString(
- bestFit[0], bestFit[1], bestFit[2],
- True, emaj, emin, epa
- );
- }
- }
- else {
- *_log << LogIO::SEVERE << "Could not deconvolve source from beam"
<< LogIO::POST;
- }
- }
- return size.str();
- }
-
- String ImageFitter::_fluxToString(uInt compNumber) const {
- Vector<String> unitPrefix(8);
- unitPrefix[0] = "T";
- unitPrefix[1] = "G";
- unitPrefix[2] = "M";
- unitPrefix[3] = "k";
- unitPrefix[4] = "";
- unitPrefix[5] = "m";
- unitPrefix[6] = "u";
- unitPrefix[7] = "n";
-
- ostringstream fluxes;
- Quantity fluxDensity = _fluxDensities[compNumber];
- Quantity fluxDensityError = _fluxDensityErrors[compNumber];
- Vector<String> polarization = _results.getStokes(compNumber);
- /*
- for (uInt i=0; i<polarization.nelements(); i++) {
- if (polarization[i] == _kludgedStokes) {
- std::complex<double> error =
_results.component(compNumber).flux().errors()[i];
- fluxDensityError.setValue(sqrt(error.real()*error.real() +
error.imag()*error.imag()));
- fluxDensityError.setUnit(fluxDensity.getUnit());
- break;
- }
- }
- */
- String unit;
- for (uInt i=0; i<unitPrefix.size(); i++) {
- unit = unitPrefix[i] + "Jy";
- if (fluxDensity.getValue(unit) > 1) {
- fluxDensity.convert(unit);
- fluxDensityError.convert(unit);
- break;
- }
- }
- Vector<Double> fd(2);
- fd[0] = fluxDensity.getValue();
- fd[1] = fluxDensityError.getValue();
- Quantity peakIntensity = _peakIntensities[compNumber];
- Quantity intensityToFluxConversion =
_image->units().getName().contains("/beam")
- ? Quantity(1.0, "beam")
- : Quantity(1.0, "pixel");
-
- Quantity tmpFlux = peakIntensity * intensityToFluxConversion;
- tmpFlux.convert("Jy");
-
- Quantity peakIntensityError = _peakIntensityErrors[compNumber];
- Quantity tmpFluxError = peakIntensityError *
intensityToFluxConversion;
-
- uInt precision = 0;
- fluxes << "Flux ---" << endl;
-
- if (! _noBeam) {
- precision = precisionForValueErrorPairs(fd, Vector<Double>());
- fluxes << std::fixed << setprecision(precision);
- fluxes << " --- Integrated: " << fluxDensity.getValue()
- << " +/- " << fluxDensityError.getValue() << " "
- << fluxDensity.getUnit() << endl;
- }
-
- for (uInt i=0; i<unitPrefix.size(); i++) {
- String unit = unitPrefix[i] + tmpFlux.getUnit();
- if (tmpFlux.getValue(unit) > 1) {
- tmpFlux.convert(unit);
- tmpFluxError.convert(unit);
- break;
- }
- }
- //String newUnit = tmpFlux.getUnit() + "/" +
intensityToFluxConversion.getUnit();
- peakIntensity = Quantity(tmpFlux.getValue(), tmpFlux.getUnit()
+ "/" + intensityToFluxConversion.getUnit());
- peakIntensityError = Quantity(tmpFluxError.getValue(),
peakIntensity.getUnit());
-
-
- Vector<Double> pi(2);
- pi[0] = peakIntensity.getValue();
- pi[1] = peakIntensityError.getValue();
- precision = precisionForValueErrorPairs(pi, Vector<Double>());
- fluxes << std::fixed << setprecision(precision);
- fluxes << " --- Peak: " << peakIntensity.getValue()
- << " +/- " << peakIntensityError.getValue() << " "
- << peakIntensity.getUnit() << endl;
- fluxes << " --- Polarization: " << _stokesString << endl;
- return fluxes.str();
- }
-
- String ImageFitter::_spectrumToString(uInt compNumber) const {
- Vector<String> unitPrefix(9);
- unitPrefix[0] = "T";
- unitPrefix[1] = "G";
- unitPrefix[2] = "M";
- unitPrefix[3] = "k";
- unitPrefix[4] = "";
- unitPrefix[5] = "c";
- unitPrefix[6] = "m";
- unitPrefix[7] = "u";
- unitPrefix[8] = "n";
- ostringstream spec;
- const SpectralModel& spectrum =
_results.component(compNumber).spectrum();
- Quantity frequency = spectrum.refFrequency().get("MHz");
- Quantity c(C::c, "m/s");
- Quantity wavelength = c/frequency;
- String prefUnit;
- for (uInt i=0; i<unitPrefix.size(); i++) {
- prefUnit = unitPrefix[i] + "Hz";
- if (frequency.getValue(prefUnit) > 1) {
- frequency.convert(prefUnit);
- break;
- }
- }
- for (uInt i=0; i<unitPrefix.size(); i++) {
- prefUnit = unitPrefix[i] + "m";
- if (wavelength.getValue(prefUnit) > 1) {
- wavelength.convert(prefUnit);
- break;
- }
- }
-
- spec << "Spectrum ---" << endl;
- spec << " --- frequency: " << frequency << " (" <<
wavelength << ")" << endl;
- return spec.str();
- }
-
- void ImageFitter::_writeLogfile(const String& output) const {
- File log(_logfileName);
- switch (File::FileWriteStatus status = log.getWriteStatus()) {
- case File::OVERWRITABLE:
- if (logfileAppend) {
- Int fd = open(_logfileName.c_str(), O_RDWR | O_APPEND);
- FiledesIO fio(fd, _logfileName.c_str());
- fio.write(output.length(), output.c_str());
- FiledesIO::close(fd);
- *_log << LogIO::NORMAL << "Appended results to file "
- << _logfileName << LogIO::POST;
- }
- // no break here to fall through to the File::CREATABLE block if
logFileAppend is false
- case File::CREATABLE:
- if (status == File::CREATABLE || ! logfileAppend) {
- // can fall throw from previous case block so status can be
File::OVERWRITABLE
- String action = (status ==
File::OVERWRITABLE) ? "Overwrote" : "Created";
- Int fd = FiledesIO::create(_logfileName.c_str());
- FiledesIO fio (fd, _logfileName.c_str());
- fio.write(output.length(), output.c_str());
- FiledesIO::close(fd);
- *_log << LogIO::NORMAL << action << " file "
- << _logfileName << " with new log file"
- << LogIO::POST;
- }
- break;
- default:
- // checks to see if the log file is not creatable or not writeable
should have already been
- // done and if so _logFileName set to the empty string so this
method wouldn't be called in
- // those cases.
- *_log << "Programming logic error. This block should never be
reached" << LogIO::EXCEPTION;
- }
- }
-
- void ImageFitter::_writeNewEstimatesFile() const {
- ostringstream out;
- for (uInt i=0; i<_results.nelements(); i++) {
-
- MDirection mdir = _results.getRefDirection(i);
- Quantity lat = mdir.getValue().getLat("rad");
- Quantity longitude = mdir.getValue().getLong("rad");
- Vector<Double> world(4,0), pixel(4,0);
- _image->coordinates().toWorld(world, pixel);
- world[0] = longitude.getValue();
- world[1] = lat.getValue();
- if (_image->coordinates().toPixel(pixel, world)) {
- out << _peakIntensities[i].getValue() << ", "
- << pixel[0] << ", " << pixel[1] << ", "
- << _majorAxes[i] << ", " << _minorAxes[i] << ", "
- << _positionAngles[i] << endl;
- }
- else {
- *_log << LogIO::WARN << "Unable to calculate pixel location of "
- << "component number " << i << " so cannot write new estimates"
- << "file" << LogIO::POST;
- return;
- }
- }
- String output = out.str();
- File estimates(_newEstimatesFileName);
- String action = (estimates.getWriteStatus() ==
File::OVERWRITABLE) ? "Overwrote" : "Created";
- Int fd = FiledesIO::create(_newEstimatesFileName.c_str());
- FiledesIO fio(fd, _logfileName.c_str());
- fio.write(output.length(), output.c_str());
- FiledesIO::close(fd);
- *_log << LogIO::NORMAL << action << " file "
- << _newEstimatesFileName << " with new estimates file"
- << LogIO::POST;
- }
-
- void ImageFitter::_writeCompList() {
- if (! _compListName.empty()) {
- switch (_writeControl) {
- case NO_WRITE:
- return;
- case WRITE_NO_REPLACE:
- {
- File file(_compListName);
- if (file.exists()) {
- LogOrigin logOrigin("ImageFitter", __FUNCTION__);
- *_log << logOrigin;
- *_log << LogIO::WARN << "Requested persistent component
list " << _compListName
- << " already exists and user does not wish to overwrite it
so "
- << "the component list will not be written" << LogIO::POST;
- return;
- }
- }
- // allow fall through
- case OVERWRITE: {
- Path path(_compListName);
- _results.rename(path, Table::New);
- *_log << LogIO::NORMAL << "Wrote component list table " <<
_compListName << endl;
- }
- return;
- default:
- // should never get here
- return;
- }
- }
- }
-
-}
-
=======================================
--- /trunk/images/Images/ImageFitter.h Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,242 +0,0 @@
-//# tSubImage.cc: Test program for class SubImage
-//# Copyright (C) 1998,1999,2000,2001,2003
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or (at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#ifndef IMAGES_IMAGEFITTER_H
-#define IMAGES_IMAGEFITTER_H
-
-#include <measures/Measures/Stokes.h>
-#include <lattices/LatticeMath/Fit2D.h>
-#include <casa/Logging/LogIO.h>
-#include <components/ComponentModels/ComponentList.h>
-#include <images/Images/ImageInterface.h>
-#include <images/Images/ImageInputProcessor.h>
-
-#include <components/ComponentModels/ComponentType.h>
-#include <casa/namespace.h>
-
-namespace casa {
-
- class ImageFitter {
- // <summary>
- // Top level interface to ImageAnalysis::fitsky to handle inputs,
bookkeeping etc and
- // ultimately call fitsky to do fitting
- // </summary>
-
- // <reviewed reviewer="" date="" tests="" demos="">
- // </reviewed>
-
- // <prerequisite>
- // </prerequisite>
-
- // <etymology>
- // Fits components to sources in images
(ImageSourceComponentFitter was deemed to be to long
- // of a name)
- // </etymology>
-
- // <synopsis>
- // ImageFitter is the top level interface for fitting image source
components. It handles most
- // of the inputs, bookkeeping etc. It can be instantiated and its
one public method, fit,
- // run from either a C++ app or python.
- // </synopsis>
-
- // <example>
- // <srcblock>
- // ImageFitter fitter(...)
- // fitter.fit()
- // </srcblock>
- // </example>
-
- public:
- enum CompListWriteControl {
- NO_WRITE,
- WRITE_NO_REPLACE,
- OVERWRITE
- };
-
- // constructor approprate for API calls.
- // Parameters:
- // <ul>
- // <li>imagename - the name of the input image in which to fit
the models</li>
- // <li>box - A 2-D rectangular box in which to use pixels for
the fitting, eg box=100,120,200,230
- // In cases where both box and region are specified, box, not
region, is used.</li>
- // <li>region - Named region to use for fitting</li>
- // <li>regionPtr - A pointer to a region. Note there are unfortunately
several different types of
- // region records throughout CASA. In this case, it must be a Record
produced by creating a
- // region via a RegionManager method.
- // <li>chanInp - Zero-based channel number on which to do the
fit. Only a single channel can be
- // specified.</li>
- // <li>stokes - Stokes plane on which to do the fit. Only a
single Stokes parameter can be
- // specified.</li>
- // <li> maskInp - Mask (as LEL) to use as a way to specify
which pixels to use </li>
- // <li> includepix - Pixel value range to include in the fit.
includepix and excludepix
- // cannot be specified simultaneously. </li>
- // <li> excludepix - Pixel value range to exclude from fit</li>
- // <li> residualInp - Name of residual image to save. Blank
means do not save residual image</li>
- // <li> modelInp - Name of the model image to save. Blank
means do not save model image</li>
- // </ul>
-
- // DEPRECATED, DO NOT USE FOR NEW CODE AND CHANGE OLD CODE TO USE ONE
OF THE CONSTRUCTORS BELOW
- ImageFitter(
- const String& imagename, const String& region, const
String& box="",
- const uInt chanInp=0, const String& stokes="I",
- const String& maskInp="",
- const Vector<Float>& includepix = Vector<Float>(0),
- const Vector<Float>& excludepix = Vector<Float>(0),
- const String& residualInp="", const String& modelInp="",
- const String& estiamtesFilename="", const String&
logfile="",
- const Bool& append=True, const String& newEstimatesInp="",
- const String& compListName="",
- const CompListWriteControl writeControl=NO_WRITE
- );
-
- // DEPRECATED, DO NOT USE FOR NEW CODE AND CHANGE OLD CODE TO
USE ONE OF THE CONSTRUCTORS BELOW
- ImageFitter(
- const String& imagename, const Record* regionPtr, const
String& box="",
- const uInt chanInp=0, const String& stokes="I",
- const String& maskInp="",
- const Vector<Float>& includepix = Vector<Float>(0),
- const Vector<Float>& excludepix = Vector<Float>(0),
- const String& residualInp="", const String& modelInp="",
- const String& estiamtesFilename="", const String&
logfile="",
- const Bool& append=True, const String& newEstimatesInp="",
- const String& compListName="",
- const CompListWriteControl writeControl=NO_WRITE
- );
-
- // use these constructors when you already have a pointer to a
valid ImageInterface object
-
- ImageFitter(
- const ImageInterface<Float>*& image, const String& region,
const String& box="",
- const uInt chanInp=0, const String& stokes="I",
- const String& maskInp="",
- const Vector<Float>& includepix = Vector<Float>(0),
- const Vector<Float>& excludepix = Vector<Float>(0),
- const String& residualInp="", const String& modelInp="",
- const String& estiamtesFilename="", const String&
logfile="",
- const Bool& append=True, const String& newEstimatesInp="",
- const String& compListName="",
- const CompListWriteControl writeControl=NO_WRITE
- );
-
- ImageFitter(
- const ImageInterface<Float>*& image, const Record*
regionPtr, const String& box="",
- const uInt chanInp=0, const String& stokes="I",
- const String& maskInp="",
- const Vector<Float>& includepix = Vector<Float>(0),
- const Vector<Float>& excludepix = Vector<Float>(0),
- const String& residualInp="", const String& modelInp="",
- const String& estiamtesFilename="", const String&
logfile="",
- const Bool& append=True, const String& newEstimatesInp="",
- const String& compListName="",
- const CompListWriteControl writeControl=NO_WRITE
- );
-
- // destructor
- ~ImageFitter();
-
- // Do the fit. If componentList is specified, store the fitted
components in
- // that object.
- ComponentList fit();
-
- // Did the fit converge? Throw AipsError if the fit has not
yet been done.
- Bool converged() const;
-
- private:
- LogIO *_log;
- ImageInterface<Float> *_image;
- Record _regionRecord;
- uInt _chan;
- String _stokesString, _mask, _residual, _model, _logfileName,
- regionString, estimatesString, _newEstimatesFileName, _compListName;
- Vector<Float> includePixelRange, excludePixelRange;
- ComponentList estimates, _results;
- Vector<String> fixed;
- Bool logfileAppend, _fitConverged, fitDone, _noBeam,
_deleteImageOnDestruct;
- Vector<Quantity> _peakIntensities, _peakIntensityErrors,
_fluxDensityErrors,
- _fluxDensities, _majorAxes, _majorAxisErrors, _minorAxes,
_minorAxisErrors,
- _positionAngles, _positionAngleErrors;
- Record _residStats, inputStats;
- Double chiSquared;
- String _kludgedStokes;
- CompListWriteControl _writeControl;
-
- // does the lion's share of constructing the object, ie checks
validity of
- // inputs, etc.
- void _construct(
- const String& imagename, const String& box, const String&
regionName,
- const Record* regionPtr, const String& estimatesFilename
- );
-
- void _construct(
- const ImageInterface<Float> *image, const String& box,
const String& regionName,
- const Record* regionPtr, const String& estimatesFilename
- );
-
- Vector<ImageInputProcessor::OutputStruct> _getOutputs();
-
- void _finishConstruction(const String& estimatesFilename);
-
- // summarize the results in a nicely formatted string
- String _resultsToString();
-
- //summarize the size details in a nicely formatted string
- String _sizeToString(const uInt compNumber) const;
-
- String _fluxToString(uInt compNumber) const;
-
- // String _fluxToString2(uInt compNumber) const;
-
-
- String _spectrumToString(uInt compNumber) const;
-
- // write output to log file
- void _writeLogfile(const String& output) const;
-
- // Write the estimates file using this fit.
- void _writeNewEstimatesFile() const;
-
- // Set the flux densities and peak intensities of the fitted
components.
- void _setFluxes();
-
- // Set the convolved sizes of the fitted components.
- void _setSizes();
-
- void _getStandardDeviations(Double& inputStdDev, Double& residStdDev)
const;
-
- void _getRMSs(Double& inputRMS, Double& residRMS) const;
-
- Double _getStatistic(const String& type, const Record& stats) const;
-
- String _statisticsToString() const;
-
- void setErrors(const Record& residStats);
-
- void _writeCompList();
- };
-}
-
-#endif
=======================================
--- /trunk/images/Images/ImageInputProcessor.cc Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,311 +0,0 @@
-
-//# Copyright (C) 1998,1999,2000,2001,2003
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or (at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-
-#include <images/Images/ImageInputProcessor.h>
-
-#include <casa/Containers/HashMap.h>
-#include <casa/Utilities/Sort.h>
-#include <casa/iostream.h>
-
-#include <images/Images/FITSImage.h>
-#include <images/Images/ImageMetaData.h>
-#include <images/Images/ImageUtilities.h>
-#include <images/Images/MIRIADImage.h>
-#include <images/Regions/WCBox.h>
-#include <lattices/Lattices/LCSlicer.h>
-
-#include <measures/Measures/Stokes.h>
-
-ImageInputProcessor::ImageInputProcessor()
-: _log(new LogIO()), _processHasRun(False),
- _nSelectedChannels(0)
-{}
-
-ImageInputProcessor::~ImageInputProcessor() {
- delete _log;
-}
-
-void ImageInputProcessor::process(
- ImageInterface<Float>*& image, Record& regionRecord,
- String& diagnostics, Vector<OutputStruct> *outputStruct,
- String& stokes, const String& imagename, const Record* regionPtr,
- const String& regionName, const String& box,
- const String& chans,
- const RegionManager::StokesControl& stokesControl, const Bool&
allowMultipleBoxes,
- const Vector<Coordinate::Type> *requiredCoordinateTypes
-) {
- LogOrigin origin("ImageInputProcessor", __FUNCTION__);
- *_log << origin;
- if (imagename.empty()) {
- *_log << "imagename cannot be blank" << LogIO::EXCEPTION;
- }
- // Register the functions to create a FITSImage or MIRIADImage object.
- FITSImage::registerOpenFunction();
- MIRIADImage::registerOpenFunction();
- image = 0;
- ImageUtilities::openImage(image, imagename, *_log);
- if (image == 0) {
- *_log << origin;
- *_log << "Unable to open image " << imagename << LogIO::EXCEPTION;
- }
- _process(
- regionRecord, diagnostics, outputStruct, stokes,
- image, regionPtr, regionName, box, chans, stokesControl,
- allowMultipleBoxes, requiredCoordinateTypes
- );
-}
-
-void ImageInputProcessor::process(
- Record& regionRecord, String& diagnostics,
- Vector<OutputStruct> *outputStruct, String& stokes,
- const ImageInterface<Float>*& image,
- const Record* regionPtr, const String& regionName,
- const String& box, const String& chans,
- const RegionManager::StokesControl& stokesControl, const Bool&
allowMultipleBoxes,
- const Vector<Coordinate::Type> *requiredCoordinateTypes
-) {
- _process(
- regionRecord, diagnostics, outputStruct, stokes,
- image, regionPtr, regionName, box, chans, stokesControl,
- allowMultipleBoxes, requiredCoordinateTypes
- );
-}
-
-void ImageInputProcessor::_process(
- Record& regionRecord,
- String& diagnostics, Vector<OutputStruct>* outputStruct,
- String& stokes, const ImageInterface<Float>* image,
- const Record*& regionPtr,
- const String& regionName, const String& box,
- const String& chans, const RegionManager::StokesControl& stokesControl,
- const Bool& allowMultipleBoxes,
- const Vector<Coordinate::Type>* requiredCoordinateTypes
-) {
- LogOrigin origin("ImageInputProcessor", __FUNCTION__);
- *_log << origin;
- if (outputStruct != 0) {
- checkOutputs(outputStruct, *_log);
- }
- *_log << origin;
- if (requiredCoordinateTypes) {
- for (
- Vector<Coordinate::Type>::const_iterator iter =
requiredCoordinateTypes->begin();
- iter != requiredCoordinateTypes->end(); iter++
- ) {
- if (image->coordinates().findCoordinate(*iter) < 0) {
- *_log << "Image " << image->name() << " does not have required
coordinate "
- << Coordinate::typeToString(*iter) << LogIO::EXCEPTION;
- }
- }
- }
- ImageMetaData metaData(*image);
- _nSelectedChannels = metaData.nChannels();
- // region specification order:
- // 1: process box if given explicitly
- // 2. else process region if pointer to region record specified
- // 3. else process region if region name specified
- // 4. else no box or region specified, process region as entire
- // positional plane anding with chans and stokes specs
-
- RegionManager regionMgr(image->coordinates());
- regionRecord = regionMgr.fromBCS(
- diagnostics, _nSelectedChannels, stokes,
- regionPtr, regionName, chans,
- stokesControl, box, image->shape(), image->name()
- );
- if (!allowMultipleBoxes && regionRecord.fieldNumber("regions") >= 0) {
- *_log << "Only a single n-dimensional rectangular region is
supported."
- << LogIO::EXCEPTION;
- }
- _processHasRun = True;
-}
-
-uInt ImageInputProcessor::nSelectedChannels() const {
- if (! _processHasRun) {
- *_log << LogOrigin("ImageInputProcessor", __FUNCTION__);
- *_log << "Programming logic error, ImageInputProcessor::process() must
be called "
- << "before ImageInputProcessor::" << __FUNCTION__ << "()" <<
LogIO::EXCEPTION;
- }
- return _nSelectedChannels;
-}
-
-String ImageInputProcessor::_stokesFromRecord(
- const Record& region, const CoordinateSystem& csys
-) const {
- // FIXME This implementation is incorrect for complex, recursive records
- String stokes = "";
- if(csys.hasPolarizationAxis()) {
- Int polAxis = csys.polarizationAxisNumber();
- uInt stokesBegin, stokesEnd;
- ImageRegion *imreg = ImageRegion::fromRecord(region, "");
- Array<Float> blc, trc;
- Bool oneRelAccountedFor = False;
- if (imreg->isLCSlicer()) {
- blc = imreg->asLCSlicer().blc();
- trc = imreg->asLCSlicer().trc();
- stokesBegin = (uInt)((Vector<Float>)blc)[polAxis];
- stokesEnd = (uInt)((Vector<Float>)trc)[polAxis];
- oneRelAccountedFor = True;
- }
- else if
(RegionManager::isPixelRegion(*(ImageRegion::fromRecord(region, "")))) {
- region.toArray("blc", blc);
- region.toArray("trc", trc);
- stokesBegin = (uInt)((Vector<Float>)blc)[polAxis];
- stokesEnd = (uInt)((Vector<Float>)trc)[polAxis];
- }
- else {
- Record blcRec = region.asRecord("blc");
- Record trcRec = region.asRecord("trc");
- stokesBegin = (Int)blcRec.asRecord(
- String("*" + String::toString(polAxis - 1))
- ).asDouble("value");
- stokesEnd = (Int)trcRec.asRecord(
- String("*" + String::toString(polAxis - 1))
- ).asDouble("value");
- }
-
- if (! oneRelAccountedFor && region.isDefined("oneRel") &&
region.asBool("oneRel")) {
- stokesBegin--;
- stokesEnd--;
- }
- for (uInt i=stokesBegin; i<=stokesEnd; i++) {
- stokes += csys.stokesAtPixel(i);
- }
- }
- return stokes;
-}
-
-void ImageInputProcessor::_setRegion(
- Record& regionRecord, String& diagnostics,
- const Record* regionPtr
-) const {
- // region record pointer provided
- regionRecord = *(regionPtr->clone());
- // set stokes from the region record
- diagnostics = "used provided region record";
-}
-
-void ImageInputProcessor::_setRegion(
- Record& regionRecord, String& diagnostics,
- const ImageInterface<Float> *image, const String& regionName
-) const {
- // region name provided
- File myFile(regionName);
- if (myFile.exists()) {
- Record *rec = RegionManager::readImageFile(regionName, "");
- regionRecord = *rec;
- delete rec;
- diagnostics = "Region read from file " + regionName;
- }
- else {
-
- ImageRegion imRegion;
- Regex otherImage("(.*)+:(.*)+");
- try {
- String imagename;
- if (regionName.matches(otherImage)) {
- String res[2];
- casa::split(regionName, res, 2, ":");
- imagename = res[0];
- PagedImage<Float> other(imagename);
- imRegion = other.getRegion(res[1]);
- }
- else {
- imRegion = image->getRegion(regionName);
- imagename = image->name();
- }
- regionRecord = Record(imRegion.toRecord(""));
- diagnostics = "Used region " + regionName + " from image "
- + imagename + " table description";
- }
- catch (AipsError) {
- *_log << "Unable to open region file or region table description " <<
regionName << LogIO::EXCEPTION;
- }
- }
-}
-
-String ImageInputProcessor::_pairsToString(const Vector<uInt>& pairs)
const {
- ostringstream os;
- uInt nPairs = pairs.size()/2;
- for (uInt i=0; i<nPairs; i++) {
- os << pairs[2*i] << " - " << pairs[2*i + 1];
- if (i < nPairs - 1) {
- os << "; ";
- }
- }
- return os.str();
-}
-
-void ImageInputProcessor::checkOutputs(
- Vector<OutputStruct> *output, LogIO& log
-) {
- for (
- Vector<OutputStruct>::iterator iter = output->begin();
- iter != output->end();
- iter++
- ) {
- String label = iter->label;
- String name = *(iter->outputFile);
- Bool required = iter->required;
- Bool replaceable = iter->replaceable;
- if (name.empty()) {
- if (required) {
- log << label << " cannot be blank" << LogIO::EXCEPTION;
- }
- else {
- continue;
- }
- }
- LogIO::Command logLevel = required ? LogIO::SEVERE : LogIO::WARN;
- LogIO::Command logAction = required ? LogIO::EXCEPTION : LogIO::POST;
- File f(name);
- switch (f.getWriteStatus()) {
- case File::NOT_CREATABLE:
- log << logLevel << "Requested " << label << " " << name
- << " cannot be created so will not be written" << logAction;
- *(iter->outputFile) = "";
- break;
- case File::NOT_OVERWRITABLE:
- log << logLevel << "There is already a file or directory named "
- << name << " which cannot be overwritten so the " << label
- << " will not be written" << logAction;
- *(iter->outputFile) = "";
- break;
- case File::OVERWRITABLE:
- if (! replaceable) {
- log << logLevel << "Replaceable flag is false and there is "
- << "already a file or directory named " << name
- << " so the " << label << " will not be written"
- << logAction;
- *(iter->outputFile) = "";
- }
- break;
- default:
- continue;
- }
- }
-}
-
=======================================
--- /trunk/images/Images/ImageInputProcessor.h Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,175 +0,0 @@
-//# Copyright (C) 1998,1999,2000,2001,2003
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or (at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-
-#ifndef IMAGES_IMAGEINPUTPROCESSOR_H
-#define IMAGES_IMAGEINPUTPROCESSOR_H
-
-#include <images/Images/ImageInterface.h>
-#include <images/Images/ImageMetaData.h>
-#include <images/Regions/RegionManager.h>
-
-#include <casa/namespace.h>
-
-namespace casa {
-
-class ImageInputProcessor {
- // <summary>
- // Collection of methods for processing inputs to image analysis
applications
- // </summary>
-
- // <reviewed reviewer="" date="" tests="" demos="">
- // </reviewed>
-
- // <prerequisite>
- // </prerequisite>
-
- // <etymology>
- // Processes inputs to image analysis apps.
- // </etymology>
-
- // <synopsis>
- // Collection of methods for processing inputs to image analysis
applications
- // </synopsis>
-
-public:
- // instruction if input stokes is blank
- /*
- enum StokesControl {
- USE_FIRST_STOKES,
- USE_ALL_STOKES
- };
- */
-
- // const static String ALL;
-
- // struct for checking output file writability
- struct OutputStruct {
- // label used for messages, eg "residual image", "estmates file"
- String label;
- // pointer to the output name
- String *outputFile;
- // is this file required to be written, or can the task continue if it
cannot be?
- Bool required;
- // If a file by the same name already exists, will the task allow it to
be overwritten?
- Bool replaceable;
- };
-
- //constructor
- ImageInputProcessor();
-
- //Destructor
- ~ImageInputProcessor();
-
- // Process the inputs. Output parameters are the pointer to the
- // opened <src>image</src>, the specified region as a record (<src>
- // regionRecord</src>, and a <src>diagnostics</src> String describing
- // how the region was chosen. If provided, <src>regionPtr</src> should
- // be a pointer to a record created by a RegionManager method.
- // <src>stokesControl</src> indicates default
- // stokes range to use if <src>stokes</src> is blank. In this case
<src>stokes</src>
- // will be set the the value of stokes that will be used. If
- // <src>allowMultipleBoxes</src> is False, an exception will be thrown if
- // the inputs specify multiple n-dimensional rectangles. This should
usually
- // be set to false if the caller can only deal with a single n-dimensional
- // rectangular region.
- void process(
- ImageInterface<Float>*& image, Record& regionRecord,
- String& diagnostics, Vector<OutputStruct> *outputStruct,
- String& stokes,
- const String& imagename, const Record* regionPtr,
- const String& regionName, const String& box,
- const String& chans,
- const RegionManager::StokesControl& stokesControl, const Bool&
allowMultipleBoxes,
- const Vector<Coordinate::Type> *requiredCoordinateTypes
- );
-
- // Process the inputs. Use this version if the associated image already
exists.
- // Output parameters the specified region as a record (<src>
- // regionRecord</src>, and a <src>diagnostics</src> String describing
- // how the region was chosen. If provided, <src>regionPtr</src> should
- // be a pointer to a record created by a RegionManager method.
- // <src>stokesControl</src> indicates default
- // stokes range to use if <src>stokes</src> is blank. In this case
<src>stokes</src>
- // will be set the the value of stokes that will be used. If
- // <src>allowMultipleBoxes</src> is False, an exception will be thrown if
- // the inputs specify multiple n-dimensional rectangles. This should
usually
- // be set to false if the caller can only deal with a single n-dimensional
- // rectangular region.
- void process(
- Record& regionRecord,
- String& diagnostics, Vector<OutputStruct> *outputStruct,
- String& stokes,
- const ImageInterface<Float>*& image,
- const Record* regionPtr,
- const String& regionName, const String& box,
- const String& chans,
- const RegionManager::StokesControl& stokesControl,
- const Bool& allowMultipleBoxes,
- const Vector<Coordinate::Type> *requiredCoordinateTypes
- );
-
- static void checkOutputs(Vector<OutputStruct> *output, LogIO& log);
-
- // Get the number of channels that have been selected. The process()
method must
- // be called prior to calling this method or an exception is thrown.
- uInt nSelectedChannels() const;
-
-private:
- LogIO *_log;
- Bool _processHasRun;
- uInt _nSelectedChannels;
-
- void _process(
- Record& regionRecord,
- String& diagnostics, Vector<OutputStruct>* outputStruct,
- String& stokes, const ImageInterface<Float>* image,
- const Record*& regionPtr,
- const String& regionName, const String& box,
- const String& chans, const RegionManager::StokesControl&
stokesControl,
- const Bool& allowMultipleBoxes,
- const Vector<Coordinate::Type>* requiredCoordinateTypes
- );
-
- // set region given a pointer to a region record.
- void _setRegion(
- Record& regionRecord, String& diagnostics,
- const Record* regionPtr
- ) const;
-
- void _setRegion(Record& regionRecord, String& diagnostics,
- const ImageInterface<Float> *image, const String& regionName
- ) const;
-
- String _stokesFromRecord(
- const Record& region, const CoordinateSystem& csys
- ) const;
-
- String _pairsToString(const Vector<uInt>& pairs) const;
-
-};
-
-}
-
-#endif /* IMAGES_IMAGEINPUTPROCESSOR_H */
=======================================
--- /trunk/images/Images/ImageProfileFit.cc Mon Jun 8 19:16:01 2009
+++ /dev/null
@@ -1,1050 +0,0 @@
-//# ImageProfileFit.cc: Class to fit spectra from images
-//# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This library is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU Library General Public License as published
by
-//# the Free Software Foundation; either version 2 of the License, or (at
your
-//# option) any later version.
-//#
-//# This library is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
-//# License for more details.
-//#
-//# You should have received a copy of the GNU Library General Public
License
-//# along with this library; if not, write to the Free Software Foundation,
-//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addessed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#include <images/Images/ImageProfileFit.h>
-
-#include <casa/Arrays/IPosition.h>
-#include <casa/Arrays/Vector.h>
-#include <casa/Arrays/ArrayMath.h>
-#include <casa/Containers/Record.h>
-#include <coordinates/Coordinates/CoordinateSystem.h>
-#include <coordinates/Coordinates/SpectralCoordinate.h>
-#include <coordinates/Coordinates/CoordinateUtil.h>
-#include <casa/Exceptions/Error.h>
-#include <images/Images/ImageInterface.h>
-#include <images/Regions/ImageRegion.h>
-#include <images/Images/SubImage.h>
-#include <lattices/Lattices/TiledLineStepper.h>
-#include <lattices/Lattices/MaskedLatticeIterator.h>
-#include <lattices/Lattices/LCRegion.h>
-#include <lattices/Lattices/MaskedLattice.h>
-#include <lattices/Lattices/LatticeUtilities.h>
-#include <casa/Logging/LogIO.h>
-#include <casa/BasicSL/Constants.h>
-#include <casa/BasicMath/Math.h>
-#include <casa/System/ProgressMeter.h>
-#include <casa/Quanta/Quantum.h>
-#include <casa/Utilities/Assert.h>
-#include <components/SpectralComponents/SpectralFit.h>
-#include <components/SpectralComponents/SpectralEstimate.h>
-#include <components/SpectralComponents/SpectralElement.h>
-#include <components/SpectralComponents/SpectralList.h>
-
-//#include <casa/OS/Timer.h>
-
-
-#include <casa/ostream.h>
-
-
-namespace casa { //# NAMESPACE CASA - BEGIN
-
-ImageProfileFit::ImageProfileFit()
-: itsImagePtr(0),
- itsSigmaImagePtr(0),
- itsFitDone(False),
- itsMask(0),
- itsSigma(0),
- itsSpectralFitPtr(0),
- itsProfileAxis(-1),
- itsDoppler(""),
- itsXAbs(True),
- itsFitRegion(False)
-{
- UnitMap::putUser("pix",UnitVal(1.0), "pixel units");
- itsSpectralFitPtr = new SpectralFit;
-}
-
-ImageProfileFit::ImageProfileFit(const ImageProfileFit& other)
-: itsImagePtr(0),
- itsSigmaImagePtr(0),
- itsFitDone(other.itsFitDone),
- itsMask(other.itsMask.copy()),
- itsSigma(other.itsSigma.copy()),
- itsCoords(other.itsCoords),
- itsProfileAxis(other.itsProfileAxis),
- itsDoppler(other.itsDoppler),
- itsXAbs(other.itsXAbs),
- itsFitRegion(other.itsFitRegion)
-{
- UnitMap::putUser("pix",UnitVal(1.0), "pixel units");
- itsX = other.itsX; // Get a copy of the vector
- itsY = other.itsY; // Get a copy of the vector
-//
- itsSpectralFitPtr = new SpectralFit (*(other.itsSpectralFitPtr));
-}
-
-
-ImageProfileFit::~ImageProfileFit()
-{
- delete itsSpectralFitPtr;
- itsSpectralFitPtr = 0;
- if (itsImagePtr) delete itsImagePtr;
- if (itsSigmaImagePtr) delete itsSigmaImagePtr;
-}
-
-
-ImageProfileFit& ImageProfileFit::operator=(const ImageProfileFit& other)
-{
- if (this != &other) {
- if (itsImagePtr) {
- delete itsImagePtr;
- itsImagePtr = 0;
- }
- if (other.itsImagePtr) {
- itsImagePtr = other.itsImagePtr->cloneII();
- }
-//
- if (itsSigmaImagePtr) {
- delete itsSigmaImagePtr;
- itsSigmaImagePtr = 0;
- }
- if (other.itsSigmaImagePtr) {
- itsSigmaImagePtr = other.itsSigmaImagePtr->cloneII();
- }
-//
- itsFitDone = other.itsFitDone;
-
-// Does this look after resizing of itsX, itsY vectors ?
-
- itsX = other.itsX;
- itsY = other.itsY;
-//
- itsMask.resize(0);
- itsMask = other.itsMask.copy();
- itsSigma.resize(0);
- itsSigma= other.itsSigma.copy();
-//
- itsCoords = other.itsCoords;
- itsProfileAxis = other.itsProfileAxis;
-//
- itsDoppler = other.itsDoppler;
- itsXAbs = other.itsXAbs;
- itsFitRegion = other.itsFitRegion;
-//
- delete itsSpectralFitPtr;
- itsSpectralFitPtr = 0;
- itsSpectralFitPtr = new SpectralFit (*(other.itsSpectralFitPtr));
- }
- return *this;
-}
-
-void ImageProfileFit::setData (const ImageInterface<Float>& image,
- const ImageRegion& region,
- uInt profileAxis, Bool average)
-{
- itsCoords = image.coordinates();
- itsProfileAxis = profileAxis;
- const SubImage<Float> subImage(image, region);
- const Slicer& sl = region.asLCRegion().boundingBox();
-//
- const ImageInterface<Float>* pSubSigma = 0;
-//
- setData(pSubSigma, subImage, sl, average);
-}
-
-void ImageProfileFit::setData (const ImageInterface<Float>& image,
- const ImageInterface<Float>& sigma,
- const ImageRegion& region,
- uInt profileAxis, Bool average)
-{
- itsCoords = image.coordinates();
- itsProfileAxis = profileAxis;
- const SubImage<Float> subImage(image, region);
- const Slicer& sl = region.asLCRegion().boundingBox();
-//
- const SubImage<Float> subSigma(sigma, region);
- const ImageInterface<Float>* pSubSigma = &subSigma;
-//
- setData(pSubSigma, subImage, sl, average);
-}
-
-
-void ImageProfileFit::setData (const ImageInterface<Float>& image,
- uInt profileAxis, Bool average)
-{
- itsCoords = image.coordinates();
- itsProfileAxis = profileAxis;
- IPosition start(image.ndim(),0);
- Slicer sl(start, image.shape(), Slicer::endIsLength);
-//
- const ImageInterface<Float>* pSubSigma = 0;
-//
- setData(pSubSigma, image, sl, average);
-}
-
-
-
-
-void ImageProfileFit::setData (const ImageInterface<Float>& image,
- const ImageInterface<Float>& sigma,
- uInt profileAxis, Bool average)
-{
- itsCoords = image.coordinates();
- itsProfileAxis = profileAxis;
- IPosition start(image.ndim(),0);
- Slicer sl(start, image.shape(), Slicer::endIsLength);
-//
- const ImageInterface<Float>* pSubSigma = &sigma;
-//
- setData(pSubSigma, image, sl, average);
-}
-
-void ImageProfileFit::setData (const CoordinateSystem& cSys,
- uInt profileAxis,
- const Quantum<Vector<Float> >& x,
- const Quantum<Vector<Float> >& y,
- const Vector<Float>& sigma,
- Bool isAbs, const String& doppler)
-{
- Vector<Bool> mask;
- setData (cSys, profileAxis, x, y, mask, sigma, isAbs, doppler);
-}
-
-void ImageProfileFit::setData (const CoordinateSystem& cSys,
- uInt profileAxis,
- const Quantum<Vector<Float> >& x,
- const Quantum<Vector<Float> >& y,
- const Vector<Bool>& mask,
- const Vector<Float>& sigma,
- Bool isAbs, const String& doppler)
-{
- uInt n = x.getValue().nelements();
- AlwaysAssert(n==y.getValue().nelements(),AipsError);
- AlwaysAssert(n>0,AipsError);
-//
- itsMask.resize(n);
- if (mask.nelements()!=n) {
- itsMask = True;
- } else {
- itsMask = mask;
- }
-//
- itsSigma.resize(n);
- if (sigma.nelements()!=n) {
- itsSigma = 1.0;
- } else {
- itsSigma = sigma.copy();
- }
-//
- itsCoords = cSys;
- itsProfileAxis = profileAxis;
- itsFitRegion = False;
- itsDoppler = doppler;
- itsXAbs = isAbs;
-
-// Does this take care of resizing ?
-
- itsX = x;
- itsY = y;
-}
-
-
-
-uInt ImageProfileFit::addElements (const RecordInterface& rec)
-//
-// Decode the record and add the SpectralElements to the
-// private fitter. We convert the units of the model
-// to those of the data source
-//
-{
- LogIO os(LogOrigin("ImageProfileFit", "addElements", WHERE));
- if (!rec.isDefined("xunit")) {
- os << "Record holding model is missing 'xunit' field" <<
LogIO::EXCEPTION;
- }
- if (!rec.isDefined("yunit")) {
- os << "Record holding model is missing 'yunit' field" <<
LogIO::EXCEPTION;
- }
- if (!rec.isDefined("xabs")) {
- os << "Record holding model is missing 'xabs' field" <<
LogIO::EXCEPTION;
- }
- if (!rec.isDefined("elements")) {
- os << "Record holding model is missing 'elements' field" <<
LogIO::EXCEPTION;
- }
-
-// Setup input units etc
-
- String xUnitIn = rec.asString("xunit");
- String dopplerIn;
- if (rec.isDefined("doppler")) {
- dopplerIn = rec.asString("doppler");
- }
- Bool xAbsIn = rec.asBool("xabs");
- String yUnitIn = rec.asString("yunit");
-
-// Setup output units
-
- String xUnitOut = itsX.getFullUnit().getName();
- Bool xAbsOut = itsXAbs;
- String dopplerOut = itsDoppler;
- String yUnitOut = itsY.getFullUnit().getName();
-
-// Loop over elements in record
-
- if (rec.dataType("elements") != TpRecord) {
- os << "Record holding model is invalid - 'elements' is not a record"
<< LogIO::EXCEPTION;
- }
- Record rec2 = rec.asRecord("elements");
-//
- const uInt n = rec2.nfields();
- if (n <=0) {
- os << "The 'elements' field of record is of zero length" <<
LogIO::EXCEPTION;
- }
-//
- String error;
- for (uInt i=0; i<n; i++) {
- Record rec3 = rec2.asRecord(i);
- SpectralElement se;
- if (!se.fromRecord(error, rec3)) {
- os << error << LogIO::EXCEPTION;
- }
-
-// Convert element to required units
-
- SpectralElement seOut = convertSpectralElement (se, xAbsIn, xAbsOut,
- xUnitIn, xUnitOut,
- dopplerIn,
dopplerOut,
- yUnitIn, yUnitOut);
-
-// See if we have the 'fixed' record; set fixed parameters as indicated
-
- Vector<Bool> fixed;
- if (rec3.isDefined("fixed")) fixed = rec3.asArrayBool("fixed");
-//
- if (fixed.nelements()!=0) {
- if (seOut.getType()==SpectralElement::GAUSSIAN) {
- if (fixed.nelements()==3) {
- if (fixed(0)) seOut.fixAmpl(True);
- if (fixed(1)) seOut.fixCenter(True);
- if (fixed(2)) seOut.fixFWHM(True);
- } else {
- os << "'fixed' vector must be of length 3 for a Gaussian"
<< LogIO::EXCEPTION;
- }
- } else {
- os << "element type not supported" << LogIO::EXCEPTION;
- }
- }
-//
- itsSpectralFitPtr->addFitElement(seOut);
- }
-//
- itsSpectralFitter = *itsSpectralFitPtr;
- return itsSpectralFitPtr->list().nelements() - 1;
-}
-
-Bool ImageProfileFit::getList (RecordInterface& rec,
- Bool xAbsOut,
- const String& xUnitOut,
- const String& dopplerOut)
-{
- const SpectralList& list = itsSpectralFitter.list();
- SpectralList list2 = filterList(list);
- Bool ok = getElements(rec, xAbsOut, xUnitOut, dopplerOut, list2);
- return ok;
-}
-
-SpectralList ImageProfileFit::getList () const
-{
- const SpectralList& list = itsSpectralFitter.list();
- return filterList(list);
-}
-
-
-
-Bool ImageProfileFit::getElements (RecordInterface& rec,
- Bool xAbsOut,
- const String& xUnitOut,
- const String& dopplerOut,
- const SpectralList& list)
-
-//
-// Convert internal model to a record
-//
-{
- LogIO os(LogOrigin("ImageProfileFit", "getElements", WHERE));
- String pix("pix");
-
-// Setup units
-
- Bool xAbsIn = itsXAbs;
- String xUnitIn = itsX.getFullUnit().getName();
- String dopplerIn = itsDoppler;
- String yUnitIn = itsY.getFullUnit().getName();
-//
- String yUnitOut = yUnitIn;
-//
- Record recVec;
- const uInt n = list.nelements();
- for (uInt i=0; i<n; i++) {
-
-// Convert element to desired units
-
- SpectralElement seOut = convertSpectralElement (list[i], xAbsIn,
xAbsOut,
- xUnitIn, xUnitOut,
- dopplerIn,
dopplerOut,
- yUnitIn, yUnitOut);
-
-// Fill record for this element
-
- Record recEl;
- String error;
- if (!seOut.toRecord(error, recEl)) {
- os << error << LogIO::EXCEPTION;
- }
-
-// Stick it in output vector record
-
- recVec.defineRecord(i, recEl);
- }
-
-// Put vector into field "elements"
-
- rec.defineRecord("elements", recVec);
-
-// Add extra fields
-
- rec.define("xunit", xUnitOut);
- rec.define("xabs", xAbsOut);
- rec.define("doppler", dopplerOut);
- rec.define("yunit", yUnitOut);
-//
- return True;
-}
-
-
-void ImageProfileFit::reset ()
-{
- delete itsSpectralFitPtr;
- itsSpectralFitPtr = 0;
- itsSpectralFitPtr = new SpectralFit;
- itsSpectralFitter = *itsSpectralFitPtr;
- itsFitDone = False;
-}
-
-
-Bool ImageProfileFit::estimate (uInt nMax)
-{
- delete itsSpectralFitPtr;
- itsSpectralFitPtr = 0;
- itsSpectralFitPtr = new SpectralFit;
-
-// The estimate is made with x-units of absolute 0-rel pixels
-// so we need to convert it to the data source units
-
- SpectralEstimate est;
- est.setQ(5);
- Vector<Float> der(itsY.getValue().nelements());
- const SpectralList& list = est.estimate(itsY.getValue(), &der);
-//
- uInt nEl = 0;
- if (nMax==0) {
- nEl = list.nelements();
- } else {
- nEl = min(nMax, list.nelements());
- }
- if (nEl==0) {
- return False;
- }
-
-// Setup input units
-
- String xUnitIn("pix");
- String dopplerIn = itsDoppler;
- Bool xAbsIn = True;
- String yUnitIn = itsY.getFullUnit().getName();
-
-// Setup output units
-
- String xUnitOut = itsX.getFullUnit().getName();
- Bool xAbsOut = itsXAbs;
- String dopplerOut = dopplerIn;
- String yUnitOut = yUnitIn;
-
-// Convert
-
- for (uInt i=0; i<nEl; i++) {
- SpectralElement se = convertSpectralElement (list[i], xAbsIn,
xAbsOut,
- xUnitIn, xUnitOut,
- dopplerIn, dopplerOut,
- yUnitIn, yUnitOut);
-
- itsSpectralFitPtr->addFitElement(se);
- }
- itsSpectralFitter = *itsSpectralFitPtr;
-//
- return True;
-}
-
-
-uInt ImageProfileFit::nElements ()
-{
- return itsSpectralFitPtr->list().nelements();
-}
-
-
-Bool ImageProfileFit::fit(Int order)
-{
- itsSpectralFitter = *itsSpectralFitPtr;
-//
- if (order >= 0) {
-
-// Add polynomial to fitter
-
- SpectralElement el(order);
- itsSpectralFitter.addFitElement(el);
- }
-
-// Do fit
-
- Bool ok = itsSpectralFitter.fit(itsSigma,
- itsY.getValue(),
- itsX.getValue(),
- itsMask);
- itsFitDone = True;
- return ok;
-}
-
-
-void ImageProfileFit::residual(Vector<Float>& resid,
- const Vector<Float>& x) const
-{
- resid.resize(0);
- resid = itsY.getValue();
- itsSpectralFitter.list().residual(resid, x);
-}
-
-
-void ImageProfileFit::residual(Vector<Float>& resid) const
-{
- resid.resize(0);
- resid = itsY.getValue();
- itsSpectralFitter.list().residual(resid, itsX.getValue());
-}
-
-
-void ImageProfileFit::model(Vector<Float>& model,
- const Vector<Float>& x) const
-{
- itsSpectralFitter.list().evaluate(model, x);
-}
-
-void ImageProfileFit::model(Vector<Float>& model) const
-{
- itsSpectralFitter.list().evaluate(model, itsX.getValue());
-}
-
-// Private functions
-
-void ImageProfileFit::collapse (Vector<Float>& profile, Vector<Bool>& mask,
- uInt profileAxis, const
MaskedLattice<Float>& lat) const
-{
- if (lat.ndim()==1) {
-
-// Nothing to collapse
-
- profile.resize(0);
- mask.resize(0);
- profile = lat.get();
- mask = lat.getMask();
- } else {
- IPosition axes = IPosition::otherAxes(lat.ndim(),
IPosition(1,profileAxis));
- Array<Float> p;
- Array<Bool> m;
- LatticeUtilities::collapse(p, m, axes, lat, True);
-//
- uInt n = p.shape()(0);
- profile.resize(n);
- mask.resize(n);
- profile = p.reform(IPosition(1,n));
- mask = m.reform(IPosition(1,n));
- }
- }
-
-
-
-
-void ImageProfileFit::setData (const ImageInterface<Float>*& pSigma,
- const ImageInterface<Float>& image,
- const Slicer& sl, Bool average)
-{
- itsXAbs = True;
- const uInt n = image.shape()(itsProfileAxis);
- if (average) {
-
-// Average data over region except along profile axis
-
- Vector<Float> y;
- collapse(y, itsMask, itsProfileAxis, image);
- itsY = Quantum<Vector<Float> >(y, image.units());
-//
- Vector<Float> x(y.nelements());
- indgen(x, Float(sl.start()(itsProfileAxis)));
- itsX = Quantum<Vector<Float> >(x, Unit("pix"));
-//
- if (pSigma) {
- AlwaysAssert (image.shape().conform(pSigma->shape()), AipsError);
-//
- IPosition axes = IPosition::otherAxes(pSigma->ndim(),
IPosition(1,itsProfileAxis));
- Array<Float> p;
- Array<Bool> m;
- LatticeUtilities::collapse(p, m, axes, *pSigma, True); //
Ignore mask of sigma
- itsSigma.resize(0);
- itsSigma = p;
- } else {
- itsSigma.resize(n);
- itsSigma = 1.0;
- }
-//
- itsFitRegion = False;
- } else {
-
-// We are going to fit all profiles in the region. However, the
-// estimate function might be used to fish out an initial
-// estimate and it uses itsX and itsY...
-// I think i need to rewrite this class !
-
- Vector<Float> x(n);
- indgen(x, Float(sl.start()(itsProfileAxis)));
- itsX = Quantum<Vector<Float> >(x, Unit("pix"));
-//
- IPosition blc(image.ndim(),0);
- IPosition trc = blc;
- trc(itsProfileAxis) = n-1;
- Vector<Float> y = image.getSlice(blc,trc-blc+1,True);
- itsY = Quantum<Vector<Float> >(y, image.units());
-//
- itsFitRegion = True;
- itsImagePtr = image.cloneII();
-//
- if (pSigma) {
- AlwaysAssert (image.shape().conform(pSigma->shape()), AipsError);
- itsSigmaImagePtr = pSigma->cloneII();
- }
- }
-}
-
-
-SpectralElement ImageProfileFit::convertSpectralElement (const
SpectralElement& elIn,
- Bool xAbsIn, Bool
xAbsOut,
- const String&
xUnitIn,
- const String&
xUnitOut,
- const String&
dopplerIn,
- const String&
dopplerOut,
- const String&
yUnitIn,
- const String&
yUnitOut)
-{
- LogIO os(LogOrigin("ImageProfileFit", "convertSpectralElement", WHERE));
-//
- SpectralElement elOut(elIn);
- if (elOut.getType()==SpectralElement::GAUSSIAN) {
-
-// Brightness
-
- if (yUnitIn != yUnitOut) {
- Unit tIn(yUnitIn);
- Quantum<Double> amp(elOut.getAmpl(), tIn);
- Quantum<Double> ampErr(elOut.getAmplErr(), tIn);
-//
- Unit tOut(yUnitOut);
- elOut.setAmpl(amp.getValue(tOut));
- Vector<Double> errors;
- elOut.getError(errors);
- errors(0) = ampErr.getValue(tOut);
- elOut.setError(errors);
- }
-//
- convertGaussianElementX (elOut, itsCoords, itsProfileAxis,
- xAbsIn, xAbsOut,
- xUnitIn, xUnitOut,
- dopplerIn, dopplerOut);
- } else {
- os << "Elements of type " << elOut.getType() << " are not yet
supported" << LogIO::EXCEPTION;
- }
-//
- return elOut;
-}
-
-
-
-
-void ImageProfileFit::convertGaussianElementX (SpectralElement& el,
- CoordinateSystem& cSys,
- uInt profileAxis,
- Bool absIn, Bool absOut,
- const String& unitIn,
- const String& unitOut,
- const String& dopplerIn,
- const String& dopplerOut)
-//
-// The data source is an image
-//
-{
- LogIO os(LogOrigin("ImageProfileFit", "convertGaussianElementX",
WHERE));
-//
- Unit velUnit(String("m/s"));
- String pix("pix");
-
-// See if we have something to do. One of the doppler strings
-// might be empty so only check if doing velocity.
-
- if (unitIn==unitOut && absIn==absOut) {
- Unit tIn(unitIn);
- Unit tOut(unitOut);
- if (tIn==velUnit && tOut==velUnit) {
- if (dopplerIn==dopplerOut) return;
- }
- }
-
-// Assumes world and pixel axes the same.
-// Should really check order as well...
-
- AlwaysAssert (cSys.nWorldAxes() == cSys.nPixelAxes(), AipsError);
- const uInt n = cSys.nWorldAxes();
-
-// Get dopplers sorted out
-
- MDoppler::Types dopplerTypeIn, dopplerTypeOut;
- if (Unit(unitIn)==velUnit) {
- if (dopplerIn.empty()) {
- os << "Input doppler is not specified" << LogIO::EXCEPTION;
- }
- if (!MDoppler::getType(dopplerTypeIn, dopplerIn)) {
- os << "Input doppler type is invalid" << LogIO::EXCEPTION;
- }
- }
- if (Unit(unitOut)==velUnit) {
- if (dopplerOut.empty()) {
- os << "Output doppler is not specified" << LogIO::EXCEPTION;
- }
- if (!MDoppler::getType(dopplerTypeOut, dopplerOut)) {
- os << "Output doppler type is invalid" << LogIO::EXCEPTION;
- }
- }
-
-// Find input fractional error
-
- Vector<Double> errors, pars;
- el.getError(errors);
- el.get(pars);
- errors /= pars;
-
-// Declare conversion vectors
-
- Vector<Double> coordIn(n), coordOut(n);
- Vector<Bool> absIns(n), absOuts(n);
- Vector<String> unitsIn(n), unitsOut(n);
-
-// Input; first make a vector at the reference of the correct abs/rel
-
- if (unitIn==pix) {
- coordIn = cSys.referencePixel();
- unitsIn = pix;
- if (!absIn) cSys.makePixelRelative(coordIn);
- } else {
- coordIn = cSys.referenceValue();
- unitsIn = cSys.worldAxisUnits();
- if (!absIn) cSys.makeWorldRelative(coordIn);
- }
- unitsIn(profileAxis) = unitIn;
- absIns = absIn;
-
-// Output; Just convert all except profile axis to pixels to make it easy.
-
- absOuts = absOut;
- unitsOut = pix;
- unitsOut(profileAxis) = unitOut;
-
-// Setup pixel offsets (always 0-rel).
-
- Double offsetIn = 0.0;
- Double offsetOut = 0.0;
-//
-// Get center and set values
-//
- Double centerValue = el.getCenter();
- coordIn(profileAxis) = centerValue;
-
-// Convert to 0-rel absolute pixels
-
- if (!cSys.convert (coordOut, coordIn, absIns, unitsIn, dopplerTypeIn,
- absOuts, unitsOut, dopplerTypeOut, offsetIn,
offsetOut)) {
- os << cSys.errorMessage() << LogIO::EXCEPTION;
- }
-
-// Fish out error then set new center value in element (sets error to 0)
-
- Double centerValueErr = el.getCenterErr();
- Double newCenterValue = coordOut(profileAxis);
- el.setCenter(newCenterValue);
-
-// Now compute new error
-
- coordIn(profileAxis) = centerValue + centerValueErr;
- if (!cSys.convert (coordOut, coordIn, absIns, unitsIn, dopplerTypeIn,
- absOuts, unitsOut, dopplerTypeOut, offsetIn,
offsetOut)) {
- throw (AipsError(cSys.errorMessage()));
- }
- Double newCenterValueErr = abs(newCenterValue - coordOut(profileAxis));
-//
-// Convert width
-//
- if (unitIn==pix) {
- coordIn = cSys.referencePixel();
- unitsIn = pix;
- cSys.makePixelRelative(coordIn);
- } else {
- coordIn = cSys.referenceValue();
- unitsIn = cSys.worldAxisUnits();
- cSys.makeWorldRelative(coordIn);
- }
- unitsIn(profileAxis) = unitIn;
- absIns = False;
- absOuts = False;
- coordIn(profileAxis) = el.getFWHM();
- if (!cSys.convert (coordOut, coordIn, absIns, unitsIn, dopplerTypeIn,
- absOuts, unitsOut, dopplerTypeOut, 0.0, 0.0)) {
- os << cSys.errorMessage() << LogIO::EXCEPTION;
- }
-
-// Set new width
-
- el.setFWHM(abs(coordOut(profileAxis)));
-
-// Compute fractional errors
-
- el.get(pars);
- errors *= pars;
-
-// Overwrite position error
-
- errors(1) = newCenterValueErr;
-
-// Set errors
-
- el.setError(abs(errors));
-}
-
-
-void ImageProfileFit::fit (Bool fillRecord, RecordInterface& rec,
- Bool xAbsOut,
- const String& xUnitOut,
- const String& dopplerOut,
- ImageInterface<Float>* pFit,
- ImageInterface<Float>* pResid,
- Int order)
-
-{
- LogIO os(LogOrigin("ImageProfileFit", "fit", WHERE));
- if (!itsFitRegion) {
- os << "You cannot call this function as you are averaging all
profiles" << LogIO::EXCEPTION;
- }
-//
- IPosition inShape = itsImagePtr->shape();
- if (pFit!=0) {
- AlwaysAssert(inShape.isEqual(pFit->shape()), AipsError);
- }
- if (pResid!=0) {
- AlwaysAssert(inShape.isEqual(pResid->shape()), AipsError);
- }
-//
- IPosition inTileShape = itsImagePtr->niceCursorShape();
- TiledLineStepper stepper (itsImagePtr->shape(), inTileShape,
itsProfileAxis);
- RO_MaskedLatticeIterator<Float> inIter(*itsImagePtr, stepper);
-//
- PtrHolder<LatticeIterator<Float> > fitIter;
- PtrHolder<LatticeIterator<Bool> > fitMaskIter;
- PtrHolder<LatticeIterator<Float> > residIter;
- PtrHolder<LatticeIterator<Bool> > residMaskIter;
-//
- LatticeIterator<Float>* pFitIter = 0;
- LatticeIterator<Bool>* pFitMaskIter = 0;
- LatticeIterator<Float>* pResidIter = 0;
- LatticeIterator<Bool>* pResidMaskIter = 0;
- if (pFit) {
- fitIter.set(new LatticeIterator<Float>(*pFit, stepper));
- pFitIter = fitIter.ptr();
- if (pFit->hasPixelMask()) {
- fitMaskIter.set(new LatticeIterator<Bool>(pFit->pixelMask(),
stepper));
- pFitMaskIter = fitMaskIter.ptr();
- }
- }
- if (pResid) {
- residIter.set(new LatticeIterator<Float>(*pResid, stepper));
- pResidIter = residIter.ptr();
- if (pResid->hasPixelMask()) {
- residMaskIter.set(new LatticeIterator<Bool>(pResid->pixelMask(),
stepper));
- pResidMaskIter = residMaskIter.ptr();
- }
- }
-//
- Int nProfiles =
itsImagePtr->shape().product()/inIter.vectorCursor().nelements();
- ProgressMeter clock(0.0, Double(nProfiles), "Profile
fitting", "Profiles fitted",
- "", "", True, max(1,Int(nProfiles/20)));
- Double meterValue = 0.0;
-//
- Vector<Float> x(inShape(itsProfileAxis));
- Vector<Float> y(inShape(itsProfileAxis));
- for (uInt i=0; i<x.nelements(); i++) x[i] = i;
-
-// Fish out initial estimate (defines number of components to fit)
-// We just use it to get the number of components
-
- const SpectralList& l = itsSpectralFitPtr->list();
- uInt nEl = max(l.nelements(),uInt(1));
-
-// Make an auto-estimator, which we use for each spectrum
-// after the initial one
-
- SpectralEstimate estimator(nEl);
-
-// Make fitter and poly element if desired
-
- SpectralFit fitter;
- SpectralElement poly(order);
- SpectralList estimate;
-//
- Vector<Bool> inMask;
- Vector<Float> inSigma;
- Bool ok = False;
- uInt nFail = 0;
-//
- while (!inIter.atEnd()) {
-
-// Get data and mask (reflects pixelMask and region mask of SubImage)
-
- const Vector<Float>& data = inIter.vectorCursor();
- inMask = inIter.getMask(True);
-
-// Make estimate
-
- estimate = estimator.estimate(data);
-
-// Set fitter
-
- fitter.clear();
- fitter.addFitElement(estimate);
- if (order >= 0) fitter.addFitElement(poly);
-//
- if (itsSigmaImagePtr) {
- inSigma = itsSigmaImagePtr->getSlice(inIter.position(),
- inIter.cursorShape(), True);
- try {
- ok = fitter.fit(inSigma, data, x, inMask);
- } catch (AipsError x) {
- ok = False;
- }
-
- } else {
- try {
- ok = fitter.fit(data, x, inMask);
- } catch (AipsError x) {
- ok = False;
- }
- }
-
-// Evaluate
-
- SpectralList list(fitter.list());
- if (ok) {
- if (pFit) {
- list.evaluate(pFitIter->rwVectorCursor());
- }
- if (pFitMaskIter) {
- pFitMaskIter->rwVectorCursor() = inMask;
- }
- if (pResid) {
- pResidIter->rwVectorCursor() = data;
- list.residual(pResidIter->rwVectorCursor());
- }
- if (pResidMaskIter) {
- pResidMaskIter->rwVectorCursor() = inMask;
- }
- } else {
- nFail++;
- if (pFit) {
- pFitIter->rwVectorCursor() = 0.0;
- }
- if (pFitMaskIter) {
- pFitMaskIter->rwVectorCursor() = False;
- }
- if (pResid) {
- pResidIter->rwVectorCursor() = 0.0;
- }
- if (pResidMaskIter) {
- pResidMaskIter->rwVectorCursor() = False;
- }
***The diff for this file has been truncated for email.***
=======================================
--- /trunk/images/Images/ImageProfileFit.h Wed Apr 2 22:56:44 2008
+++ /dev/null
@@ -1,305 +0,0 @@
-//# ImageProfileFit.h: Class to fit profiles in images
-//# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This library is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU Library General Public License as published
by
-//# the Free Software Foundation; either version 2 of the License, or (at
your
-//# option) any later version.
-//#
-//# This library is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
-//# License for more details.
-//#
-//# You should have received a copy of the GNU Library General Public
License
-//# along with this library; if not, write to the Free Software Foundation,
-//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#ifndef IMAGES_IMAGEPROFILEFIT_H
-#define IMAGES_IMAGEPROFILEFIT_H
-
-//# Includes
-#include <casa/aips.h>
-#include <casa/Arrays/Vector.h>
-#include <coordinates/Coordinates/CoordinateSystem.h>
-#include <casa/Quanta/Quantum.h>
-#include <casa/Quanta/Unit.h>
-#include <measures/Measures/MDoppler.h>
-
-#include <components/SpectralComponents/SpectralElement.h>
-#include <components/SpectralComponents/SpectralFit.h>
-#include <casa/Utilities/PtrHolder.h>
-
-
-namespace casa { //# NAMESPACE CASA - BEGIN
-
-//# Forward declarations
-template<class T> class ImageInterface;
-template<class T> class MaskedLattice;
-class GlishRecord;
-class ImageRegion;
-class IPosition;
-class Slicer;
-class LogIO;
-
-// <summary>
-// Class to fit profiles in images
-// </summary>
-
-// <use visibility=export>
-
-// <reviewed reviewer="" date="" tests="">
-// </reviewed>
-
-// <prerequisite>
-// <li> <linkto class=></linkto>
-// </prerequisite>
-
-// <synopsis>
-// The Record used to contain models/fits is as follows. Other fields
-// will be ignored.
-//
-// <srcblock>
-// Field Type Description
-//-----------------------------------------------------------
-// xabs Bool Are the x-values absolute or
-// relative to the reference
pixel
-// xunit String The x unit
-// yunit String The y unit
-// doppler String doppler type
('radio', 'optical', 'true' etc)
-// Only required if x unit
consistent with m/s
-// elements Record Holds SpectralElement
descriptions.
-// elements[i] Record There will be N of these
where N is the
-// number of spectral
elements. Each elements[i]
-// holds a record description
that SpectralElement::fromRecord
-// can read or writes
-// .type String Type of SE (gaussian,
polynomial, etc)
-// .parameters Vector<Double> The parameters of the
element.
-// .errors Vector<Double> The errors of parameters of
the element.
-// .fixed Vector<Bool> Says whether parameter is
going to be held fixed when fitting
-// This field is not used by
SpectralElement. If its
-// not there, all parameters
are fitted for.
-//
-// </srcblock>
-// </synopsis>
-// <example>
-// <srcblock>
-// </srcblock>
-// </example>
-
-// <todo asof="2001/03/01">
-// </todo>
-
-class ImageProfileFit
-{
-public:
- // Constructor
- explicit ImageProfileFit();
-
- // Destructor
- ~ImageProfileFit();
-
- // Copy constructor. Uses copy semantics.
- ImageProfileFit(const ImageProfileFit& other);
-
- // Assignment operator. Uses copy semantics.
- ImageProfileFit& operator=(const ImageProfileFit& other);
-
- // Set data via an image. <src>profileAxis</src> specifies the profile
pixel
- // axis. You can either average the data over all other axes in the
- // image (<src>average=True</src>) or fit all profiles in the
- // image.
- //<group>
- void setData (const ImageInterface<Float>& image,
- const ImageRegion& region,
- uInt profileAxis, Bool average=True);
- void setData (const ImageInterface<Float>& image,
- const ImageInterface<Float>& sigma,
- const ImageRegion& region,
- uInt profileAxis, Bool average=True);
-//
- void setData (const ImageInterface<Float>& image,
- uInt profileAxis, Bool average=True);
- void setData (const ImageInterface<Float>& image,
- const ImageInterface<Float>& sigma,
- uInt profileAxis, Bool average=True);
-
- //</group>
-
- // Set the data directly, and provide a coordinate system
- // and specify the profile axis in the coordinate system.
- // The x-units can be 'pix'. If absolute
- // they must be 0-rel pixels. <src>isAbs</src> specifies whether
- // the coordinates are absolute or relative to the reference pixel.
- // If the weights vector, <src>sigma</src> is of zero length,
- // it is treated as all unity.
- //<group>
- void setData (const CoordinateSystem& cSys,
- uInt ProfileAxis,
- const Quantum<Vector<Float> >& x,
- const Quantum<Vector<Float> >& y,
- const Vector<Float>& sigma,
- Bool isAbs, const String& doppler);
-
- void setData (const CoordinateSystem& cSys,
- uInt ProfileAxis,
- const Quantum<Vector<Float> >& x,
- const Quantum<Vector<Float> >& y,
- const Vector<Bool>& mask,
- const Vector<Float>& sigma,
- Bool isAbs, const String& doppler);
- //</group>
-
- // Makes an estimate from the set data. This means it generates
SpectralElements
- // holding the estimate, and replaces all elements you might have
- // put in place with function <src>addElement</src>. Returns
- // False if it can't find any elements.
- Bool estimate (uInt nMax = 0);
-
- // Decode the Glish record holding SpectralElements and add
- // them to the fitter. Absolute pixel coordinate units are assumed
to be
- // 1-rel on input. Return the number of the element added.
- uInt addElements (const RecordInterface& estimate);
-
- // Gets the internal SpectralElements (either estimate or fit
- // depending on what function you called last) into a record.
- // Only returns False if the field is already defined. Absolute pixel
- // coordinate units are 1-rel on output.
- Bool getList (RecordInterface& rec,
- Bool xAbsOut,
- const String& xUnitOut,
- const String& dopplerOut);
-
- // Gets the internal SpectralElements (either estimate or fit
- // depending on what function you called last) into a SpectralList
- // Only returns False if the field is already defined. Absolute pixel
- // coordinate units are 1-rel on output.
- SpectralList getList () const;
-
- // Reset the internal list of SpectralElements to null
- void reset ();
-
- // Return number of SpectralElements set
- uInt nElements ();
-
- // Do the fit of the averaged profile. Specify the
- // order of the baseline you would also like to fit for.
- //<group>
- Bool fit(Int order=-1);
- //</group>
-
- // Fit all profiles with shapes + optional polynomial in the region
and write out images.
- // If fillRecord is True, the output record is filled with the the
parameters of
- // every fit. This can get VERY large so use with care. If the
output images
- // have a writable mask, the input mask is transferred to the output.
- //<group>
- void fit (Bool fillRecord, RecordInterface& rec,
- Bool xAbsRec,
- const String& xUnitRec,
- const String& dopplerRec,
- ImageInterface<Float>* pFit,
- ImageInterface<Float>* pResid,
- Int order=-1);
- //</group>
-
- // Find the residuals (fit or estimate) of the averaged profile
- //<group>
- void residual(Vector<Float>& resid, const Vector<Float>& x) const;
- void residual(Vector<Float>& resid) const;
- //</group>
-
- // Find the model (fit or estimate) of the averaged profile
- //<group>
- void model (Vector<Float>& model, const Vector<Float>& x) const;
- void model (Vector<Float>& model) const;
- //</group>
-
-private:
-
-// Images holding data and weights. Will be set only if fitting all
profiles in region
-
- ImageInterface<Float>* itsImagePtr;
- ImageInterface<Float>* itsSigmaImagePtr;
-//
- Bool itsFitDone;
-
-// Holds the abcissa, ordinate and mask. x-units will be pixels
-// if data source is an image, else as specified in setData
-
- Quantum<Vector<Float> > itsX, itsY;
- Vector<Bool> itsMask;
- Vector<Float> itsSigma;
-
-// The fitters. The first one does not have a polynomial in it
-// The second one may.
- SpectralFit* itsSpectralFitPtr;
- SpectralFit itsSpectralFitter;
-
-// The coordinate system if the data source was an image
-// itsProfileAxis specified the profile axis in the image
-// Will be -1 if data source was just vectors
-
- CoordinateSystem itsCoords;
- Int itsProfileAxis;
-
-// If the data source was an image, these give the doppler type
-// (if any) and x-unit IF an estimate was specified by the user via
-// function addElements. They are used in function getElements
-// so that the output record has the same units.
-
- String itsDoppler; // The doppler of the data source
- Bool itsXAbs; // Is the data source in absolute
coordinates ?
- Bool itsFitRegion; // True to fit all profiles in region
-//
- void collapse (Vector<Float>& profile, Vector<Bool>& mask,
- uInt profileAxis, const MaskedLattice<Float>&
lat) const;
-
-// Convert SE
- SpectralElement convertSpectralElement (const SpectralElement& elIn,
- Bool xAbsIn, Bool xAbsOut,
- const String& xUnitIn,
- const String& xUnitOut,
- const String& dopplerIn,
- const String& dopplerOut,
- const String& yUnitIn,
- const String& yUnitOut);
-
-// Convert Gaussian model x-units when data source is an image
- void convertGaussianElementX (SpectralElement& el,
- CoordinateSystem& cSys,
- uInt profileAxis,
- Bool absIn, Bool absOut,
- const String& unitIn,
- const String& unitOut,
- const String& dopplerIn,
- const String& dopplerOut);
-
-//
- SpectralList filterList (const SpectralList& listIn) const;
-
-//
- Bool getElements (RecordInterface& rec,
- Bool xAbsOut,
- const String& xUnitOut,
- const String& dopplerOut,
- const SpectralList& list);
-//
- void setData (const ImageInterface<Float>*& pSigma,
- const ImageInterface<Float>& image,
- const Slicer& sl, Bool average);
-};
-
-
-} //# NAMESPACE CASA - END
-
-#endif
=======================================
--- /trunk/images/Images/ImageProfileFitter.cc Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,756 +0,0 @@
-//# tSubImage.cc: Test program for class SubImage
-//# Copyright (C) 1998,1999,2000,2001,2003
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or (at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#include <images/Images/ImageProfileFitter.h>
-
-#include <casa/Quanta/MVAngle.h>
-#include <casa/Quanta/MVTime.h>
-#include <casa/Utilities/Precision.h>
-#include <images/Images/ImageAnalysis.h>
-#include <images/Images/ImageCollapser.h>
-#include <images/Images/PagedImage.h>
-#include <images/Images/TempImage.h>
-#include <scimath/Mathematics/Combinatorics.h>
-
-namespace casa {
-
-ImageProfileFitter::ImageProfileFitter(
- const String& imagename, const String& region, const String& box,
- const String& chans, const String& stokes,
- const String& mask, const Int axis, const Bool multiFit,
- const String& residual, const String& model, const uInt ngauss,
- const Int polyOrder, const String& ampName,
- const String& ampErrName, const String& centerName,
- const String& centerErrName, const String& fwhmName,
- const String& fwhmErrName
-) : _log(new LogIO()), _image(0),
- _regionName(region), _box(box), _chans(chans), _stokes(stokes),
- _mask(mask), _residual(residual),
- _model(model), _regionString(""), _xUnit(""),
- _centerName(centerName), _centerErrName(centerErrName),
- _fwhmName(fwhmName), _fwhmErrName(fwhmErrName),
- _ampName(ampName), _ampErrName(ampErrName),
- _multiFit(multiFit), _deleteImageOnDestruct(True),
- _polyOrder(polyOrder), _fitAxis(axis), _ngauss(ngauss),
- _results(Record()) {
- _construct(imagename);
-}
-
-ImageProfileFitter::ImageProfileFitter(
- const ImageInterface<Float> * const image, const String& region,
- const String& box, const String& chans, const String& stokes,
- const String& mask, const Int axis, const Bool multiFit,
- const String& residual, const String& model, const uInt ngauss,
- const Int polyOrder, const String& ampName,
- const String& ampErrName, const String& centerName,
- const String& centerErrName, const String& fwhmName,
- const String& fwhmErrName
-) : _log(new LogIO()), _image(image->cloneII()),
- _regionName(region), _box(box), _chans(chans), _stokes(stokes),
- _mask(mask), _residual(residual),
- _model(model), _regionString(""), _xUnit(""),
- _centerName(centerName), _centerErrName(centerErrName),
- _fwhmName(fwhmName), _fwhmErrName(fwhmErrName),
- _ampName(ampName), _ampErrName(ampErrName),
- _multiFit(multiFit), _deleteImageOnDestruct(False),
- _polyOrder(polyOrder), _fitAxis(axis), _ngauss(ngauss),
- _results(Record()) {
- _construct(_image);
-}
-
-ImageProfileFitter::~ImageProfileFitter() {
- delete _log;
- delete _image;
-}
-
-Record ImageProfileFitter::fit() {
- LogOrigin logOrigin("ImageProfileFitter", __FUNCTION__);
- *_log << logOrigin;
-
- ImageAnalysis ia(_image);
- Record estimate;
-
- String weightsImageName = "";
- try {
- if (_multiFit) {
- // FIXME need to be able to specify the weights image
- _fitters = ia.fitallprofiles(
- _regionRecord, _subImage, _xUnit, _fitAxis, _mask,
- _ngauss, _polyOrder, weightsImageName,
- _model, _residual
- );
- }
- else {
- ImageFit1D<Float> fitter = ia.fitprofile(
- _regionRecord, _subImage, _xUnit, _fitAxis, _mask,
- estimate, _ngauss, _polyOrder, _model,
- _residual
- );
- Vector<uInt> axes(1, _fitAxis);
- ImageCollapser collapser(
- &_subImage, axes, True,
- ImageCollapser::MEAN, "", True
- );
- ImageInterface<Float> *x =
static_cast<SubImage<Float>*>(collapser.collapse(True));
- _subImage = (SubImage<Float>)(*x);
- delete x;
- _fitters.resize(1);
- _fitters[0] = fitter;
- }
- }
- catch (AipsError exc) {
- *_log << "Exception during fit: " << exc.getMesg()
- << LogIO::EXCEPTION;
- }
- _setResults();
- *_log << LogIO::NORMAL << _resultsToString() << LogIO::POST;
- return _results;
-}
-
-Record ImageProfileFitter::getResults() const {
- return _results;
-}
-
-void ImageProfileFitter::_construct(const String& imagename) {
- LogOrigin logOrigin("ImageProfileFitter", __FUNCTION__);
- *_log << logOrigin;
-
- _checkNGaussAndPolyOrder();
- ImageInputProcessor inputProcessor;
- Vector<ImageInputProcessor::OutputStruct> outputStruct;
- _getOutputStruct(outputStruct);
- Vector<ImageInputProcessor::OutputStruct>* outputPtr =
outputStruct.size() == 0
- ? 0
- : &outputStruct;
- String diagnostics;
- inputProcessor.process(
- _image, _regionRecord, diagnostics,
- outputPtr, _stokes, imagename, 0,
- _regionName, _box, _chans,
- RegionManager::USE_FIRST_STOKES,
- False, 0
- );
- _finishConstruction();
-}
-
-void ImageProfileFitter::_construct(const ImageInterface<Float>* image) {
- LogOrigin logOrigin("ImageProfileFitter", __FUNCTION__);
- *_log << logOrigin;
- _checkNGaussAndPolyOrder();
- ImageInputProcessor inputProcessor;
- Vector<ImageInputProcessor::OutputStruct> outputStruct;
- _getOutputStruct(outputStruct);
- Vector<ImageInputProcessor::OutputStruct>* outputPtr =
outputStruct.size() == 0
- ? 0
- : &outputStruct;
- String diagnostics;
- inputProcessor.process(
- _regionRecord, diagnostics,
- outputPtr, _stokes, image, 0,
- _regionName, _box, _chans,
- RegionManager::USE_FIRST_STOKES,
- False, 0
- );
- _finishConstruction();
-}
-
-void ImageProfileFitter::_getOutputStruct(
- Vector<ImageInputProcessor::OutputStruct>& outputs
-) {
- outputs.resize(0);
- if (! _model.empty()) {
- ImageInputProcessor::OutputStruct modelImage;
- modelImage.label = "model image";
- modelImage.outputFile = &_model;
- modelImage.required = True;
- modelImage.replaceable = False;
- outputs.resize(1);
- outputs[0] = modelImage;
- }
- if (! _residual.empty()) {
- ImageInputProcessor::OutputStruct residImage;
- residImage.label = "residual image";
- residImage.outputFile = &_residual;
- residImage.required = True;
- residImage.replaceable = False;
- outputs.resize(outputs.size() + 1, True);
- outputs[outputs.size() - 1] = residImage;
- }
-}
-
-void ImageProfileFitter::_checkNGaussAndPolyOrder() const {
- if (_ngauss == 0 && _polyOrder < 0) {
- *_log << "Number of gaussians is 0 and polynomial order is less than
zero. "
- << "According to these inputs there is nothing to fit."
- << LogIO::EXCEPTION;
- }
-}
-
-void ImageProfileFitter::_finishConstruction() {
- if (! _box.empty() && ! _regionName.empty()) {
- // for output later
- _regionName = "";
- }
- if (_fitAxis >= (Int)_image->ndim()) {
- *_log << "Specified fit axis " << _fitAxis
- << " must be less than the number of image axes ("
- << _image->ndim() << ")" << LogIO::EXCEPTION;
- }
- if (_fitAxis < 0) {
- Int specCoord =
_image->coordinates().findCoordinate(Coordinate::SPECTRAL);
- if (specCoord < 0) {
- _fitAxis = 0;
- *_log << LogIO::WARN << "No spectral coordinate found in image, "
- << "using axis 0 as fit axis" << LogIO::POST;
- }
- else {
- _fitAxis = _image->coordinates().pixelAxes(specCoord)[0];
- *_log << LogIO::NORMAL << "Using spectral axis (axis " << _fitAxis
- << ") as fit axis" << LogIO::POST;
- }
- }
-}
-
-void ImageProfileFitter::_setResults() {
- uInt nComps = _polyOrder < 0 ? _ngauss : _ngauss + 1;
- Vector<Bool> convergedArr(_fitters.size());
- Vector<Int> niterArr(_fitters.size(), 0);
- Matrix<Double> centerMat(_fitters.size(), nComps, -1);
- Matrix<Double> fwhmMat(_fitters.size(), nComps, -1);
- Matrix<Double> ampMat(_fitters.size(), nComps, -1);
- Matrix<Double> centerErrMat(_fitters.size(), nComps, -1);
- Matrix<Double> fwhmErrMat(_fitters.size(), nComps, -1);
- Matrix<Double> ampErrMat(_fitters.size(), nComps, -1);
- Matrix<String> typeMat(_fitters.size(), nComps, "");
-
- Vector<Int> nCompArr(_fitters.size(), 0);
-
- IPosition inTileShape = _subImage.niceCursorShape();
- TiledLineStepper stepper (_subImage.shape(), inTileShape, _fitAxis);
- RO_MaskedLatticeIterator<Float> inIter(_subImage, stepper);
- Vector<ImageFit1D<Float> >::const_iterator fitter = _fitters.begin();
- SpectralList solutions;
-
- uInt count = 0;
- for (
- inIter.reset();
- ! inIter.atEnd() && fitter != _fitters.end();
- inIter++, fitter++, count++
- ) {
- convergedArr[count] = fitter->converged();
- if (fitter->converged()) {
- niterArr[count] = (Int)fitter->getNumberIterations();
- nCompArr[count] = (Int)fitter->getList().nelements();
- solutions = fitter->getList();
- for (uInt i=0; i<solutions.nelements(); i++) {
- typeMat(count, i) = SpectralElement::fromType(solutions[i].getType());
- if (solutions[i].getType() == SpectralElement::GAUSSIAN) {
- centerMat(count, i) = solutions[i].getCenter();
- fwhmMat(count, i) = solutions[i].getFWHM();
- ampMat(count, i) = solutions[i].getFWHM();
- ampMat(count, i) = solutions[i].getAmpl();
- centerErrMat(count, i) = solutions[i].getCenterErr();
- fwhmErrMat(count, i) = solutions[i].getFWHMErr();
- ampErrMat(count, i) = solutions[i].getAmplErr();
- }
- }
- }
- }
- _results.define("converged", convergedArr);
- _results.define("niter", niterArr);
- _results.define("ncomps", nCompArr);
- _results.define("xUnit", _xUnit);
- _results.define("yUnit", _image->units().getName());
-
- String key;
- TempImage<Float> *tmp = static_cast<TempImage<Float>*
>(_image->cloneII());
- Vector<uInt> axes(1, _fitAxis);
- ImageCollapser collapser(
- tmp, axes, False, ImageCollapser::ZERO, String(""), False
- );
- TempImage<Float> *myTemplate = static_cast<TempImage<Float>*
>(collapser.collapse(True));
- delete tmp;
- IPosition shape = myTemplate->shape();
- CoordinateSystem csys = myTemplate->coordinates();
- delete myTemplate;
- uInt gaussCount = 0;
- if (
- ! _multiFit && (
- ! _centerName.empty() || ! _centerErrName.empty()
- || ! _fwhmName.empty() || ! _fwhmErrName.empty()
- || ! _ampName.empty() || ! _ampErrName.empty()
- )
- ) {
- *_log << LogIO::WARN << "This was not a multi-pixel fit request so
solution "
- << "images will not be written" << LogIO::POST;
- }
-
- for (uInt i=0; i<nComps; i++) {
- String num = String::toString(i);
- if (_multiFit && solutions[i].getType() == SpectralElement::GAUSSIAN) {
- String gnum = String::toString(gaussCount);
- String mUnit = _xUnit;
- if (!_centerName.empty()) {
- _makeSolutionImage(
- _centerName + "_" + gnum, shape,
- csys, centerMat.column(i), mUnit
- );
- }
- if (!_centerErrName.empty()) {
- _makeSolutionImage(
- _centerErrName + "_" + gnum, shape,
- csys, centerErrMat.column(i), mUnit
- );
- }
- if (!_fwhmName.empty()) {
- _makeSolutionImage(
- _fwhmName + "_" + gnum, shape,
- csys, fwhmMat.column(i), mUnit
- );
- }
- if (!_fwhmErrName.empty()) {
- _makeSolutionImage(
- _fwhmErrName + "_" + gnum, shape,
- csys, fwhmErrMat.column(i), mUnit
- );
- }
- mUnit = _image->units().getName();
- if (!_ampName.empty()) {
- _makeSolutionImage(
- _ampName + "_" + gnum, shape,
- csys, ampMat.column(i), mUnit
- );
- }
- if (!_ampErrName.empty()) {
- _makeSolutionImage(
- _ampErrName + "_" + gnum, shape,
- csys, ampErrMat.column(i), mUnit
- );
- }
- gaussCount++;
- }
- key = "center" + num;
- _results.define(key, centerMat.column(i));
- key = "fwhm" + num;
- _results.define(key, fwhmMat.column(i));
- key = "amp" + num;
- _results.define(key, ampMat.column(i));
- key = "centerErr" + num;
- _results.define(key, centerErrMat.column(i));
- key = "fwhmErr" + num;
- _results.define(key, fwhmErrMat.column(i));
- key = "ampErr" + num;
- _results.define(key, ampErrMat.column(i));
- key = "type" + num;
- _results.define(key, typeMat.column(i));
- }
-}
-
-String ImageProfileFitter::_radToRa(Float ras) const{
-
- Int h, m;
- Float rah = ras * 12 / C::pi;
- h = (int)floor(rah);
- Float ram = (rah - h) * 60;
- m = (int)floor(ram);
- ras = (ram - m) * 60;
- ras = (int)(1000 * ras) / 1000.;
-
- String raStr = (h < 10) ? "0" : "";
- raStr.append(String::toString(h)).append(String(":"))
- .append(String((m < 10) ? "0" : ""))
- .append(String::toString(m)).append(String(":"))
- .append(String((ras < 10) ? "0" : ""))
- .append(String::toString(ras));
-
- return raStr;
-
-}
-
-
-String ImageProfileFitter::_resultsToString() const {
- ostringstream summary;
- summary << "****** Fit performed at " << Time().toString() << "******" <<
endl << endl;
- summary << "Input parameters ---" << endl;
- summary << " --- imagename: " << _image->name() << endl;
- summary << " --- region: " << _regionName << endl;
- summary << " --- box: " << _box << endl;
- summary << " --- channels: " << _chans << endl;
- summary << " --- stokes: " << _stokes << endl;
- summary << " --- mask: " << _mask << endl;
- summary << " --- polynomial order: " << _polyOrder << endl;
- summary << " --- number of Gaussians: " << _ngauss << endl;
-
- if (_multiFit) {
- summary << " --- Multiple spectra fit, one per pixel over selected
region" << endl;
- }
- else {
- summary << " --- One spectrum fit, averaged over several pixels if
necessary" << endl;
- }
-
- IPosition inTileShape = _subImage.niceCursorShape();
- TiledLineStepper stepper (_subImage.shape(), inTileShape, _fitAxis);
- RO_MaskedLatticeIterator<Float> inIter(_subImage, stepper);
- CoordinateSystem csysSub = _subImage.coordinates();
- CoordinateSystem csys = _image->coordinates();
- Vector<Double> worldStart;
- if (! csysSub.toWorld(worldStart, inIter.position())) {
- *_log << csysSub.errorMessage() << LogIO::EXCEPTION;
- }
- CoordinateSystem csysIm = _image->coordinates();
- Vector<Double> pixStart;
- if (! csysIm.toPixel(pixStart, worldStart)) {
- *_log << csysIm.errorMessage() << LogIO::EXCEPTION;
- }
- if (_multiFit) {
- for (uInt i=0; i<pixStart.size(); i++) {
- pixStart[i] = (Int)std::floor( pixStart[i] + 0.5);
- }
- }
- Vector<ImageFit1D<Float> >::const_iterator fitter = _fitters.begin();
- Vector<String> axesNames = csysSub.worldAxisNames();
- Vector<Double> imPix(pixStart.size());
- Vector<Double> world;
- IPosition subimPos;
- SpectralList solutions;
- String axisUnit = csysIm.worldAxisUnits()[_fitAxis];
- Int pixPrecision = _multiFit ? 0 : 3;
- Int basePrecision = summary.precision(1);
- summary.precision(basePrecision);
- for (
- Vector<String>::iterator iter = axesNames.begin();
- iter != axesNames.end(); iter++
- ) {
- iter->upcase();
- }
- for (
- inIter.reset();
- ! inIter.atEnd() && fitter != _fitters.end();
- inIter++, fitter++
- ) {
- subimPos = inIter.position();
- if (csysSub.toWorld(world, subimPos)) {
- summary << "Fit centered at:" << endl;
- for (uInt i=0; i<world.size(); i++) {
- if ((Int)i != _fitAxis) {
- //summary << "ax=" << axesNames[i]
- // << " world=" << world[i]
- // << endl;
- if (axesNames[i].startsWith("RIG")) {
- // right ascension
- summary << " RA : "
- // << MVTime(world[i]).string(MVTime::TIME, 9) << endl;
- << _radToRa(world[i]) << endl;
- }
- else if (axesNames[i].startsWith("DEC")) {
- // declination
- summary << " Dec : "
- << MVAngle(world[i]).string(MVAngle::ANGLE_CLEAN, 8) << endl;
- }
- else if (axesNames[i].startsWith("FREQ")) {
- // frequency
- summary << " Freq : "
- << world[i]
- << csysSub.spectralCoordinate(i).formatUnit() << endl;
- }
- else if (axesNames[i].startsWith("STO")) {
- // stokes
- summary << " Stokes : "
- << Stokes::name(Stokes::type((Int)world[i])) << endl;
- }
- }
- }
- }
- else {
- *_log << csysSub.errorMessage() << LogIO::EXCEPTION;
- }
- for (uInt i=0; i<pixStart.size(); i++) {
- imPix[i] = pixStart[i] + subimPos[i];
- }
-
- summary.setf(ios::fixed);
- summary << setprecision(pixPrecision) << " Pixel : [";
- for (uInt i=0; i<imPix.size(); i++) {
- if (i == (uInt)_fitAxis) {
- summary << " *";
- }
- else {
- summary << imPix[i];
- }
- if (i != imPix.size()-1) {
- summary << ", ";
- }
- }
- summary << "]" << setprecision(basePrecision) << endl;
- summary.unsetf(ios::fixed);
- String converged = fitter->converged() ? "YES" : "NO";
- summary << " Converged : " << converged << endl;
- if (fitter->converged()) {
- solutions = fitter->getList();
- summary << " Iterations : " << fitter->getNumberIterations() << endl;
- for (uInt i=0; i<solutions.nelements(); i++) {
- summary << " Results for component " << i << ":" << endl;
- if (solutions[i].getType() == SpectralElement::GAUSSIAN) {
- summary << _gaussianToString(
- solutions[i], csys, world.copy()
- );
- }
- else if (solutions[i].getType() == SpectralElement::POLYNOMIAL) {
- summary << _polynomialToString(solutions[i], csys, imPix, world);
- }
- }
- }
- summary << endl;
- }
- return summary.str();
-}
-
-String ImageProfileFitter::_elementToString(
- const Double value, const Double error,
- const String& unit
-) const {
- ostringstream out;
-
- Unit myUnit(unit);
- Vector<String> unitPrefix;
- String outUnit;
- Quantity qVal(value, unit);
- Quantity qErr(error, unit);
-
- if (myUnit.getValue() == UnitVal::ANGLE) {
- Vector<String> angUnits(5);
- angUnits[0] = "deg";
- angUnits[1] = "arcmin";
- angUnits[2] = "arcsec";
- angUnits[3] = "marcsec";
- angUnits[4] = "uarcsec";
- for (uInt i=0; i<angUnits.size(); i++) {
- outUnit = angUnits[i];
- if (fabs(qVal.getValue(outUnit)) > 1) {
- qVal.convert(outUnit);
- qErr.convert(outUnit);
- break;
- }
- }
- }
- else if (unit.empty() || Quantity(1,
myUnit).isConform(Quantity(1, "m/s"))) {
- // do nothing
- }
- else {
- Vector<String> unitPrefix(10);
- unitPrefix[0] = "T";
- unitPrefix[1] = "G";
- unitPrefix[2] = "M";
- unitPrefix[3] = "k";
- unitPrefix[4] = "";
- unitPrefix[5] = "m";
- unitPrefix[6] = "u";
- unitPrefix[7] = "n";
- unitPrefix[8] = "p";
- unitPrefix[9] = "f";
-
- for (uInt i=0; i<unitPrefix.size(); i++) {
- outUnit = unitPrefix[i] + unit;
- if (fabs(qVal.getValue(outUnit)) > 1) {
- qVal.convert(outUnit);
- qErr.convert(outUnit);
- break;
- }
- }
- }
- Vector<Double> valErr(2);
- valErr[0] = qVal.getValue();
- valErr[1] = qErr.getValue();
-
- uInt precision = precisionForValueErrorPairs(valErr, Vector<Double>());
- out << std::fixed << setprecision(precision);
- out << qVal.getValue() << " +/- " << qErr.getValue()
- << " " << qVal.getUnit();
- return out.str();
-}
-
-String ImageProfileFitter::_gaussianToString(
- const SpectralElement& gauss, const CoordinateSystem& csys,
- const Vector<Double> world
-) const {
- Vector<Double> myWorld = world;
- String yUnit = _image->units().getName();
- ostringstream summary;
- summary << " Type : GAUSSIAN" << endl;
- summary << " Peak : "
- << _elementToString(
- gauss.getAmpl(), gauss.getAmplErr(), yUnit
- )
- << endl;
- Double center = gauss.getCenter();
- Double centerErr = gauss.getCenterErr();
- Double fwhm = gauss.getFWHM();
- Double fwhmErr = gauss.getFWHMErr();
-
- Double pCenter = 0;
- Double pCenterErr = 0;
- Double pFWHM = 0;
- Double pFWHMErr = 0;
- Int specCoordIndex = csys.findCoordinate(Coordinate::SPECTRAL);
- Bool convertedCenterToPix = True;
- Bool convertedFWHMToPix = True;
-
- if (
- specCoordIndex >= 0
- && _fitAxis == csys.pixelAxes(specCoordIndex)[0]
- && ! csys.spectralCoordinate(specCoordIndex).velocityUnit().empty()
- ) {;
- if (csys.spectralCoordinate(specCoordIndex).velocityToPixel(pCenter,
center)) {
- Double nextVel;
- csys.spectralCoordinate(specCoordIndex).pixelToVelocity(nextVel,
pCenter+1);
- Double velInc = fabs(center - nextVel);
- pCenterErr = centerErr/velInc;
- pFWHM = fwhm/velInc;
- pFWHMErr = fwhmErr/velInc;
- }
- else {
- convertedCenterToPix = False;
- convertedFWHMToPix = False;
- }
- }
- else {
- Vector<Double> pixel(myWorld.size());
- myWorld[_fitAxis] = center;
- Double delta = csys.increment()[_fitAxis];
- if (csys.toPixel(pixel, myWorld)) {
- pCenter = pixel[_fitAxis];
- pCenterErr = centerErr/delta;
- }
- else {
- convertedCenterToPix = False;
- }
- pFWHM = fwhm/delta;
- pFWHMErr = fwhmErr/delta;
- }
- summary << " Center : "
- << _elementToString(
- center, centerErr, _xUnit
- )
- << endl;
- if (convertedCenterToPix) {
- summary << " "
- << _elementToString(
- pCenter, pCenterErr, "pixel"
- )
- << endl;
- }
- else {
- summary << " Could not convert world to pixel for center"
<< endl;
- }
- summary << " FWHM : "
- << _elementToString(
- fwhm, fwhmErr, _xUnit
- )
- << endl;
- if (convertedFWHMToPix) {
- summary << " " << _elementToString(
- pFWHM, pFWHMErr, "pixel"
- )
- << endl;
- }
- else {
- summary << " Could not convert FWHM to pixel" << endl;
- }
- return summary.str();
-}
-
-String ImageProfileFitter::_polynomialToString(
- const SpectralElement& poly, const CoordinateSystem& csys,
- const Vector<Double> imPix, const Vector<Double> world
-) const {
- ostringstream summary;
- summary << " Type: POLYNOMIAL" << endl;
- Vector<Double> parms, errs;
- poly.get(parms);
- poly.getError(errs);
- for (uInt j=0; j<parms.size(); j++) {
- String unit = _image->units().getName();
- if (j > 0) {
- String denom = _xUnit.find("/") != String::npos
- ? "(" + _xUnit + ")"
- : _xUnit;
- unit = unit + "/" + denom + "^" + String::toString(j);
- }
- summary << " c" << j << " : "
- << _elementToString(parms[j], errs[j], unit) << endl;
- }
- // coefficents in pixels
- Double x0, deltaX;
- if (Quantity(1,_xUnit).isConform(Quantity(1, "m/s"))) {
- Double x1;
-
csys.spectralCoordinate(csys.findCoordinate(Coordinate::SPECTRAL)).pixelToVelocity(x0,
0);
-
csys.spectralCoordinate(csys.findCoordinate(Coordinate::SPECTRAL)).pixelToVelocity(x1,
1);
- deltaX = x1 - x0;
- }
- else {
- Vector<Double> p0 = imPix;
- p0[_fitAxis] = 0;
- Vector<Double> world0 = world;
- csys.toWorld(world0, p0);
- x0 = world0[_fitAxis];
- deltaX = csys.increment()[_fitAxis];
- }
- Vector<Double> pCoeff(_polyOrder + 1, 0);
- Vector<Double> pCoeffErr(_polyOrder + 1, 0);
- for (Int j=0; j<=_polyOrder; j++) {
- Double sumsq = 0;
- for (Int k=j; k<=_polyOrder; k++) {
- Double multiplier = Combinatorics::choose(k, j)
- * casa::pow(x0, Float(k-j))
- * casa::pow(deltaX, Float(j));
-
- pCoeff[j] += multiplier * parms[k];
- Double errCoeff = multiplier * errs[k];
- sumsq += errCoeff * errCoeff;
- }
- pCoeffErr[j] = casa::sqrt(sumsq);
- summary << " c" << j << " : ";
- String unit = _image->units().getName() + "/(pixel)^" +
String::toString(j);
- summary << _elementToString(pCoeff[j], pCoeffErr[j], unit) << endl;
- }
- return summary.str();
-}
-
-void ImageProfileFitter::_makeSolutionImage(
- const String& name, const IPosition& shape, const CoordinateSystem& csys,
- const Vector<Double>& values, const String& unit
-) {
- Vector<Float> tmpVec(values.size());
- for (uInt i=0; i<tmpVec.size(); i++) {
- tmpVec[i] = values[i];
- }
- PagedImage<Float> image( shape, csys, name);
- image.put(tmpVec.reform(shape));
- image.setUnits(Unit(unit));
-}
-}
-
=======================================
--- /trunk/images/Images/ImageProfileFitter.h Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,176 +0,0 @@
-//# tSubImage.cc: Test program for class SubImage
-//# Copyright (C) 1998,1999,2000,2001,2003
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or (at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#ifndef IMAGES_IMAGEPROFILEFITTER_H
-#define IMAGES_IMAGEPROFILEFITTER_H
-
-#include <casa/Logging/LogIO.h>
-#include <images/Images/ImageFit1D.h>
-#include <images/Images/ImageInterface.h>
-#include <images/Images/ImageInputProcessor.h>
-
-#include <casa/namespace.h>
-
-namespace casa {
-
-class ImageProfileFitter {
- // <summary>
- // Top level interface for one-dimensional profile fits.
- // </summary>
-
- // <reviewed reviewer="" date="" tests="" demos="">
- // </reviewed>
-
- // <prerequisite>
- // </prerequisite>
-
- // <etymology>
- // Fits one-dimensional profiles.
- // </etymology>
-
- // <synopsis>
- // Top level interface for one-dimensional profile fits.
- // </synopsis>
-
- // <example>
- // <srcblock>
- // ImageProfileFitter fitter(...)
- // fitter.fit()
- // </srcblock>
- // </example>
-
-public:
- // constructor appropriate for API calls.
- // Parameters:
- // <src>imagename</src> - the name of the input image in which to fit the
models
- // <src>box</src> - A 2-D rectangular box in direction space in which to
use pixels for the fitting, eg box=100,120,200,230
- // In cases where both box and region are specified, box, not region, is
used.
- // <src>region</src> - Named region to use for fitting
- // <src>chans</src> - Zero-based channel range on which to do the fit.
- // <src>stokes</src> - Stokes plane on which to do the fit. Only a single
Stokes parameter can be
- // specified.
- // <src>mask</src> - Mask (as LEL) to use as a way to specify which
pixels to use
- // <src>axis</src> - axis along which to do the fits. If <0, use spectral
axis, and if no spectral
- // axis, use zeroth axis.
- // <src>multiFit</src> True means do fit on each pixel in the specified
region or box. False means
- // average together the pixels in the specified box or region and do a
single fit to that average profile.
- // <src>residual</src> - Name of residual image to save. Blank means do
not save residual image.
- // <src>model</src> - Name of the model image to save. Blank means do not
save model image.
- // The output solution images are only written if multiFit is true.
-
- ImageProfileFitter(
- const String& imagename, const String& region, const String& box,
- const String& chans, const String& stokes, const String& mask,
- const Int axis, const Bool multiFit, const String& residual,
- const String& model, const uInt ngauss, const Int polyOrder,
- const String& ampName = "", const String& ampErrName = "",
- const String& centerName = "", const String& centerErrName = "",
- const String& fwhmName = "", const String& fwhmErrName = ""
- );
-
- ImageProfileFitter(
- const ImageInterface<Float> * const image, const String& region,
- const String& box, const String& chans, const String& stokes,
- const String& mask, const Int axis, const Bool multiFit,
- const String& residual, const String& model, const uInt ngauss,
- const Int polyOrder, const String& ampName = "",
- const String& ampErrName = "", const String& centerName = "",
- const String& centerErrName = "", const String& fwhmName = "",
- const String& fwhmErrName = ""
- );
-
- // destructor
- ~ImageProfileFitter();
-
- // Do the fit.
- Record fit();
-
- // get the fit results
- Record getResults() const;
-
-
-private:
- LogIO *_log;
- ImageInterface<Float> *_image;
- Record _regionRecord;
- String _regionName, _box, _chans, _stokes, _mask,
- _residual, _model, _regionString, _xUnit,
- _centerName, _centerErrName, _fwhmName,
- _fwhmErrName, _ampName, _ampErrName;
- Bool _logfileAppend, _fitConverged, _fitDone, _multiFit,
_deleteImageOnDestruct;
- Int _polyOrder, _fitAxis;
- uInt _ngauss;
- Vector<ImageFit1D<Float> > _fitters;
- SubImage<Float> _subImage;
- Record _results;
-
- // disallow default constructor
- ImageProfileFitter();
-
- // does the lion's share of constructing the object, ie checks validity of
- // inputs, etc.
-
- void _construct(const String& imagename);
-
- void _construct(const ImageInterface<Float>* image);
-
- void _getOutputStruct(
- Vector<ImageInputProcessor::OutputStruct>& outputs
- );
-
- void _checkNGaussAndPolyOrder() const;
-
- void _finishConstruction();
-
- void _setResults();
-
- String _radToRa(const Float ras) const;
- String _resultsToString() const;
-
- String _elementToString(
- const Double value, const Double error,
- const String& unit
- ) const;
-
- String _gaussianToString(
- const SpectralElement& gauss, const CoordinateSystem& csys,
- const Vector<Double> world
- ) const;
-
- String _polynomialToString(
- const SpectralElement& poly, const CoordinateSystem& csys,
- const Vector<Double> imPix, const Vector<Double> world
- ) const;
-
- static void _makeSolutionImage(
- const String& name, const IPosition& shape, const CoordinateSystem&
csys,
- const Vector<Double>& values, const String& unit
- );
-};
-}
-
-#endif
=======================================
--- /trunk/images/Images/test/dImageProfileFit.cc Sun Jul 5 22:04:32 2009
+++ /dev/null
@@ -1,130 +0,0 @@
-//# dImageProfileFit.cc:
-//# Copyright (C) 1996,1997,1999,2000,2001,2002,2003,2004
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This library is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU Library General Public License as published
by
-//# the Free Software Foundation; either version 2 of the License, or (at
your
-//# option) any later version.
-//#
-//# This library is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
-//# License for more details.
-//#
-//# You should have received a copy of the GNU Library General Public
License
-//# along with this library; if not, write to the Free Software Foundation,
-//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id:
-//
-#include <images/Images/ImageProfileFit.h>
-
-#include <casa/aips.h>
-#include <casa/Exceptions/Error.h>
-#include <casa/Inputs/Input.h>
-#include <casa/Logging.h>
-#include <tables/Tables/Table.h>
-#include <casa/BasicSL/String.h>
-#include <casa/Containers/Record.h>
-#include <casa/Utilities/Assert.h>
-
-#include <images/Images/PagedImage.h>
-#include <casa/iostream.h>
-
-
-
-#include <casa/namespace.h>
-int main (int argc, const char* argv[])
-{
-
-try {
-
- Input inputs(1);
-
-// Get inputs
-
- String name = "test_image.im";
- inputs.create("in", name, "Input file name");
- inputs.create("out", "dImageProfileFit_tmp", "Output root image name");
- inputs.create("axis", "2", "axis to fit (default is 2)");
- inputs.readArguments(argc, argv);
-//
- const String in = inputs.getString("in");
- const String out = inputs.getString("out");
- Int axis = inputs.getInt("axis");
-//
- if (in.empty()) {
- cout << "You must give an input image name" << endl;
- return 1;
- }
- if (out.empty()) {
- cout << "You must give an output image root name" << endl;
- return 1;
- }
-//
- DataType imageType = imagePixelType(in);
- if (imageType!=TpFloat) {
- cout << "The image must be of type Float" << endl;
- return 1;
- }
-
-// Construct images
-
- PagedImage<Float> inImage(in);
- if (axis > Int(inImage.ndim())) {
- cout << "Axis must be <=" << inImage.ndim() << endl;
- return 1;
- }
-//
- String outFitName(out+String(".fit"));
- PagedImage<Float>* pOutFit = new
PagedImage<Float>(TiledShape(inImage.shape()),
- inImage.coordinates(), outFitName);
-
- String outResidName(out+String(".residual"));
- PagedImage<Float>* pOutResid = new
PagedImage<Float>(TiledShape(inImage.shape()),
- inImage.coordinates(), outResidName);
-// Make fitter
-
- ImageProfileFit fitter;
-
-// Set Data
-
- fitter.setData (inImage, axis, False);
-
-// Initial estimate
-
- Int nMax = 1;
- Bool ok = fitter.estimate(nMax);
- AlwaysAssert(ok, AipsError);
-// Fit
-
- Record rec;
- Bool xAbs = True;
- String xUnit ("pix");
- String doppler("RADIO");
- // For the time being use fillRecord=False, otherwise an exception is
- // thrown by ImageProfileFit::getElements that it cannot handle
polynomials.
- //// fitter.fit (True, rec, xAbs, xUnit, doppler, pOutFit, pOutResid,
2);
- fitter.fit (False, rec, xAbs, xUnit, doppler, pOutFit, pOutResid, 2);
-//
- delete pOutFit;
- delete pOutResid;
- Table::deleteTable (outFitName, True);
- Table::deleteTable (outResidName, True);
-
-} catch (AipsError x) {
- cerr << "aipserror: error " << x.getMesg() << endl;
- return 1;
-}
-
- return 0;
-}
-
=======================================
--- /trunk/images/Images/test/dImageProfileFit.run Tue Nov 1 07:01:00 2005
+++ /dev/null
@@ -1,6 +0,0 @@
-#!/bin/sh
-
-# dImageProfileFit often fails or hangs for unknown reason.
-# So skip this test until resolved.
- exit 3
-
=======================================
--- /trunk/images/Images/test/tImageCollapser.cc Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,255 +0,0 @@
-//# tImageFitter.cc: test the PagedImage class
-//# Copyright (C) 1994,1995,1998,1999,2000,2001,2002
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or(at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#include <images/Images/ImageCollapser.h>
-
-#include <images/Images/ImageInputProcessor.h>
-#include <images/Images/ImageUtilities.h>
-#include <images/Images/FITSImage.h>
-#include <images/Images/PagedImage.h>
-
-#include <casa/OS/Directory.h>
-#include <casa/namespace.h>
-
-#include <sys/types.h>
-
-#include <unistd.h>
-#include <iomanip>
-
-
-
-uInt testNumber = 0;
-
-void writeTestString(const String& test) {
- cout << "\n" << "*** " << testNumber << ": "
- << test << " ***" << endl;
- testNumber++;
-}
-
-String dirName;
-
-String outname() {
- return dirName + "/" + "test_" + String::toString(testNumber) + ".im";
-
-}
-
-void checkImage(
- const ImageInterface<Float> *gotImage, const String& expectedName
-) {
- FITSImage expectedImage(expectedName);
- AlwaysAssert(gotImage->shape() == expectedImage.shape(), AipsError);
- Array<Float> diffData = gotImage->get() - expectedImage.get();
- AlwaysAssert(max(abs(diffData)) == 0, AipsError);
- CoordinateSystem gotCsys = gotImage->coordinates();
- CoordinateSystem expectedCsys = expectedImage.coordinates();
- Array<Double> diffPixels = gotCsys.referencePixel() -
expectedCsys.referencePixel();
- AlwaysAssert(max(abs(diffPixels)) == 0, AipsError);
- Array<Double> fracDiffRef = (
- gotCsys.referenceValue() - expectedCsys.referenceValue()
- )/expectedCsys.referenceValue();
- AlwaysAssert(max(abs(fracDiffRef)) <= 1.5e-6, AipsError);
-}
-
-void checkImage(
- const String& gotName, const String& expectedName
-) {
- PagedImage<Float> gotImage(gotName);
- checkImage(&gotImage, expectedName);
-}
-
-
-void testException(
- const String& test, const String& aggString,
- const String& imagename, const String& region,
- const String& box, const String& chans,
- const String& stokes, const String& mask,
- const uInt compressionAxis
-) {
- writeTestString(test);
- Bool exceptionThrown = true;
- try {
- ImageCollapser collapser(
- aggString, imagename, region, box,
- chans, stokes, mask, compressionAxis,
- outname(), False
- );
- // should not get here, fail if we do.
- exceptionThrown = false;
- AlwaysAssert(false, AipsError);
- }
- catch (AipsError x) {
- AlwaysAssert(exceptionThrown, AipsError);
- }
-}
-
-int main() {
- pid_t pid = getpid();
- ostringstream os;
- os << "tImageCollapser_tmp_" << pid;
- dirName = os.str();
- Directory workdir(dirName);
- String goodImage("collapse_in.fits");
- const String ALL = RegionManager::ALL;
- workdir.create();
- uInt retVal = 0;
- try {
- testException(
- "Exception if no image name given", "mean",
- "", "", "", "", "", "", 0
- );
- testException(
- "Exception if bogus image name given", "mean",
- "mybogus.im", "", "", "", "", "", 0
- );
- testException(
- "Exception if no aggregate string given", "",
- goodImage, "", "", "", "", "", 0
- );
- testException(
- "Exception if bogus aggregate string given", "bogus function",
- goodImage, "", "", "", "", "", 0
- );
- testException(
- "Exception if bogus region string given", "mean",
- goodImage, "bogus_region", "", "", "", "", 0
- );
- testException(
- "Exception if bogus box string given #1", "mean",
- goodImage, "", "abc", "", "", "", 0
- );
- testException(
- "Exception if bogus box string given #2", "mean",
- goodImage, "", "0,0,1000,1000", "", "", "", 0
- );
- {
- writeTestString("average full image collapse along axis 0");
- ImageCollapser collapser(
- "mean", goodImage, "", "", ALL,
- ALL, "", 0, outname(), False
- );
- collapser.collapse(False);
- checkImage(outname(), "collapse_avg_0.fits");
- }
- {
- writeTestString("average full image collapse along axis 2");
- ImageCollapser collapser(
- "mean", goodImage, "", "", ALL,
- ALL, "", 2, outname(), False
- );
- collapser.collapse(False);
- checkImage(outname(), "collapse_avg_2.fits");
- }
- {
- writeTestString("sum subimage collapse along axis 1");
- ImageCollapser *collapser = new ImageCollapser(
- "sum", goodImage, "", "1,1,2,2", "1~2",
- "qu", "", 2, outname(), False
- );
- collapser->collapse(False);
- delete collapser;
- // and check that we can overwrite the previous output
- collapser = new ImageCollapser(
- "sum", goodImage, "", "1,1,2,2", "1~2",
- "qu", "", 1, outname(), True
- );
- collapser->collapse(False);
- delete collapser;
- checkImage(outname(), "collapse_sum_1.fits");
- }
- {
- writeTestString("Check not specifying out file is ok");
- ImageCollapser collapser(
- "mean", goodImage, "", "", ALL,
- ALL, "", 2, "", False
- );
- ImageInterface<Float> *collapsed = collapser.collapse(True);
- checkImage(collapsed, "collapse_avg_2.fits");
- delete collapsed;
- }
- {
- writeTestString("Check not wanting return pointer results in a NULL
pointer being returned");
- ImageCollapser collapser(
- "mean", goodImage, "", "", ALL,
- ALL, "", 2, "", False
- );
- ImageInterface<Float> *collapsed = collapser.collapse(False);
- AlwaysAssert(collapsed == NULL, AipsError);
- }
- {
- writeTestString("average full image collapse along all axes but 0");
- Vector<uInt> axes(3);
- axes[0] = 1;
- axes[1] = 2;
- axes[2] = 3;
-
- ImageCollapser collapser(
- "max", goodImage, "", "", ALL,
- ALL, "", axes, outname(), False
- );
- collapser.collapse(False);
- checkImage(outname(), "collapse_max_0_a.fits");
- }
- {
- writeTestString("average full temporary image collapse along axis
0");
- ImageInterface<Float> *pIm;
- LogIO log;
- ImageUtilities::openImage(pIm, goodImage, log);
- IPosition shape = pIm->shape();
- CoordinateSystem csys = pIm->coordinates();
- Array<Float> vals = pIm->get();
- delete pIm;
- TempImage<Float> tIm(shape, csys);
- tIm.put(vals);
- ImageCollapser collapser(
- "mean", &tIm, "", "", ALL,
- ALL, "", 0, outname(), False
- );
- collapser.collapse(False);
- checkImage(outname(), "collapse_avg_0.fits");
- }
- {
- writeTestString("full image collapse along axes 0, 1");
- Vector<uInt> axes(2, 0);
- axes[1] = 1;
- ImageCollapser collapser(
- "mean", goodImage, "", "", ALL,
- ALL, "", axes, outname(), False
- );
- collapser.collapse(False);
- checkImage(outname(), "collapse_avg_0_1.fits");
- }
- cout << "ok" << endl;
- }
- catch (AipsError x) {
- cerr << "Exception caught: " << x.getMesg() << endl;
- retVal = 1;
- }
- workdir.removeRecursive();
-
- return retVal;
-}
-
=======================================
--- /trunk/images/Images/test/tImageFitter.cc Tue Oct 18 00:39:05 2011
+++ /dev/null
@@ -1,663 +0,0 @@
-//# tImageFitter.cc: test the PagedImage class
-//# Copyright (C) 1994,1995,1998,1999,2000,2001,2002
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or(at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-
-#include <casa/BasicMath/Math.h>
-#include <components/ComponentModels/ComponentList.h>
-#include <components/ComponentModels/Flux.h>
-#include <images/Images/ImageFitter.h>
-#include <measures/Measures/MDirection.h>
-#include <components/ComponentModels/SpectralModel.h>
-#include <components/ComponentModels/ComponentShape.h>
-#include <components/ComponentModels/GaussianShape.h>
-#include <images/Images/ImageAnalysis.h>
-#include <images/Images/FITSImage.h>
-#include <images/Images/ImageMetaData.h>
-#include <images/Regions/RegionManager.h>
-#include <casa/BasicSL/Constants.h>
-#include <casa/OS/Directory.h>
-#include <casa/namespace.h>
-
-#include <sys/types.h>
-#include <unistd.h>
-
-
-void writeTestString(const String& test) {
- cout << "\n" << "*** " << test << " ***" << endl;
-}
-
-void checkImage(
- const String& gotImage, const String& expectedImage,
- const String& differenceImage
- ) {
- ImageAnalysis ia;
- ia.open(gotImage);
- String expr = "\"" + gotImage + "\" - \"" + expectedImage + "\"";
- ia.imagecalc(differenceImage, expr, True);
- ia.open(differenceImage);
- Record stats;
- Vector<Int> axes(2);
- axes[0] = 0;
- axes[1] = 1;
- Record region;
- Vector<String> plotstats(0);
- ia.statistics(stats, axes, region, "", plotstats, Vector<Float>(0),
Vector<Float>(0));
-
- Array<Double> minArray = stats.asArrayDouble("min");
- Array<Double> maxArray = stats.asArrayDouble("max");
-
- vector<Double> min, max;
- minArray.tovector(min);
- maxArray.tovector(max);
-
- AlwaysAssert(min[0] == 0 && max[0] == 0, AipsError);
-}
-
-int main() {
- pid_t pid = getpid();
- ostringstream os;
- os << "tImageFitter_tmp_" << pid;
- String dirName = os.str();
- Directory workdir(dirName);
- workdir.create(True);
- const Double DEGREES_PER_RADIAN = 180/C::pi;
- Double arcsecsPerRadian = DEGREES_PER_RADIAN*3600;
- String test;
- const ImageInterface<Float> *gaussianModel = new
FITSImage("gaussian_model.fits");
- const ImageInterface<Float> *noisyImage = new
FITSImage("gaussian_model_with_noise.fits");
- const ImageInterface<Float> *convolvedModel = new
FITSImage("gaussian_convolved.fits");
- const ImageInterface<Float> *jykms = new
FITSImage("jyperbeamkmpersec.fits");
- const ImageInterface<Float> *gaussNoPol = new
FITSImage("gauss_no_pol.fits");
- const ImageInterface<Float> *twoGauss = new
FITSImage("two_gaussian_model.fits");
- const Path compTable(dirName + "/myCompList.cl");
- // Directory compTableDir(compTable.absoluteName());
- int returnValue = 0;
-
- try {
- {
- writeTestString(
- "test fitter using all available image pixels with model
with no noise"
- );
- ImageFitter fitter = ImageFitter(gaussianModel, "", "");
- // test to ensure exception is thrown if convergence is
checked for before fit is done
- try {
- fitter.converged();
- // should never get there
- AlwaysAssert(false, AipsError);
- }
- catch (AipsError) {
- // got here, just continue
- }
- ComponentList compList = fitter.fit();
- AlwaysAssert(fitter.converged(), AipsError);
- Vector<Quantity> flux;
- compList.getFlux(flux,0);
- // I stokes flux test
- AlwaysAssert(near(flux(0).getValue(), 60318.5801, 1e-4),
AipsError);
- // Q stokes flux test
- AlwaysAssert(flux(1).getValue() == 0, AipsError);
- MDirection direction = compList.getRefDirection(0);
-
AlwaysAssert(near(direction.getValue().getLong("rad").getValue(),
0.000213318, 1e-5), AipsError);
-
AlwaysAssert(near(direction.getValue().getLat("rad").getValue(),
1.939254e-5, 1e-5), AipsError);
-
- Vector<Double> parameters = compList.getShape(0)->parameters();
-
- Double majorAxis = arcsecsPerRadian*parameters(0);
- AlwaysAssert(near(majorAxis, 23.548201, 1e-7), AipsError);
-
- Double minorAxis = arcsecsPerRadian*parameters(1);
- AlwaysAssert(near(minorAxis, 18.838560, 1e-7), AipsError);
-
- Double positionAngle = DEGREES_PER_RADIAN*parameters(2);
- AlwaysAssert(near(positionAngle, 120.0, 1e-7), AipsError);
- }
- {
- writeTestString(
- "test fitter using all available image pixels with model
with noise added"
- );
- ImageFitter fitter = ImageFitter(noisyImage, 0, "");
- ComponentList compList = fitter.fit();
- AlwaysAssert(fitter.converged(), AipsError);
- Vector<Quantity> flux;
- compList.getFlux(flux,0);
- // I stokes flux test
- cout << flux(0).getValue() << endl;
- AlwaysAssert(near(flux(0).getValue(), 60291.80, 1e-5),
AipsError);
- // Q stokes flux test
- AlwaysAssert(flux(1).getValue() == 0, AipsError);
- MDirection direction = compList.getRefDirection(0);
-
AlwaysAssert(nearAbs(direction.getValue().getLong("rad").getValue(),
0.000213379, 1e-5), AipsError);
-
AlwaysAssert(nearAbs(direction.getValue().getLat("rad").getValue(),
1.9358247e-5, 1e-5), AipsError);
-
- Vector<Double> parameters = compList.getShape(0)->parameters();
-
- Double majorAxis = arcsecsPerRadian*parameters(0);
- AlwaysAssert(near(majorAxis, 23.53002154, 1e-7), AipsError);
-
- Double minorAxis = arcsecsPerRadian*parameters(1);
- AlwaysAssert(near(minorAxis, 18.86212502, 1e-7), AipsError);
-
- Double positionAngle = DEGREES_PER_RADIAN*parameters(2);
- AlwaysAssert(nearAbs(positionAngle, 119.881851057, 1e-7),
AipsError);
- }
- {
- writeTestString(
- "test fitter using a box region with model with noise
added"
- );
- ImageFitter fitter(noisyImage, 0, "130,89,170,129");
- ComponentList compList = fitter.fit();
- AlwaysAssert(fitter.converged(), AipsError);
- Vector<Quantity> flux;
- compList.getFlux(flux,0);
- // I stokes flux test
- AlwaysAssert(near(flux(0).getValue(), 60319.860, 1e-5),
AipsError);
- // Q stokes flux test
- AlwaysAssert(flux(1).getValue() == 0, AipsError);
- MDirection direction = compList.getRefDirection(0);
-
AlwaysAssert(nearAbs(direction.getValue().getLong("rad").getValue(),
0.000213372, 1e-5), AipsError);
-
AlwaysAssert(nearAbs(direction.getValue().getLat("rad").getValue(),
1.9359058e-5, 1e-5), AipsError);
-
- Vector<Double> parameters = compList.getShape(0)->parameters();
-
- Double majorAxis = arcsecsPerRadian*parameters(0);
- AlwaysAssert(near(majorAxis, 23.545212, 1e-7), AipsError);
-
- Double minorAxis = arcsecsPerRadian*parameters(1);
- AlwaysAssert(near(minorAxis, 18.864505, 1e-7), AipsError);
-
- Double positionAngle = DEGREES_PER_RADIAN*parameters(2);
- AlwaysAssert(nearAbs(positionAngle, 119.81297, 1e-5),
AipsError);
- }
- {
- writeTestString(
- "test fitter using a region record with model with noise added"
- );
- IPosition imShape = noisyImage->shape();
- Vector<Double> blc(imShape.nelements(), 0);
- Vector<Double> trc(imShape.nelements(), 0);
- Vector<Double> inc(imShape.nelements(), 1);
-
- Vector<Int> dirNums =
ImageMetaData(*noisyImage).directionAxesNumbers();
- blc[dirNums[0]] = 130;
- blc[dirNums[1]] = 89;
- trc[dirNums[0]] = 170;
- trc[dirNums[1]] = 129;
-
- RegionManager rm;
- Record *box = rm.box(blc, trc, inc, "abs", False);
- ImageFitter fitter(noisyImage, box);
- ComponentList compList = fitter.fit();
- AlwaysAssert(fitter.converged(), AipsError);
-
- Vector<Quantity> flux;
- compList.getFlux(flux,0);
- // I stokes flux test
- AlwaysAssert(near(flux(0).getValue(), 60319.8604, 1e-5),
AipsError);
- // Q stokes flux test
- AlwaysAssert(flux(1).getValue() == 0, AipsError);
- MDirection direction = compList.getRefDirection(0);
-
AlwaysAssert(nearAbs(direction.getValue().getLong("rad").getValue(),
0.000213372, 1e-5), AipsError);
-
AlwaysAssert(nearAbs(direction.getValue().getLat("rad").getValue(),
1.9359058e-5, 1e-5), AipsError);
-
- Vector<Double> parameters = compList.getShape(0)->parameters();
-
- Double majorAxis = arcsecsPerRadian*parameters(0);
- AlwaysAssert(near(majorAxis, 23.545212, 1e-7), AipsError);
-
- Double minorAxis = arcsecsPerRadian*parameters(1);
- AlwaysAssert(near(minorAxis, 18.864505, 1e-7), AipsError);
-
- Double positionAngle = DEGREES_PER_RADIAN*parameters(2);
- AlwaysAssert(near(positionAngle, 119.81297, 1e-5), AipsError);
- }
- {
- // test fitter using an includepix (i=0) and excludepix (i=1)
range with model with noise
- ImageAnalysis ia;
- String outname = dirName + "/myout.im";
- ImageInterface<Float>* outIm = ia.newimagefromfits(outname,
noisyImage->name(), 0, 0, False, True);
- ImageAnalysis ia2(outIm);
- String goodMask = "\"" + outname + "\">40";
- Record r;
- ia2.calcmask(goodMask, r, "mymask", True);
- // it appears this call is explicitly needed even though the
previous statement should have made
- // the new mask the default mask
- outIm->setDefaultMask("mymask");
-
- for (uInt i=0; i<4; i++) {
- String mask;
- Vector<Float> includepix, excludepix;
- Vector<const ImageInterface<Float>* > images(4,
noisyImage);
- images[3] = outIm;
- switch (i) {
- case 0:
- writeTestString("test using includepix range");
- includepix.resize(2);
- includepix(0) = 40;
- includepix(1) = 121;
- mask = "";
- break;
- case 1:
- writeTestString("test using excludepix range");
- includepix.resize(0);
- excludepix.resize(2);
- excludepix(0) = -10;
- excludepix(1) = 40;
- mask = "";
- break;
- case 2:
- includepix.resize(0);
- excludepix.resize(0);
- mask = "\"gaussian_model_with_noise.fits\">40";
- writeTestString("test using LEL mask " + mask);
- break;
- case 3:
- includepix.resize(0);
- excludepix.resize(0);
- mask = "";
- writeTestString("test using pixel mask " +
images[i]->name() + "/mymask");
- break;
- }
- ImageFitter fitter(
- images[i], "", "", 0, "I", mask, includepix, excludepix
- );
- ComponentList compList = fitter.fit();
-
- //cout << "*** def mask " << outIm->getMask() << endl;
- cout << "*** image " << images[i]->name() << endl;
-
- AlwaysAssert(fitter.converged(), AipsError);
- Vector<Quantity> flux;
- compList.getFlux(flux,0);
- // I stokes flux test
- cout << "flux " << flux(0).getValue() << endl;
- AlwaysAssert(near(flux(0).getValue(), 60354.3232, 1e-5),
AipsError);
- // Q stokes flux test
- AlwaysAssert(flux(1).getValue() == 0, AipsError);
- MDirection direction = compList.getRefDirection(0);
-
AlwaysAssert(near(direction.getValue().getLong("rad").getValue(),
0.000213391, 1e-5), AipsError);
-
AlwaysAssert(near(direction.getValue().getLat("rad").getValue(),
1.93449e-05, 1e-5), AipsError);
-
- Vector<Double> parameters =
compList.getShape(0)->parameters();
-
- Double majorAxis = arcsecsPerRadian*parameters(0);
- AlwaysAssert(near(majorAxis, 23.541712, 1e-7), AipsError);
-
- Double minorAxis = arcsecsPerRadian*parameters(1);
- AlwaysAssert(near(minorAxis, 18.882029, 1e-7), AipsError);
-
- Double positionAngle = DEGREES_PER_RADIAN*parameters(2);
- AlwaysAssert(near(positionAngle, 119.769648, 1e-7),
AipsError);
- }
- }
- {
- writeTestString("test writing of residual and mdoel images");
-
- String residImage = dirName + "/residualImage";
- String modelImage = dirName + "/modelImage";
- String residDiff = dirName + "/residualImage.diff";
- String modelDiff = dirName + "/modelImage.diff";
- ImageFitter fitter(
- noisyImage, "", "100,100,200,200", 0, "I", "",
- Vector<Float>(0), Vector<Float>(0), residImage,
- modelImage
- );
- fitter.fit();
- AlwaysAssert(fitter.converged(), AipsError);
- writeTestString("test residual image correctness");
- checkImage(
- residImage, "gaussian_model_with_noise_resid.fits",
- residDiff
- );
- writeTestString("test model image correctness");
- checkImage(
- modelImage, "gaussian_model_with_noise_model.fits",
- modelDiff
- );
-
- writeTestString("test fit succeeds when model and residual
cannot be written");
- residImage = "/residualImage";
- modelImage = "/modelImage";
-
- ImageFitter fitter2(
- noisyImage, "", "100,100,200,200", 0, "I", "",
- Vector<Float>(0), Vector<Float>(0), residImage,
- modelImage
- );
- fitter2.fit();
- AlwaysAssert(fitter2.converged(), AipsError);
- }
- {
- writeTestString("test fitting model gaussian that has been
convolved with a beam");
- ImageFitter fitter(convolvedModel, "", "");
- ComponentList compList = fitter.fit();
- AlwaysAssert(fitter.converged(), AipsError);
- Vector<Quantity> flux;
-
- compList.getFlux(flux,0);
- // I stokes flux test
- AlwaysAssert(near(flux(0).getValue(), 60318.6, 1e-5), AipsError);
- // Q stokes flux test
- AlwaysAssert(flux(1).getValue() == 0, AipsError);
- MDirection direction = compList.getRefDirection(0);
- AlwaysAssert(near(direction.getValue().getLong("rad").getValue(),
0.000213318, 1e-5), AipsError);
- AlwaysAssert(near(direction.getValue().getLat("rad").getValue(),
1.939254e-5, 1e-5), AipsError);
-
- Vector<Double> parameters = compList.getShape(0)->parameters();
-
- Double majorAxis = arcsecsPerRadian*parameters(0);
- AlwaysAssert(near(majorAxis, 26.50461508, 1e-7), AipsError);
-
- Double minorAxis = arcsecsPerRadian*parameters(1);
- AlwaysAssert(near(minorAxis, 23.99821851, 1e-7), AipsError);
-
- Double positionAngle = DEGREES_PER_RADIAN*parameters(2);
- AlwaysAssert(near(positionAngle, 126.3211060, 1e-7), AipsError);
-
- Quantity xPosError =
compList.getShape(0)->refDirectionErrorLong();
- AlwaysAssert(near(xPosError.getValue(), 8.8704e-08, 1e-4),
AipsError);
-
- Quantity yPosError = compList.getShape(0)->refDirectionErrorLat();
- AlwaysAssert(near(yPosError.getValue(), 1.05117e-07, 1e-4),
AipsError);
-
- Vector<Double> compErrors = compList.getShape(0)->errors();
- AlwaysAssert(near(compErrors[0], 8.57789e-09, 1e-4), AipsError);
- AlwaysAssert(near(compErrors[1], 7.11494e-09, 1e-4), AipsError);
- AlwaysAssert(near(compErrors[2], 3.40562e-05, 1e-4), AipsError);
-
- GaussianShape *gauss = dynamic_cast<GaussianShape
*>(compList.getShape(0)->clone());
- AlwaysAssert(gauss->majorAxisError().getUnit() ==
gauss->majorAxis().getUnit(), AipsError);
- AlwaysAssert(gauss->minorAxisError().getUnit() ==
gauss->minorAxis().getUnit(), AipsError);
- AlwaysAssert(gauss->positionAngleError().getUnit() ==
gauss->positionAngle().getUnit(), AipsError);
- AlwaysAssert(near(gauss->majorAxisError().getValue(), 0.00176932,
1e-5), AipsError);
- AlwaysAssert(near(gauss->minorAxisError().getValue(), 0.00146756,
1e-5), AipsError);
- AlwaysAssert(near(gauss->positionAngleError().getValue(),
0.00195128, 1e-5), AipsError);
-
-
- }
- {
- writeTestString(
- String("test fitting model gaussian that has been convolved with
a beam and fix ")
- + String("the peak intensity to be artificially low")
- );
- ImageFitter fitter(
- convolvedModel, "", "", 0, "I", "",
- Vector<Float>(0), Vector<Float>(0), "",
- "", "estimates_convolved.txt"
- );
- ComponentList compList = fitter.fit();
- AlwaysAssert(fitter.converged(), AipsError);
- Vector<Quantity> flux;
-
- compList.getFlux(flux,0);
- // I stokes flux test
- AlwaysAssert(near(flux(0).getValue(), 60082.6, 1e-5), AipsError);
- // Q stokes flux test
- AlwaysAssert(flux(1).getValue() == 0, AipsError);
- MDirection direction = compList.getRefDirection(0);
-
AlwaysAssert(nearAbs(direction.getValue().getLong("rad").getValue(),
0.000213318, 1e-5), AipsError);
-
AlwaysAssert(nearAbs(direction.getValue().getLat("rad").getValue(),
1.939254e-5, 1e-5), AipsError);
-
- Vector<Double> parameters = compList.getShape(0)->parameters();
-
- Double majorAxis = arcsecsPerRadian*parameters(0);
- AlwaysAssert(near(majorAxis, 28.21859344, 1e-7), AipsError);
-
- Double minorAxis = arcsecsPerRadian*parameters(1);
-
- AlwaysAssert(near(minorAxis, 25.55011520, 1e-7), AipsError);
-
- Double positionAngle = DEGREES_PER_RADIAN*parameters(2);
- AlwaysAssert(nearAbs(positionAngle, 126.3211050, 1e-7),
AipsError);
- }
- {
- writeTestString("Fit two gaussians");
- ImageFitter fitter(
- twoGauss, "", "", 0, "I", "",
- Vector<Float>(0), Vector<Float>(0), "",
- "", "estimates_2gauss.txt"
- );
- ComponentList compList = fitter.fit();
- AlwaysAssert(fitter.converged(), AipsError);
- Vector<Quantity> flux;
- MDirection direction;
- Vector<Double> parameters;
- AlwaysAssert(compList.nelements() == 2, AipsError);
- Vector<Double> expectedFlux(2);
- expectedFlux[0] = 60318.5820312;
- expectedFlux[1] = 112174.6953125;
- Vector<Double> expectedLong(2);
- expectedLong[0] = 2.1331802e-04;
- expectedLong[1] = -2.2301344e-04;
- Vector<Double> expectedLat(2);
- expectedLat[0] = 1.9392547e-05;
- expectedLat[1] = 4.5572321e-04;
- Vector<Double> expectedMajorAxis(2);
- expectedMajorAxis[0] = 23.548201;
- expectedMajorAxis[1] = 46.582182;
- Vector<Double> expectedMinorAxis(2);
- expectedMinorAxis[0] = 18.838561;
- expectedMinorAxis[1] = 23.613296;
- Vector<Double> expectedPositionAngle(2);
- expectedPositionAngle[0] = 120.0;
- expectedPositionAngle[1] = 140.07385;
-
- for (uInt i = 0; i < compList.nelements(); i++) {
- compList.getFlux(flux,i);
- // I stokes flux test
- AlwaysAssert(near(flux(0).getValue(), expectedFlux[i], 1e-7),
AipsError);
- // Q stokes flux test
- AlwaysAssert(flux(1).getValue() == 0, AipsError);
- direction = compList.getRefDirection(i);
-
-
AlwaysAssert(nearAbs(direction.getValue().getLong("rad").getValue(),
expectedLong[i], 1e-7), AipsError);
-
AlwaysAssert(nearAbs(direction.getValue().getLat("rad").getValue(),
expectedLat[i], 1e-7), AipsError);
- Vector<Double> parameters =
compList.getShape(i)->parameters();
- Double majorAxis = arcsecsPerRadian*parameters(0);
- AlwaysAssert(near(majorAxis, expectedMajorAxis[i], 1e-7),
AipsError);
-
- Double minorAxis = arcsecsPerRadian*parameters(1);
- AlwaysAssert(near(minorAxis, expectedMinorAxis[i], 1e-7),
AipsError);
-
- Double positionAngle = DEGREES_PER_RADIAN*parameters(2);
- AlwaysAssert(nearAbs(positionAngle,
expectedPositionAngle[i], 5e-6), AipsError);
- }
- }
- {
- writeTestString("Test of nonconvergence");
- ImageFitter fitter(noisyImage, "", "0,0,20,20");
- fitter.fit();
- AlwaysAssert(! fitter.converged(), AipsError);
- }
- {
- writeTestString("Test of fitting in a multi-polarization image");
- Vector<String> stokes(4);
- stokes[0] = "I";
- stokes[1] = "Q";
- stokes[2] = "U";
- stokes[3] = "V";
-
- Vector<Double> expectedFlux(stokes.size());
- expectedFlux[0] = 133.60641;
- expectedFlux[1] = 400.81921;
- expectedFlux[2] = 375.76801;
- expectedFlux[3] = -1157.92212;
-
- Vector<Double> expectedRA(stokes.size());
- expectedRA[0] = 1.2479113396;
- expectedRA[1] = 1.2479113694;
- expectedRA[2] = 1.2478908580;
- expectedRA[3] = 1.2478908284;
-
- Vector<Double> expectedDec(stokes.size());
- expectedDec[0] = 0.782579122;
- expectedDec[1] = 0.782593666;
- expectedDec[2] = 0.782593687;
- expectedDec[3] = 0.782579143;
-
- Vector<Double> expectedMajorAxis(stokes.size());
- expectedMajorAxis[0] = 7.992524398;
- expectedMajorAxis[1] = 11.988806751;
- expectedMajorAxis[2] = 8.991589959;
- expectedMajorAxis[3] = 12.987878913;
-
- Vector<Double> expectedMinorAxis(stokes.size());
- expectedMinorAxis[0] = 5.994405977;
- expectedMinorAxis[1] = 5.994395540;
- expectedMinorAxis[2] = 4.995338093;
- expectedMinorAxis[3] = 7.992524265;
-
- Vector<Double> expectedPositionAngle(stokes.size());
- expectedPositionAngle[0] = 40.083248;
- expectedPositionAngle[1] = 160.083213;
- expectedPositionAngle[2] = 50.082442;
- expectedPositionAngle[3] = 135.08243;
-
- for (uInt i=0; i<stokes.size(); i++) {
- ImageFitter fitter("imfit_stokes.fits", "", "", 0, stokes[i]);
- ComponentList compList = fitter.fit();
- AlwaysAssert(fitter.converged(), AipsError);
- Vector<Quantity> flux;
- MDirection direction = compList.getRefDirection(0);
- compList.getFlux(flux,0);
- AlwaysAssert(compList.nelements() == 1, AipsError);
- AlwaysAssert(near(flux(i).getValue(), expectedFlux[i], 1e-5),
AipsError);
-
AlwaysAssert(nearAbs(direction.getValue().getLong("rad").getValue(),
expectedRA[i], 1e-8), AipsError);
-
AlwaysAssert(nearAbs(direction.getValue().getLat("rad").getValue(),
expectedDec[i], 1e-8), AipsError);
- Vector<Double> parameters =
compList.getShape(0)->parameters();
- Double majorAxis = arcsecsPerRadian*parameters(0);
- AlwaysAssert(near(majorAxis, expectedMajorAxis[i], 1e-7),
AipsError);
- Double minorAxis = arcsecsPerRadian*parameters(1);
- AlwaysAssert(near(minorAxis, expectedMinorAxis[i], 1e-7),
AipsError);
- Double positionAngle = DEGREES_PER_RADIAN*parameters(2);
- AlwaysAssert(nearAbs(positionAngle,
expectedPositionAngle[i], 5e-6), AipsError);
- }
- }
- {
- writeTestString("Test of CAS-2318 fix");
-
- ImageFitter fitter(
- gaussNoPol, "", "", 0, "", "",
- Vector<Float>(0), Vector<Float>(0), "",
- "", ""
- );
- ComponentList compList = fitter.fit();
- // Just the fact that an exception isn't thrown verifies the
fix
- Vector<Quantity> flux;
- compList.getFlux(flux,0);
- AlwaysAssert(near(flux(0).getValue(), 394312.65593496, 1e-5),
AipsError);
- }
- {
- writeTestString("test fitting image with units of Jy km/s
(CAS-1233");
- ImageFitter fitter(jykms, "", "");
- ComponentList compList = fitter.fit();
- AlwaysAssert(fitter.converged(), AipsError);
- Vector<Quantity> flux;
-
- compList.getFlux(flux,0);
- // I stokes flux test
- AlwaysAssert(near(flux(0).getValue(), 60318.6, 1e-5), AipsError);
- AlwaysAssert(flux(0).getUnit() == "Jy.km/s", AipsError);
- // Q stokes flux test
- AlwaysAssert(flux(1).getValue() == 0, AipsError);
-
- Vector<std::complex<double> > fluxErrors =
compList.component(0).flux().errors();
- AlwaysAssert(near(fluxErrors[0].real(), 0.000863785, 1e-5),
AipsError);
- AlwaysAssert(fluxErrors[0].imag() == 0, AipsError);
- AlwaysAssert(fluxErrors[1].real() == 0, AipsError);
- AlwaysAssert(fluxErrors[1].imag() == 0, AipsError);
-
- MDirection direction = compList.getRefDirection(0);
- AlwaysAssert(near(direction.getValue().getLong("rad").getValue(),
0.000213318, 1e-5), AipsError);
- AlwaysAssert(near(direction.getValue().getLat("rad").getValue(),
1.939254e-5, 1e-5), AipsError);
-
- Vector<Double> parameters = compList.getShape(0)->parameters();
-
- Double majorAxis = arcsecsPerRadian*parameters(0);
- AlwaysAssert(near(majorAxis, 26.50461508, 1e-7), AipsError);
-
- Double minorAxis = arcsecsPerRadian*parameters(1);
- AlwaysAssert(near(minorAxis, 23.99821851, 1e-7), AipsError);
-
- Double positionAngle = DEGREES_PER_RADIAN*parameters(2);
- AlwaysAssert(near(positionAngle, 126.3211060, 1e-7), AipsError);
-
- // CAS-2633
- Double refFreq =
compList.component(0).spectrum().refFrequency().getValue();
- AlwaysAssert(near(refFreq, 1.415e9, 1e-7), AipsError);
- }
- {
- writeTestString("test writing component list (CAS-2595");
- {
- ImageFitter fitter(
- noisyImage, "", "", 0, "I", "", Vector<Float>(0),
Vector<Float>(0),
- "", "", "", "", True, "", compTable.absoluteName(),
ImageFitter::WRITE_NO_REPLACE
- );
- fitter.fit();
- ComponentList c1(compTable);
- AlwaysAssert(c1.nelements() == 1, AipsError);
- }
- {
- ImageFitter fitter(
- twoGauss, "", "", 0, "I", "", Vector<Float>(0),
Vector<Float>(0),
- "", "", "estimates_2gauss.txt", "", True, "",
compTable.absoluteName(),
- ImageFitter::WRITE_NO_REPLACE
- );
- fitter.fit();
- ComponentList c1(compTable);
- AlwaysAssert(c1.nelements() == 1, AipsError);
- }
- {
- ImageFitter fitter(
- twoGauss, "", "", 0, "I", "", Vector<Float>(0),
Vector<Float>(0),
- "", "", "estimates_2gauss.txt", "", True, "",
compTable.absoluteName(),
- ImageFitter::OVERWRITE
- );
- fitter.fit();
- ComponentList c1(compTable);
- AlwaysAssert(c1.nelements() == 2, AipsError);
- }
- }
- cout << "ok" << endl;
- }
- catch (AipsError x) {
- cerr << "Exception caught: " << x.getMesg() << endl;
- returnValue = 1;
- }
- /*
- if(workdir.exists()) {
- workdir.removeRecursive();
- }
- */
- delete gaussianModel;
- delete noisyImage;
- delete convolvedModel;
- delete jykms;
- delete gaussNoPol;
-
- return returnValue;
-}
-
=======================================
--- /trunk/images/Images/test/tImageInputProcessor.cc Tue Oct 18 00:39:05
2011
+++ /dev/null
@@ -1,567 +0,0 @@
-//# tImageFitter.cc: test the PagedImage clas
-//# Copyright (C) 1994,1995,1998,1999,2000,2001,2002
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or(at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#include <images/Images/ImageInputProcessor.h>
-#include <images/Images/FITSImage.h>
-#include <casa/namespace.h>
-#include <casa/iomanip.h>
-
-
-void writeTestString(const String& test) {
- cout << "\n" << "*** " << test << " ***" << endl;
-}
-
-void _checkCorner(const Record& gotRecord, const Vector<Double>& expected)
{
- for (uInt i=0; i<expected.size(); i++) {
- uInt fieldNumber = i+1;
- Double got = gotRecord.asRecord(RecordFieldId(
- "*" + String::toString(fieldNumber))
- ).asDouble(RecordFieldId("value"));
- /*
- cout << setprecision(10);
- cout << "got " << got << endl;
- cout << "expected " << expected << endl;
- */
- AlwaysAssert(fabs((got-expected[i])/expected[i]) < 3e-9, AipsError);
- }
-}
-
-void testException(
- const String& test,
- Vector<ImageInputProcessor::OutputStruct> *outputStruct,
- const String& imagename, const Record *regionPtr,
- const String& regionName, const String& box,
- const String& chans, String& stokes,
- const RegionManager::StokesControl& stokesControl,
- const Bool allowMultipleBoxes
-) {
- ImageInputProcessor processor;
- ImageInterface<Float> *image = 0;
- Record region;
- String diagnostics;
- Bool fail = True;
- try {
- writeTestString(test);
- processor.process(
- image, region, diagnostics, outputStruct,
- stokes, imagename, regionPtr, regionName,
- box, chans, stokesControl,
- allowMultipleBoxes, 0
- );
- // should not get here
- fail = False;
- delete image;
- AlwaysAssert(false, AipsError);
- }
- catch (AipsError) {
- // should get here with fail = true
- AlwaysAssert(fail, AipsError);
- }
-}
-
-void testSuccess(
- const String& test,
- Vector<ImageInputProcessor::OutputStruct> *outputStruct,
- const String& imagename, const Record *regionPtr,
- const String& regionName, const String& box,
- const String& chans, String& stokes,
- const RegionManager::StokesControl& stokesControl,
- const Bool allowMultipleBoxes,
- const Vector<Double>& expectedBlc,
- const Vector<Double>& expectedTrc
-) {
- ImageInputProcessor processor;
- ImageInterface<Float> *image = 0;
- Record region;
- String diagnostics;
- writeTestString(test);
- processor.process(
- image, region, diagnostics, outputStruct, stokes,
- imagename, regionPtr, regionName, box, chans,
- stokesControl, allowMultipleBoxes, 0
- );
- _checkCorner(region.asRecord(RecordFieldId("blc")), expectedBlc);
- _checkCorner(region.asRecord(RecordFieldId("trc")), expectedTrc);
- AlwaysAssert(processor.nSelectedChannels() == 1, AipsError);
-}
-
-void runProcess(
- const String& test,
- Vector<ImageInputProcessor::OutputStruct> *outputStruct,
- const String& imagename, const Record *regionPtr,
- const String& regionName, const String& box,
- const String& chans, String& stokes,
- const RegionManager::StokesControl& stokesControl,
- const Bool allowMultipleBoxes
-) {
- ImageInputProcessor processor;
- ImageInterface<Float> *image = 0;
- Record region;
- String diagnostics;
- writeTestString(test);
- processor.process(
- image, region, diagnostics, outputStruct, stokes,
- imagename, regionPtr, regionName, box, chans,
- stokesControl, allowMultipleBoxes, 0
- );
-}
-
-
-int main() {
- try {
- String goodImage = "image_input_processor.im";
- ImageInputProcessor processor;
- ImageInterface<Float> *image = 0;
- Record region;
- String diagnostics;
- Bool fail = True;
- String none = "";
- testException("Bad image name throws exception",
- 0, "bogus_image", 0, "", "", "", none,
- RegionManager::USE_ALL_STOKES, True
- );
- testException("Bad region name throws exception", 0, goodImage,
0, "bogus.rgn",
- "", "", none, RegionManager::USE_ALL_STOKES, True
- );
- testException("Bad region name in another image throws exception",
- 0, goodImage, 0, "bogus.im:bogus.rgn",
- "", "", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Bad box spec #1 throws exception",
- 0, goodImage, 0, "", "-1,0,10,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Bad box spec #2 throws exception",
- 0, goodImage, 0, "", "0,-1,10,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- testException("Bad box spec #3 throws exception",
- 0, goodImage, 0, "", "0,0,100 ,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Bad box spec #4 throws exception",
- 0, goodImage, 0, "", "0, 0,10 ,100",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Bad box spec #5 throws exception",
- 0, goodImage, 0, "", "5, 0, 0,10 ,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Bad box spec #6 throws exception",
- 0, goodImage, 0, "", "a, 0,10 ,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- testException("Bad box spec #7 throws exception",
- 0, goodImage, 0, "", "1a, 0,10 ,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Valid box spec with invalid channel spec #1 throws exception",
- 0, goodImage, 0, "", "0, 0,10 ,10",
- "1", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Valid box spec with invalid channel spec #2 throws exception",
- 0, goodImage, 0, "", "0, 0,10 ,10",
- "a", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Valid box spec with invalid channel spec #3 throws exception",
- 0, goodImage, 0, "", "0, 0,10 ,10",
- "a-b", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Valid box spec with invalid channel spec #4 throws exception",
- 0, goodImage, 0, "", "0, 0,10 ,10",
- "0-b", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Valid box spec with invalid channel spec #5 throws exception",
- 0, goodImage, 0, "", "0, 0,10 ,10",
- ">0", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Valid box spec with invalid channel spec #6 throws exception",
- 0, goodImage, 0, "", "0, 0,10 ,10",
- "-1", none, RegionManager::USE_ALL_STOKES, True
- );
- testException(
- "Valid box spec with invalid channel spec #7 throws exception",
- 0, goodImage, 0, "", "0, 0,10 ,10",
- "<5", none, RegionManager::USE_ALL_STOKES, True
- );
- String stokes = "b";
- testException(
- "Valid box spec with invalid stokes spec #1 throws exception",
- 0, goodImage, 0, "", "0, 0,10 ,10",
- "", stokes, RegionManager::USE_ALL_STOKES, True
- );
- stokes = "yy";
- testException(
- "Valid box spec with invalid stokes spec #2 throws exception",
- 0, goodImage, 0, "", "0, 0,10 ,10",
- "", stokes, RegionManager::USE_ALL_STOKES, True
- );
- try {
- writeTestString("Calling nSelectedChannels() before process() throws
an exception");
- ImageInputProcessor processor;
- processor.nSelectedChannels();
- // should not get here
- fail = False;
- AlwaysAssert(False, AipsError);
- }
- catch (AipsError) {
- // should get here with fail = true
- AlwaysAssert(fail, AipsError);
- }
- {
- ImageInputProcessor::OutputStruct output;
- String out = "/cannot_create";
- output.label = "file";
- output.outputFile = &out;
- output.required = True;
- output.replaceable = True;
- Vector<ImageInputProcessor::OutputStruct> outs(1);
- outs[0] = output;
- testException(
- "Non-createable output file throws exception",
- &outs, goodImage, 0, "", "0, 0, 10,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- }
- {
- writeTestString("Non-overwriteable output file throws exception");
- ImageInputProcessor::OutputStruct output;
- String out = "/usr";
- output.label = "file";
- output.outputFile = &out;
- output.required = True;
- output.replaceable = True;
- Vector<ImageInputProcessor::OutputStruct> outs(1);
- outs[0] = output;
- testException(
- "Non-overwriteable output file throws exception",
- &outs, goodImage, 0, "", "0, 0, 10,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- }
-
- {
- writeTestString("Non-replaceable output file throws exception");
- ImageInputProcessor::OutputStruct output;
- String out = "input_image_processor_dont_replace_me";
- output.label = "file";
- output.outputFile = &out;
- output.required = True;
- output.replaceable = False;
- Vector<ImageInputProcessor::OutputStruct> outs(1);
- outs[0] = output;
- testException(
- "Non-replaceable output file throws exception",
- &outs, goodImage, 0, "", "0, 0, 10,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- }
- testException(
- "Multiple boxes with allowMultipleRegions = False throws
exception",
- 0, goodImage, 0, "", "0, 0,10 ,10, 20,20,30,30",
- "", none, RegionManager::USE_ALL_STOKES, False
- );
- stokes = "iu";
- testException(
- "Multiple stokes ranges with allowMultipleRegions = False throws
exception",
- 0, goodImage, 0, "", "0, 0,10 ,10",
- "", stokes, RegionManager::USE_ALL_STOKES, False
- );
- stokes = "";
- testException(
- "Bad region name throws exception",
- 0, goodImage, 0, "mybox", "",
- "", stokes, RegionManager::USE_ALL_STOKES, False
- );
- {
- stokes = "iu";
- writeTestString("Multiple stokes ranges with allowMultipleRegions
= True succeeds");
- processor.process(
- image, region, diagnostics, 0, stokes, goodImage, 0, "", "0,
0,10 ,10",
- none, RegionManager::USE_ALL_STOKES, True, 0
- );
- delete image;
- // FIXME just checks that no excpetion is thrown at this point,
need to
- // do region checking.
- }
- Vector<Double> expectedBlc(4);
- expectedBlc[0] = 1.24795230e+00;
- expectedBlc[1] = 7.82549990e-01;
- expectedBlc[2] = 4.73510000e+09;
- expectedBlc[3] = 1;
- Vector<Double> expectedTrc(4);
- expectedTrc[0] = 1.24784989e+00;
- expectedTrc[1] = 7.82622817e-01;
- expectedTrc[2] = 4.73510000e+09;
- expectedTrc[3] = 4;
-
- testSuccess(
- "Nothing specified gives entire image as region",
- 0, goodImage, 0, "", "", "", none,
- RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- expectedTrc[0] = 1.24793182e+00;
- expectedTrc[1] = 7.82564556e-01;
- expectedTrc[2] = 4.73510000e+09;
- expectedTrc[3] = 4;
- testSuccess(
- "Valid box specification succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10",
- "", none, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- testSuccess(
- "Valid box specification with valid channel specification #1
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10",
- "0~0", none, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- testSuccess(
- "Valid box specification with valid channel specification #2
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10",
- "0", none, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- testSuccess(
- "Valid box specification with valid channel specification #3
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10",
- "0,0,0", none, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- testSuccess(
- "Valid box specification with valid channel specification #4
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10",
- "0;0;0", none, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- testSuccess(
- "Valid box specification with valid channel specification #5
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10",
- "0,0;0", none, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- stokes = "QVIU";
- testSuccess(
- "Valid box specification with valid stokes specification #1
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10",
- "", stokes, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- {
- expectedTrc[3] = 3;
- stokes = "QIU";
- testSuccess(
- "Valid box specification with valid stokes specification #2
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10", "",
- stokes, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- }
- {
- stokes = "Q";
- expectedBlc[3] = expectedTrc[3] = 2;
- testSuccess(
- "Valid box specification with valid stokes specification #3
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10",
- "", stokes, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- }
- {
- expectedBlc[3] = 1;
- expectedTrc[3] = 3;
- stokes = "Q,I,U";
- testSuccess(
- "Valid box specification with valid stokes specification #4
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10", "",
- stokes, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- }
- {
- expectedBlc[3] = 1;
- expectedTrc[3] = 3;
- stokes = "Q;I;U";
- testSuccess(
- "Valid box specification with valid stokes specification #5
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10", "",
- stokes, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- }
- {
- expectedBlc[3] = 1;
- expectedTrc[3] = 3;
- stokes = "Q,I;U";
- testSuccess(
- "Valid box specification with valid stokes specification #6
succeeds",
- 0, goodImage, 0, "", "0, 0, 10,10", "",
- stokes, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- }
- {
- expectedBlc[3] = 1;
- expectedTrc[3] = 4;
- testSuccess(
- "Valid box specification using all polarizations for blank
stokes",
- 0, goodImage, 0, "", "0, 0, 10,10",
- "", none, RegionManager::USE_ALL_STOKES, True,
- expectedBlc, expectedTrc
- );
- }
- {
- expectedTrc[3] = 1;
- expectedBlc[3] = 1;
- stokes = " ";
- testSuccess(
- "Valid box specification using first polarizations for blank
stokes",
- 0, goodImage, 0, "", "0, 0, 10,10",
- "", stokes, RegionManager::USE_FIRST_STOKES, True,
- expectedBlc, expectedTrc
- );
- }
- {
- // these are one relative
- expectedBlc.set(1.0);
- expectedTrc[0] = 21;
- expectedTrc[1] = 21;
- expectedTrc[2] = 1;
- expectedTrc[3] = 3;
- stokes = "";
- String regionFile = goodImage + "/mybox.rgn";
- testSuccess(
- "Valid region file",
- 0, goodImage, 0, regionFile, "",
- "", stokes, RegionManager::USE_FIRST_STOKES, True,
- expectedBlc, expectedTrc
- );
- }
- {
- // these are one relative
- expectedBlc.set(1.0);
- expectedTrc[0] = 21;
- expectedTrc[1] = 21;
- expectedTrc[2] = 1;
- expectedTrc[3] = 3;
- stokes = "";
- String region = goodImage + ":" + "mybox2";
- testSuccess(
- "Valid region description from image table",
- 0, goodImage, 0, region, "",
- "", stokes, RegionManager::USE_FIRST_STOKES, True,
- expectedBlc, expectedTrc
- );
- }
- {
- writeTestString("Non-required, non-overwriteable output file is
set to blank");
- ImageInputProcessor::OutputStruct output;
- String out = "/usr";
- output.label = "file";
- output.outputFile = &out;
- output.required = False;
- output.replaceable = True;
- Vector<ImageInputProcessor::OutputStruct> outs(1);
- outs[0] = output;
- runProcess(
- "Non-required, non-overwriteable output file is set to blank",
- &outs, goodImage, 0, "", "0, 0, 10,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- AlwaysAssert(out.empty(), AipsError);
- }
- {
- ImageInputProcessor::OutputStruct output;
- String out = "/cannot_write_me";
- output.label = "file";
- output.outputFile = &out;
- output.required = False;
- output.replaceable = True;
- Vector<ImageInputProcessor::OutputStruct> outs(1);
- outs[0] = output;
- runProcess(
- "Non-required, non-createable output file is set to blank",
- &outs, goodImage, 0, "", "0, 0, 10,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- AlwaysAssert(out.empty(), AipsError);
- }
- {
- ImageInputProcessor::OutputStruct output;
- String out = "input_image_processor_dont_replace_me";
- output.required = False;
- output.label = "file";
- output.outputFile = &out;
- output.required = False;
- output.replaceable = False;
- Vector<ImageInputProcessor::OutputStruct> outs(1);
- outs[0] = output;
- runProcess(
- "Non-required, non-replaceable output file is set to blank",
- &outs, goodImage, 0, "", "0, 0, 10,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- AlwaysAssert(out.empty(), AipsError);
- }
- {
- ImageInputProcessor::OutputStruct output;
- output.label = "file";
- String name = "youcanwritemedddslsl";
- String save = name;
- output.outputFile = &name;
- output.required = True;
- output.replaceable = False;
- Vector<ImageInputProcessor::OutputStruct> outs(1);
- outs[0] = output;
- runProcess(
- "Writeable file is not reset",
- &outs, goodImage, 0, "", "0, 0, 10,10",
- "", none, RegionManager::USE_ALL_STOKES, True
- );
- AlwaysAssert(name == save, AipsError);
- }
- }
- catch (AipsError x) {
- cerr << "Exception caught: " << x.getMesg() << endl;
- return 1;
- }
- cout << "ok" << endl;
- return 0;
-}
-
=======================================
--- /trunk/images/Images/test/tImageProfileFitter.cc Tue Oct 18 00:39:05
2011
+++ /dev/null
@@ -1,331 +0,0 @@
-//# tImageFitter.cc: test the PagedImage class
-//# Copyright (C) 1994,1995,1998,1999,2000,2001,2002
-//# Associated Universities, Inc. Washington DC, USA.
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the
Free
-//# Software Foundation; either version 2 of the License, or(at your
option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
WITHOUT
-//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for
-//# more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning AIPS++ should be addressed as follows:
-//# Internet email: aips2-...@nrao.edu.
-//# Postal address: AIPS++ Project Office
-//# National Radio Astronomy Observatory
-//# 520 Edgemont Road
-//# Charlottesville, VA 22903-2475 USA
-//#
-//# $Id$
-
-#include <images/Images/ImageProfileFitter.h>
-
-#include <images/Images/FITSImage.h>
-#include <images/Images/ImageUtilities.h>
-
-#include <casa/OS/Directory.h>
-#include <casa/namespace.h>
-
-#include <sys/types.h>
-
-#include <unistd.h>
-#include <iomanip>
-
-
-
-uInt testNumber = 0;
-
-void writeTestString(const String& test) {
- cout << "\n" << "*** " << testNumber << ": "
- << test << " ***" << endl;
- testNumber++;
-}
-
-String dirName;
-
-
-void checkImage(
- const ImageInterface<Float> *gotImage, const String& expectedName
-) {
- Float ftol = 1.0e-12; // 3.0e-13 is too small.
-
- FITSImage expectedImage(expectedName);
- AlwaysAssert(gotImage->shape() == expectedImage.shape(), AipsError);
- Array<Float> diffData = gotImage->get() - expectedImage.get();
- Float maxdiff = max(abs(diffData));
- if(maxdiff > ftol){
- cerr << "For expectedImage " << expectedName << ":" << endl;
- cerr << "\tmaxdiff = " << maxdiff << endl;
- cerr << "\t ftol = " << ftol << endl;
- }
- AlwaysAssert(max(abs(diffData)) <= ftol, AipsError);
- CoordinateSystem gotCsys = gotImage->coordinates();
- CoordinateSystem expectedCsys = expectedImage.coordinates();
- Array<Double> diffPixels = gotCsys.referencePixel() -
expectedCsys.referencePixel();
- AlwaysAssert(max(abs(diffPixels)) == 0, AipsError);
- Array<Double> fracDiffRef = (
- gotCsys.referenceValue() - expectedCsys.referenceValue()
- )/expectedCsys.referenceValue();
- AlwaysAssert(max(abs(fracDiffRef)) <= ftol, AipsError);
-}
-
-void checkImage(
- const String& gotName, const String& expectedName
-) {
- PagedImage<Float> gotImage(gotName);
- checkImage(&gotImage, expectedName);
-}
-
-
-void testException(
- const String& test,
- const String& imagename, const String& region,
- const String& box, const String& chans,
- const String& stokes, const String& mask,
- const uInt axis, const Bool multiFit,
- const String& residual, const String& model,
- const uInt ngauss, const Int polyOrder
-) {
- writeTestString(test);
- Bool exceptionThrown = true;
- try {
- ImageProfileFitter fitter(
- imagename, region, box,
- chans, stokes,
- mask, axis, multiFit,
- residual, model, ngauss,
- polyOrder
- );
-
- // should not get here, fail if we do.
- exceptionThrown = false;
- AlwaysAssert(false, AipsError);
- }
- catch (AipsError x) {
- AlwaysAssert(exceptionThrown, AipsError);
- }
-}
-
-int main() {
- pid_t pid = getpid();
- ostringstream os;
- os << "tImageProfileFitter_tmp_" << pid;
- dirName = os.str();
- Directory workdir(dirName);
- String goodImage = "specfit_multipix_2gauss.fits";
- String goodPolyImage = "specfit_multipix_poly_2gauss.fits";
- workdir.create();
- uInt retVal = 0;
- try {
- testException(
- "Exception if no image name given", "", "",
- "", "", "", "", 2, False, "", "", 1, -1
- );
- testException(
- "Exception if bogus image name given", "my bad", "",
- "", "", "", "", 2, False, "", "", 1, -1
- );
- testException(
- "Exception if given axis is out of range", goodImage, "",
- "", "", "", "", 5, False, "", "", 1, -1
- );
- testException(
- "Exception if given axis is out of range", goodImage, "bogus",
- "", "", "", "", 2, False, "", "", 1, -1
- );
- testException(
- "Exception if bogus box string given #1", goodImage, "",
- "abc", "", "", "", 2, False, "", "", 1, -1
- );
- testException(
- "Exception if bogus box string given #2", goodImage, "",
- "0,0,1000,1000", "", "", "", 2, False, "", "", 1, -1
- );
- testException(
- "Exception if bogus chans string given #1", goodImage, "",
- "", "abc", "", "", 2, False, "", "", 1, -1
- );
- testException(
- "Exception if bogus chans string given #2", goodImage, "",
- "", "0-200", "", "", 2, False, "", "", 1, -1
- );
- testException(
- "Exception if bogus stokes string given #1", goodImage, "",
- "", "", "abc", "", 2, False, "", "", 1, -1
- );
- testException(
- "Exception if bogus stokes string given #2", goodImage, "",
- "", "", "v", "", 2, False, "", "", 1, -1
- );
- testException(
- "Exception if no gaussians and no polynomial specified",
goodImage, "",
- "", "", "", "", 2, False, "", "", 0, -1
- );
- {
- writeTestString("test results of non-multi-pixel two gaussian fit");
- ImageProfileFitter *fitter;
- ImageInterface<Float> *image;
- LogIO log;
- ImageUtilities::openImage(image, goodImage, log);
- for (uInt i=0; i<1; i++) {
- if (i == 0) {
- fitter = new ImageProfileFitter(
- goodImage, "", "", "", "", "", 2, False,
- "", "", 2, -1
- );
- }
- else {
- new ImageProfileFitter(
- image, "", "", "", "", "", 2, False,
- "", "", 2, -1
- );
- }
-
- }
- Record results = fitter->fit();
- delete fitter;
- Vector<Bool> converged = results.asArrayBool("converged");
-
- writeTestString(" -- Results arrays have one member");
- AlwaysAssert(converged.size() == 1, AipsError);
-
- writeTestString(" -- Fit converged");
- AlwaysAssert(converged[0], AipsError);
-
- writeTestString(" -- Only one gaussian fit");
- AlwaysAssert(((Vector<Int>)results.asArrayInt("ncomps"))[0] == 1,
AipsError);
-
- writeTestString(" -- It is a gaussian and not something else");
- AlwaysAssert(((Vector<String>)results.asArrayString("type0"))[0]
== "GAUSSIAN", AipsError);
-
- writeTestString(" -- Various tests of the fit values");
- AlwaysAssert((((Vector<Double>)results.asArrayDouble("amp0"))[0] -
(49.7) < 0.1), AipsError);
- AlwaysAssert((((Vector<Double>)results.asArrayDouble("ampErr0"))[0]
- (4.0) < 0.1), AipsError);
- AlwaysAssert((((Vector<Double>)results.asArrayDouble("center0"))[0]
- (-237.7) < 0.1), AipsError);
-
AlwaysAssert((((Vector<Double>)results.asArrayDouble("centerErr0"))[0] -
(1.7) < 0.1), AipsError);
- AlwaysAssert((((Vector<Double>)results.asArrayDouble("fwhm0"))[0] -
(42.4) < 0.1), AipsError);
- AlwaysAssert((((Vector<Double>)results.asArrayDouble("fwhmErr0"))[0]
- (4.0) < 0.1), AipsError);
-
- writeTestString(" -- Test of fit units");
- AlwaysAssert(results.asString("xUnit") == "km/s", AipsError);
- AlwaysAssert(results.asString("yUnit") == "Jy", AipsError);
- }
- {
- writeTestString("test results of multi-pixel two gaussian fit");
- ImageProfileFitter fitter(
- goodImage, "", "", "", "", "", 2, True,
- "", "", 2, -1
- );
- Record results = fitter.fit();
-
- writeTestString("t -- test correct number of fits performed");
- Vector<Bool> converged = results.asArrayBool("converged");
- AlwaysAssert(converged.size() == 81, AipsError);
-
- writeTestString(" -- test all fits converged");
- for (uInt i=0; i<converged.size(); i++) {
- AlwaysAssert(converged[i], AipsError);
- }
- writeTestString(
- " -- test all fits were for 2 gaussians except for the first which
was for 1"
- );
-
- Vector<Int> nComps = results.asArrayInt("ncomps");
- for (uInt i=0; i<nComps.size(); i++) {
- Int expected = i > 0 ? 2 : 1;
- AlwaysAssert(nComps[i] == expected, AipsError);
- }
-
- writeTestString(" -- Test fit component types");
- for(uInt i=0; i<2; i++) {
- String num = String::toString(i);
- Vector<String> types = results.asArrayString("type" + num);
- for (uInt j=0; j<types.size(); j++) {
- String expected = (i == 1 && j == 0) ? "" : "GAUSSIAN";
- AlwaysAssert(types[j] == expected, AipsError);
- }
- }
- writeTestString(" -- Test of fit units");
- AlwaysAssert(results.asString("xUnit") == "km/s", AipsError);
- AlwaysAssert(results.asString("yUnit") == "Jy", AipsError);
- }
- {
- writeTestString("test writing result images of multi-pixel two
gaussian fit");
- String center = dirName + "/center";
- String centerErr = dirName + "/centerErr";
- String fwhm = dirName + "/fwhm";
- String fwhmErr = dirName + "/fwhmErr";
- String amp = dirName + "/amp";
- String ampErr = dirName + "/ampErr";
-
- ImageProfileFitter fitter(
- goodImage, "", "", "", "", "", 2, True,
- "", "", 2, -1, amp, ampErr, center, centerErr,
- fwhm, fwhmErr
- );
- Record results = fitter.fit();
-
- Vector<String> names(12);
- names[0] = "center_0";
- names[1] = "centerErr_0";
- names[2] = "center_1";
- names[3] = "centerErr_1";
- names[4] = "amp_0";
- names[5] = "ampErr_0";
- names[6] = "amp_1";
- names[7] = "ampErr_1";
- names[8] = "fwhm_0";
- names[9] = "fwhmErr_0";
- names[10] = "fwhm_1";
- names[11] = "fwhmErr_1";
-
- for (uInt i=0; i<names.size(); i++) {
- checkImage(dirName + "/" + names[i], names[i] + ".fits");
- }
- }
-
- {
- writeTestString("test results of multi-pixel two gaussian, order 3
polynomial fit");
- ImageProfileFitter fitter(
- goodPolyImage, "", "", "", "", "", 2, True,
- "", "", 2, 3
- );
- Record results = fitter.fit();
-
- writeTestString("t -- test correct number of fits attempted");
- Vector<Bool> converged = results.asArrayBool("converged");
- AlwaysAssert(converged.size() == 81, AipsError);
-
- writeTestString(" -- test all but one fits converged");
- for (uInt i=0; i<converged.size(); i++) {
- if (i == 72) {
- AlwaysAssert(! converged[i], AipsError);
- }
- else {
- AlwaysAssert(converged[i], AipsError);
- }
- }
- writeTestString(" -- Test of fit units");
- AlwaysAssert(results.asString("xUnit") == "km/s", AipsError);
- AlwaysAssert(results.asString("yUnit") == "Jy", AipsError);
- }
-
- cout << endl << "All " << testNumber << " tests succeeded" << endl;
- cout << "ok" << endl;
- }
- catch (AipsError x) {
- cerr << "Exception caught: " << x.getMesg() << endl;
- retVal = 1;
- }
- workdir.removeRecursive();
-
- return retVal;
-}
-
=======================================
--- /trunk/images/CMakeLists.txt Mon Mar 26 23:10:00 2012
+++ /trunk/images/CMakeLists.txt Mon Apr 2 05:08:24 2012
@@ -1,6 +1,7 @@
#
# CASA Images
#
+
BISON_TARGET (ImageExprGram Images/ImageExprGram.yy
${CMAKE_CURRENT_BINARY_DIR}/ImageExprGram.ycc COMPILE_FLAGS "-p
ImageExprGram")
FLEX_TARGET (ImageExprGram Images/ImageExprGram.ll
${CMAKE_CURRENT_BINARY_DIR}/ImageExprGram.lcc
COMPILE_FLAGS "-PImageExprGram")

@@ -34,7 +35,6 @@
Images/ImageExprGram.cc
Images/MaskSpecifier.cc
Images/ImageMomentsProgress.cc
-Images/ImageAnalysis.cc
Images/FITSImage.cc
Images/FITSImgParser.cc
Images/PagedImage2.cc
@@ -45,7 +45,6 @@
Images/ImageAttrHandlerCasa.cc
Images/ImageAttrHandlerHDF5.cc
Images/ImagePolProxy.cc
-Images/ImageProfileFit.cc
Images/ImageExprParse.cc
Images/ImageOpener.cc
Images/LELImageCoord.cc
@@ -56,12 +55,7 @@
Images/HDF5Image2.cc
Images/ImageProxy.cc
Images/ComponentImager.cc
-Images/ImageCollapser.cc
-Images/ImageInputProcessor.cc
Images/ImageReorderer.cc
-Images/ImageFitter.cc
-Images/ImageProfileFitter.cc
-IO/FitterEstimatesFileParser.cc
Wnbt/ComponentUpdate.cc
${BISON_ImageExprGram_OUTPUTS}
${FLEX_ImageExprGram_OUTPUTS}
@@ -98,7 +92,6 @@
Images/HDF5Image.tcc
Images/Image2DConvolver.h
Images/Image2DConvolver.tcc
-Images/ImageAnalysis.h
Images/ImageAttrGroup.h
Images/ImageAttrGroupCasa.h
Images/ImageAttrGroupHDF5.h
@@ -109,7 +102,6 @@
Images/ImageConcat.tcc
Images/ImageConvolver.h
Images/ImageConvolver.tcc
-Images/ImageCollapser.h
Images/ImageDecomposer.h
Images/ImageDecomposer.tcc
Images/ImageExpr.h
@@ -119,13 +111,11 @@
Images/ImageFFT.h
Images/ImageFITSConverter.h
Images/ImageFITSConverter.tcc
-Images/ImageFitter.h
Images/ImageFit1D.h
Images/ImageFit1D.tcc
Images/ImageHistograms.h
Images/ImageHistograms.tcc
Images/ImageInfo.h
-Images/ImageInputProcessor.h
Images/ImageInterface.h
Images/ImageInterface.tcc
Images/ImageMetaData.h
@@ -135,8 +125,6 @@
Images/ImageOpener.h
Images/ImagePolProxy.h
Images/ImagePolarimetry.h
-Images/ImageProfileFit.h
-Images/ImageProfileFitter.h
Images/ImageProxy.h
Images/ImageRegrid.h
Images/ImageRegrid.tcc
@@ -195,11 +183,6 @@
DESTINATION include/casacore/images/Regions
)

-install (FILES
-IO/FitterEstimatesFileParser.h
-DESTINATION include/casacore/images/IO
-)
-
install (FILES
Regions.h
Images.h
@@ -211,3 +194,8 @@
add_subdirectory (Wnbt/test)
add_subdirectory (Regions/test)
endif (BUILD_TESTING)
+
+The MSSS HBA observations should only use the inner 24 tiles of the remote
stations, while the core stations should be used in dual mode. Thus need a
new mode HBA_DUAL_INNER needs to be supported.
+At the same time future modes HBA_JOINED_INNER, HBA_ZERO_INNER, and
HBA_ONE_INNER could be added as well.
+
+Add support for HBA_DUAL_INNER and other inner modes
=======================================
--- /trunk/images/Images/ImageUtilities2.tcc Wed Mar 21 05:55:04 2012
+++ /trunk/images/Images/ImageUtilities2.tcc Mon Apr 2 05:08:24 2012
@@ -43,7 +43,7 @@
#include <lattices/Lattices/TempLattice.h>
#include <lattices/Lattices/TiledLineStepper.h>
#include <lattices/Lattices/MaskedLatticeIterator.h>
-#include <components/SpectralComponents/SpectralElement.h>
+#include <components/SpectralComponents/PolynomialSpectralElement.h>
#include <casa/System/ProgressMeter.h>
#include <casa/Logging/LogIO.h>
#include <casa/Quanta/Unit.h>
@@ -237,7 +237,7 @@
throw AipsError(errMsg);
}

- SpectralElement polyEl(poly);
+ PolynomialSpectralElement polyEl(poly);

IPosition inTileShape = inImage.niceCursorShape();
TiledLineStepper stepper (inImage.shape(), inTileShape, axis);
=======================================
--- /trunk/images/Images/test/CMakeLists.txt Wed Mar 21 05:55:04 2012
+++ /trunk/images/Images/test/CMakeLists.txt Mon Apr 2 05:08:24 2012
@@ -74,7 +74,6 @@
dImageHistograms
dImageInterface
dImageMoments
-dImageProfileFit
dImageStatistics
dImageSummary
dPagedImage
@@ -87,7 +86,6 @@
tFITSExtImageII
tHDF5Image
tImageAttrHandler
-tImageCollapser
tImageConcat
tImageConvolver
tImageDecomposer
@@ -97,13 +95,10 @@
tImageExprGram
tImageExprParse_addDir
tImageFFT
-tImageFit1D
-tImageFitter
+#tImageFit1D
tImageInfo
-tImageInputProcessor
tImageMetaData
tImagePolarimetry
-tImageProfileFitter
tImageRegrid
tImageReorderer
tImageSourceFinder
=======================================
--- /trunk/images/apps/CMakeLists.txt Tue Oct 18 00:39:05 2011
+++ /trunk/images/apps/CMakeLists.txt Mon Apr 2 05:08:24 2012
@@ -1,5 +1,5 @@
-foreach(prog image2fits imagecalc imageregrid imageslice imcollapse imfit
imreorder)
+foreach(prog image2fits imagecalc imageregrid imageslice imreorder)
add_executable (${prog} ${prog}.cc)
target_link_libraries (${prog} casa_images)
install(TARGETS ${prog} DESTINATION bin)
-endforeach(prog image2fits imagecalc imageregrid imageslice imcollapse
imfit imreorder)
+endforeach(prog image2fits imagecalc imageregrid imageslice imreorder)