ImagePolarimetry.h

void rotationMeasure(ImageInterface<Float>*& rmPtr, ImageInterface<Float>*& rmErrPtr, ImageInterface<Float>*& pa0Ptr, ImageInterface<Float>*& pa0ErrPtr, ImageInterface<Float>*& nTurns, ImageInterface<Float>*& rChiSqPtr, PGPlotter& plotter, Int spectralAxis, Float rmMax, Float maxPaErr=1.0e30, Float sigma=-1.0, Float rmFg=0.0, Bool showProgress=False)

static ImageExpr<Float> depolarizationRatio (const ImageInterface<Float>& im1, const ImageInterface<Float>& im2, Bool debias, Float clip=10.0, Float sigma=-1.0)

static ImageExpr<Float> sigmaDepolarizationRatio (const ImageInterface<Float>& im1, const ImageInterface<Float>& im2, Bool debias, Float clip=10.0, Float sigma=-1.0)

Private Members

void cleanup()

void copyDataAndMask(ImageInterface<Float>& out, ImageInterface<Float>& in) const

Bool dealWithMask (Lattice<Bool>*& pMask, ImageInterface<Float>*& pIm, LogIO& os, const LogIO& type) const

void fiddleStokesCoordinate(ImageInterface<Float>& ie, Stokes::StokesTypes type) const

void fiddleStokesCoordinate(CoordinateSystem& cSys, Stokes::StokesTypes type) const

void fiddleStokesCoordinate(ImageInterface<Complex>& ie, Complex::StokesTypes type) const

void fiddleTimeCoordinate(ImageInterface<Complex>& ie, const Complex<Double>& f, Int coord) const

Quantum<Double> findCentralFrequency(const Coordinate& coord, Int shape) const

Bool findRotationMeasure (Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted, Float& pa0ErrFitted, Float& rChiSqFitted, Float& nTurns, const Vector<uInt>& sortidx, const Vector<Float>& wsq, const Vector<Float>& pa, const Array<Bool>& paMask, const Array<Float>& paerr, Float rmfg, Float rmmax, Float paErrMax, PGPlotter& plotter, const String& posString)

void findStokes()

Int findSpectralCoordinate(const CoordinateSystem& cSys, LogIO& os, Bool fail) const

void findFrequencyAxis (Int& spectralCoord, Int& fAxis, const CoordinateSystem& cSys, Int spectralAxis) const

void hasQU () const

LatticeExprNode makePolIntNode(LogIO& os, Bool debias, Float clip, Float sigma, Bool doLin, Bool doCirc)

ImageExpr<Float> makeStokesExpr(ImageInterface<Float>* imPtr, Stokes::StokesTypes type, const String& name) const

ImageInterface<Float>* makeSubImage (IPosition& blc, IPosition& trc, Int axis, Int pix) const

Bool rmLsqFit (Vector<Float>& pars, const Vector<Float>& wsq, const Vector<Float> pa, const Vector<Float>& paerr) const

Bool rmPrimaryFit (Float& nTurns, Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted, Float& pa0ErrFitted, Float& rChiSqFitted, const Vector<Float>& wsq, const Vector<Float>& pa, const Vector<Float>& paerr, Float rmmax, PGPlotter& plotter, const String& posString)

Bool rmSupplementaryFit (Float& nTurns, Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted, Float& pa0ErrFitted, Float& rChiSqFitted, const Vector<Float>& wsq, const Vector<Float>& pa, const Vector<Float>& paerr)

String stokesName (ImagePolarimetry::StokesTypes index) const

Stokes::StokesTypes stokesType (ImagePolarimetry::StokesTypes index) const

Float sigma (ImagePolarimetry::StokesTypes index, Float clip)

void subtractProfileMean (ImageInterface<Float>& im, uInt axis) const

Description

Review Status

Date Reviewed:: yyyy/mm/dd

Prerequisite

Etymology

Polarimetric analysis of Images

Synopsis

This class provides polarimetric image analysis capability. It takes an image with a Stokes axis (some combination of IQUV is needed) as its input.

Many functions return ImageExpr objects. These are read-only images.

Sometimes the standard deviation of the noise is needed. This is for debiasing polarized intensity images or working out error images. By default it is worked out for you with a clipped mean algorithm. However, you can provide sigma if you know it accurately. It should be the standard deviation of the noise in the absence of signal. You won't measure that very well from Stokes I if it is dynamic range limited. Better to get it from V or Q or U. When this class needs the standard deviation of the noise, it will try and get it from V or Q and U and finally I.

