00001 //# ImagePolarimetry.h: Polarimetric analysis of images 00002 //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002 00003 //# Associated Universities, Inc. Washington DC, USA. 00004 //# 00005 //# This library is free software; you can redistribute it and/or modify it 00006 //# under the terms of the GNU Library General Public License as published by 00007 //# the Free Software Foundation; either version 2 of the License, or (at your 00008 //# option) any later version. 00009 //# 00010 //# This library is distributed in the hope that it will be useful, but WITHOUT 00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public 00013 //# License for more details. 00014 //# 00015 //# You should have received a copy of the GNU Library General Public License 00016 //# along with this library; if not, write to the Free Software Foundation, 00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. 00018 //# 00019 //# Correspondence concerning AIPS++ should be addressed as follows: 00020 //# Internet email: aips2-request@nrao.edu. 00021 //# Postal address: AIPS++ Project Office 00022 //# National Radio Astronomy Observatory 00023 //# 520 Edgemont Road 00024 //# Charlottesville, VA 22903-2475 USA 00025 //# 00026 //# $Id: ImagePolarimetry.h 19817 2006-12-22 05:24:00Z gvandiep $ 00027 00028 #ifndef IMAGES_IMAGEPOLARIMETRY_H 00029 #define IMAGES_IMAGEPOLARIMETRY_H 00030 00031 00032 //# Includes 00033 #include <casa/aips.h> 00034 #include <casa/Arrays/Vector.h> 00035 #include <casa/Containers/Block.h> 00036 #include <measures/Measures/Stokes.h> 00037 #include <casa/BasicSL/Complex.h> 00038 #include <images/Images/ImageInterface.h> 00039 #include <scimath/Fitting/LinearFitSVD.h> 00040 00041 00042 namespace casa { //# NAMESPACE CASA - BEGIN 00043 00044 //# Forward Declarations 00045 template <class T> class SubImage; 00046 template <class T> class ImageExpr; 00047 template <class T> class Quantum; 00048 template <class T> class LatticeStatistics; 00049 // 00050 class CoordinateSystem; 00051 class IPosition; 00052 class LatticeExprNode; 00053 class LCBox; 00054 class LogIO; 00055 00056 // <summary> 00057 // Polarimetric analysis of images 00058 // </summary> 00059 00060 // <use visibility=export> 00061 00062 // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos=""> 00063 // </reviewed> 00064 00065 // <prerequisite> 00066 // <li> <linkto class=ImageExpr>ImageExpr</linkto> 00067 // <li> <linkto class=ImageInterface>ImageInterface</linkto> 00068 // </prerequisite> 00069 00070 // <etymology> 00071 // Polarimetric analysis of Images 00072 // </etymology> 00073 00074 // <synopsis> 00075 // This class provides polarimetric image analysis capability. 00076 // It takes an image with a Stokes axis (some combination 00077 // of IQUV is needed) as its input. 00078 // 00079 // Many functions return ImageExpr objects. These are 00080 // read-only images. 00081 // 00082 // Sometimes the standard deviation of the noise is needed. 00083 // This is for debiasing polarized intensity images or working out 00084 // error images. By default it is worked out for you with a 00085 // clipped mean algorithm. However, you can provide sigma if you 00086 // know it accurately. It should be the standard deviation of the noise in 00087 // the absence of signal. You won't measure that very well from 00088 // Stokes I if it is dynamic range limited. Better to get it from 00089 // V or Q or U. When this class needs the standard deviation of 00090 // the noise, it will try and get it from V or Q and U and finally I. 00091 // 00092 // However, note that the functions sigmaStokes{I,Q,U,V} DO return the standard 00093 // deviation of the noise for that specific Stokes type. 00094 // 00095 // The ImageExpr objects returned have the brightness units and ImageInfo 00096 // set. The MiscInfo (a permanent record) and logSink are not set. 00097 // 00098 // </synopsis> 00099 // 00100 // <motivation> 00101 // Basic image analysis capability 00102 // </motivation> 00103 00104 // <todo asof="1999/11/01"> 00105 // <li> plotting for function rotationMeasure 00106 // <li> some assessment of the curvature of pa-l**2 00107 // </todo> 00108 00109 00110 class ImagePolarimetry 00111 { 00112 public: 00113 00114 // Stokes types 00115 enum StokesTypes {I, Q, U, V}; 00116 00117 // Constructor. The input image must have a Stokes 00118 // axis with some subset of I,Q,U, and V 00119 ImagePolarimetry (const ImageInterface<Float>& image); 00120 00121 // Copy constructor (reference semantics) 00122 ImagePolarimetry(const ImagePolarimetry& other); 00123 00124 // Destructor 00125 virtual ~ImagePolarimetry (); 00126 00127 // Assignment operator (reference semantics) 00128 ImagePolarimetry& operator=(const ImagePolarimetry& other); 00129 00130 // Summary. Just invokes the ImageSummary list function 00131 // to summarize the header of the construction image 00132 void summary(LogIO& os) const; 00133 00134 // Get the ImageInterface pointer of the construction image 00135 // Don't delete it ! 00136 const ImageInterface<Float>* imageInterface() const {return itsInImagePtr;}; 00137 00138 // Get the CoordinateSystem of the construction image 00139 CoordinateSystem coordinates() const {return itsInImagePtr->coordinates();}; 00140 00141 // Get the shape of the construction image 00142 IPosition shape() const {return itsInImagePtr->shape();}; 00143 00144 // Is the construction image masked ? 00145 Bool isMasked() const {return itsInImagePtr->isMasked();}; 00146 00147 // Get the shape and CoordinateSystem of an image for a single Stokes pixel 00148 // Thus, if the construction image shape was [10,10,4,20] where 00149 // axis 2 (shape 4) is the Stokes axis, this function 00150 // would return [10,10,1,20] Specify the type of Stokes pixel 00151 // you want. 00152 IPosition singleStokesShape(CoordinateSystem& cSys, Stokes::StokesTypes type) const; 00153 00154 // Complex linear polarization 00155 ImageExpr<Complex> complexLinearPolarization (); 00156 00157 // Complex fractional linear polarization 00158 ImageExpr<Complex> complexFractionalLinearPolarization (); 00159 00160 // Get the Stokes I image and the standard deviation of the 00161 // I image. This is worked out by first clipping 00162 // outliers from the mean at the specified level. 00163 // <group> 00164 ImageExpr<Float> stokesI() const; 00165 Float sigmaStokesI (Float clip=10.0); 00166 // </group> 00167 00168 // Get the Stokes Q image and the standard deviation 00169 // of the Q image. This is worked out by first clipping 00170 // outliers from the mean at the specified level. 00171 // <group> 00172 ImageExpr<Float> stokesQ() const; 00173 Float sigmaStokesQ (Float clip=10.0); 00174 // </group> 00175 00176 // Get the Stokes U image and the standard deviation 00177 // of the U image. This is worked out by first clipping 00178 // outliers from the mean at the specified level. 00179 // <group> 00180 ImageExpr<Float> stokesU() const; 00181 Float sigmaStokesU (Float clip=10.0); 00182 // </group> 00183 00184 // Get the Stokes V image and the standard deviation 00185 // of the V image. This is worked out by first clipping 00186 // outliers from the mean at the specified level. 00187 // <group> 00188 ImageExpr<Float> stokesV() const; 00189 Float sigmaStokesV (Float clip=10.0); 00190 // </group> 00191 00192 // Get the specified Stokes image and the standard deviation 00193 // of the image. This is worked out by first clipping 00194 // outliers from the mean at the specified level. 00195 // <group> 00196 ImageExpr<Float> stokes(ImagePolarimetry::StokesTypes index) const; 00197 Float sigmaStokes (ImagePolarimetry::StokesTypes index, Float clip=10.0); 00198 // </group> 00199 00200 // Get the best estimate of the statistical noise. This gives you 00201 // the standard deviation with outliers from the mean 00202 // clipped first. The idea is to not be confused by source or dynamic range issues. 00203 // Generally Stokes V is empty of sources (not always), then Q and U are generally 00204 // less bright than I. So this function first tries V, then Q and U 00205 // and lastly I to give you its noise estimate 00206 Float sigma (Float clip=10.0); 00207 00208 // Get the linearly polarized intensity image and its 00209 // standard deviation. If wish to debias the image, you 00210 // can either provide <src>sigma</src> (the standard 00211 // deviation of the termal noise ) or if <src>sigma</src> is non-positive, 00212 // it will be worked out for you with a clipped mean algorithm. 00213 // <group> 00214 ImageExpr<Float> linPolInt(Bool debias, Float clip=10.0, Float sigma=-1.0); 00215 Float sigmaLinPolInt (Float clip=10.0, Float sigma=-1.0); 00216 // </group> 00217 00218 // Get the total polarized intensity (from whatever combination 00219 // of Q, U, and V the construction image has) image and its error 00220 // (standard deviation). If wish to debias the image, you 00221 // can either provide <src>sigma</src> (the standard deviation 00222 // of the thermal noise) or if <src>sigma</src> is 00223 // non-positive, it will be worked out for you with a 00224 // clipped mean algorithm. 00225 // <group> 00226 ImageExpr<Float> totPolInt(Bool debias, Float clip=10.0, Float sigma=-1.0); 00227 Float sigmaTotPolInt (Float clip=10.0, Float sigma=-1.0); 00228 // </group> 00229 00230 // Get linearly polarized position angle (degrees or radians) image 00231 // and error (standard deviation). If you provide 00232 // <src>sigma</src> it is the standard deviation of 00233 // the termal noise. If <src>sigma</src> is non-positive, it will be 00234 // worked out for you with a clipped mean algorithm. 00235 // <group> 00236 ImageExpr<Float> linPolPosAng(Bool radians) const; 00237 ImageExpr<Float> sigmaLinPolPosAng (Bool radians, Float clip=10.0, Float sigma=-1.0); 00238 // </group> 00239 00240 // Get fractional linear polarization image 00241 // and error (standard deviation). If wish to debias the image, you 00242 // can either provide <src>sigma</src> (the standard 00243 // deviation of the termal noise) or if <src>sigma</src> is non-positive, 00244 // it will be worked out for you with a clipped mean algorithm. 00245 // <group> 00246 ImageExpr<Float> fracLinPol(Bool debias, Float clip=10.0, Float sigma=-1.0); 00247 ImageExpr<Float> sigmaFracLinPol (Float clip=10.0, Float sigma=-1.0); 00248 // </group> 00249 00250 // Get Fractional total polarization and error (standard deviation) 00251 // <src>var</src> is the standard deviation of the thermal noise. 00252 // If <src>sigma</src> is non-positive, 00253 // it will be worked out for you with a clipped mean algorithm. 00254 // <group> 00255 ImageExpr<Float> fracTotPol(Bool debias, Float clip=10.0, Float sigma=-1.0); 00256 ImageExpr<Float> sigmaFracTotPol (Float clip=10.0, Float sigma=-1.0); 00257 // </group> 00258 00259 // Fourier Rotation Measure. The output image is the complex polarization 00260 // (Q + iU) with the spectral axis replaced by a RotationMeasure axis. 00261 // The appropriate shape and CoordinateSystem must be obtained with 00262 // function singleStokesShape (ask for type STokes::Plinear). 00263 // Howeverm this function will replace the SpectralCoordinate 00264 // by a LinearCoordinate describing the Rotation Measure. 00265 // ImageInfo, and Units are copied to the output. MiscInfo and 00266 // history are not. If the output has a mask, 00267 // and the input is masked, the mask is copied. If the output 00268 // has a mask, it should already have been initialized to True 00269 void fourierRotationMeasure(ImageInterface<Complex>& pol, 00270 Bool zeroZeroLag); 00271 00272 // This function is used in concert with the rotationMeasure function. 00273 // It tells you what the shape of the output RM image should be, and 00274 // gives you its CoordinateSystem. Because the ImagePolarimetry 00275 // construction image may house the frequencies coordinate description 00276 // in a Spectral, Tabular or Linear coordinate, you may explicitly 00277 // specify which axis is the Spectral axis (spectralAxis). By default, 00278 // it tries to find the SpectralCoordinate. If there is none, it will 00279 // look for Tabular or Linear coordinates with a "FREQ" label. 00280 // It returns to you the frequencyAxis (i.e. the one it is concluded 00281 // houses the frequency spectrum) and the stokesAxis that it 00282 // finds. 00283 IPosition rotationMeasureShape(CoordinateSystem& cSys, 00284 Int& frequencyAxis, Int& stokesAxis, 00285 LogIO& os, Int spectralAxis=-1) const; 00286 00287 // This function is used in concert with the rotationMeasure function. 00288 // It tells you what the shape of the output Position Angle image should be, and 00289 // gives you its CoordinateSystem. Because the ImagePolarimetry 00290 // construction image may house the frequencies coordinate description 00291 // in a Spectral, Tabular or Linear coordinate, you may explicitly 00292 // specify which axis is the Spectral axis (spectralAxis). By default, 00293 // it tries to find the SpectralCoordinate. If there is none, it will 00294 // look for Tabular or Linear coordinates with a "FREQ" label. 00295 IPosition positionAngleShape(CoordinateSystem& cSys, 00296 Int& frequencyAxis, Int& stokesAxis, 00297 LogIO& os, Int spectralAxis=-1) const; 00298 00299 // This function applies a traditional (i.e. non-Fourier) Rotation Measure 00300 // algorithm (Leahy et al, A&A, 156, 234) approach. For the RM images 00301 // you must get the shape and CoordinateSYstem from function 00302 // rotationMeasureShape. For the position angle images, use function 00303 // singleStokesShape. Null pointer ImageInterface objects are 00304 // not accessed so you can select which output images you want. 00305 // Any output images not masked will be given a mask. 00306 // The position angles are all in degrees. The RM images in rad/m/m. 00307 // ImageInfo and Units, are copied to the output. MiscInfo and history are not. 00308 // You specify which axis houses the frequencies, the noise level of 00309 // Q and U if you know it (by default it is worked out for you) 00310 // for error images, the value of the foreground RM if you know it 00311 // (helps for unwinding ambiguity), the absolute maximum RM it should 00312 // solve for, and the maximum error in the position angle that should 00313 // be allowed. The state of the plotter should be set by the caller 00314 // (e.g. character size, number of plots in x and y etc). 00315 void rotationMeasure(ImageInterface<Float>*& rmPtr, ImageInterface<Float>*& rmErrPtr, 00316 ImageInterface<Float>*& pa0Ptr, ImageInterface<Float>*& pa0ErrPtr, 00317 ImageInterface<Float>*& nTurns, ImageInterface<Float>*& rChiSqPtr, 00318 Int spectralAxis, Float rmMax, Float maxPaErr=1.0e30, 00319 Float sigma=-1.0, Float rmFg=0.0, Bool showProgress=False); 00320 00321 // Depolarization ratio image and error. Requires two images hence static 00322 // functions. 00323 // <group> 00324 static ImageExpr<Float> depolarizationRatio (const ImageInterface<Float>& im1, 00325 const ImageInterface<Float>& im2, 00326 Bool debias, Float clip=10.0, Float sigma=-1.0); 00327 00328 static ImageExpr<Float> sigmaDepolarizationRatio (const ImageInterface<Float>& im1, 00329 const ImageInterface<Float>& im2, 00330 Bool debias, Float clip=10.0, Float sigma=-1.0); 00331 // </group> 00332 00333 private: 00334 const ImageInterface<Float>* itsInImagePtr; 00335 LinearFitSVD<Float>* itsFitterPtr; 00336 Float itsOldClip; 00337 00338 // These blocks are always size 4, with IQUV in slots 0,1,2,3 00339 // If your image is IV only, they still use slots 0 and 3 00340 00341 PtrBlock<ImageInterface<Float>* > itsStokesPtr; 00342 PtrBlock<LatticeStatistics<Float>* > itsStokesStatsPtr; 00343 00344 Matrix<Bool> _beamsEqMat; 00345 00346 00347 // Delete all private pointers 00348 void cleanup(); 00349 00350 // Copy data and mask 00351 void copyDataAndMask(ImageInterface<Float>& out, 00352 ImageInterface<Float>& in) const; 00353 00354 // For traiditional RM approach, give output a mask if possible 00355 Bool dealWithMask (Lattice<Bool>*& pMask, ImageInterface<Float>*& pIm, LogIO& os, const String& type) const; 00356 00357 // Change the Stokes Coordinate for the given float image to be of the specified Stokes type 00358 void fiddleStokesCoordinate(ImageInterface<Float>& ie, Stokes::StokesTypes type) const; 00359 void fiddleStokesCoordinate(CoordinateSystem& cSys, Stokes::StokesTypes type) const; 00360 00361 // Change the Stokes Coordinate for the given complex image to be of the specified Stokes type 00362 void fiddleStokesCoordinate(ImageInterface<Complex>& ie, Stokes::StokesTypes type) const; 00363 00364 // Change the time coordinate to be rotation measure 00365 void fiddleTimeCoordinate(ImageInterface<Complex>& ie, const Quantum<Double>& f, Int coord) const; 00366 00367 // Find the central frequency from the given spectral coordinate 00368 Quantum<Double> findCentralFrequency(const Coordinate& coord, Int shape) const; 00369 00370 // Fit the spectrum of position angles to find the rotation measure via Leahy algorithm 00371 Bool findRotationMeasure (Float& rmFitted, Float& rmErrFitted, 00372 Float& pa0Fitted, Float& pa0ErrFitted, Float& rChiSqFitted, 00373 Float& nTurns, 00374 const Vector<uInt>& sortidx, const Vector<Float>& wsq, 00375 const Vector<Float>& pa, 00376 const Array<Bool>& paMask, 00377 const Array<Float>& paerr, 00378 Float rmfg, Float rmmax, Float paErrMax, 00379 const String& posString); 00380 00381 // Find the Stokes in the construction image and assign pointers 00382 void findStokes(); 00383 00384 // Find the spectral coordinate. 00385 Int findSpectralCoordinate(const CoordinateSystem& cSys, LogIO& os, Bool fail) const; 00386 00387 // FInd frequency axis 00388 void findFrequencyAxis (Int& spectralCoord, Int& fAxis, 00389 const CoordinateSystem& cSys, Int spectralAxis) const; 00390 // So we have Q and U ? Excpetion if not 00391 void hasQU () const; 00392 00393 // Make a LEN for the give types of polarized intensity 00394 LatticeExprNode makePolIntNode(LogIO& os, Bool debias, Float clip, Float sigma, 00395 Bool doLin, Bool doCirc); 00396 00397 // Make an IE for the specified Stokes 00398 ImageExpr<Float> makeStokesExpr(ImageInterface<Float>* imPtr, 00399 Stokes::StokesTypes type, const String& name) const; 00400 00401 // Make a SubImage from the construction image for the specified pixel 00402 // along the specified pixel axis 00403 ImageInterface<Float>* makeSubImage (IPosition& blc, IPosition& trc, 00404 Int axis, Int pix) const; 00405 00406 // Least squares fit to find RM from position angles 00407 Bool rmLsqFit (Vector<Float>& pars, const Vector<Float>& wsq, 00408 const Vector<Float> pa, const Vector<Float>& paerr) const; 00409 00410 // Fit the spectrum of position angles to find the rotation measure via Leahy algorithm 00411 // for primary (n>2) points 00412 Bool rmPrimaryFit (Float& nTurns, Float& rmFitted, Float& rmErrFitted, 00413 Float& pa0Fitted, Float& pa0ErrFitted, 00414 Float& rChiSqFitted, const Vector<Float>& wsq, 00415 const Vector<Float>& pa, const Vector<Float>& paerr, 00416 Float rmmax, const String& posString); 00417 00418 // Fit the spectrum of position angles to find the rotation measure via Leahy algorithm 00419 // for supplementary (n==2) points 00420 Bool rmSupplementaryFit (Float& nTurns, Float& rmFitted, Float& rmErrFitted, 00421 Float& pa0Fitted, Float& pa0ErrFitted, 00422 Float& rChiSqFitted, const Vector<Float>& wsq, 00423 const Vector<Float>& pa, const Vector<Float>& paerr); 00424 00425 // Return I, Q, U or V for specified integer index (0-3) 00426 String stokesName (ImagePolarimetry::StokesTypes index) const; 00427 00428 // Return I, Q, U or V for specified integer index (0-3) 00429 Stokes::StokesTypes stokesType (ImagePolarimetry::StokesTypes index) const; 00430 00431 // Find the standard deviation for the Stokes image specified by the integer index 00432 Float sigma (ImagePolarimetry::StokesTypes index, Float clip); 00433 00434 // Subtract profile mean from image 00435 void subtractProfileMean (ImageInterface<Float>& im, uInt axis) const; 00436 00437 void _createBeamsEqMat(); 00438 00439 Bool _checkBeams( 00440 const Vector<StokesTypes>& stokes, 00441 const Bool requireChannelEquality, 00442 const Bool throws=True 00443 ) const; 00444 00445 Bool _checkIQUBeams( 00446 const Bool requireChannelEquality, 00447 const Bool throws=True 00448 ) const; 00449 00450 Bool _checkIVBeams( 00451 const Bool requireChannelEquality, 00452 const Bool throws=True 00453 ) const; 00454 00455 Bool _checkQUBeams( 00456 const Bool requireChannelEquality, 00457 const Bool throws=True 00458 ) const; 00459 00460 static void _checkBeams( 00461 const ImagePolarimetry& im1, const ImagePolarimetry& im2, 00462 const Vector<StokesTypes>& stokes 00463 ); 00464 00465 void _setInfo(ImageInterface<Complex>& im, const StokesTypes stokes) const; 00466 00467 void _setInfo(ImageInterface<Float>& im, const StokesTypes stokes) const; 00468 00469 void _setDoLinDoCirc(Bool& doLin, Bool& doCirc) const; 00470 00471 }; 00472 00473 00474 } //# NAMESPACE CASA - END 00475 00476 #endif