Line data Source code
1 : //# GridFT.h: Definition for GridFT
2 : //# Copyright (C) 1996-2012
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_rGRIDFT_H
30 : #define SYNTHESIS_rGRIDFT_H
31 :
32 : #include <synthesis/TransformMachines/FTMachine.h>
33 : #include <casacore/casa/Arrays/Matrix.h>
34 : #include <casacore/scimath/Mathematics/FFTServer.h>
35 : #include <msvis/MSVis/VisBuffer.h>
36 : #include <casacore/images/Images/ImageInterface.h>
37 : #include <casacore/images/Images/ImageInterface.h>
38 : #include <casacore/casa/Containers/Block.h>
39 : #include <casacore/casa/Arrays/Array.h>
40 : #include <casacore/casa/Arrays/Vector.h>
41 : #include <casacore/casa/Arrays/Matrix.h>
42 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
43 : #include <casacore/lattices/Lattices/LatticeCache.h>
44 : #include <casacore/lattices/Lattices/ArrayLattice.h>
45 : //#include <synthesis/MeasurementComponents/SynthesisPeek.h>
46 :
47 :
48 : namespace casacore{
49 :
50 : class UVWMachine;
51 : }
52 :
53 : namespace casa { //# NAMESPACE CASA - BEGIN
54 :
55 : // <summary> An FTMachine for Gridded Fourier transforms </summary>
56 :
57 : // <use visibility=export>
58 :
59 : // <reviewed reviewer="" date="" tests="" demos="">
60 :
61 : // <prerequisite>
62 : // <li> <linkto class=FTMachine>FTMachine</linkto> module
63 : // <li> <linkto class=SkyEquation>SkyEquation</linkto> module
64 : // <li> <linkto class=VisBuffer>VisBuffer</linkto> module
65 : // </prerequisite>
66 : //
67 : // <etymology>
68 : // FTMachine is a Machine for Fourier Transforms. GridFT does
69 : // Grid-based Fourier transforms.
70 : // </etymology>
71 : //
72 : // <synopsis>
73 : // The <linkto class=SkyEquation>SkyEquation</linkto> needs to be able
74 : // to perform Fourier transforms on visibility data. GridFT
75 : // allows efficient Fourier Transform processing using a
76 : // <linkto class=VisBuffer>VisBuffer</linkto> which encapsulates
77 : // a chunk of visibility (typically all baselines for one time)
78 : // together with all the information needed for processing
79 : // (e.g. UVW coordinates).
80 : //
81 : // Gridding and degridding in GridFT are performed using a
82 : // novel sort-less algorithm. In this approach, the gridded plane is
83 : // divided into small patches, a cache of which is maintained in memory
84 : // using a general-purpose <linkto class=casacore::LatticeCache>LatticeCache</linkto> class. As the (time-sorted)
85 : // visibility data move around slowly in the Fourier plane, patches are
86 : // swapped in and out as necessary. Thus, optimally, one would keep at
87 : // least one patch per baseline.
88 : //
89 : // A grid cache is defined on construction. If the gridded uv plane is smaller
90 : // than this, it is kept entirely in memory and all gridding and
91 : // degridding is done entirely in memory. Otherwise a cache of tiles is
92 : // kept an paged in and out as necessary. Optimally the cache should be
93 : // big enough to hold all polarizations and frequencies for all
94 : // baselines. The paging rate will then be small. As the cache size is
95 : // reduced below this critical value, paging increases. The algorithm will
96 : // work for only one patch but it will be very slow!
97 : //
98 : // This scheme works well for arrays having a moderate number of
99 : // antennas since the saving in space goes as the ratio of
100 : // baselines to image size. For the ATCA, VLBA and WSRT, this ratio is
101 : // quite favorable. For the VLA, one requires images of greater than
102 : // about 200 pixels on a side to make it worthwhile.
103 : //
104 : // The FFT step is done plane by plane for images having less than
105 : // 1024 * 1024 pixels on each plane, and line by line otherwise.
106 : //
107 : // The gridding and degridding steps are implemented in Fortran
108 : // for speed. In gridding, the visibilities are added onto the
109 : // grid points in the neighborhood using a weighting function.
110 : // In degridding, the value is derived by a weight summ of the
111 : // same points, using the same weighting function.
112 : // </synopsis>
113 : //
114 : // <example>
115 : // See the example for <linkto class=SkyModel>SkyModel</linkto>.
116 : // </example>
117 : //
118 : // <motivation>
119 : // Define an interface to allow efficient processing of chunks of
120 : // visibility data
121 : // </motivation>
122 : //
123 : // <todo asof="97/10/01">
124 : // <ul> Deal with large VLA spectral line case
125 : // </todo>
126 :
127 : class rGridFT : public FTMachine {
128 : public:
129 :
130 : // Constructor: cachesize is the size of the cache in words
131 : // (e.g. a few million is a good number), tilesize is the
132 : // size of the tile used in gridding (cannot be less than
133 : // 12, 16 works in most cases), and convType is the type of
134 : // gridding used (SF is prolate spheriodal wavefunction,
135 : // and BOX is plain box-car summation). mLocation is
136 : // the position to be used in some phase rotations. If
137 : // mTangent is specified then the uvw rotation is done for
138 : // that location iso the image center.
139 : // <group>
140 : rGridFT(casacore::Long cachesize, casacore::Int tilesize, casacore::String convType="SF",
141 : casacore::Float padding=1.0, casacore::Bool usezero=true, casacore::Bool useDoublePrec=false);
142 : rGridFT(casacore::Long cachesize, casacore::Int tilesize, casacore::String convType,
143 : casacore::MPosition mLocation, casacore::Float padding=1.0, casacore::Bool usezero=true,
144 : casacore::Bool useDoublePrec=false);
145 : rGridFT(casacore::Long cachesize, casacore::Int tilesize, casacore::String convType,
146 : casacore::MDirection mTangent, casacore::Float padding=1.0, casacore::Bool usezero=true,
147 : casacore::Bool useDoublePrec=false);
148 : rGridFT(casacore::Long cachesize, casacore::Int tilesize, casacore::String convType,
149 : casacore::MPosition mLocation, casacore::MDirection mTangent, casacore::Float passing=1.0,
150 : casacore::Bool usezero=true, casacore::Bool useDoublePrec=false);
151 : // </group>
152 :
153 : // Construct from a casacore::Record containing the GridFT state
154 : rGridFT(const casacore::RecordInterface& stateRec);
155 :
156 : // Copy constructor
157 : rGridFT(const rGridFT &other);
158 :
159 : // Assignment operator
160 : rGridFT &operator=(const rGridFT &other);
161 :
162 : ~rGridFT();
163 :
164 : // Initialize transform to Visibility plane using the image
165 : // as a template. The image is loaded and Fourier transformed.
166 : void initializeToVis(casacore::ImageInterface<casacore::Complex>& image,
167 : const VisBuffer& vb);
168 :
169 : // Finalize transform to Visibility plane: flushes the image
170 : // cache and shows statistics if it is being used.
171 : void finalizeToVis();
172 :
173 : // Initialize transform to Sky plane: initializes the image
174 : void initializeToSky(casacore::ImageInterface<casacore::Complex>& image, casacore::Matrix<casacore::Float>& weight,
175 : const VisBuffer& vb);
176 :
177 :
178 : // Finalize transform to Sky plane: flushes the image
179 : // cache and shows statistics if it is being used. DOES NOT
180 : // DO THE FINAL TRANSFORM!
181 : void finalizeToSky();
182 :
183 :
184 : // Get actual coherence from grid by degridding
185 : void get(VisBuffer& vb, casacore::Int row=-1);
186 :
187 :
188 : // Put coherence to grid by gridding.
189 : void put(const VisBuffer& vb, casacore::Int row=-1, casacore::Bool dopsf=false,
190 : FTMachine::Type type=FTMachine::OBSERVED);
191 :
192 : // Make the entire image
193 : void makeImage(FTMachine::Type type,
194 : VisSet& vs,
195 : casacore::ImageInterface<casacore::Complex>& image,
196 : casacore::Matrix<casacore::Float>& weight);
197 :
198 : // Get the final image: do the Fourier transform and
199 : // grid-correct, then optionally normalize by the summed weights
200 : casacore::ImageInterface<casacore::Complex>& getImage(casacore::Matrix<casacore::Float>&, casacore::Bool normalize=true);
201 0 : virtual void normalizeImage(casacore::Lattice<casacore::Complex>& /*skyImage*/,
202 : const casacore::Matrix<casacore::Double>& /*sumOfWts*/,
203 : casacore::Lattice<casacore::Float>& /*sensitivityImage*/,
204 : casacore::Bool /*fftNorm*/)
205 0 : {throw(casacore::AipsError("GridFT::normalizeImage() called"));}
206 :
207 : // Get the final weights image
208 : void getWeightImage(casacore::ImageInterface<casacore::Float>&, casacore::Matrix<casacore::Float>&);
209 :
210 : // Save and restore the GridFT to and from a record
211 : virtual casacore::Bool toRecord(casacore::String& error, casacore::RecordInterface& outRec,
212 : casacore::Bool withImage=false, const casacore::String diskimage="");
213 : virtual casacore::Bool fromRecord(casacore::String& error, const casacore::RecordInterface& inRec);
214 :
215 : // Can this FTMachine be represented by Fourier convolutions?
216 0 : virtual casacore::Bool isFourier() {return true;}
217 :
218 0 : virtual void setNoPadding(casacore::Bool nopad){noPadding_p=nopad;};
219 :
220 : virtual casacore::String name() const;
221 0 : virtual void setMiscInfo(const casacore::Int qualifier){(void)qualifier;};
222 0 : virtual void ComputeResiduals(VisBuffer&/*vb*/, casacore::Bool /*useCorrected*/) {};
223 :
224 : protected:
225 :
226 :
227 : // Padding in FFT
228 : casacore::Float padding_p;
229 :
230 : // Get the appropriate data pointer
231 : casacore::Array<casacore::Complex>* getDataPointer(const casacore::IPosition&, casacore::Bool);
232 :
233 : void ok();
234 :
235 : void init();
236 :
237 : //Prepare the grid for degridding
238 : void prepGridForDegrid();
239 :
240 : // Is this record on Grid? check both ends. This assumes that the
241 : // ends bracket the middle
242 : casacore::Bool recordOnGrid(const VisBuffer& vb, casacore::Int rownr) const;
243 :
244 :
245 : // Image cache
246 : casacore::LatticeCache<casacore::Complex> * imageCache;
247 :
248 : // Sizes
249 : casacore::Long cachesize;
250 : casacore::Int tilesize;
251 :
252 : // Gridder
253 : casacore::ConvolveGridder<casacore::Double, casacore::Complex>* gridder;
254 :
255 : // Is this tiled?
256 : casacore::Bool isTiled;
257 :
258 : // casacore::Array lattice
259 : casacore::CountedPtr<casacore::Lattice<casacore::Complex> > arrayLattice;
260 :
261 : // Lattice. For non-tiled gridding, this will point to arrayLattice,
262 : // whereas for tiled gridding, this points to the image
263 : casacore::CountedPtr<casacore::Lattice<casacore::Complex> > lattice;
264 :
265 : casacore::String convType;
266 :
267 : casacore::Float maxAbsData;
268 :
269 : // Useful IPositions
270 : casacore::IPosition centerLoc, offsetLoc;
271 :
272 : // Image Scaling and offset
273 : casacore::Vector<casacore::Double> uvScale, uvOffset;
274 :
275 :
276 : casacore::Int priorCacheSize;
277 :
278 : // Grid/degrid zero spacing points?
279 :
280 : casacore::Bool usezero_p;
281 :
282 : //force no padding
283 : casacore::Bool noPadding_p;
284 :
285 : //Check if using put that avoids non-necessary reads
286 : casacore::Bool usePut2_p;
287 :
288 : //machine name
289 : casacore::String machineName_p;
290 :
291 : casacore::Double timemass_p, timegrid_p;
292 : // casa::async::SynthesisAsyncPeek *peek;
293 : };
294 :
295 : } //# NAMESPACE CASA - END
296 :
297 : #endif
|
