Line data    Source code

       1             : // -*- C++ -*-
       2             : //# AWProjectWBFT.cc: Implementation of AWProjectWBFT class
       3             : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
       4             : //# Associated Universities, Inc. Washington DC, USA.
       5             : //#
       6             : //# This library is free software; you can redistribute it and/or modify it
       7             : //# under the terms of the GNU Library General Public License as published by
       8             : //# the Free Software Foundation; either version 2 of the License, or (at your
       9             : //# option) any later version.
      10             : //#
      11             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      12             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      13             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      14             : //# License for more details.
      15             : //#
      16             : //# You should have received a copy of the GNU Library General Public License
      17             : //# along with this library; if not, write to the Free Software Foundation,
      18             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      19             : //#
      20             : //# Correspondence concerning AIPS++ should be addressed as follows:
      21             : //#        Internet email: aips2-request@nrao.edu.
      22             : //#        Postal address: AIPS++ Project Office
      23             : //#                        National Radio Astronomy Observatory
      24             : //#                        520 Edgemont Road
      25             : //#                        Charlottesville, VA 22903-2475 USA
      26             : //#
      27             : //# $Id$
      28             : 
      29             : #include <synthesis/TransformMachines/AWProjectWBFTNew.h>
      30             : 
      31             : #include <synthesis/TransformMachines/AWVisResampler.h>
      32             : #include <synthesis/TransformMachines/StokesImageUtil.h>
      33             : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
      34             : #include <casacore/scimath/Mathematics/FFTServer.h>
      35             : #include <casacore/scimath/Mathematics/Convolver.h>
      36             : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
      37             : #include <casacore/images/Images/ImageInterface.h>
      38             : #include <casacore/images/Images/PagedImage.h>
      39             : #include <msvis/MSVis/VisBuffer.h>
      40             : #include <casacore/casa/Arrays/Vector.h>
      41             : #include <casacore/casa/Arrays/Slice.h>
      42             : #include <casacore/casa/Arrays/ArrayMath.h>
      43             : #include <casacore/casa/Arrays/Array.h>
      44             : #include <casacore/casa/OS/HostInfo.h>
      45             : #include <sstream>
      46             : #include <casacore/casa/Utilities/CompositeNumber.h>
      47             : 
      48             : using namespace casacore;
      49             : namespace casa { //# NAMESPACE CASA - BEGIN
      50             :   //---------------------------------------------------------------
      51             :   // 
      52             : 
      53             : //-----------------------------------------------------------------------
      54           0 :   FTMachine* AWProjectWBFTNew::cloneFTM(){
      55           0 :     return new AWProjectWBFTNew(*this);
      56             :   }
      57             : 
      58             : 
      59           0 :   void AWProjectWBFTNew::ftWeightImage(Lattice<Complex>& wtImage, 
      60             :                                     const Matrix<Float>& sumWt,
      61             :                                     const Bool& doFFTNorm)
      62             :   {
      63           0 :     LogIO log_l(LogOrigin("AWProjectWBFTNew", "ftWeightImage[R&D]"));
      64           0 :     if (wtImageFTDone_p) return;
      65             : 
      66             :     // Bool doSumWtNorm=true;
      67             :     // if (sumWt.shape().nelements()==0) doSumWtNorm=false;
      68             : 
      69           0 :     if ((sumWt.shape().nelements() < 2) || 
      70           0 :         (sumWt.shape()(0) != wtImage.shape()(2)) || 
      71           0 :         (sumWt.shape()(1) != wtImage.shape()(3)))
      72           0 :       log_l << "Sum of weights per poln and chan is required" << LogIO::EXCEPTION;
      73             :     //
      74             :     // Use only the amplitude of the gridded weights.  LatticeExpr
      75             :     // classes while useful, appear to be too strict about types (the
      76             :     // code below will not compile if the abs(wtImage) is not
      77             :     // converted back to a complex type).
      78             : 
      79             :     // LatticeExpr<Complex> le(abs(wtImage)*Complex(1,0));
      80             :     // wtImage.copyData(le);
      81             : 
      82             :     // {
      83             :     //   String name("wtimg.im");
      84             :     //   storeArrayAsImage(name,image->coordinates(),wtImage.get());
      85             :     // }
      86           0 :     LatticeFFT::cfft2d(wtImage,false);
      87             :     // {
      88             :     //   String name("ftwtimg.im");
      89             :     //   storeArrayAsImage(name,image->coordinates(),wtImage.get());
      90             :     // }
      91           0 :     wtImageFTDone_p=true;
      92             : 
      93           0 :     Int sizeX=wtImage.shape()(0), sizeY=wtImage.shape()(1);
      94             : 
      95           0 :     Array<Complex> wtBuf; wtImage.get(wtBuf,false);
      96           0 :     ArrayLattice<Complex> wtLat(wtBuf,true);
      97             :     //
      98             :     // Copy one 2D plane at a time, normalizing by the sum of weights
      99             :     // and possibly 2D FFT.
     100             :     //
     101             :     // Set up Lattice iterators on wtImage and sensitivityImage
     102             :     //
     103           0 :     IPosition axisPath(4, 0, 1, 2, 3);
     104           0 :     IPosition cursorShape(4, sizeX, sizeY, 1, 1);
     105             : 
     106           0 :     LatticeStepper wtImStepper(wtImage.shape(), cursorShape, axisPath);
     107           0 :     LatticeIterator<Complex> wtImIter(wtImage, wtImStepper);
     108             :     //
     109             :     // Iterate over channel and polarization axis
     110             :     //
     111           0 :     if (!doFFTNorm) sizeX=sizeY=1;
     112             :     //
     113             :     // Normalize each frequency and polarization plane of the complex
     114             :     // sensitivity pattern
     115             :     //
     116             :     // getSumOfCFWeights() returns the sum-of-weight appropriate for
     117             :     // computing sensitivity patterns (which could be different from
     118             :     // the data-sum-of-weights.
     119             :     //
     120             : 
     121             :     // USEFUL DEBUG MESSAGE
     122             :     //cerr << "SumCFWt: " << getSumOfCFWeights() << " " << max(wtBuf) << " " << sensitivityPatternQualifier_p << endl;
     123           0 :     for(wtImIter.reset(); !wtImIter.atEnd(); wtImIter++)
     124             :       {
     125             :     // These variables no longe appear to be used
     126             :         //Int pol_l=wtImIter.position()(2), chan_l=wtImIter.position()(3);
     127             :         // Lets write some mildly obfuscated code ~[8-)
     128             :         //if ((sensitivityPatternQualifier_p == -1) && (doSumWtNorm))
     129             :         //  sumwt_l = ((sumwt_l = getSumOfCFWeights()(pol_l,chan_l))==0)?1.0:sumwt_l;
     130             : 
     131             :         //sumwt_l = getSumOfCFWeights()(pol_l,chan_l);
     132             : 
     133           0 :         wtImIter.rwCursor() = (wtImIter.rwCursor()
     134           0 :                                *Float(sizeX)*Float(sizeY)
     135             :                                //U                             /sumwt_l
     136           0 :                                );
     137             : 
     138             :         ////////////////////    wtImIter.rwCursor() = sqrt( fabs(wtImIter.rwCursor()) );
     139             : 
     140             :         //Double maxval = fabs( max( wtImIter.rwCursor() ) );
     141             :         //cout << "sumwt from WBAWPNew::ftWeightImage : " << sumwt_l << "  max val in wtimg : " << maxval << endl;
     142             : 
     143             : 
     144             : 
     145             :         //      sumwt_l = getSumOfCFWeights()(pol_l,chan_l);
     146             :         //      weightRatio_p = maxval * Float(sizeX)*Float(sizeY) / sumwt_l;
     147             :       }
     148             : 
     149             :   }
     150             : 
     151             :  
     152             : } //# NAMESPACE CASA - END