Line data Source code
1 : //# MosaicFT.h: Definition for MosaicFT
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2002
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: 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 : //#
27 : //# $Id$
28 :
29 : #ifndef SYNTHESIS_MOSAICFT_H
30 : #define SYNTHESIS_MOSAICFT_H
31 :
32 : #include <synthesis/TransformMachines/FTMachine.h>
33 : #include <synthesis/TransformMachines/SkyJones.h>
34 : #include <casacore/casa/Arrays/Matrix.h>
35 : #include <casacore/scimath/Mathematics/FFTServer.h>
36 : #include <msvis/MSVis/VisBuffer.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 : class SimplePBConvFunc;
133 :
134 :
135 : class MosaicFT : public FTMachine {
136 : public:
137 :
138 : // Constructor: cachesize is the size of the cache in words
139 : // (e.g. a few million is a good number), tilesize is the
140 : // size of the tile used in gridding (cannot be less than
141 : // 12, 16 works in most cases).
142 : // <group>
143 : // Note: this takes ownership of the sj pointer
144 : MosaicFT(SkyJones* sj, casacore::MPosition mloc, casacore::String stokes,
145 : casacore::Long cachesize, casacore::Int tilesize=16,
146 : casacore::Bool usezero=true, casacore::Bool useDoublePrec=false);
147 : // </group>
148 :
149 : // Construct from a casacore::Record containing the MosaicFT state
150 : MosaicFT(const casacore::RecordInterface& stateRec);
151 :
152 : // Copy constructor--convolution function is referenced
153 : MosaicFT(const MosaicFT &other);
154 :
155 : // Assignment operator -- convolution function is referenced
156 : MosaicFT &operator=(const MosaicFT &other);
157 :
158 : ~MosaicFT();
159 :
160 : // Initialize transform to Visibility plane using the image
161 : // as a template. The image is loaded and Fourier transformed.
162 : void initializeToVis(casacore::ImageInterface<casacore::Complex>& image,
163 : const VisBuffer& vb);
164 :
165 : // Finalize transform to Visibility plane: flushes the image
166 : // cache and shows statistics if it is being used.
167 : void finalizeToVis();
168 :
169 : // Initialize transform to Sky plane: initializes the image
170 : void initializeToSky(casacore::ImageInterface<casacore::Complex>& image, casacore::Matrix<casacore::Float>& weight,
171 : const VisBuffer& vb);
172 : //////
173 :
174 : // Finalize transform to Sky plane: flushes the image
175 : // cache and shows statistics if it is being used. DOES NOT
176 : // DO THE FINAL TRANSFORM!
177 : void finalizeToSky();
178 :
179 : // Get actual coherence from grid by degridding
180 : void get(VisBuffer& vb, casacore::Int row=-1);
181 :
182 :
183 : // Put coherence to grid by gridding.
184 : void put(const VisBuffer& vb, casacore::Int row=-1, casacore::Bool dopsf=false,
185 : FTMachine::Type type=FTMachine::OBSERVED);
186 :
187 : // Make the entire image
188 : void makeImage(FTMachine::Type type,
189 : VisSet& vs,
190 : casacore::ImageInterface<casacore::Complex>& image,
191 : casacore::Matrix<casacore::Float>& weight);
192 :
193 : // Get the final image: do the Fourier transform and
194 : // grid-correct, then optionally normalize by the summed weights
195 : casacore::ImageInterface<casacore::Complex>& getImage(casacore::Matrix<casacore::Float>&, casacore::Bool normalize=true);
196 0 : virtual void normalizeImage(casacore::Lattice<casacore::Complex>& /*skyImage*/,
197 : const casacore::Matrix<casacore::Double>& /*sumOfWts*/,
198 : casacore::Lattice<casacore::Float>& /*sensitivityImage*/,
199 : casacore::Bool /*fftNorm*/)
200 0 : {throw(casacore::AipsError("MosaicFT::normalizeImage() called"));}
201 :
202 :
203 : // Get the final weights image
204 : void getWeightImage(casacore::ImageInterface<casacore::Float>&, casacore::Matrix<casacore::Float>&);
205 :
206 : // Get a flux (divide by this to get a flux density correct image)
207 : // image if there is one
208 : virtual void getFluxImage(casacore::ImageInterface<casacore::Float>& image);
209 :
210 : // Save and restore the MosaicFT to and from a record
211 :
212 : casacore::Bool toRecord(casacore::String& error, casacore::RecordInterface& outRec,
213 : casacore::Bool withImage=false, const casacore::String diskimage="");
214 : casacore::Bool fromRecord(casacore::String& error, const casacore::RecordInterface& inRec);
215 :
216 : // Can this FTMachine be represented by Fourier convolutions?
217 0 : casacore::Bool isFourier() {return true;}
218 :
219 : // Return name of this machine
220 :
221 : virtual casacore::String name() const;
222 0 : virtual casacore::Bool useWeightImage(){return true;};
223 :
224 :
225 : // Copy convolution function etc to another FT machine
226 : // necessary if ft and ift are distinct but can share convfunctions
227 :
228 : void setConvFunc(casacore::CountedPtr<SimplePBConvFunc>& pbconvFunc);
229 : casacore::CountedPtr<SimplePBConvFunc>& getConvFunc();
230 :
231 : casacore::CountedPtr<casacore::TempImage<casacore::Float> >& getConvWeightImage();
232 :
233 : //reset weight image
234 : virtual void reset();
235 0 : virtual void setMiscInfo(const casacore::Int qualifier){(void)qualifier;};
236 0 : virtual void ComputeResiduals(VisBuffer&/*vb*/, casacore::Bool /*useCorrected*/) {};
237 :
238 : protected:
239 :
240 : casacore::Int nint(casacore::Double val) {return casacore::Int(floor(val+0.5));};
241 :
242 : // Find the convolution function
243 : void findConvFunction(const casacore::ImageInterface<casacore::Complex>& image,
244 : const VisBuffer& vb);
245 :
246 : // void girarUVW(casacore::Matrix<casacore::Double>& uvw, casacore::Vector<casacore::Double>& dphase,
247 : // const VisBuffer& vb);
248 :
249 : void addBeamCoverage(casacore::ImageInterface<casacore::Complex>& image);
250 : void prepGridForDegrid();
251 :
252 : // shared: this is assigned in the copy constructor
253 : std::shared_ptr<SkyJones> sj_p;
254 :
255 :
256 : // Get the appropriate data pointer
257 : casacore::Array<casacore::Complex>* getDataPointer(const casacore::IPosition&, casacore::Bool);
258 :
259 : void ok();
260 :
261 : void init();
262 :
263 : // Is this record on Grid? check both ends. This assumes that the
264 : // ends bracket the middle
265 : casacore::Bool recordOnGrid(const VisBuffer& vb, casacore::Int rownr) const;
266 :
267 :
268 : // Image cache
269 : casacore::LatticeCache<casacore::Complex> * imageCache;
270 :
271 : // Sizes
272 : casacore::Long cachesize;
273 : casacore::Int tilesize;
274 :
275 : // Gridder
276 : std::unique_ptr<casacore::ConvolveGridder<casacore::Double, casacore::Complex>> gridder;
277 :
278 : // Is this tiled?
279 : casacore::Bool isTiled;
280 :
281 : // casacore::Array lattice
282 : casacore::CountedPtr<casacore::Lattice<casacore::Complex> > arrayLattice;
283 :
284 : // Lattice. For non-tiled gridding, this will point to arrayLattice,
285 : // whereas for tiled gridding, this points to the image
286 : casacore::CountedPtr<casacore::Lattice<casacore::Complex> > lattice;
287 : casacore::CountedPtr<casacore::Lattice<casacore::Complex> > weightLattice;
288 :
289 : casacore::Float maxAbsData;
290 :
291 : // Useful IPositions
292 : casacore::IPosition centerLoc, offsetLoc;
293 :
294 : // Image Scaling and offset
295 : casacore::Vector<casacore::Double> uvScale, uvOffset;
296 :
297 : // casacore::Array for non-tiled gridding
298 : casacore::Array<casacore::Complex> griddedWeight;
299 : casacore::Array<casacore::DComplex> griddedWeight2;
300 : // Pointing columns
301 : casacore::MSPointingColumns* mspc;
302 :
303 : // Antenna columns
304 : casacore::MSAntennaColumns* msac;
305 :
306 : casacore::DirectionCoordinate directionCoord;
307 :
308 : casacore::MDirection::Convert* pointingToImage;
309 :
310 : casacore::Vector<casacore::Double> xyPos;
311 :
312 : casacore::MDirection worldPosMeas;
313 :
314 : casacore::Int priorCacheSize;
315 :
316 : // Grid/degrid zero spacing points?
317 : casacore::Bool usezero_p;
318 :
319 : casacore::Array<casacore::Complex> convFunc;
320 : casacore::Array<casacore::Complex> weightConvFunc_p;
321 : casacore::Int convSampling;
322 : casacore::Int convSize;
323 : casacore::Int convSupport;
324 : casacore::Vector<casacore::Int> convSupportPlanes_p;
325 : casacore::Vector<casacore::Int> convSizePlanes_p;
326 : casacore::Vector<casacore::Int> convRowMap_p;
327 : casacore::Vector<casacore::Int> convChanMap_p;
328 : casacore::Vector<casacore::Int> convPolMap_p;
329 :
330 : casacore::Int wConvSize;
331 :
332 : casacore::Int lastIndex_p;
333 :
334 : casacore::Int getIndex(const casacore::MSPointingColumns& mspc, const casacore::Double& time,
335 : const casacore::Double& interval);
336 :
337 : casacore::Bool getXYPos(const VisBuffer& vb, casacore::Int row);
338 :
339 :
340 : casacore::CountedPtr<casacore::TempImage<casacore::Float> >skyCoverage_p;
341 : casacore::TempImage<casacore::Complex>* convWeightImage_p;
342 : casacore::CountedPtr<SimplePBConvFunc> pbConvFunc_p;
343 : //Later this
344 : casacore::String machineName_p;
345 : casacore::Bool doneWeightImage_p;
346 : casacore::String stokes_p;
347 :
348 :
349 : };
350 :
351 : } //# NAMESPACE CASA - END
352 :
353 : #endif
|
