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