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