However, note that the functions sigmaStokes{I,Q,U,V} DO return the standard deviation of the noise for that specific Stokes type.

The ImageExpr objects returned have the brightness units and ImageInfo set. The MiscInfo (a permanent record) and logSink are not set.

Motivation

Basic image analysis capability

To Do

plotting for function rotationMeasure
some assessment of the curvature of pa-l**2

Member Description

enum StokesTypes

Stokes types

ImagePolarimetry (const ImageInterface<Float>& image)

Constructor. The input image must have a Stokes axis with some subset of I,Q,U, and V

ImagePolarimetry(const ImagePolarimetry& other)

Copy constructor (reference semantics)

virtual ~ImagePolarimetry ()

Destructor

ImagePolarimetry& operator=(const ImagePolarimetry& other)

Assignment operator (reference semantics)

void summary(LogIO& os) const

Summary. Just invokes the ImageSummary list function to summarize the header of the construction image

const ImageInterface<Float>* imageInterface() const

Get the ImageInterface pointer of the construction image Don't delete it !

CoordinateSystem coordinates() const

Get the CoordinateSystem of the construction image

IPosition shape() const

Get the shape of the construction image

Bool isMasked() const

Is the construction image masked ?

IPosition singleStokesShape(CoordinateSystem& cSys, Stokes::StokesTypes type) const

Get the shape and CoordinateSystem of an image for a single Stokes pixel Thus, if the construction image shape was [10,10,4,20] where axis 2 (shape 4) is the Stokes axis, this function would return [10,10,1,20] Specify the type of Stokes pixel you want.

ImageExpr<Complex> complexLinearPolarization ()

Complex linear polarization

ImageExpr<Complex> complexFractionalLinearPolarization ()

Complex fractional linear polarization

ImageExpr<Float> stokesI() const
Float sigmaStokesI (Float clip=10.0)

Get the Stokes I image and the standard deviation of the I image. This is worked out by first clipping outliers from the mean at the specified level.

ImageExpr<Float> stokesQ() const
Float sigmaStokesQ (Float clip=10.0)

Get the Stokes Q image and the standard deviation of the Q image. This is worked out by first clipping outliers from the mean at the specified level.

ImageExpr<Float> stokesU() const
Float sigmaStokesU (Float clip=10.0)

Get the Stokes U image and the standard deviation of the U image. This is worked out by first clipping outliers from the mean at the specified level.

ImageExpr<Float> stokesV() const
Float sigmaStokesV (Float clip=10.0)

Get the Stokes V image and the standard deviation of the V image. This is worked out by first clipping outliers from the mean at the specified level.

ImageExpr<Float> stokes(ImagePolarimetry::StokesTypes index) const
Float sigmaStokes (ImagePolarimetry::StokesTypes index, Float clip=10.0)

Get the specified Stokes image and the standard deviation of the image. This is worked out by first clipping outliers from the mean at the specified level.

Float sigma (Float clip=10.0)

Get the best estimate of the statistical noise. This gives you the standard deviation with outliers from the mean clipped first. The idea is to not be confused by source or dynamic range issues. Generally Stokes V is empty of sources (not always), then Q and U are generally less bright than I. So this function first tries V, then Q and U and lastly I to give you its noise estimate

ImageExpr<Float> linPolInt(Bool debias, Float clip=10.0, Float sigma=-1.0)
Float sigmaLinPolInt (Float clip=10.0, Float sigma=-1.0)

Get the linearly polarized intensity image and its standard deviation. If wish to debias the image, you can either provide sigma<\src> (the standard deviation of the termal noise ) or if <src>sigma is non-positive, it will be worked out for you with a clipped mean algorithm.

ImageExpr<Float> totPolInt(Bool debias, Float clip=10.0, Float sigma=-1.0)
Float sigmaTotPolInt (Float clip=10.0, Float sigma=-1.0)

Get the total polarized intensity (from whatever combination of Q, U, and V the construction image has) image and its error (standard deviation). If wish to debias the image, you can either provide sigma<\src> (the standard deviation of the thermal noise) or if <src>sigma is non-positive, it will be worked out for you with a clipped mean algorithm.

ImageExpr<Float> linPolPosAng(Bool radians) const
ImageExpr<Float> sigmaLinPolPosAng (Bool radians, Float clip=10.0, Float sigma=-1.0)

