Line data    Source code

       1             : //# Mosaicft.cc: Implementation of MosaicFTNew class
       2             : //# Copyright (C) 2016
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: aips2-request@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //# $Id$
      27             : 
      28             : #include <synthesis/TransformMachines/MosaicFTNew.h>
      29             : 
      30             : 
      31             : #include <msvis/MSVis/VisibilityIterator.h>
      32             : #include <casacore/casa/Quanta/UnitMap.h>
      33             : #include <casacore/casa/Quanta/MVTime.h>
      34             : #include <casacore/casa/Quanta/UnitVal.h>
      35             : #include <casacore/measures/Measures/Stokes.h>
      36             : #include <casacore/measures/Measures/UVWMachine.h>
      37             : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
      38             : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
      39             : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
      40             : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
      41             : #include <casacore/coordinates/Coordinates/Projection.h>
      42             : #include <casacore/ms/MeasurementSets/MSColumns.h>
      43             : #include <casacore/casa/BasicSL/Constants.h>
      44             : #include <casacore/scimath/Mathematics/FFTServer.h>
      45             : #include <synthesis/TransformMachines/SimplePBConvFunc.h>
      46             : #include <synthesis/TransformMachines/HetArrayConvFunc.h>
      47             : #include <synthesis/TransformMachines/PBMath.h>
      48             : #include <synthesis/TransformMachines/VPSkyJones.h>
      49             : #include <casacore/scimath/Mathematics/RigidVector.h>
      50             : #include <msvis/MSVis/StokesVector.h>
      51             : #include <synthesis/TransformMachines/StokesImageUtil.h>
      52             : #include <msvis/MSVis/VisBuffer.h>
      53             : #include <msvis/MSVis/VisSet.h>
      54             : #include <casacore/images/Images/ImageInterface.h>
      55             : #include <casacore/images/Images/PagedImage.h>
      56             : #include <casacore/images/Images/SubImage.h>
      57             : #include <casacore/images/Regions/ImageRegion.h>
      58             : #include <casacore/images/Regions/WCBox.h>
      59             : #include <casacore/casa/Containers/Block.h>
      60             : #include <casacore/casa/Containers/Record.h>
      61             : #include <casacore/casa/Arrays/ArrayLogical.h>
      62             : #include <casacore/casa/Arrays/ArrayMath.h>
      63             : #include <casacore/casa/Arrays/Array.h>
      64             : #include <casacore/casa/Arrays/MaskedArray.h>
      65             : #include <casacore/casa/Arrays/Vector.h>
      66             : #include <casacore/casa/Arrays/Slice.h>
      67             : #include <casacore/casa/Arrays/Matrix.h>
      68             : #include <casacore/casa/Arrays/Cube.h>
      69             : #include <casacore/casa/Arrays/MatrixIter.h>
      70             : #include <casacore/casa/BasicSL/String.h>
      71             : #include <casacore/casa/Utilities/Assert.h>
      72             : #include <casacore/casa/Exceptions/Error.h>
      73             : #include <casacore/lattices/Lattices/ArrayLattice.h>
      74             : #include <casacore/lattices/Lattices/SubLattice.h>
      75             : #include <casacore/lattices/LRegions/LCBox.h>
      76             : #include <casacore/lattices/LEL/LatticeExpr.h>
      77             : #include <casacore/lattices/Lattices/LatticeCache.h>
      78             : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
      79             : #include <casacore/lattices/Lattices/LatticeIterator.h>
      80             : #include <casacore/lattices/Lattices/LatticeStepper.h>
      81             : #include <casacore/casa/Utilities/CompositeNumber.h>
      82             : #include <casacore/casa/OS/Timer.h>
      83             : #include <casacore/casa/OS/HostInfo.h>
      84             : #include <sstream>
      85             : #ifdef _OPENMP
      86             : #include <omp.h>
      87             : #endif
      88             : using namespace casacore;
      89             : namespace casa { //# NAMESPACE CASA - BEGIN
      90             : 
      91           0 :   FTMachine* MosaicFTNew::cloneFTM(){
      92           0 :     Record rec;
      93           0 :     String err;
      94           0 :     if(!(this->toRecord(err, rec)))
      95           0 :        throw(AipsError("Error in copying FTMachine"));
      96           0 :     MosaicFTNew *clone=new MosaicFTNew(rec);
      97           0 :     clone->reset();
      98           0 :     return clone;
      99             :   }
     100             : 
     101             : // Finalize the FFT to the Sky. Here we actually do the FFT and
     102             : // return the resulting image
     103           0 : ImageInterface<Complex>& MosaicFTNew::getImage(Matrix<Float>& weights,
     104             :                                             Bool normalize) 
     105             : {
     106             :   //AlwaysAssert(lattice, AipsError);
     107           0 :   AlwaysAssert(image, AipsError);
     108             : 
     109           0 :    if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
     110           0 :     throw(AipsError("Programmer error ...request for image without right order of calls"));
     111             :   }
     112             : 
     113           0 :   logIO() << LogOrigin("MosaicFTNew", "getImage") << LogIO::NORMAL;
     114             :   
     115           0 :   weights.resize(sumWeight.shape());
     116             :   
     117           0 :   convertArray(weights, sumWeight);
     118             :   
     119             :   // If the weights are all zero then we cannot normalize
     120             :   // otherwise we don't care.
     121           0 :   if(max(weights)==0.0) {
     122           0 :     if(normalize) {
     123             :       logIO() << LogIO::SEVERE << "No useful data in MosaicFTNew: weights all zero"
     124           0 :               << LogIO::POST;
     125             :     }
     126             :     else {
     127             :       logIO() << LogIO::WARN << "No useful data in MosaicFTNew: weights all zero"
     128           0 :               << LogIO::POST;
     129             :     }
     130             :   }
     131             :   else {
     132             :     
     133             :     //const IPosition latticeShape = lattice->shape();
     134             :     
     135             :     logIO() << LogIO::DEBUGGING
     136           0 :             << "Starting FFT and scaling of image" << LogIO::POST;
     137           0 :     if(useDoubleGrid_p){
     138           0 :       ArrayLattice<DComplex> darrayLattice(griddedData2);
     139           0 :       LatticeFFT::cfft2d(darrayLattice,false);
     140           0 :       griddedData.resize(griddedData2.shape());
     141           0 :       convertArray(griddedData, griddedData2);
     142             :       
     143             :       //Don't need the double-prec grid anymore...
     144           0 :       griddedData2.resize();
     145           0 :       arrayLattice = new ArrayLattice<Complex>(griddedData);
     146           0 :       lattice=arrayLattice;
     147             : 
     148             :     }
     149             :     else{
     150           0 :       arrayLattice = new ArrayLattice<Complex>(griddedData);
     151           0 :       lattice=arrayLattice;
     152           0 :       LatticeFFT::cfft2d(*lattice,false);
     153             :     }
     154             :     
     155             :     {
     156           0 :       Int inx = lattice->shape()(0);
     157           0 :       Int iny = lattice->shape()(1);
     158           0 :       Vector<Complex> correction(inx);
     159           0 :       correction=Complex(1.0, 0.0);
     160           0 :       Vector<Float> sincConvX(nx);
     161           0 :       for (Int ix=0;ix<nx;ix++) {
     162           0 :         Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
     163           0 :         if(ix==nx/2) {
     164           0 :           sincConvX(ix)=1.0;
     165             :         }
     166             :         else {
     167           0 :           sincConvX(ix)=sin(x)/x;
     168             :         }
     169             :       }
     170           0 :       Vector<Float> sincConvY(ny);
     171           0 :       for (Int ix=0;ix<ny;ix++) {
     172           0 :         Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
     173           0 :         if(ix==ny/2) {
     174           0 :           sincConvY(ix)=1.0;
     175             :         }
     176             :         else {
     177           0 :           sincConvY(ix)=sin(x)/x;
     178             :         }
     179             :       }
     180             : 
     181             : 
     182           0 :       IPosition cursorShape(4, inx, 1, 1, 1);
     183           0 :       IPosition axisPath(4, 0, 1, 2, 3);
     184           0 :       LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
     185           0 :       LatticeIterator<Complex> lix(*lattice, lsx);
     186           0 :       for(lix.reset();!lix.atEnd();lix++) 
     187             :         {
     188           0 :         Int pol=lix.position()(2);
     189           0 :         Int chan=lix.position()(3);
     190           0 :         Int iy=lix.position()(1);
     191           0 :         for (Int ix=0;ix<nx;ix++) {
     192           0 :             correction(ix)=1.0/(sincConvX(ix)*sincConvY(iy));
     193             :         }
     194           0 :         if(normalize) 
     195             :           {
     196           0 :             if(weights(pol, chan)!=0.0) 
     197             :               {
     198           0 :                 Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
     199           0 :                 lix.rwVectorCursor()*=rnorm*correction;
     200             :               }
     201             :             else {
     202           0 :               lix.woCursor()=0.0;
     203             :             }
     204             :         }         else {
     205           0 :           Complex rnorm(Float(inx)*Float(iny));
     206           0 :           lix.rwVectorCursor()*=rnorm*correction;
     207             :         }
     208             :       }
     209             :     }
     210             : 
     211             :     //if(!isTiled) 
     212             :     {
     213             :       // Check the section from the image BEFORE converting to a lattice 
     214           0 :       IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
     215           0 :                     (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
     216           0 :       IPosition stride(4, 1);
     217           0 :       IPosition trc(blc+image->shape()-stride);
     218             :       
     219             :       // Do the copy
     220           0 :       IPosition start(4, 0);
     221           0 :       image->put(griddedData(blc, trc));
     222             :     }
     223             :   }
     224             :   
     225           0 :   if(!arrayLattice.null()) arrayLattice=0;
     226           0 :   if(!lattice.null()) lattice=0;
     227           0 :   griddedData.resize();
     228           0 :   image->clearCache();
     229           0 :   return *image;
     230             : }
     231             : 
     232             : // Get weight image
     233           0 : void MosaicFTNew::getWeightImage(ImageInterface<Float>& weightImage,
     234             :                               Matrix<Float>& weights) 
     235             : {
     236             :   
     237           0 :   logIO() << LogOrigin("MosaicFTNew", "getWeightImage") << LogIO::NORMAL;
     238             :   
     239           0 :   weights.resize(sumWeight.shape());
     240           0 :   convertArray(weights,sumWeight);
     241             : 
     242           0 :   Float inx = skyCoverage_p->shape()(0);
     243           0 :   Float iny = skyCoverage_p->shape()(1);
     244             : 
     245           0 :   weightImage.copyData( (LatticeExpr<Float>) ( (*skyCoverage_p)*inx*iny ) );
     246             : 
     247           0 :   skyCoverage_p->tempClose();
     248             : 
     249           0 : }
     250             : 
     251             : } //# NAMESPACE CASA - END
     252             : 
     253             :