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