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