Line data Source code
1 : // -*- C++ -*-
2 : //# CFBuffer.h: Definition of the CFBuffer class
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : #ifndef SYNTHESIS_TRANSFORM2_CFBUFFER_H
29 : #define SYNTHESIS_TRANSFORM2_CFBUFFER_H
30 : #include <synthesis/TransformMachines/SynthesisError.h>
31 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
32 : #include <synthesis/TransformMachines2/CFDefs.h>
33 : #include <synthesis/TransformMachines/CFCell.h>
34 : #include <synthesis/TransformMachines2/Utils.h>
35 : #include <casacore/images/Images/ImageInterface.h>
36 : #include <casacore/casa/Utilities/CountedPtr.h>
37 : #include <casacore/casa/Utilities/Sort.h>
38 : #include <casacore/casa/Logging/LogOrigin.h>
39 : #include <msvis/MSVis/VisBuffer2.h>
40 : #include <casacore/casa/Logging/LogSink.h>
41 : #include <casacore/casa/Logging/LogIO.h>
42 :
43 : //
44 : // <summary> defines interface for the storage for convolution functions </summary>
45 : // <use visibility=export>
46 : //
47 : // <prerequisite>
48 : // </prerequisite>
49 : //
50 : // <etymology>
51 : //
52 : // CFBuffer is the basic in-memory storage for convolution functions
53 : // as a function of polarization, W-value and frequency at a
54 : // particular value of Parallactic Angle and baseline type.
55 : //
56 : //</etymology>
57 : //
58 : // <synopsis>
59 : //
60 : // The CFBuffer class encapsulates the storage and associated
61 : // auxillary information required for the convolution functions. The
62 : // <linto class=CFStore>CFStore</linkto> class is a collection of
63 : // CFBuffer objects. A collection of CFStore objects is held and
64 : // managed by the <linto class=FTMachine>FTMachine</linkto> and the
65 : // appropriate one, depending on the casacore::Time/PA value, polarization and
66 : // frequency of the data in the <linkto
67 : // class=VisBuffer2>VisBuffer2</linkto>, is supplied to the <linkto
68 : // class=VisibilityResampler>VisibilityResampler</linkto> object for
69 : // re-sampling the data onto a grid (or vice versa).
70 : //
71 : // Conceptually, this object holds the convolution functions
72 : // parameterized by the properties of the electromagnetic radiation
73 : // (polarization state, frequency and the w-value which is the Fresnel
74 : // term and implicitly a function of the frequency). The <linkto
75 : // class=CFStore>CFStore</linkto> object holds a list of this object
76 : // index by the telescope related parameters (antenna1,
77 : // antenna2, Parallactic Angle or casacore::Time, etc.).
78 : //
79 : //</synopsis>
80 : //
81 : // <example>
82 : // </example>
83 : //
84 : // <motivation>
85 : //
86 : // To factor out the details of effecient in-memory storage of the
87 : // convolution functions into a separate class. This class can then
88 : // be optmized by specializations for various types of convolution
89 : // functions and imaging algorithms without the need to change the
90 : // imaging framework.
91 : //
92 : // </motivation>
93 : //
94 :
95 : using namespace std;
96 : namespace casacore{
97 :
98 : // extern template class Array<casacore::CountedPtr<CFBuffer> >;
99 : }
100 :
101 : namespace casa {
102 :
103 : using namespace casa::vi;
104 : // template <class T>
105 : typedef casacore::Complex TT;
106 : namespace refim{
107 : using namespace CFDefs;
108 :
109 : struct CFBStruct {
110 : CFCStruct *CFBStorage;
111 : int shape[3];
112 : casacore::Double *freqValues, *wValues, *pointingOffset;
113 : casacore::Double fIncr, wIncr;
114 : casacore::Int **muellerElementsIndex, **muellerElements,
115 : **conjMuellerElementsIndex, **conjMuellerElements,
116 : **conjFreqNdxMap, **freqNdxMap;
117 : casacore::Int nMueller;
118 :
119 :
120 : CFCStruct* getCFB(int i, int j, int k)
121 : { return &(CFBStorage[i + (shape[1]-1)*j + (shape[1]-1)*(shape[2]-1)*k]);}
122 : // { return &(CFBStorage[i + (shape[1]-1)*j + (shape[2]-1)*k]);}
123 : } ;
124 :
125 : class CFBuffer
126 : {
127 : public:
128 : //
129 : //========================= Administrative Parts ==========================
130 : //------------------------------------------------------------------
131 : //
132 184 : CFBuffer(): wValues_p(), maxXSupport_p(-1), maxYSupport_p(-1), pointingOffset_p(), cfHitsStats(),
133 184 : freqNdxMapsReady_p(false), freqNdxMap_p(), conjFreqNdxMap_p(), cfCacheDirName_p(), maxCFSize_p(-1)
134 184 : {};
135 :
136 : CFBuffer(casacore::Int maxXSup, casacore::Int maxYSup):
137 : wValues_p(), maxXSupport_p(maxXSup), maxYSupport_p(maxYSup), pointingOffset_p(), cfHitsStats(),
138 : freqNdxMapsReady_p(false), freqNdxMap_p(), conjFreqNdxMap_p(), cfCacheDirName_p()
139 : {
140 : // storage_p.resize(1,1,1);
141 : // storage_p(0,0,0) = new casacore::Array<TT>(dataPtr);
142 : // coordSys_p.resize(1,1,1);
143 : // coordSys_p(0,0,0) = cs;
144 : };
145 :
146 184 : ~CFBuffer()
147 184 : {
148 : //cerr << "############### " << "~CFBuffer() called" << endl;
149 : // casacore::LogIO log_l(casacore::LogOrigin("CFBuffer","~CFBuffer[R&D]"));
150 : // log_l << "CF Hits stats gathered: " << cfHitsStats << endl;
151 184 : };
152 :
153 : casacore::CountedPtr<CFBuffer> clone();
154 : void allocCells(const casacore::Cube<casacore::CountedPtr<CFCell> >& cells);
155 : void setParams(const CFBuffer& other);
156 : //
157 : //============================= casacore::Functional Parts ============================
158 : //------------------------------------------------------------------
159 : //
160 : // CFBuffer& operator=(const CFBuffer& other);
161 : //
162 : // Get the single convolution function as an casacore::Array<T> for the
163 : // supplied value of the frequency and the muellerElement.
164 : // Mueller element is essentially the polarization product, but
165 : // can be any of the of 16 elements of the outer product.
166 : //
167 : //-------------------------------------------------------------------------
168 : //
169 : inline casacore::Int nChan() {return nChan_p;}
170 : inline casacore::Int nW() {return nW_p;}
171 : inline casacore::Int nMuellerElements() {return nPol_p;}
172 44 : inline casacore::IPosition shape() {casacore::IPosition shp(3,nChan_p, nW_p, nPol_p); return shp;}
173 :
174 : inline casacore::Vector<casacore::Double> getFreqList() {return freqValues_p;};
175 : inline casacore::Vector<casacore::Double> getWList() {return wValues_p;};
176 :
177 : CFCell& getCFCell(const casacore::Double& freqVal, const casacore::Double& wValue,
178 : const casacore::Int & muellerElement);
179 : // muellerElement: (i,j) of the Mueller Matrix
180 : casacore::CountedPtr<CFCell>& getCFCellPtr(const casacore::Double& freqVal, const casacore::Double& wValue,
181 : const casacore::Int & muellerElement);
182 36 : CFCell& operator()(const casacore::Int& i, const casacore::Int& j, const casacore::Int& k) {return *cfCells_p(i,j,k);}
183 : CFCell& getCFCell(const casacore::Int& i, const casacore::Int& j, const casacore::Int& k);
184 :
185 : casacore::CountedPtr<CFCell >& getCFCellPtr(const casacore::Int& i, const casacore::Int& j, const casacore::Int& k);
186 :
187 : //=========================================================================
188 : casacore::Array<TT>& getCF(const casacore::Double& freqVal, const casacore::Double& wValue,
189 : const casacore::Int & muellerElement)
190 : {return *(getCFCell(freqVal, wValue, muellerElement).storage_p);}
191 : // muellerElement: (i,j) of the Mueller Matrix
192 :
193 : casacore::CountedPtr<casacore::Array<TT> >& getCFPtr(const casacore::Double& freqVal, const casacore::Double& wValue,
194 : const casacore::Int & muellerElement)
195 : {return getCFCellPtr(freqVal, wValue, muellerElement)->storage_p;}
196 :
197 : casacore::Array<TT>& getCF(const casacore::Int& i, const casacore::Int& j, const casacore::Int& k)
198 : {return *(getCFCell(i,j,k).storage_p);}
199 :
200 : casacore::CountedPtr<casacore::Array<TT> >& getCFPtr(const casacore::Int& i, const casacore::Int& j, const casacore::Int& k)
201 : {return getCFCellPtr(i,j,k)->storage_p;}
202 :
203 :
204 : //
205 : // Get the parameters of a the CFs indexed by values. The version
206 : // which returns also the casacore::Coordinate System associated with the
207 : // CFs are slow (casacore::CoordinateSystem::operator=() is surprisingly
208 : // expensive!). So do not use this in tight loops. If it is
209 : // required, use the version without the co-ordinate system below.
210 : //
211 : void getParams(casacore::CoordinateSystem& cs, casacore::Float& sampling,
212 : casacore::Int& xSupport, casacore::Int& ySupport, casacore::String& bandName,
213 : const casacore::Double& freqVal, const casacore::Double& wValue,
214 : const casacore::Int& muellerElement);
215 : //-------------------------------------------------------------------------
216 : // Get CF by directly indexing in the list of CFs (data vector)
217 248 : inline void getParams(casacore::CoordinateSystem& cs, casacore::Float& sampling,
218 : casacore::Int& xSupport, casacore::Int& ySupport, casacore::String& bandName,
219 : const casacore::Int& i, const casacore::Int& j, const casacore::Int& k)
220 : {
221 248 : cs = cfCells_p(i,j,k)->coordSys_p;
222 248 : sampling = cfCells_p(i,j,k)->sampling_p;
223 248 : xSupport = cfCells_p(i,j,k)->xSupport_p;
224 248 : ySupport = cfCells_p(i,j,k)->ySupport_p;
225 248 : bandName = cfCells_p(i,j,k)->bandName_p;
226 248 : }
227 3628308 : void getParams(casacore::Double& freqVal, casacore::Float& sampling,
228 : casacore::Int& xSupport, casacore::Int& ySupport,
229 : const casacore::Int& iFreq, const casacore::Int& iW, const casacore::Int& iPol)
230 : {
231 3628308 : sampling = cfCells_p(iFreq,iW,iPol)->sampling_p;
232 3628308 : xSupport = cfCells_p(iFreq,iW,iPol)->xSupport_p;
233 3628308 : ySupport = cfCells_p(iFreq,iW,iPol)->ySupport_p;
234 3628308 : freqVal = freqValues_p(iFreq);
235 3628308 : }
236 :
237 33789 : inline void getCoordList(casacore::Vector<casacore::Double>& freqValues, casacore::Vector<casacore::Double>& wValues,
238 : PolMapType& muellerElementsIndex, PolMapType& muellerElements,
239 : PolMapType& conjMuellerElementsIndex, PolMapType& conjMuellerElements,
240 : casacore::Double& fIncr, casacore::Double& wIncr)
241 : {
242 33789 : freqValues.assign(freqValues_p);wValues.assign(wValues_p);
243 33789 : muellerElements.assign(muellerElements_p); muellerElementsIndex.assign(muellerElementsIndex_p);
244 33789 : conjMuellerElements.assign(conjMuellerElements_p); conjMuellerElementsIndex.assign(conjMuellerElementsIndex_p);
245 33789 : fIncr = freqValIncr_p; wIncr = wValIncr_p;
246 33789 : }
247 :
248 : casacore::Int nearestNdx(const casacore::Double& val, const casacore::Vector<casacore::Double>& valList, const casacore::Double& incr);
249 :
250 : casacore::Int nearestFreqNdx(const casacore::Double& freqVal) ;
251 :
252 3628308 : inline casacore::Int nearestWNdx(const casacore::Double& wVal)
253 : {
254 7256616 : return max(0,
255 : min(
256 10884924 : (int)round(sqrt(wValIncr_p*abs(wVal))),
257 3628308 : (int)(wValues_p.nelements()-1)
258 : )
259 3628308 : ); // CAS-13191
260 : //return max(0,min((int)round(sqrt(abs(wVal)/wValIncr_p)),(int)(wValues_p.nelements()-1))); // CAS-13191
261 : }
262 :
263 : casacore::Double nearest(casacore::Bool& found, const casacore::Double& val, const casacore::Vector<casacore::Double>& valList, const casacore::Double& incr);
264 :
265 : inline casacore::Double nearestFreq(casacore::Bool& found, const casacore::Double& freqVal)
266 : {return nearest(found, freqVal, freqValues_p, freqValIncr_p);}
267 :
268 : inline casacore::Double nearestWVal(casacore::Bool& found, const casacore::Double& wVal)
269 : {return nearest(found, wVal, wValues_p, wValIncr_p);}
270 :
271 : //-------------------------------------------------------------------------
272 : //
273 : // Generate a map for the given frequency and Mueller element list
274 : // to the index in the internal list of CFs. This can be used in
275 : // tight loops to get get direct access to the required CF.
276 : //
277 : void makeCFBufferMap(const casacore::Vector<casacore::Double>& freqVals,
278 : const casacore::Vector<casacore::Double>& wValues,
279 : const MuellerMatrixType& muellerElements);
280 : //-------------------------------------------------------------------------
281 : //
282 : // Add a Convolution casacore::Function with associated parameters.
283 : //
284 : void addCF(casacore::Array<TT>*, //dataPtr,
285 : casacore::CoordinateSystem&,// cs,
286 : casacore::Float& ,//sampling,
287 : casacore::Int& ,//xSupport,
288 : casacore::Int& ,//ySupport,
289 : casacore::Double& ,//freqValue,
290 : casacore::Double& ,//wValue,
291 : casacore::Int& //muellerElement
292 : )
293 : {throw(casacore::AipsError("CFBuffer::addCF called"));}
294 : //-------------------------------------------------------------------------
295 : //
296 0 : void resize(const casacore::IPosition& size) {cfCells_p.resize(size);};
297 : void resize(const casacore::Double& wIncr, const casacore::Double& freqIncr,
298 : const casacore::Vector<casacore::Double>& wValues,
299 : const casacore::Vector<casacore::Double>& freqValues,
300 : const PolMapType& muellerElements,
301 : const PolMapType& muellerElementsIndex,
302 : const PolMapType& conjMuellerElements,
303 : const PolMapType& conjMuellerElementsIndex);
304 : casacore::Int noOfMuellerElements(const PolMapType& muellerElements);
305 : //-------------------------------------------------------------------------
306 : // Set only the CF parameters. Return to index of the CF that was set.
307 : //
308 : casacore::RigidVector<casacore::Int, 3> setParams(const casacore::Int& inu, const casacore::Int& iw, const casacore::Int& ipx, const casacore::Int& ipy,
309 : const casacore::Double& freqValue,
310 : const casacore::Double& wValue,
311 : const casacore::Int& muellerElement,
312 : casacore::CoordinateSystem& cs,
313 : const casacore::TableRecord& miscInfo);
314 :
315 :
316 : casacore::RigidVector<casacore::Int, 3> setParams(const casacore::Int& i, const casacore::Int& j, const casacore::Int& ipx, const casacore::Int& ipy,
317 : const casacore::Double& freqValue, const casacore::String& bandName,
318 : const casacore::Double& wValue,
319 : const casacore::Int& muellerElement,
320 : casacore::CoordinateSystem& cs,
321 : casacore::Float& sampling,
322 : casacore::Int& xSupport, casacore::Int& ySupport,
323 : const casacore::String& fileName=casacore::String(),
324 : const casacore::Double& conjFreq=0.0,
325 : const casacore::Int& conjPol=-1,
326 : const casacore::String& telescopeName=casacore::String("EVLA"),
327 : const casacore::Float& diameter=25.0);
328 : // casacore::RigidVector<casacore::Int, 3> setParams(const casacore::Int& inu, const casacore::Int& iw, const casacore::Int& muellerElement,
329 : // const casacore::TableRecord& miscInfo);
330 3684130 : void setPointingOffset(const casacore::Vector<casacore::Vector<casacore::Double> >& offset)
331 3684130 : {pointingOffset_p.assign(offset);};
332 : casacore::Vector<casacore::Vector<casacore::Double> > getPointingOffset() {return pointingOffset_p;};
333 : //
334 : // Also set the size of the CF in x and y.
335 : //
336 : void setParams(casacore::Int& nx, casacore::Int& ny, casacore::CoordinateSystem& cs, casacore::Float& sampling,
337 : casacore::Int& xSupport, casacore::Int& ySupport,
338 : const casacore::Double& freqVal, const casacore::Double& wValue,
339 : const casacore::Int& muellerElement,
340 : const casacore::String& fileName);
341 : void setPA(casacore::Float& pa);
342 366 : void setDir(const casacore::String& Dir) {cfCacheDirName_p=Dir;}
343 : void clear();
344 1606 : const casacore::String& getCFCacheDir() {return cfCacheDirName_p;};
345 :
346 : casacore::RigidVector<casacore::Int,3> getIndex(const casacore::Double& freqVal, const casacore::Double& wValue,
347 : const casacore::Int& muellerElement);
348 : //-------------------------------------------------------------------------
349 : //
350 : // Copy just the parameters from other to this.
351 : //
352 : void copyParams(const CFBuffer& other)
353 : {
354 : cfCells_p = other.cfCells_p;
355 : // coordSys_p = other.coordSys_p; sampling_p.assign(other.sampling_p);
356 : // xSupport_p.assign(other.xSupport_p); ySupport_p.assign(other.ySupport_p);
357 : maxXSupport_p=other.maxXSupport_p; maxYSupport_p=other.maxYSupport_p;
358 : }
359 : //-------------------------------------------------------------------------
360 : //
361 : // Write the description of the storage on the supplied ostream.
362 : // Used mostly for debugging, but might be useful for user
363 : // feedback/logging.
364 : //
365 : void show(const char *Mesg=NULL,ostream &os=cerr);
366 : //
367 : // Returns true if the internal storage is not yet initialized.
368 : //
369 : casacore::Bool null() {return (cfCells_p.nelements() == 0);};
370 :
371 130724 : casacore::Cube<casacore::CountedPtr<CFCell> >& getStorage() {return cfCells_p;};
372 : void makePersistent(const char *dir, const char *cfName="");
373 :
374 : void primeTheCache();
375 : void initMaps(const VisBuffer2& vb,const casacore::Matrix<casacore::Double>& freqSelection,const casacore::Double& imRefFreq);
376 : void initPolMaps(PolMapType& polMap, PolMapType& conjPolMap);
377 : //
378 : // For CUDA kernel
379 : //
380 : void getFreqNdxMaps(casacore::Vector<casacore::Vector<casacore::Int> >& freqNdx, casacore::Vector<casacore::Vector<casacore::Int> >& conjFreqNdx);
381 3628308 : inline casacore::Int nearestFreqNdx(const casacore::Int& spw, const casacore::Int& chan, const casacore::Bool conj=false)
382 : {
383 : // Single pixel along the freq. axis but chan, spw > 0 indicates
384 : // use of W-only projection on a wide-band data.
385 3628308 : if (nChan_p == 1) return 0;
386 3551571 : if (conj) return conjFreqNdxMap_p[spw][chan];
387 2437623 : else return freqNdxMap_p[spw][chan];
388 : }
389 :
390 : void getAsStruct(CFBStruct& st);
391 :
392 86 : static void initCFBStruct(CFBStruct& cfbSt)
393 : {
394 86 : cfbSt.CFBStorage=NULL;
395 86 : cfbSt.freqValues=NULL;
396 86 : cfbSt.wValues=NULL;
397 86 : cfbSt.muellerElementsIndex=NULL;
398 86 : cfbSt.muellerElements=NULL;
399 86 : cfbSt.conjMuellerElementsIndex=NULL;
400 86 : cfbSt.conjMuellerElements=NULL;
401 86 : cfbSt.shape[0]=cfbSt.shape[1]=cfbSt.shape[2]=0;
402 86 : cfbSt.fIncr=cfbSt.wIncr=0.0;
403 86 : }
404 : void fill(const casacore::Int& nx, const casacore::Int& ny,
405 : const casacore::Vector<casacore::Double>& freqValues,
406 : const casacore::Vector<casacore::Double>& wValues,
407 : const PolMapType& muellerElements);
408 :
409 0 : casacore::IPosition getShape() {return cfCells_p.shape();}
410 :
411 : int getMaxCFSize();
412 : bool finitePointingOffsets();
413 : //
414 : //============================= Protected Parts ============================
415 : //------------------------------------------------------------------
416 : //
417 : protected:
418 : //
419 : // The storage buffer for the pixel values in CFCell is casacore::Array<T>
420 : // rather than casacore::Matrix<T> to accomodate rotationally symmetric CFs
421 : // (like the Prolate Spheroidal) which can be held as a casacore::Vector of
422 : // values.
423 : //
424 : casacore::Cube<casacore::CountedPtr<CFCell> > cfCells_p;// freqValues x wValues x muellerElements
425 : casacore::Vector<casacore::Double> wValues_p, freqValues_p;
426 : PolMapType muellerElements_p, muellerElementsIndex_p,conjMuellerElements_p,conjMuellerElementsIndex_p;
427 : casacore::Double wValIncr_p, freqValIncr_p;
428 : MuellerMatrixType muellerMask_p;
429 :
430 : casacore::Int nPol_p, nChan_p, nW_p, maxXSupport_p, maxYSupport_p;
431 : casacore::Vector<casacore::Vector<casacore::Double> > pointingOffset_p;
432 : casacore::Cube<casacore::Int> cfHitsStats;
433 : casacore::Bool freqNdxMapsReady_p;
434 : casacore::Vector<casacore::Vector<casacore::Int> > freqNdxMap_p, conjFreqNdxMap_p;
435 : void ASSIGNVVofI(casacore::Int** &target,casacore::Vector<casacore::Vector<casacore::Int> >& source, casacore::Bool& doAlloc);
436 : casacore::String cfCacheDirName_p;
437 :
438 : int maxCFSize_p;
439 : };
440 :
441 : // declare a commonly used template extern
442 :
443 : } // end namesapce refim
444 : } // end namespace casa
445 :
446 : #endif
|
