Line data    Source code

       1             : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
       2             : //# Associated Universities, Inc. Washington DC, USA.
       3             : //#
       4             : //# This library is free software; you can redistribute it and/or modify it
       5             : //# under the terms of the GNU Library General Public License as published by
       6             : //# the Free Software Foundation; either version 2 of the License, or (at your
       7             : //# option) any later version.
       8             : //#
       9             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      10             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      11             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      12             : //# License for more details.
      13             : //#
      14             : //# You should have received a copy of the GNU Library General Public License
      15             : //# along with this library; if not, write to the Free Software Foundation,
      16             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      17             : //#
      18             : //# Correspondence concerning AIPS++ should be addressed as follows:
      19             : //#        Internet email: aips2-request@nrao.edu.
      20             : //#        Postal address: AIPS++ Project Office
      21             : //#                        National Radio Astronomy Observatory
      22             : //#                        520 Edgemont Road
      23             : //#                        Charlottesville, VA 22903-2475 USA
      24             : //#
      25             : 
      26             : #include <imageanalysis/ImageAnalysis/ImagePolTask.h>
      27             : 
      28             : #include <casacore/lattices/LRegions/LCSlicer.h>
      29             : 
      30             : using namespace casacore;
      31             : namespace casa {
      32             : 
      33             : const String ImagePolTask::CLASS_NAME = "ImagePolTask";
      34             : 
      35          22 : ImagePolTask::ImagePolTask(
      36             :     const SPCIIF image, const String& outname, Bool overwrite
      37          22 : ) : ImageTask<Float>(image, nullptr, "", outname, overwrite) {
      38          22 :     _findStokes();
      39          22 :     _createBeamsEqMat();
      40          22 : }
      41             : 
      42          22 : ImagePolTask::~ImagePolTask() {}
      43             : 
      44           0 : String ImagePolTask::getClass() const {
      45           0 :     return CLASS_NAME;
      46             : }
      47             : 
      48           0 : Float ImagePolTask::sigmaLinPolInt(Float clip, Float sigma) {
      49           0 :     *_getLog() <<LogOrigin(CLASS_NAME, __func__, WHERE);
      50           0 :     ThrowIf(
      51             :         ! _stokesImage[Q] && ! _stokesImage[U]==0,
      52             :         "This image does not have Stokes Q and U so cannot provide linear polarization"
      53             :     );
      54           0 :     _checkQUBeams(False);
      55           0 :     Float sigma2 = 0.0;
      56           0 :     if (sigma > 0) {
      57           0 :         sigma2 = sigma;
      58             :     }
      59             :     else {
      60           0 :        Float sq = _sigma(Q, clip);
      61           0 :        Float su = _sigma(U, clip);
      62           0 :        sigma2 = (sq+su)/2.0;
      63           0 :        *_getLog() << LogIO::NORMAL << "Determined noise from Q&U images to be "
      64           0 :            << sigma2 << LogIO::POST;
      65             :    }
      66           0 :    return sigma2;
      67             : }
      68             : 
      69          26 : Bool ImagePolTask::_checkBeams(
      70             :     const std::vector<StokesTypes>& stokes,
      71             :     Bool requireChannelEquality, Bool throws
      72             : ) const {
      73          84 :     for (const auto i: stokes) {
      74         192 :         for (const auto j: stokes) {
      75         134 :             if (i == j) {
      76          58 :                 continue;
      77             :             }
      78          76 :             if (! _beamsEqMat(i, j)) {
      79           0 :                 ThrowIf(
      80             :                     throws,
      81             :                     "Input image has multiple beams and the corresponding "
      82             :                     "beams for the stokes planes necessary for this computation "
      83             :                     "are not equal."
      84             :                 );
      85           0 :                 return False;
      86             :             }
      87             :         }
      88             :     }
      89          52 :     auto image0 = _stokesImage[stokes[0]];
      90          26 :     if (
      91             :         requireChannelEquality
      92           0 :         && image0->coordinates().hasSpectralAxis()
      93          26 :         && image0->imageInfo().hasMultipleBeams()
      94             :     ) {
      95           0 :         const auto& beamSet = image0->imageInfo().getBeamSet().getBeams();
      96           0 :         auto firstBeam = *(beamSet.begin());
      97           0 :         for (const auto& beam: beamSet) {
      98           0 :             if (beam != firstBeam) {
      99           0 :                 ThrowIf(
     100             :                     throws, "At least one beam in this image is not equal "
     101             :                     "to all the others along its spectral axis so this "
     102             :                     "computation cannot be performed"
     103             :                 );
     104           0 :                 return False;
     105             :             }
     106             :         }
     107             :     }
     108          26 :     return True;
     109             : }
     110             : 
     111          22 : void ImagePolTask::_createBeamsEqMat() {
     112          22 :     const auto hasMultiBeams = _getImage()->imageInfo().hasMultipleBeams();
     113         110 :     for (uInt i=0; i<4; ++i) {
     114         308 :         for (uInt j=i; j<4; ++j) {
     115         220 :             if (! _stokesImage[i] || ! _stokesImage[j]) {
     116          94 :                 _beamsEqMat(i, j) = False;
     117             :             }
     118         126 :             else if (i == j) {
     119          62 :                 _beamsEqMat(i, j) = True;
     120             :             }
     121          64 :             else if (hasMultiBeams) {
     122           0 :                 _beamsEqMat(i, j) = (
     123           0 :                     _stokesImage[i]->imageInfo().getBeamSet()
     124           0 :                     == _stokesImage[j]->imageInfo().getBeamSet()
     125             :                 );
     126             :             }
     127             :             else {
     128          64 :                 _beamsEqMat(i, j) = True;
     129             :             }
     130         220 :             _beamsEqMat(j, i) = _beamsEqMat(i, j);
     131             :         }
     132             :     }
     133          22 : }
     134             : 
     135           6 : casacore::Bool ImagePolTask::_checkIQUBeams(
     136             :     Bool requireChannelEquality, Bool throws
     137             : ) const {
     138          12 :     std::vector<StokesTypes> types { I, Q, U };
     139          12 :     return _checkBeams(types, requireChannelEquality, throws);
     140             : }
     141             : 
     142           4 : Bool ImagePolTask::_checkIVBeams(
     143             :     Bool requireChannelEquality, Bool throws
     144             : ) const {
     145           8 :     std::vector<StokesTypes> types {I, V};
     146           8 :     return _checkBeams(types, requireChannelEquality, throws);
     147             : }
     148             : 
     149          16 : Bool ImagePolTask::_checkQUBeams(
     150             :     Bool requireChannelEquality, Bool throws
     151             : ) const {
     152          32 :     std::vector<StokesTypes> types {Q, U};
     153          32 :     return _checkBeams(types, requireChannelEquality, throws);
     154             : }
     155             : 
     156          22 : void ImagePolTask::_fiddleStokesCoordinate(
     157             :     ImageInterface<Float>& im, Stokes::StokesTypes type
     158             : ) const {
     159          44 :    CoordinateSystem cSys = im.coordinates();
     160          22 :    _fiddleStokesCoordinate(cSys, type);
     161          22 :    im.setCoordinateInfo(cSys);
     162          22 : }
     163             : 
     164          22 : void ImagePolTask::_fiddleStokesCoordinate(
     165             :     CoordinateSystem& cSys, Stokes::StokesTypes type
     166             : ) const {
     167          22 :     Int afterCoord = -1;
     168          22 :     Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord);
     169          44 :     Vector<Int> which(1);
     170          22 :     which(0) = Int(type);
     171          44 :     StokesCoordinate stokes(which);
     172          22 :     cSys.replaceCoordinate(stokes, iStokes);
     173          22 : }
     174             : 
     175             : 
     176          22 : void ImagePolTask::_findStokes() {
     177          22 :     *_getLog() << LogOrigin(CLASS_NAME, __func__, WHERE);
     178             :     // Do we have any Stokes ?
     179          22 :     const auto& csys = _getImage()->coordinates();
     180          22 :     ThrowIf(
     181             :         ! csys.hasPolarizationCoordinate(),
     182             :         "There is no Stokes Coordinate in this image"
     183             :     );
     184             :     // Find the pixel axis of the image which is Stokes
     185          22 :     auto stokesAxis = csys.polarizationAxisNumber();
     186             :     // Make the regions
     187          22 :     const auto& stokes = csys.stokesCoordinate();
     188          22 :     const auto ndim = _getImage()->ndim();
     189          44 :     const auto shape = _getImage()->shape();
     190          44 :     IPosition blc(ndim, 0);
     191          44 :     auto trc = shape - 1;
     192             :     Int pix;
     193          22 :     if (stokes.toPixel(pix, Stokes::I)) {
     194          10 :         _stokesImage[I] = _makeSubImage(blc, trc, stokesAxis, pix);
     195             :     }
     196          22 :     if (stokes.toPixel(pix, Stokes::Q)) {
     197          22 :         _stokesImage[Q] = _makeSubImage(blc, trc, stokesAxis, pix);
     198             :     }
     199          22 :     if (stokes.toPixel(pix, Stokes::U)) {
     200          22 :         _stokesImage[U] = _makeSubImage(blc, trc, stokesAxis, pix);
     201             :     }
     202          22 :     if (stokes.toPixel(pix, Stokes::V)) {
     203           8 :         _stokesImage[V] = _makeSubImage(blc, trc, stokesAxis, pix);
     204             :     }
     205          22 :     ThrowIf (
     206             :         (_stokesImage[Q] && ! _stokesImage[U])
     207             :         || (! _stokesImage[Q] && _stokesImage[U]),
     208             : 
     209             :         "This Stokes coordinate has only one of Q and U. This is not useful"
     210             :     );
     211          22 :     ThrowIf(
     212             :         ! _stokesImage[Q] && ! _stokesImage[U] && ! _stokesImage[V],
     213             :         "This image has no Stokes Q, U, or V.  This is not useful"
     214             :     );
     215          22 : }
     216             : 
     217          32 : SPCIIF ImagePolTask::_getStokesImage(StokesTypes type) const {
     218          32 :     return _stokesImage[type];
     219             : }
     220             : 
     221          15 : LatticeExprNode ImagePolTask::_makePolIntNode(
     222             :     Bool debias, Float clip, Float sigma,
     223             :     Bool doLin, Bool doCirc
     224             : ) {
     225          30 :     LatticeExprNode linNode, circNode;
     226          15 :     if (doLin) {
     227          30 :         linNode = LatticeExprNode(
     228          30 :             pow(*_stokesImage[U],2) + pow(*_stokesImage[Q],2)
     229          15 :         );
     230             :     }
     231          15 :     if (doCirc) {
     232           6 :         circNode = LatticeExprNode(pow(*_stokesImage[V],2));
     233             :     }
     234          15 :     auto node = (doLin && doCirc) ? linNode + circNode
     235          30 :         : doLin ? linNode : circNode;
     236          15 :     if (debias) {
     237           1 :         auto sigma2 = sigma > 0 ? sigma : _sigma(clip);
     238           1 :         node = node - LatticeExprNode(sigma2*sigma2);
     239             :         // node = iif(node >= 0, node, 0);
     240           2 :         *_getLog() << LogIO::NORMAL << "Debiasing with sigma = "
     241           1 :             << sigma2 << LogIO::POST;
     242             :     }
     243          30 :     return sqrt(node);
     244             : }
     245             : 
     246          62 : SPIIF ImagePolTask::_makeSubImage(
     247             :     IPosition& blc, IPosition& trc, Int axis, Int pix
     248             : ) const {
     249          62 :     blc(axis) = pix;
     250          62 :     trc(axis) = pix;
     251         124 :     LCSlicer slicer(blc, trc, RegionType::Abs);
     252          62 :     ImageRegion region(slicer);
     253         124 :     return SPIIF(new SubImage<Float>(*_getImage(), region));
     254             : }
     255             : 
     256           1 : void ImagePolTask::_maskAndZeroNaNs(SPIIF out) {
     257           2 :     auto isnan = isNaN(*out);
     258           1 :     if (any(isnan).getBool()) {
     259           0 :         LatticeExpr<Bool> mask(! isnan);
     260           0 :         if (! out->hasPixelMask()) {
     261           0 :             String x;
     262           0 :             ImageMaskAttacher::makeMask(*out, x, True, True, *_getLog(), False);
     263             :         }
     264           0 :         auto& pixelMask = out->pixelMask();
     265           0 :         LatticeExpr<Bool> newMask(mask && LatticeExpr<Bool>(out->pixelMask()));
     266           0 :         pixelMask.copyData(newMask);
     267           0 :         out->copyData(LatticeExpr<Float>(iif(isnan, 0, *out)));
     268             :     }
     269           1 : }
     270             : 
     271          12 : void ImagePolTask::_setDoLinDoCirc(Bool& doLin, Bool& doCirc, Bool requireI) const {
     272          12 :     *_getLog() << LogOrigin(CLASS_NAME, __func__, WHERE);
     273          12 :     doLin = _stokesImage[Q] && _stokesImage[U];
     274          12 :     doCirc = Bool(_stokesImage[V]);
     275          12 :     AlwaysAssert((doLin||doCirc), AipsError);    // Should never happen
     276          12 :     if (requireI && ! _stokesImage[I]) {
     277           0 :         *_getLog() << "This image does not have Stokes I so this calculation cannot be carried out"
     278           0 :             << LogIO::EXCEPTION;
     279             :     }
     280          12 :     if (doLin) {
     281          12 :         if (_stokesImage[I]) {
     282           6 :             if (! _checkIQUBeams(False, False)) {
     283           0 :                 *_getLog() << LogIO::WARN
     284             :                     << "I, Q, and U beams are not the same, cannot do linear portion"
     285           0 :                     << LogIO::POST;
     286           0 :                 doLin = False;
     287             :             }
     288             :         }
     289             :         else {
     290           6 :             if (! _checkQUBeams(False, False)) {
     291           0 :                 *_getLog() << LogIO::WARN
     292             :                     << "Q, and U beams are not the same, cannot do linear portion"
     293           0 :                     << LogIO::POST;
     294           0 :                 doLin = False;
     295             :             }
     296             :         }
     297             :     }
     298          12 :     if (doCirc) {
     299           6 :         if (_stokesImage[I] && ! _checkIVBeams(False, False)) {
     300           0 :             *_getLog() << LogIO::WARN
     301             :                 << "I and V beams are not the same, cannot do circular portion"
     302           0 :                 << LogIO::POST;
     303           0 :             doCirc = False;
     304             :         }
     305             :     }
     306          12 :     if (! doLin && ! doCirc) {
     307           0 :         throw AipsError("Can do neither linear nor circular portions");
     308             :     }
     309          12 : }
     310             : 
     311          15 : void ImagePolTask::_setInfo(ImageInterface<Float>& im, const StokesTypes stokes) const {
     312          30 :     ImageInfo info = _getImage()->imageInfo();
     313          15 :     if (info.hasMultipleBeams()) {
     314           0 :         info.setBeams(_stokesImage[stokes]->imageInfo().getBeamSet());
     315             :     }
     316          15 :     im.setImageInfo(info);
     317          15 : }
     318             : 
     319           0 : Float ImagePolTask::_sigma(Float clip) {
     320           0 :     *_getLog() << LogOrigin(CLASS_NAME, __func__, WHERE);
     321           0 :     Float sigma2 = 0.0;
     322           0 :     if (_stokesImage[V]) {
     323           0 :         *_getLog() << LogIO::NORMAL << "Determined noise from V image to be ";
     324           0 :         sigma2 = _sigma(V, clip);
     325             :     }
     326           0 :     else if (
     327           0 :         _stokesImage[Q] && _stokesImage[U] && _checkQUBeams(False, False)
     328             :     ) {
     329           0 :         sigma2 = sigmaLinPolInt(clip);
     330             :     }
     331           0 :     else if (_stokesImage[Q]) {
     332           0 :         *_getLog() << LogIO::NORMAL << "Determined noise from Q image to be " << LogIO::POST;
     333           0 :         sigma2 = _sigma(Q, clip);
     334             :     }
     335           0 :     else if (_stokesImage[U]) {
     336           0 :          *_getLog() << LogIO::NORMAL << "Determined noise from U image to be " << LogIO::POST;
     337           0 :          sigma2 = _sigma(U, clip);
     338             :     }
     339           0 :     else if (_stokesImage[I]) {
     340           0 :         *_getLog() << LogIO::NORMAL << "Determined noise from I image to be " << LogIO::POST;
     341           0 :         sigma2 = _sigma(I, clip);
     342             :     }
     343           0 :     *_getLog() << sigma2 << LogIO::POST;
     344           0 :     return sigma2;
     345             : }
     346             : 
     347           0 : Float ImagePolTask::_sigma (StokesTypes index, Float clip) {
     348           0 :     auto clip2 = abs(clip);
     349           0 :     if (clip2==0.0) {
     350           0 :         clip2 = 10.0;
     351             :     }
     352           0 :     if (clip2 != _oldClip) {
     353           0 :         _stokesStats[index].reset();
     354             :     }
     355           0 :     if (! _stokesStats[index]) {
     356             :         // Find sigma for all points inside +/- clip-sigma of the mean
     357             :         // More joys of LEL
     358           0 :         const auto p = _stokesImage[index];
     359           0 :         LatticeExprNode n1 (*p);
     360           0 :         LatticeExprNode n2 (n1[abs(n1-mean(n1)) < clip2*stddev(n1)]);
     361           0 :         LatticeExpr<Float> le(n2);
     362           0 :         _stokesStats[index].reset(new LatticeStatistics<Float>(le, false, false));
     363             :     }
     364           0 :     Array<Float> sigmaA;
     365           0 :     _stokesStats[index]->getConvertedStatistic(sigmaA, LatticeStatsBase::SIGMA);
     366           0 :     ThrowIf(
     367             :         sigmaA.empty(), "No good points in clipped determination of the noise "
     368             :         "for the Stokes " + _stokesName(index) + " image"
     369             :     );
     370           0 :     _oldClip = clip2;
     371           0 :     return sigmaA(IPosition(1,0));
     372             : }
     373             : 
     374           0 : String ImagePolTask::_stokesName (StokesTypes index) const {
     375           0 :     switch(index) {
     376           0 :     case I: return "I";
     377           0 :     case Q: return "Q";
     378           0 :     case U: return "U";
     379           0 :     case V: return "V";
     380           0 :     default:
     381           0 :         ThrowCc("Unsupported stokes index " + String::toString(index));
     382             :     }
     383             : }
     384             : 
     385             : }
     386             :