Line data Source code
1 : //# MosaicFT.h: Definition for MosaicFT
2 : //# Copyright (C) 2003-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 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 :
29 : #ifndef SYNTHESIS_TRANSFORM2_MOSAICFT_H
30 : #define SYNTHESIS_TRANSFORM2_MOSAICFT_H
31 :
32 : #include <synthesis/TransformMachines2/FTMachine.h>
33 : #include <synthesis/TransformMachines2/SkyJones.h>
34 : #include <casacore/casa/Arrays/Matrix.h>
35 : #include <casacore/scimath/Mathematics/FFTServer.h>
36 : #include <msvis/MSVis/VisBuffer2.h>
37 : #include <casacore/images/Images/ImageInterface.h>
38 : #include <casacore/images/Images/ImageInterface.h>
39 : #include <casacore/casa/Containers/Block.h>
40 : #include <casacore/casa/Arrays/Array.h>
41 : #include <casacore/casa/Arrays/Vector.h>
42 : #include <casacore/casa/Utilities/CountedPtr.h>
43 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
44 : #include <casacore/lattices/Lattices/LatticeCache.h>
45 : #include <casacore/lattices/Lattices/ArrayLattice.h>
46 : #include <casacore/ms/MeasurementSets/MSColumns.h>
47 : #include <casacore/measures/Measures/Measure.h>
48 : #include <casacore/measures/Measures/MDirection.h>
49 : #include <casacore/measures/Measures/MPosition.h>
50 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
51 :
52 : namespace casacore{
53 :
54 : class MPosition;
55 : class UVWMachine;
56 : }
57 :
58 : namespace casa { //# NAMESPACE CASA - BEGIN
59 :
60 : // <summary> An FTMachine for Gridded Fourier transforms </summary>
61 :
62 : // <use visibility=export>
63 :
64 : // <reviewed reviewer="" date="" tests="" demos="">
65 :
66 : // <prerequisite>
67 : // <li> <linkto class=FTMachine>FTMachine</linkto> module
68 : // <li> <linkto class=SkyEquation>SkyEquation</linkto> module
69 : // <li> <linkto class=VisBuffer>VisBuffer</linkto> module
70 : // </prerequisite>
71 : //
72 : // <etymology>
73 : // FTMachine is a Machine for Fourier Transforms. MosaicFT does
74 : // Grid-based Fourier transforms.
75 : // </etymology>
76 : //
77 : // <synopsis>
78 : // The <linkto class=SkyEquation>SkyEquation</linkto> needs to be able
79 : // to perform Fourier transforms on visibility data. MosaicFT
80 : // allows efficient Fourier Transform processing using a
81 : // <linkto class=VisBuffer>VisBuffer</linkto> which encapsulates
82 : // a chunk of visibility (typically all baselines for one time)
83 : // together with all the information needed for processing
84 : // (e.g. UVW coordinates).
85 : //
86 : // Gridding and degridding in MosaicFT are performed using a
87 : // novel sort-less algorithm. In this approach, the gridded plane is
88 : // divided into small patches, a cache of which is maintained in memory
89 : // using a general-purpose <linkto class=casacore::LatticeCache>LatticeCache</linkto> class. As the (time-sorted)
90 : // visibility data move around slowly in the Fourier plane, patches are
91 : // swapped in and out as necessary. Thus, optimally, one would keep at
92 : // least one patch per baseline.
93 : //
94 : // A grid cache is defined on construction. If the gridded uv plane is smaller
95 : // than this, it is kept entirely in memory and all gridding and
96 : // degridding is done entirely in memory. Otherwise a cache of tiles is
97 : // kept an paged in and out as necessary. Optimally the cache should be
98 : // big enough to hold all polarizations and frequencies for all
99 : // baselines. The paging rate will then be small. As the cache size is
100 : // reduced below this critical value, paging increases. The algorithm will
101 : // work for only one patch but it will be very slow!
102 : //
103 : // This scheme works well for arrays having a moderate number of
104 : // antennas since the saving in space goes as the ratio of
105 : // baselines to image size. For the ATCA, VLBA and WSRT, this ratio is
106 : // quite favorable. For the VLA, one requires images of greater than
107 : // about 200 pixels on a side to make it worthwhile.
108 : //
109 : // The FFT step is done plane by plane for images having less than
110 : // 1024 * 1024 pixels on each plane, and line by line otherwise.
111 : //
112 : // The gridding and degridding steps are implemented in Fortran
113 : // for speed. In gridding, the visibilities are added onto the
114 : // grid points in the neighborhood using a weighting function.
115 : // In degridding, the value is derived by a weight summ of the
116 : // same points, using the same weighting function.
117 : // </synopsis>
118 : //
119 : // <example>
120 : // See the example for <linkto class=SkyModel>SkyModel</linkto>.
121 : // </example>
122 : //
123 : // <motivation>
124 : // Define an interface to allow efficient processing of chunks of
125 : // visibility data
126 : // </motivation>
127 : //
128 : // <todo asof="97/10/01">
129 : // <ul> Deal with large VLA spectral line case
130 : // </todo>
131 : class MosaicFT;
132 :
133 :
134 : namespace refim{
135 : class SimplePBConvFunc;
136 :
137 : class MosaicFT : public FTMachine {
138 : public:
139 :
140 : // Constructor: cachesize is the size of the cache in words
141 : // (e.g. a few million is a good number), tilesize is the
142 : // size of the tile used in gridding (cannot be less than
143 : // 12, 16 works in most cases).
144 : // <group>
145 : // Note: this takes ownership of the sj pointer
146 : MosaicFT(SkyJones* sj, casacore::MPosition mloc, casacore::String stokes,
147 : casacore::Long cachesize, casacore::Int tilesize=16,
148 : casacore::Bool usezero=true, casacore::Bool useDoublePrec=false, casacore::Bool useConjConvFunc=false, casacore::Bool usePointingTable=false);
149 : // </group>
150 :
151 : // Construct from a casacore::Record containing the MosaicFT state
152 : MosaicFT(const casacore::RecordInterface& stateRec);
153 :
154 : // Copy constructor
155 : MosaicFT(const MosaicFT &other);
156 :
157 : // Assignment operator
158 : MosaicFT &operator=(const MosaicFT &other);
159 :
160 : ~MosaicFT();
161 :
162 : // Initialize transform to Visibility plane using the image
163 : // as a template. The image is loaded and Fourier transformed.
164 : void initializeToVis(casacore::ImageInterface<casacore::Complex>& image,
165 : const vi::VisBuffer2& vb);
166 :
167 : // Finalize transform to Visibility plane: flushes the image
168 : // cache and shows statistics if it is being used.
169 : void finalizeToVis();
170 :
171 : // Initialize transform to Sky plane: initializes the image
172 : void initializeToSky(casacore::ImageInterface<casacore::Complex>& image, casacore::Matrix<casacore::Float>& weight,
173 : const vi::VisBuffer2& vb);
174 : //////
175 :
176 : // Finalize transform to Sky plane: flushes the image
177 : // cache and shows statistics if it is being used. DOES NOT
178 : // DO THE FINAL TRANSFORM!
179 : void finalizeToSky();
180 :
181 : // Get actual coherence from grid by degridding
182 : virtual void get(vi::VisBuffer2& vb, casacore::Int row=-1);
183 :
184 :
185 : // Put coherence to grid by gridding.
186 : virtual void put(const vi::VisBuffer2& vb, casacore::Int row=-1, casacore::Bool dopsf=false,
187 : FTMachine::Type type=FTMachine::OBSERVED);
188 :
189 : virtual void gridImgWeights(const vi::VisBuffer2& vb);
190 :
191 : // Make the entire image
192 : void makeImage(FTMachine::Type type,
193 : vi::VisibilityIterator2& vs,
194 : casacore::ImageInterface<casacore::Complex>& image,
195 : casacore::Matrix<casacore::Float>& weight);
196 :
197 : // Get the final image: do the Fourier transform and
198 : // grid-correct, then optionally normalize by the summed weights
199 : casacore::ImageInterface<casacore::Complex>& getImage(casacore::Matrix<casacore::Float>&, casacore::Bool normalize=true);
200 0 : virtual void normalizeImage(casacore::Lattice<casacore::Complex>& /*skyImage*/,
201 : const casacore::Matrix<casacore::Double>& /*sumOfWts*/,
202 : casacore::Lattice<casacore::Float>& /*sensitivityImage*/,
203 : casacore::Bool /*fftNorm*/)
204 0 : {throw(casacore::AipsError("MosaicFT::normalizeImage() called"));}
205 :
206 :
207 : // Get the final weights image
208 : void getWeightImage(casacore::ImageInterface<casacore::Float>&, casacore::Matrix<casacore::Float>&);
209 :
210 : // Get a flux (divide by this to get a flux density correct image)
211 : // image if there is one
212 : virtual void getFluxImage(casacore::ImageInterface<casacore::Float>& image);
213 :
214 : // Save and restore the MosaicFT to and from a record
215 : virtual casacore::Bool toRecord(casacore::String& error, casacore::RecordInterface& outRec,
216 : casacore::Bool withImage=false, const casacore::String diskimage="");
217 : virtual casacore::Bool fromRecord(casacore::String& error, const casacore::RecordInterface& inRec);
218 :
219 : // Can this FTMachine be represented by Fourier convolutions?
220 0 : casacore::Bool isFourier() {return true;}
221 :
222 : // Return name of this machine
223 :
224 : virtual casacore::String name() const;
225 495 : virtual casacore::Bool useWeightImage(){return true;};
226 :
227 :
228 : // Copy convolution function etc to another FT machine
229 : // necessary if ft and ift are distinct but can share convfunctions
230 :
231 : void setConvFunc(casacore::CountedPtr<SimplePBConvFunc>& pbconvFunc);
232 : casacore::CountedPtr<SimplePBConvFunc>& getConvFunc();
233 :
234 : casacore::CountedPtr<casacore::TempImage<casacore::Float> >& getConvWeightImage();
235 :
236 : //reset weight image
237 : virtual void reset();
238 117 : virtual void setMiscInfo(const casacore::Int qualifier){(void)qualifier;};
239 0 : virtual void ComputeResiduals(vi::VisBuffer2&/*vb*/, casacore::Bool /*useCorrected*/) {};
240 :
241 : //Sometimes weightimage is already calculated ...just use it
242 : virtual void setWeightImage(casacore::CountedPtr<casacore::ImageInterface<casacore::Float> >& wgtimage);
243 : protected:
244 :
245 : casacore::Int nint(casacore::Double val) {return casacore::Int(floor(val+0.5));};
246 :
247 : // Find the convolution function
248 : virtual void findConvFunction(const casacore::ImageInterface<casacore::Complex>& image,
249 : const vi::VisBuffer2& vb, const casacore::Matrix<casacore::Double>& uvw);
250 :
251 :
252 :
253 : void addBeamCoverage(casacore::ImageInterface<casacore::Complex>& image);
254 : void prepGridForDegrid();
255 :
256 : // shared: this is assigned in the copy constructor
257 : std::shared_ptr<SkyJones> sj_p;
258 :
259 :
260 : // Get the appropriate data pointer
261 : casacore::Array<casacore::Complex>* getDataPointer(const casacore::IPosition&, casacore::Bool);
262 :
263 : void ok();
264 :
265 : virtual void init(const vi::VisBuffer2& vb);
266 :
267 : // Is this record on Grid? check both ends. This assumes that the
268 : // ends bracket the middle
269 : casacore::Bool recordOnGrid(const vi::VisBuffer2& vb, casacore::Int rownr) const;
270 :
271 :
272 : // Image cache
273 : casacore::LatticeCache<casacore::Complex> * imageCache;
274 :
275 : // Sizes
276 : casacore::Long cachesize;
277 : casacore::Int tilesize;
278 :
279 : // Gridder
280 : std::unique_ptr<casacore::ConvolveGridder<casacore::Double, casacore::Complex>> gridder;
281 :
282 : // Is this tiled?
283 : casacore::Bool isTiled;
284 :
285 : // casacore::Array lattice
286 : casacore::CountedPtr<casacore::Lattice<casacore::Complex> > arrayLattice;
287 :
288 : // Lattice. For non-tiled gridding, this will point to arrayLattice,
289 : // whereas for tiled gridding, this points to the image
290 : casacore::CountedPtr<casacore::Lattice<casacore::Complex> > lattice;
291 : casacore::CountedPtr<casacore::Lattice<casacore::Complex> > weightLattice;
292 :
293 : casacore::Float maxAbsData;
294 :
295 : // Useful IPositions
296 : casacore::IPosition centerLoc, offsetLoc;
297 :
298 : // Image Scaling and offset
299 : casacore::Vector<casacore::Double> uvScale, uvOffset;
300 :
301 : // casacore::Array for non-tiled gridding
302 : casacore::Array<casacore::Complex> griddedWeight;
303 : casacore::Array<casacore::DComplex> griddedWeight2;
304 : // Pointing columns
305 : casacore::MSPointingColumns* mspc;
306 :
307 : // Antenna columns
308 : casacore::MSAntennaColumns* msac;
309 :
310 : casacore::DirectionCoordinate directionCoord;
311 :
312 : casacore::MDirection::Convert* pointingToImage;
313 :
314 : casacore::Vector<casacore::Double> xyPos;
315 :
316 : casacore::MDirection worldPosMeas;
317 :
318 : casacore::Int priorCacheSize;
319 :
320 : // Grid/degrid zero spacing points?
321 : casacore::Bool usezero_p;
322 :
323 : casacore::Array<casacore::Complex> convFunc;
324 : casacore::Array<casacore::Complex> weightConvFunc_p;
325 : casacore::Int convSampling;
326 : casacore::Int convSize;
327 : casacore::Int convSupport;
328 : casacore::Vector<casacore::Int> convSupportPlanes_p;
329 : casacore::Vector<casacore::Int> convSizePlanes_p;
330 : casacore::Vector<casacore::Int> convRowMap_p;
331 : casacore::Vector<casacore::Int> convChanMap_p;
332 : casacore::Vector<casacore::Int> convPolMap_p;
333 :
334 : casacore::Int wConvSize;
335 :
336 : casacore::Int lastIndex_p;
337 :
338 : casacore::Int getIndex(const casacore::MSPointingColumns& mspc, const casacore::Double& time,
339 : const casacore::Double& interval);
340 :
341 : casacore::Bool getXYPos(const vi::VisBuffer2& vb, casacore::Int row);
342 :
343 : casacore::CountedPtr<casacore::TempImage<casacore::Float> >skyCoverage_p;
344 : casacore::TempImage<casacore::Complex>* convWeightImage_p;
345 : casacore::CountedPtr<SimplePBConvFunc> pbConvFunc_p;
346 :
347 : //Later this
348 : casacore::String machineName_p;
349 : casacore::Bool doneWeightImage_p;
350 : casacore::String stokes_p;
351 : casacore::Bool useConjConvFunc_p;
352 : casacore::Bool usePointingTable_p;
353 : casacore::Double timemass_p, timegrid_p, timedegrid_p;
354 :
355 : };
356 :
357 : }// namespace refim ends
358 :
359 : } //# NAMESPACE CASA - END
360 :
361 : #endif
|