Get linearly polarized position angle (degrees or radians) image and error (standard deviation). If you provide sigma<\src> it is the standard deviation of the termal noise. If <src>sigma is non-positive, it will be worked out for you with a clipped mean algorithm.

ImageExpr<Float> fracLinPol(Bool debias, Float clip=10.0, Float sigma=-1.0)
ImageExpr<Float> sigmaFracLinPol (Float clip=10.0, Float sigma=-1.0)

Get fractional linear polarization image and error (standard deviation). If wish to debias the image, you can either provide sigma<\src> (the standard deviation of the termal noise) or if <src>sigma is non-positive, it will be worked out for you with a clipped mean algorithm.

ImageExpr<Float> fracTotPol(Bool debias, Float clip=10.0, Float sigma=-1.0)
ImageExpr<Float> sigmaFracTotPol (Float clip=10.0, Float sigma=-1.0)

Get Fractional total polarization and error (standard deviation) var<\src> is the standard deviation of the thermal noise. If <src>sigma is non-positive, it will be worked out for you with a clipped mean algorithm.

void fourierRotationMeasure(ImageInterface<Complex>& pol, Bool zeroZeroLag)

Fourier Rotation Measure. The output image is the complex polarization (Q + iU) with the spectral axis replaced by a RotationMeasure axis. The appropriate shape and CoordinateSystem must be obtained with function singleStokesShape (ask for type STokes::Plinear). Howeverm this function will replace the SpectralCoordinate by a LinearCoordinate describing the Rotation Measure. ImageInfo, and Units are copied to the output. MiscInfo and history are not. If the output has a mask, and the input is masked, the mask is copied. If the output has a mask, it should already have been initialized to True

IPosition rotationMeasureShape(CoordinateSystem& cSys, Int& frequencyAxis, Int& stokesAxis, LogIO& os, Int spectralAxis=-1) const

This function is used in concert with the rotationMeasure function. It tells you what the shape of the output RM image should be, and gives you its CoordinateSystem. Because the ImagePolarimetry construction image may house the frequencies coordinate description in a Spectral, Tabular or Linear coordinate, you may explicitly specify which axis is the Spectral axis (spectralAxis). By default, it tries to find the SpectralCoordinate. If there is none, it will look for Tabular or Linear coordinates with a "FREQ" label. It returns to you the frequencyAxis (i.e. the one it is concluded houses the frequency spectrum) and the stokesAxis that it finds.

IPosition positionAngleShape(CoordinateSystem& cSys, Int& frequencyAxis, Int& stokesAxis, LogIO& os, Int spectralAxis=-1) const

This function is used in concert with the rotationMeasure function. It tells you what the shape of the output Position Angle image should be, and gives you its CoordinateSystem. Because the ImagePolarimetry construction image may house the frequencies coordinate description in a Spectral, Tabular or Linear coordinate, you may explicitly specify which axis is the Spectral axis (spectralAxis). By default, it tries to find the SpectralCoordinate. If there is none, it will look for Tabular or Linear coordinates with a "FREQ" label.

void rotationMeasure(ImageInterface<Float>& rmPtr, ImageInterface<Float>& rmErrPtr, ImageInterface<Float>& pa0Ptr, ImageInterface<Float>& pa0ErrPtr, ImageInterface<Float>& nTurns, ImageInterface<Float>& rChiSqPtr, PGPlotter& plotter, Int spectralAxis, Float rmMax, Float maxPaErr=1.0e30, Float sigma=-1.0, Float rmFg=0.0, Bool showProgress=False)

This function applies a traditional (i.e. non-Fourier) Rotation Measure algorithm (Leahy et al, A&A, 156, 234) approach. For the RM images you must get the shape and CoordinateSYstem from function rotationMeasureShape. For the position angle images, use function singleStokesShape. Null pointer ImageInterface objects are not accessed so you can select which output images you want. Any output images not masked will be given a mask. The position angles are all in degrees. The RM images in rad/m/m. ImageInfo and Units, are copied to the output. MiscInfo and history are not. You specify which axis houses the frequencies, the noise level of Q and U if you know it (by default it is worked out for you) for error images, the value of the foreground RM if you know it (helps for unwinding ambiguity), the absolute maximum RM it should solve for, and the maximum error in the position angle that should be allowed. The state of the plotter should be set by the caller (e.g. character size, number of plots in x and y etc).

