Line data    Source code

       1             : //# PSTerm.cc: implementation of PSTerm
       2             : //# Copyright (C) 2007
       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 Library 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 Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library 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 adressed as follows:
      20             : //#        Internet email: casa-feedback@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             : //#
      27             : //# $Id$
      28             : #include <synthesis/TransformMachines2/PSTerm.h>
      29             : #include <synthesis/TransformMachines2/Utils.h>
      30             : #include <synthesis/TransformMachines/SynthesisError.h>
      31             : #ifdef _OPENMP
      32             : #include <omp.h>
      33             : #endif
      34             : 
      35             : using namespace casacore;
      36             : namespace casa { //# NAMESPACE CASA - BEGIN
      37             :   using namespace refim;
      38             :   //
      39             :   //-------------------------------------------------------------------------
      40             :   //           
      41           0 :   void PSTerm::init(const IPosition shape, 
      42             :                     const Vector<Double>& uvScale,
      43             :                     const Vector<Double>& uvOffset,
      44             :                     const Double& psScale)
      45             :   {
      46           0 :     shape_p = shape;
      47           0 :     uvScale_p=uvScale;
      48           0 :     uvOffset_p = uvOffset;
      49           0 :     psScale_p = psScale;
      50             :     
      51           0 :     psCtor_p   = new ConvolveGridder<Double, Complex>(shape_p, uvScale_p, uvOffset_p, "SF");
      52           0 :     sampling_p = (*psCtor_p).cSupport();
      53           0 :     support_p  = (*psCtor_p).cSampling();
      54           0 :   }
      55           0 :   void PSTerm::reinit(const IPosition shape, 
      56             :                       const Vector<Double>& uvScale,
      57             :                       const Vector<Double>& uvOffset,
      58             :                       const Double& psScale)
      59             :   {
      60           0 :     if (shape.nelements() > 0) shape_p = shape;
      61           0 :     if (uvScale.nelements() > 0) uvScale_p=uvScale;
      62           0 :     if (uvOffset.nelements() > 0) uvOffset_p = uvOffset;
      63           0 :     if (psScale > 0) psScale_p = psScale;
      64             :     
      65           0 :     psCtor_p   = new ConvolveGridder<Double, Complex>(shape_p, uvScale_p, uvOffset_p, "SF");
      66           0 :     sampling_p = (*psCtor_p).cSupport();
      67           0 :     support_p  = (*psCtor_p).cSampling();
      68           0 :   }
      69             :   //
      70             :   //-------------------------------------------------------------------------
      71             :   //           
      72           0 :   Matrix<Complex>& PSTerm::operator=(Matrix<Complex>& buf)
      73             :   {
      74           0 :     applySky(buf);
      75           0 :     return buf;
      76             :   }
      77             :   //
      78             :   //-------------------------------------------------------------------------
      79             :   //           
      80           0 :   Matrix<Complex>& PSTerm::operator*=(Matrix<Complex>& buf)
      81             :   {
      82           0 :     applySky(buf,true);
      83           0 :     return buf;
      84             :   }
      85             :   //
      86             :   //-------------------------------------------------------------------------
      87             :   //           
      88           0 :   void PSTerm::applySky(Matrix<Complex>& screen, Bool multiply)
      89             :   {
      90             :     //    AlwaysAssert(psCtor_p.null()==false, SynthesisFTMachineError);
      91           0 :     Int nx=screen.shape()(0), ny=screen.shape()(1);
      92           0 :     Int convOrig=nx/2;
      93           0 :     Float xpart, psScale_local=psScale_p;
      94             : #ifdef _OPENMP
      95           0 :     Int Nth=max(omp_get_max_threads()-2,1);
      96             : #endif
      97             : 
      98           0 :     if (!isNoOp())
      99             :       {
     100           0 :         for (Int i=0; i<nx;i++)
     101             :           {
     102           0 :             xpart = square(i-convOrig);
     103             : #ifdef _OPENMP
     104             : //#pragma omp parallel default(none) firstprivate(xpart,convOrig, i) shared(screen,psScale_local,ny,multiply) num_threads(Nth)
     105           0 : #pragma omp parallel firstprivate(xpart,convOrig, i) shared(psScale_local,ny,multiply) num_threads(Nth)
     106             : #endif
     107             :    {
     108             : #ifdef _OPENMP
     109             : #pragma omp for
     110             : #endif
     111             :             for (Int j=0;j<ny;j++)
     112             :               {
     113             :                 Float ypart;
     114             :                 ypart = sqrt(xpart + square(j-convOrig))*psScale_local;
     115             :                 Float val = SynthesisUtils::libreSpheroidal(ypart)*(1.0-ypart*ypart);;
     116             :                 if (val > 0)
     117             :                   {
     118             :                     if (multiply)  screen(i, j) *= val;
     119             :                     else           screen(i, j)  = val;
     120             :                   }
     121             :               }
     122             :     }
     123             :           }
     124             :       }
     125           0 :   }
     126             : 
     127           0 :   void PSTerm::applySky(Matrix<Complex>& screen, 
     128             :                         const Vector<Double>& sampling,
     129             :                         const Int inner)
     130             :   {
     131             :     //GG: convsize: 2048 inner: 512 uvscale: [-0.0397159, 0.0397159, 2.33547e-310] uvoffset:[1024, 1024, 0]
     132             :     // Vector<Double> uvScale(3);uvScale[0]=-0.0397159;uvScale[1]=0.0397159;uvScale[2]=0.0;
     133             :     // Vector<Double> uvOffset(3,0.0); uvOffset[0]=uvOffset[1]=1024;
     134             :     // psCtor_p=new ConvolveGridder<Double, Complex>(IPosition(2, inner, inner),
     135             :     //                                            uvScale, uvOffset,
     136             :     //                                            "SF");
     137             : 
     138             : 
     139           0 :     AlwaysAssert(psCtor_p.null()==false, SynthesisFTMachineError);
     140             : 
     141           0 :     Int convSize = screen.shape()[0];
     142             :     //    Vector<Complex> correction(inner);
     143           0 :     Vector<Complex> correction;
     144           0 :     if (!isNoOp())
     145             :       {
     146           0 :         for (Int iy=-inner/2;iy<inner/2;iy++) 
     147             :           {
     148           0 :             psCtor_p->correctX1D(correction, iy+inner/2);
     149           0 :             Int k=0;
     150           0 :             for (Int ix=-inner/2;ix<inner/2;ix++) 
     151             :               {
     152           0 :                 screen(ix+convSize/2,iy+convSize/2)=correction(k++);//correction(ix+inner/2);
     153             :               }
     154             :           }
     155             :       }
     156             :     (void)sampling; 
     157           0 :   }
     158             :   //
     159             :   //-------------------------------------------------------------------------
     160             :   //           
     161           0 :   void PSTerm::normalizeImage(Lattice<Complex>& skyImage,
     162             :                               const Matrix<Float>& weights)
     163             :   {
     164           0 :     Int inx = skyImage.shape()(0);
     165             :     //    Int iny = skyImage.shape()(1);
     166           0 :     Vector<Complex> correction(inx);
     167           0 :     correction=Complex(1.0, 0.0);
     168           0 :     IPosition cursorShape(4, inx, 1, 1, 1);
     169           0 :     IPosition axisPath(4, 0, 1, 2, 3);
     170           0 :     LatticeStepper lsx(skyImage.shape(), cursorShape, axisPath);
     171           0 :     LatticeIterator<Complex> lix(skyImage, lsx);
     172             : 
     173           0 :     for(lix.reset();!lix.atEnd();lix++) 
     174             :       {
     175           0 :         Int pol=lix.position()(2), chan=lix.position()(3);
     176             : 
     177           0 :         if(weights(pol, chan)!=0.0) 
     178             :           {
     179           0 :             Int iy=lix.position()(1);
     180           0 :             psCtor_p->correctX1D(correction,iy);
     181             :             // for (Int ix=0;ix<inx;ix++) 
     182             :             //   correction(ix)*=sincConv(ix)*sincConv(iy);
     183             : 
     184           0 :             lix.rwVectorCursor()/=correction;
     185             :           }
     186             :         else 
     187           0 :           lix.woCursor()=0.0;
     188             :       }
     189           0 :   }
     190             : };