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