Line data Source code
1 : //# SDGrid.h: Definition for SDGrid
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
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: 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_SDGRID_H
30 : #define SYNTHESIS_SDGRID_H
31 :
32 : #include <casacore/casa/Arrays/Array.h>
33 : #include <casacore/casa/Arrays/Matrix.h>
34 : #include <casacore/casa/Arrays/Vector.h>
35 : #include <casacore/casa/Containers/Block.h>
36 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
37 : #include <casacore/images/Images/ImageInterface.h>
38 : #include <casacore/lattices/Lattices/ArrayLattice.h>
39 : #include <casacore/lattices/Lattices/LatticeCache.h>
40 : #include <casacore/measures/Measures/Measure.h>
41 : #include <casacore/measures/Measures/MDirection.h>
42 : #include <casacore/measures/Measures/MPosition.h>
43 : #include <casacore/ms/MeasurementSets/MSColumns.h>
44 : #include <msvis/MSVis/VisBuffer2.h>
45 : #include <msvis/MSVis/VisibilityIterator2.h>
46 : #include <casacore/scimath/Mathematics/FFTServer.h>
47 : #include <synthesis/TransformMachines2/FTMachine.h>
48 : #include <synthesis/TransformMachines2/SkyJones.h>
49 : #include <synthesis/Utilities/SDPosInterpolator.h>
50 :
51 : namespace casa { //# NAMESPACE CASA - BEGIN
52 :
53 : namespace refim { //# namespace for imaging refactor
54 :
55 : // <summary> An FTMachine for Gridding Single Dish data
56 : // </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. SDGrid does
70 : // Single Dish gridding in a similar way
71 : // </etymology>
72 : //
73 : // <synopsis>
74 : // The <linkto class=SkyEquation>SkyEquation</linkto> needs to be able
75 : // to perform Fourier transforms on visibility data and to grid
76 : // single dish data.
77 : // SDGrid allows efficient Single Dish processing using a
78 : // <linkto class=VisBuffer>VisBuffer</linkto> which encapsulates
79 : // a chunk of visibility (typically all baselines for one time)
80 : // together with all the information needed for processing
81 : // (e.g. direction coordinates).
82 : //
83 : // Gridding and degridding in SDGrid are performed using a
84 : // novel sort-less algorithm. In this approach, the gridded plane is
85 : // divided into small patches, a cache of which is maintained in memory
86 : // using a general-purpose <linkto class=casacore::LatticeCache>LatticeCache</linkto> class. As the (time-sorted)
87 : // visibility data move around slowly in the image plane, patches are
88 : // swapped in and out as necessary. Thus, optimally, one would keep at
89 : // least one patch per scan line of data.
90 : //
91 : // A grid cache is defined on construction. If the gridded image plane is smaller
92 : // than this, it is kept entirely in memory and all gridding and
93 : // degridding is done entirely in memory. Otherwise a cache of tiles is
94 : // kept an paged in and out as necessary. Optimally the cache should be
95 : // big enough to hold all polarizations and frequencies for one
96 : // complete scan line.
97 : // The paging rate will then be small. As the cache size is
98 : // reduced below this critical value, paging increases. The algorithm will
99 : // work for only one patch but it will be very slow!
100 : //
101 : // The gridding and degridding steps are implemented in Fortran
102 : // for speed. In gridding, the visibilities are added onto the
103 : // grid points in the neighborhood using a weighting function.
104 : // In degridding, the value is derived by a weight summ of the
105 : // same points, using the same weighting function.
106 : // </synopsis>
107 : //
108 : // <example>
109 : // See the example for <linkto class=SkyModel>SkyModel</linkto>.
110 : // </example>
111 : //
112 : // <motivation>
113 : // Define an interface to allow efficient processing of chunks of
114 : // visibility data
115 : // </motivation>
116 : //
117 : // <todo asof="97/10/01">
118 : // <ul> Deal with large VLA spectral line case
119 : // </todo>
120 :
121 : class SDGrid final : public FTMachine {
122 : public:
123 :
124 : // Constructor: cachesize is the size of the cache in words
125 : // (e.g. a few million is a good number), tilesize is the
126 : // size of the tile used in gridding (cannot be less than
127 : // 12, 16 works in most cases), and convType is the type of
128 : // gridding used (SF is prolate spheriodal wavefunction,
129 : // and BOX is plain box-car summation). mLocation is
130 : // the position to be used in some phase rotations. If
131 : // mTangent is specified then the uvw rotation is done for
132 : // that location iso the image center. userSupport is to allow
133 : // larger support for the convolution if the user wants it ..-1 will
134 : // use the default i.e 1 for BOX and 3 for others
135 : // USEIMAGINGWEIGHT
136 : // The parameter useImagingWeight in the constructors is to explicitly
137 : // use vb.imagingweight while gridding,
138 : // When doing just SD imaging then setting it to false is fine (in fact recommended as vb.imagingweight
139 : // is set to zero if any pol is flagged this may change later .....today being 2014/08/06)
140 : // when using it in conjuction with interferometer gridding then set useImagingWeight to true
141 : // this is to allow for proper non natural weighting scheme while imaging
142 : // <group>
143 : SDGrid(SkyJones& sj, casacore::Int cachesize, casacore::Int tilesize,
144 : casacore::String convType="BOX", casacore::Int userSupport=-1, casacore::Bool useImagingWeight=false);
145 : SDGrid(casacore::MPosition& ml, SkyJones& sj, casacore::Int cachesize,
146 : casacore::Int tilesize, casacore::String convType="BOX", casacore::Int userSupport=-1,
147 : casacore::Float minweight=0., casacore::Bool clipminmax=false, casacore::Bool useImagingWeight=false);
148 : SDGrid(casacore::Int cachesize, casacore::Int tilesize,
149 : casacore::String convType="BOX", casacore::Int userSupport=-1, casacore::Bool useImagingWeight=false);
150 : SDGrid(casacore::MPosition& ml, casacore::Int cachesize, casacore::Int tilesize,
151 : casacore::String convType="BOX", casacore::Int userSupport=-1, casacore::Float minweight=0., casacore::Bool clipminmax=false,
152 : casacore::Bool useImagingWeight=false);
153 : SDGrid(casacore::MPosition& ml, casacore::Int cachesize, casacore::Int tilesize,
154 : casacore::String convType="TGAUSS", casacore::Float truncate=-1.0,
155 : casacore::Float gwidth=0.0, casacore::Float jwidth=0.0, casacore::Float minweight=0., casacore::Bool clipminmax=false,
156 : casacore::Bool useImagingWeight=false);
157 : // </group>
158 :
159 : // Copy constructor
160 : SDGrid(const SDGrid &other);
161 :
162 : // Assignment operator
163 : SDGrid &operator=(const SDGrid &other);
164 :
165 : ~SDGrid();
166 :
167 : // Initialize transform to Visibility plane using the image
168 : // as a template. The image is loaded and Fourier transformed.
169 : void initializeToVis(casacore::ImageInterface<casacore::Complex>& image,
170 : const vi::VisBuffer2& vb);
171 :
172 : // Finalize transform to Visibility plane: flushes the image
173 : // cache and shows statistics if it is being used.
174 : void finalizeToVis();
175 :
176 : // Initialize transform to Sky plane: initializes the image
177 : void initializeToSky(casacore::ImageInterface<casacore::Complex>& image, casacore::Matrix<casacore::Float>& weight,
178 : const vi::VisBuffer2& vb);
179 :
180 : // Finalize transform to Sky plane: flushes the image
181 : // cache and shows statistics if it is being used. DOES NOT
182 : // DO THE FINAL TRANSFORM!
183 : void finalizeToSky();
184 :
185 : // Get actual coherence from grid by degridding
186 : void get(vi::VisBuffer2& vb, casacore::Int row=-1);
187 :
188 : // Put coherence to grid by gridding.
189 : void put(const vi::VisBuffer2& vb, casacore::Int row=-1, casacore::Bool dopsf=false,
190 : FTMachine::Type type=FTMachine::OBSERVED);
191 :
192 : // Make the entire image using a ROVisIter...
193 : // This is an overload for FTMachine version as
194 : //SDGrid now does everything in memory
195 : // so for large cube ..proceed by slices that fit in memory here.
196 : virtual 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:
202 : // 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("SDGrid::normalizeImage() called"));}
209 :
210 : // SDGrid needs to fill weightimage
211 0 : virtual casacore::Bool useWeightImage(){return true;};
212 :
213 : // Get the final weights image
214 : void getWeightImage(casacore::ImageInterface<casacore::Float>&, casacore::Matrix<casacore::Float>&);
215 :
216 : // Has this operator changed since the last application?
217 : virtual casacore::Bool changed(const vi::VisBuffer2& vb);
218 0 : virtual void setMiscInfo(const casacore::Int qualifier){(void)qualifier;};
219 0 : virtual void ComputeResiduals(vi::VisBuffer2& /*vb*/, casacore::Bool /*useCorrected*/) {};
220 :
221 : // Interpolation-Conversion processing scheme
222 : enum class ConvertFirst {
223 : NEVER = 0,
224 : ALWAYS = 1,
225 : AUTO = 2
226 : };
227 : static const casacore::String & toString(const ConvertFirst convertFirst);
228 : static ConvertFirst fromString(const casacore::String & name);
229 : void setConvertFirst(const casacore::String &convertFirst);
230 :
231 : virtual casacore::String name() const;
232 :
233 : private:
234 0 : casacore::Bool isSD() const override {return true;}
235 :
236 : virtual void initUVWMachine(const vi::VisBuffer2& vb) override;
237 :
238 : // Find the Primary beam and convert it into a convolution buffer
239 : void findPBAsConvFunction(const casacore::ImageInterface<casacore::Complex>& image,
240 : const vi::VisBuffer2& vb);
241 :
242 : SkyJones* sj_p;
243 :
244 : // Get the appropriate data pointer
245 : casacore::Array<casacore::Complex>* getDataPointer(const casacore::IPosition&, casacore::Bool);
246 : casacore::Array<casacore::Float>* getWDataPointer(const casacore::IPosition&, casacore::Bool);
247 :
248 : void ok();
249 :
250 : void init();
251 :
252 : // Image cache
253 : casacore::LatticeCache<casacore::Complex> * imageCache;
254 : casacore::LatticeCache<casacore::Float> * wImageCache;
255 :
256 : // Sizes
257 : casacore::Int cachesize, tilesize;
258 :
259 : // Is this tiled?
260 : casacore::Bool isTiled;
261 :
262 : // Storage for weights
263 : casacore::ImageInterface<casacore::Float>* wImage;
264 :
265 : // casacore::Array lattice
266 : casacore::Lattice<casacore::Complex> * arrayLattice;
267 : casacore::Lattice<casacore::Float> * wArrayLattice;
268 :
269 : // Lattice. For non-tiled gridding, this will point to arrayLattice,
270 : // whereas for tiled gridding, this points to the image
271 : casacore::Lattice<casacore::Complex>* lattice;
272 : casacore::Lattice<casacore::Float>* wLattice;
273 :
274 : casacore::String convType;
275 :
276 : // Useful IPositions
277 : casacore::IPosition centerLoc, offsetLoc;
278 :
279 : // casacore::Array for non-tiled gridding
280 : casacore::Array<casacore::Float> wGriddedData;
281 :
282 :
283 : casacore::DirectionCoordinate directionCoord;
284 :
285 : casacore::MDirection::Convert* pointingToImage;
286 :
287 : casacore::Vector<casacore::Double> xyPos;
288 : //Original xypos of moving source
289 : casacore::Vector<casacore::Double> xyPosMovingOrig_p;
290 :
291 : casacore::MDirection worldPosMeas;
292 :
293 : casacore::Cube<casacore::Int> flags;
294 :
295 : casacore::Vector<casacore::Float> convFunc;
296 : casacore::Int convSampling;
297 : casacore::Int convSize;
298 : casacore::Int convSupport;
299 : casacore::Int userSetSupport_p;
300 :
301 : casacore::Float truncate_p;
302 : casacore::Float gwidth_p;
303 : casacore::Float jwidth_p;
304 :
305 : casacore::Float minWeight_p;
306 :
307 : casacore::Int lastIndex_p;
308 : casacore::Vector<casacore::Int> lastIndexPerAnt_p;
309 : casacore::Bool useImagingWeight_p;
310 : casacore::Int lastAntID_p;
311 : casacore::Int msId_p;
312 :
313 : casacore::Bool isSplineInterpolationReady;
314 : SDPosInterpolator* interpolator;
315 :
316 : // for minmax clipping
317 : casacore::Bool clipminmax_;
318 : casacore::Array<casacore::Complex> gmin_;
319 : casacore::Array<casacore::Complex> gmax_;
320 : casacore::Array<casacore::Float> wmin_;
321 : casacore::Array<casacore::Float> wmax_;
322 : casacore::Array<casacore::Int> npoints_;
323 : void clipMinMax();
324 :
325 : // Interpolation-Conversion processing scheme
326 : ConvertFirst convertFirst;
327 : casacore::MSPointing ramPointingTable;
328 : std::shared_ptr<casacore::MSPointingColumns> ramPointingColumnsPtr;
329 : void convertPointingColumn(
330 : const MeasurementSet &ms,
331 : const MSPointingEnums::PredefinedColumns columnToConvert,
332 : const MDirection::Types directionRef
333 : );
334 : casacore::Bool mustConvertPointingColumn(const casacore::MeasurementSet &ms);
335 : void handleNewMs(const MeasurementSet &ms, ImageInterface<Complex>& image);
336 : void handleNewMs(const casacore::MeasurementSet & ms,
337 : casacore::CountedPtr<SIImageStore> imstore);
338 :
339 : casacore::Int getIndex(const casacore::MSPointingColumns& mspc, const casacore::Double& time,
340 : const casacore::Double& interval=-1.0, const casacore::Int& antid=-1);
341 :
342 : casacore::Bool getXYPos(const vi::VisBuffer2& vb, casacore::Int row);
343 :
344 : //get the casacore::MDirection from a chosen column of pointing table
345 : casacore::MDirection directionMeas(const casacore::MSPointingColumns& mspc, const casacore::Int& index);
346 : casacore::MDirection directionMeas(const casacore::MSPointingColumns& mspc, const casacore::Int& index, const casacore::Double& time);
347 : casacore::MDirection interpolateDirectionMeas(const casacore::MSPointingColumns& mspc, const casacore::Double& time,
348 : const casacore::Int& index, const casacore::Int& index1, const casacore::Int& index2);
349 :
350 : void pickWeights(const vi::VisBuffer2&vb, casacore::Matrix<casacore::Float>& weight);
351 :
352 : //for debugging
353 : //FILE *pfile;
354 : };
355 :
356 : } //End of namespace refim
357 : } //# NAMESPACE CASA - END
358 :
359 : #endif
|