static ImageExpr<Float> depolarizationRatio (const ImageInterface<Float>& im1, const ImageInterface<Float>& im2, Bool debias, Float clip=10.0, Float sigma=-1.0)
static ImageExpr<Float> sigmaDepolarizationRatio (const ImageInterface<Float>& im1, const ImageInterface<Float>& im2, Bool debias, Float clip=10.0, Float sigma=-1.0)

Depolarization ratio image and error. Requires two images hence static functions.

void cleanup()

Delete all private pointers

void copyDataAndMask(ImageInterface<Float>& out, ImageInterface<Float>& in) const

Copy data and mask

Bool dealWithMask (Lattice<Bool>& pMask, ImageInterface<Float>& pIm, LogIO& os, const LogIO& type) const

For traiditional RM approach, give output a mask if possible

void fiddleStokesCoordinate(ImageInterface<Float>& ie, Stokes::StokesTypes type) const

Change the Stokes Coordinate for the given float image to be of the specified Stokes type

void fiddleStokesCoordinate(CoordinateSystem& cSys, Stokes::StokesTypes type) const

void fiddleStokesCoordinate(ImageInterface<Complex>& ie, Complex::StokesTypes type) const

Change the Stokes Coordinate for the given complex image to be of the specified Stokes type

void fiddleTimeCoordinate(ImageInterface<Complex>& ie, const Complex<Double>& f, Int coord) const

Change the time coordinate to be rotation measure

Quantum<Double> findCentralFrequency(const Coordinate& coord, Int shape) const

Find the central frequency from the given spectral coordinate

Bool findRotationMeasure (Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted, Float& pa0ErrFitted, Float& rChiSqFitted, Float& nTurns, const Vector<uInt>& sortidx, const Vector<Float>& wsq, const Vector<Float>& pa, const Array<Bool>& paMask, const Array<Float>& paerr, Float rmfg, Float rmmax, Float paErrMax, PGPlotter& plotter, const String& posString)

Fit the spectrum of position angles to find the rotation measure via Leahy algorithm

void findStokes()

Find the Stokes in the construction image and assign pointers

Int findSpectralCoordinate(const CoordinateSystem& cSys, LogIO& os, Bool fail) const

Find the spectral coordinate.

void findFrequencyAxis (Int& spectralCoord, Int& fAxis, const CoordinateSystem& cSys, Int spectralAxis) const

FInd frequency axis

void hasQU () const

So we have Q and U ? Excpetion if not

LatticeExprNode makePolIntNode(LogIO& os, Bool debias, Float clip, Float sigma, Bool doLin, Bool doCirc)

Make a LEN for the give types of polarized intensity

ImageExpr<Float> makeStokesExpr(ImageInterface<Float>* imPtr, Stokes::StokesTypes type, const String& name) const

Make an IE for the specified Stokes

ImageInterface<Float>* makeSubImage (IPosition& blc, IPosition& trc, Int axis, Int pix) const

Make a SubImage from the construction image for the specified pixel along the specified pixel axis

Bool rmLsqFit (Vector<Float>& pars, const Vector<Float>& wsq, const Vector<Float> pa, const Vector<Float>& paerr) const

Least squares fit to find RM from position angles

Bool rmPrimaryFit (Float& nTurns, Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted, Float& pa0ErrFitted, Float& rChiSqFitted, const Vector<Float>& wsq, const Vector<Float>& pa, const Vector<Float>& paerr, Float rmmax, PGPlotter& plotter, const String& posString)

Fit the spectrum of position angles to find the rotation measure via Leahy algorithm for primary (n>2) points

Bool rmSupplementaryFit (Float& nTurns, Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted, Float& pa0ErrFitted, Float& rChiSqFitted, const Vector<Float>& wsq, const Vector<Float>& pa, const Vector<Float>& paerr)

Fit the spectrum of position angles to find the rotation measure via Leahy algorithm for supplementary (n==2) points

String stokesName (ImagePolarimetry::StokesTypes index) const

Return I, Q, U or V for specified integer index (0-3)

Stokes::StokesTypes stokesType (ImagePolarimetry::StokesTypes index) const

Return I, Q, U or V for specified integer index (0-3)

Float sigma (ImagePolarimetry::StokesTypes index, Float clip)

Find the standard deviation for the Stokes image specified by the integer index

void subtractProfileMean (ImageInterface<Float>& im, uInt axis) const

Subtract profile mean from image