casa
$Rev:20696$
|
00001 //# AWProjectFT.h: Definition for AWProjectFT 00002 //# Copyright (C) 1996,1997,1998,1999,2000,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 adressed 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 //# 00027 //# $Id$ 00028 00029 #ifndef SYNTHESIS_AWPROJECTFT_H 00030 #define SYNTHESIS_AWPROJECTFT_H 00031 00032 #include <synthesis/TransformMachines/VLACalcIlluminationConvFunc.h> 00033 #include <synthesis/TransformMachines/VLAIlluminationConvFunc.h> 00034 #include <synthesis/TransformMachines/AWVisResampler.h> 00035 //#include <synthesis/MeasurementComponents/ConvolutionFunction.h> 00036 #include <synthesis/TransformMachines/EVLAConvFunc.h> 00037 #include <synthesis/MeasurementComponents/SolvableVisCal.h> 00038 #include <synthesis/TransformMachines/VPSkyJones.h> 00039 #include <synthesis/TransformMachines/FTMachine.h> 00040 //#include <synthesis/MeasurementComponents/CFCache.h> 00041 #include <synthesis/TransformMachines/Utils.h> 00042 00043 #include <scimath/Mathematics/FFTServer.h> 00044 #include <synthesis/MSVis/VisBuffer.h> 00045 00046 #include <casa/Containers/Block.h> 00047 #include <casa/Arrays/Array.h> 00048 #include <casa/Arrays/Vector.h> 00049 #include <casa/Arrays/Matrix.h> 00050 00051 #include <scimath/Mathematics/ConvolveGridder.h> 00052 #include <lattices/Lattices/LatticeCache.h> 00053 #include <lattices/Lattices/ArrayLattice.h> 00054 #include <ms/MeasurementSets/MSColumns.h> 00055 #include <measures/Measures/Measure.h> 00056 #include <measures/Measures/MDirection.h> 00057 #include <measures/Measures/MPosition.h> 00058 #include <images/Images/ImageInterface.h> 00059 #include <coordinates/Coordinates/DirectionCoordinate.h> 00060 00061 #include <casa/OS/Timer.h> 00062 00063 namespace casa { //# NAMESPACE CASA - BEGIN 00064 00065 // <summary> An FTMachine for Gridded Fourier transforms including effects of primary beam and pointing offsets and the w-term</summary> 00066 00067 // <use visibility=export> 00068 00069 // <reviewed reviewer="" date="" tests="" demos=""> 00070 00071 // <prerequisite> 00072 // <li> <linkto class=FTMachine>FTMachine</linkto> module 00073 // <li> <linkto class=SkyEquation>SkyEquation</linkto> module 00074 // <li> <linkto class=VisBuffer>VisBuffer</linkto> module 00075 // <li> <linto class=EPJones>EPJones</linkto> module 00076 // </prerequisite> 00077 // 00078 // <etymology> 00079 // FTMachine is a Machine for Fourier Transforms. Like 00080 // WProjectFT, AWProjectFT does Grid-based Fourier transforms but 00081 // also includes the effects of primary beam and antenna pointing 00082 // offsets. 00083 // </etymology> 00084 // 00085 // <synopsis> 00086 // 00087 // The <linkto class=SkyEquation>SkyEquation</linkto> needs to be 00088 // able to perform Fourier transforms on visibility 00089 // data. AWProjectFT allows efficient handling of direction 00090 // dependent effects due to the primary beam and antenna pointing 00091 // offsets using a <linkto class=VisBuffer>VisBuffer</linkto> which 00092 // encapsulates a chunk of visibility (typically all baselines for 00093 // one time) together with all the information needed for processing 00094 // (e.g. UVW coordinates). 00095 // 00096 // Using this FTMachine, errors due antenna pointing offsets can be 00097 // corrected during deconvolution. One form of antenna pointing 00098 // error which is known a-priori is the VLA polarization squint 00099 // (about 6% of the Primary beam width at any frequency). For 00100 // Stokes imaging, using this FTMachine, the VLA polarization squint 00101 // and beam polarization can also be corrected. Also since the 00102 // effects of antenna pointing errors is strongest in the range of 00103 // 1-2GHz band (where the sky is not quite empty while the beams are 00104 // not too large either), this FTMachine can also be setup to 00105 // correct for the w-term. 00106 // 00107 // Switches are provided in the get() method to compute the 00108 // derivatives with respect to the parameters of the primary beam 00109 // (only pointing offsets for now). This is used in the pointing 00110 // offset solver. 00111 // 00112 // See the documentation of other FTMachines for details about the 00113 // design of the FTMachines in general. 00114 // 00115 // </synopsis> 00116 // 00117 // <example> 00118 // See the example for <linkto class=SkyModel>SkyModel</linkto>. 00119 // </example> 00120 // 00121 // <motivation> 00122 // 00123 // Encapsulate the correction of direction dependent effects via 00124 // visibility plane convolutions with a potentially different 00125 // convolution function for each baseline. 00126 // 00127 // </motivation> 00128 // 00129 // <todo asof="2005/07/21"> 00130 // 00131 // <ul> Include the antenna based complex gain term as well since 00132 // that can interfere with the effects of pointing offsets. 00133 // 00134 // <ul> Factor out the actual convolution functions as a separate 00135 // class making FTMachines for various direction dependent effects 00136 // generic. 00137 // 00138 // </todo> 00139 00140 // class EPJones; 00141 class SolvableVisJones; 00142 class AWProjectFT : public FTMachine { 00143 public: 00144 AWProjectFT(); 00145 00146 // Constructor: cachesize is the size of the cache in words 00147 // (e.g. a few million is a good number), tilesize is the 00148 // size of the tile used in gridding (cannot be less than 00149 // 12, 16 works in most cases). 00150 // <group> 00151 AWProjectFT(Int nFacets, Long cachesize, 00152 CountedPtr<CFCache>& cfcache, 00153 CountedPtr<ConvolutionFunction>& cf, 00154 CountedPtr<VisibilityResamplerBase>& visResampler, 00155 Bool applyPointingOffset=True, 00156 Bool doPBCorr=True, 00157 Int tilesize=16, 00158 Float pbLimit=5e-4, 00159 Bool usezero=False, 00160 Bool conjBeams_p=True, 00161 Bool doublePrecGrid=False); 00162 // </group> 00163 00164 // Construct from a Record containing the AWProjectFT state 00165 AWProjectFT(const RecordInterface& stateRec); 00166 00167 // Copy constructor 00168 AWProjectFT(const AWProjectFT &other); 00169 00170 // Assignment operator 00171 AWProjectFT &operator=(const AWProjectFT &other); 00172 00173 ~AWProjectFT(); 00174 00175 // void setEPJones(EPJones* ep_j) {epJ = ep_j;} 00176 void setEPJones(SolvableVisJones* ep_j) {epJ_p = ep_j;} 00177 00178 virtual void setDOPBCorrection(Bool doit=True) {doPBCorrection=doit;}; 00179 virtual Bool getDOPBCorrection() {return doPBCorrection;}; 00180 virtual void setConjBeams(Bool useit=True) {conjBeams_p=useit;}; 00181 virtual Bool getConjBeams() {return conjBeams_p;}; 00182 00183 virtual Float getPBLimit() {return pbLimit_p;}; 00184 // Initialize transform to Visibility plane using the image 00185 // as a template. The image is loaded and Fourier transformed. 00186 00187 void setObservatoryLocation(const MPosition& mLocation) {mLocation_p=mLocation;}; 00188 00189 // Vectorized version of initializeToVis. Required since 00190 // MultiTermFTM needs vectorized version. And this is implemented 00191 // in FTMachine (the baseclass). And one relies on static_cast in 00192 // Imager::createFTMachine() (or is it in CSE?) to cast the 00193 // pointer to specific types such that this methods gets called 00194 // when the FTMachine pointer is of type AWProjectFT. 00195 virtual void initializeToVis(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec, 00196 PtrBlock<SubImage<Float> *> & modelImageVec, 00197 PtrBlock<SubImage<Float> *>& weightImageVec, 00198 PtrBlock<SubImage<Float> *>& fluxScaleVec, 00199 Block<Matrix<Float> >& weightsVec, 00200 const VisBuffer& vb); 00201 00202 virtual void initializeToVis(ImageInterface<Complex>& image, 00203 const VisBuffer& vb); 00204 // This version returns the gridded vis...should be used in conjunction 00205 // with the version of 'get' that needs the gridded visdata 00206 virtual void initializeToVis(ImageInterface<Complex>& image, 00207 const VisBuffer& vb, Array<Complex>& griddedVis, 00208 Vector<Double>& uvscale); 00209 00210 // Finalize transform to Visibility plane: flushes the image 00211 // cache and shows statistics if it is being used. 00212 virtual void finalizeToVis(); 00213 00214 // Initialize transform to Sky plane: initializes the image 00215 virtual void initializeToSky(ImageInterface<Complex>& image, Matrix<Float>& weight, 00216 const VisBuffer& vb); 00217 00218 // Finalize transform to Sky plane: flushes the image 00219 // cache and shows statistics if it is being used. DOES NOT 00220 // DO THE FINAL TRANSFORM! 00221 virtual void finalizeToSky(); 00222 00223 virtual void initVisBuffer(VisBuffer& vb, Type whichVBColumn); 00224 void initVisBuffer(VisBuffer& vb, Type whichVBColumn, Int row); 00225 00226 // Get actual coherence from grid by degridding 00227 void get(VisBuffer& vb, Int row=-1); 00228 00229 // Get the coherence from grid return it in the degrid 00230 // is used especially when scratch columns are not 00231 // present in ms. 00232 void get(VisBuffer& vb, Cube<Complex>& degrid, 00233 Array<Complex>& griddedVis, Vector<Double>& scale, 00234 Int row=-1); 00235 00236 void get(VisBuffer& vb, Cube<Float>& pointingOffsets, Int row=-1, 00237 Type whichVBColumn=FTMachine::MODEL,Int Conj=0) 00238 { 00239 get(vb,vb,vb,pointingOffsets,row,whichVBColumn,whichVBColumn,Conj,0); 00240 } 00241 00242 void get(VisBuffer& vb, VisBuffer& gradAzVB,VisBuffer& gradElVB, 00243 Cube<Float>& pointingOffsets,Int row=-1, 00244 Type whichVBColumn=FTMachine::MODEL, 00245 Type whichGradVBColumn=FTMachine::MODEL, 00246 Int Conj=0, Int doGrad=1) ; 00247 void nget(VisBuffer& vb, 00248 // These offsets should be appropriate for the VB 00249 Array<Float>& l_off, Array<Float>& m_off, 00250 Cube<Complex>& Mout, 00251 Cube<Complex>& dMout1, 00252 Cube<Complex>& dMout2, 00253 Int Conj=0, Int doGrad=1); 00254 // Get the coherence from grid return it in the degrid 00255 // is used especially when scratch columns are not 00256 // present in ms. 00257 void get(VisBuffer& vb, Cube<Complex>& degrid, 00258 Array<Complex>& griddedVis, Vector<Double>& scale, 00259 Cube<Float>& pointingOffsets,Int row=-1); 00260 00261 00262 // Put coherence to grid by gridding. 00263 void put(const VisBuffer&, 00264 TempImage<Complex>&, Vector<Double>&, int, 00265 UVWMachine*, Bool) 00266 { 00267 // throw(AipsError("AWProjectFT::put is not implemented")); 00268 } 00269 void put(const VisBuffer& vb, Int row=-1, Bool dopsf=False, 00270 FTMachine::Type type=FTMachine::OBSERVED); 00271 00272 // Make the entire image using a ROVisIter 00273 virtual void makeImage(FTMachine::Type, 00274 ROVisibilityIterator&, 00275 ImageInterface<Complex>&, 00276 Matrix<Float>&) {}; 00277 00278 // Make the entire image 00279 void makeImage(FTMachine::Type type, 00280 VisSet& vs, 00281 ImageInterface<Complex>& image, 00282 Matrix<Float>& weight); 00283 00284 // Get the final image: do the Fourier transform and 00285 // grid-correct, then optionally normalize by the summed weights 00286 virtual ImageInterface<Complex>& getImage(Matrix<Float>&, Bool normalize=True); 00287 00288 // Get the final weights image 00289 void getWeightImage(ImageInterface<Float>&, Matrix<Float>&); 00290 00291 // Save and restore the AWProjectFT to and from a record 00292 Bool toRecord(RecordInterface& outRec, Bool withImage=False); 00293 Bool fromRecord(const RecordInterface& inRec); 00294 00295 // Can this FTMachine be represented by Fourier convolutions? 00296 Bool isFourier() {return True;} 00297 00298 // Bool changed(const VisBuffer& vb) {return vpSJ->changed(vb,1);}; 00299 // Bool changed(const VisBuffer& vb) {return False;} 00300 00301 virtual Int findPointingOffsets(const VisBuffer&, Array<Float>&, Array<Float>&, 00302 Bool Evaluate=True); 00303 virtual Int findPointingOffsets(const VisBuffer&, Cube<Float>&, 00304 Array<Float>&, Array<Float>&, 00305 Bool Evaluate=True); 00306 virtual Double getVBPA(const VisBuffer& vb) 00307 { 00308 // if (!rotateApertureOTP_p) return currentCFPA; 00309 // else return getPA(vb); 00310 return getPA(vb); 00311 }; 00312 MDirection::Convert makeCoordinateMachine(const VisBuffer&, 00313 const MDirection::Types&, 00314 const MDirection::Types&, 00315 MEpoch& last); 00316 // 00317 // Make a sensitivity image (sensitivityImage), given the gridded 00318 // weights (wtImage). These are related to each other by a 00319 // Fourier transform and normalization by the sum-of-weights 00320 // (sumWt) and normalization by the product of the 2D FFT size 00321 // along each axis. If doFFTNorm=False, normalization by the FFT 00322 // size is not done. If sumWt is not provided, normalization by 00323 // the sum of weights is also not done. 00324 // 00325 virtual void makeSensitivityImage(Lattice<Complex>&,// wtImage, 00326 ImageInterface<Float>&,// sensitivityImage, 00327 const Matrix<Float>&,// sumWt=Matrix<Float>(), 00328 const Bool& ) {}; 00329 virtual void makeSensitivityImage(const VisBuffer& vb, const ImageInterface<Complex>& imageTemplate, 00330 ImageInterface<Float>& sensitivityImage); 00331 00332 // 00333 // Given the sky image (Fourier transform of the visibilities), 00334 // sum of weights and the sensitivity image, this method replaces 00335 // the skyImage with the normalized image of the sky. 00336 // 00337 // These are now implement in the FTMachine base class. 00338 // 00339 // These can't be left here since now the design depends on 00340 // getting the correct version called due to static type casting 00341 // of FTM pointer and not on inheretance and and specialization 00342 // via overloading. 00343 00344 // virtual void normalizeImage(Lattice<Complex>& skyImage, 00345 // const Matrix<Double>& sumOfWts, 00346 // Lattice<Float>& sensitivityImage, 00347 // Bool fftNorm=True); 00348 // virtual void normalizeImage(Lattice<Complex>& skyImage, 00349 // const Matrix<Double>& sumOfWts, 00350 // Lattice<Float>& sensitivityImage, 00351 // Lattice<Complex>& sensitivitySqImage, 00352 // Bool fftNorm=True); 00353 00354 virtual ImageInterface<Float>& getSensitivityImage() {return *avgPB_p;} 00355 virtual Matrix<Double>& getSumOfWeights() {return sumWeight;}; 00356 virtual Matrix<Double>& getSumOfCFWeights() {return sumCFWeight;}; 00357 00358 void makeConjPolMap(const VisBuffer& vb, const Vector<Int> cfPolMap, Vector<Int>& conjPolMap); 00359 // Vector<Int> makeConjPolMap(const VisBuffer& vb); 00360 void makeCFPolMap(const VisBuffer& vb, const Vector<Int>& cfstokes, Vector<Int>& polM); 00361 // void reset() {vpSJ->reset();} 00362 void reset() {paChangeDetector.reset();} 00363 00364 void setPAIncrement(const Quantity &computePAIncr, const Quantity &rotateOTFPAIncr); 00365 00366 Vector<Int>& getPolMap() {return polMap;}; 00367 virtual String name() const { return "AWProjectFT";}; 00368 virtual Bool verifyAvgPB(ImageInterface<Float>& pb, ImageInterface<Float>& sky) 00369 {return verifyShapes(pb.shape(),sky.shape());} 00370 00371 virtual Bool verifyAvgPB(ImageInterface<Float>& pb, ImageInterface<Complex>& sky) 00372 {return verifyShapes(pb.shape(),sky.shape());} 00373 00374 virtual Bool verifyShapes(IPosition shape0, IPosition shape1); 00375 00376 inline virtual Float pbFunc(const Float& a, const Float& limit) 00377 // {Float tt=sqrt(a);return (abs(tt) >= limit)?tt:1.0;}; 00378 {if (abs(a) >= limit) return (a);else return 1.0;}; 00379 inline virtual Complex pbFunc(const Complex& a, const Float& limit) 00380 {if (abs(a)>=limit) return (a); else return Complex(1.0,0.0);}; 00381 00382 virtual void setMiscInfo(const Int qualifier) 00383 { 00384 sensitivityPatternQualifier_p=qualifier; 00385 sensitivityPatternQualifierStr_p = "_tt"+String::toString(sensitivityPatternQualifier_p); 00386 } 00387 virtual void ComputeResiduals(VisBuffer&vb, Bool useCorrected); 00388 void makeWBCFWt(CFStore2& cfs,const Double imRefFreq); 00389 00390 protected: 00391 00392 Int nint(Double val) {return Int(floor(val+0.5));}; 00393 // Locate convolution functions on the disk 00394 // Int locateConvFunction(const Int Nw, const Float pa); 00395 // void cacheConvFunction(Int which, Array<Complex>& cf, CoordinateSystem& coord); 00396 // Find the convolution function 00397 void findConvFunction(const ImageInterface<Complex>& image, 00398 const VisBuffer& vb); 00399 00400 // Get the appropriate data pointer 00401 Array<Complex>* getDataPointer(const IPosition&, Bool); 00402 00403 void ok(); 00404 00405 void init(); 00406 // virtual void initPolInfo(const VisBuffer& vb); 00407 // Is this record on Grid? check both ends. This assumes that the 00408 // ends bracket the middle 00409 Bool recordOnGrid(const VisBuffer& vb, Int rownr) const; 00410 00411 // Padding in FFT 00412 Float padding_p; 00413 00414 Int nWPlanes_p; 00415 // Image cache 00416 LatticeCache<Complex> * imageCache; 00417 00418 // Sizes 00419 Long cachesize; 00420 Int tilesize; 00421 00422 // Gridder 00423 ConvolveGridder<Double, Complex>* gridder; 00424 00425 // Is this tiled? 00426 Bool isTiled; 00427 00428 // Array lattice 00429 CountedPtr<Lattice<Complex> > arrayLattice; 00430 00431 // Lattice. For non-tiled gridding, this will point to arrayLattice, 00432 // whereas for tiled gridding, this points to the image 00433 CountedPtr<Lattice<Complex> > lattice; 00434 00435 Float maxAbsData; 00436 00437 // Useful IPositions 00438 IPosition centerLoc, offsetLoc; 00439 00440 // Image Scaling and offset 00441 Vector<Double> uvScale, uvOffset; 00442 00443 // Array for non-tiled gridding 00444 Array<Complex> griddedData; 00445 Array<DComplex> griddedData2; 00446 00447 // DirectionCoordinate directionCoord; 00448 MDirection::Convert* pointingToImage; 00449 00450 // Grid/degrid zero spacing points? 00451 Bool usezero_p; 00452 00453 // CountedPtr<ConvolutionFunction> telescopeConvFunc_p; 00454 // CFStore cfs_p, cfwts_p; 00455 // Array<Complex> convFunc_p, convWeights_p; 00456 // 00457 // Vector to hold the support size info. for the convolution 00458 // functions pointed to by the elements of convFunctions_p. The 00459 // co-ordinates of this array are (W-term, Poln, PA). 00460 // 00461 Int convSize, convSampling, wConvSize, lastIndex_p; 00462 00463 // 00464 // The average PB for sky image normalization 00465 // 00466 CountedPtr<ImageInterface<Float> > avgPB_p; 00467 CountedPtr<ImageInterface<Complex> > avgPBSq_p; 00468 // 00469 // No. of vis. polarization planes used in making the user defined 00470 // Stokes images 00471 // 00472 Int maxConvSupport; 00473 // 00474 // Percentage of the peak of the PB after which the image is set 00475 // to zero. 00476 // 00477 CountedPtr<SolvableVisJones> epJ_p; 00478 Double sigma; 00479 Int Nant_p, doPointing; 00480 Bool doPBCorrection, makingPSF, conjBeams_p; 00481 00482 // CountedPtr<CFCache> cfCache_p; 00483 ParAngleChangeDetector paChangeDetector; 00484 Double rotateOTFPAIncr_p, computePAIncr_p; 00485 00486 Unit Second, Radian, Day; 00487 Array<Float> l_offsets,m_offsets; 00488 Vector<Float> pbPeaks, paList; 00489 00490 Double currentCFPA, cfRefFreq_p, imRefFreq_p; 00491 Float lastPAUsedForWtImg; 00492 Bool pbNormalized_p; 00493 Vector<Bool> paNdxProcessed_p; 00494 // 00495 //---------------------------------------------------------------------- 00496 // 00497 virtual void normalizeAvgPB(); 00498 virtual void normalizeAvgPB(ImageInterface<Complex>& inImage, 00499 ImageInterface<Float>& outImage); 00500 virtual void resampleDataToGrid(Array<Complex>& griddedData, VBStore& vbs, 00501 const VisBuffer& vb, Bool& dopsf); 00502 virtual void resampleDataToGrid(Array<DComplex>& griddedData, VBStore& vbs, 00503 const VisBuffer& vb, Bool& dopsf); 00504 virtual void resampleGridToData(VBStore& vbs, Array<Complex>& griddedData, 00505 const VisBuffer& vb); 00506 00507 virtual void setupVBStore(VBStore& vbs, 00508 const VisBuffer& vb, 00509 const Matrix<Float>& imagingweight, 00510 const Cube<Complex>& visData, 00511 const Matrix<Double>& uvw, 00512 const Cube<Int>& flagCube, 00513 const Vector<Double>& dphase); 00514 00515 // AWVisResampler visResampler_p; 00516 CountedPtr<VisibilityResamplerBase> visResampler_p; 00517 Int sensitivityPatternQualifier_p; 00518 String sensitivityPatternQualifierStr_p; 00519 CFStore rotatedConvFunc_p; 00520 CountedPtr<CFStore2> cfs2_p, cfwts2_p; 00521 Vector<Int> ConjCFMap_p, CFMap_p; 00522 00523 Timer timer_p; 00524 Double runTime1_p; 00525 00526 #include "AWProjectFT.FORTRANSTUFF.INC" 00527 }; 00528 } //# NAMESPACE CASA - END 00529 00530 #endif