Line data Source code
1 : // -*- C++ -*-
2 : //# CFBuffer.cc: Implementation 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 : #include <synthesis/TransformMachines2/CFBuffer.h>
29 : #include <synthesis/TransformMachines2/Utils.h>
30 : #include <casacore/casa/Utilities/BinarySearch.h>
31 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
32 : #include <casacore/casa/Utilities/Assert.h>
33 :
34 : using namespace casacore;
35 : using namespace casa::vi;
36 : using namespace casacore;
37 :
38 : namespace casa{
39 : using namespace refim;
40 : // Instantiate a commonly used extern template
41 :
42 : //
43 : //---------------------------------------------------------------
44 : //
45 : //
46 0 : void CFBuffer::setParams(const CFBuffer& other)
47 : {
48 0 : wValues_p.assign(other.wValues_p);
49 0 : freqValues_p.assign(other.freqValues_p);
50 0 : muellerElements_p.assign(other.muellerElements_p);
51 0 : muellerElementsIndex_p.assign(other.muellerElementsIndex_p);
52 0 : conjMuellerElements_p.assign(other.conjMuellerElements_p);
53 0 : conjMuellerElementsIndex_p.assign(other.conjMuellerElementsIndex_p);
54 0 : wValIncr_p = other.wValIncr_p;
55 0 : freqValIncr_p = other.freqValIncr_p;
56 0 : muellerMask_p.assign(other.muellerMask_p);
57 0 : pointingOffset_p.assign(other.pointingOffset_p);
58 0 : freqNdxMapsReady_p = other.freqNdxMapsReady_p;
59 :
60 0 : freqNdxMap_p.assign(other.freqNdxMap_p);
61 0 : for(uInt i=0;i<freqNdxMap_p.nelements();i++) freqNdxMap_p[i].assign(other.freqNdxMap_p[i]);
62 0 : conjFreqNdxMap_p.assign(other.conjFreqNdxMap_p);
63 0 : for(uInt i=0;i<conjFreqNdxMap_p.nelements();i++) conjFreqNdxMap_p[i].assign(other.conjFreqNdxMap_p[i]);
64 0 : cfCacheDirName_p = other.cfCacheDirName_p;
65 0 : }
66 : //---------------------------------------------------------------
67 : //
68 : //
69 0 : CountedPtr<CFBuffer> CFBuffer::clone()
70 : {
71 0 : CountedPtr<CFBuffer> clone=new CFBuffer();
72 0 : clone->setParams(*this);
73 :
74 0 : IPosition shp(cfCells_p.shape());
75 0 : clone->resize(shp);
76 0 : clone->allocCells(cfCells_p);
77 : // clone->show("####CLONE: ",cerr);
78 0 : return clone;
79 : }
80 0 : void CFBuffer::allocCells(const Cube<CountedPtr<CFCell> >& cells)
81 : {
82 0 : IPosition shp(cells.shape());
83 0 : for (Int i=0;i < shp[0]; i++)
84 0 : for (Int j=0; j < shp[1]; j++)
85 0 : for (Int k=0; k < shp[2]; k++)
86 : {
87 0 : cfCells_p(i,j,k) = cells(i,j,k)->clone();
88 : }
89 0 : }
90 : //---------------------------------------------------------------
91 : //
92 : // template<class T> Int CFBuffer<T>
93 196 : Int CFBuffer::noOfMuellerElements(const PolMapType& muellerElements)
94 : {
95 196 : Int n=0,nrows = muellerElements.nelements();
96 588 : for (Int i=0;i<nrows;i++)
97 392 : n+=muellerElements(i).nelements();
98 196 : return n;
99 : // Int n=0;
100 : // for(Int i=0;i<muellerElements.nelements();i++)
101 : // for(Int j=0;j<muellerElements(i).nelements();j++)
102 : // if (muellerElements(i)(j) > n) n=muellerElements(i)(j);
103 : // return n;
104 : }
105 : //
106 : //---------------------------------------------------------------
107 : //
108 439 : void CFBuffer::clear()
109 : {
110 878 : IPosition shp(cfCells_p.shape());
111 1744 : for (Int i=0;i < shp[0]; i++)
112 2655 : for (Int j=0; j < shp[1]; j++)
113 4050 : for (Int k=0; k < shp[2]; k++)
114 : {
115 2700 : cfCells_p(i,j,k)->clear();
116 : }
117 439 : }
118 : //
119 : //---------------------------------------------------------------
120 : //
121 : // template<class T> void CFBuffer<T>
122 196 : void CFBuffer::resize(const Double& wIncr, const Double& freqIncr,
123 : const Vector<Double>& wValues,
124 : const Vector<Double>& freqValues,
125 : const PolMapType& muellerElements,
126 : const PolMapType& muellerElementsIndex,
127 : const PolMapType& conjMuellerElements,
128 : const PolMapType& conjMuellerElementsIndex
129 : )
130 : {
131 196 : wValues_p.assign(wValues);
132 196 : freqValues_p.assign(freqValues);
133 196 : wValIncr_p = wIncr;
134 196 : freqValIncr_p = freqIncr;
135 : // muellerMask_p.assign(muellerElements);
136 :
137 : // nPol_p=noOfMuellerElements(muellerMask_p);
138 196 : nPol_p=noOfMuellerElements(muellerElementsIndex);
139 196 : nChan_p = freqValues_p.nelements();
140 196 : nW_p = wValues_p.nelements();
141 :
142 : // muellerElements_p.resize(nPol_p);
143 196 : muellerElements_p.assign(muellerElements);
144 196 : muellerElementsIndex_p.assign(muellerElementsIndex);
145 196 : conjMuellerElements_p.assign(conjMuellerElements);
146 196 : conjMuellerElementsIndex_p.assign(conjMuellerElementsIndex);
147 : // Resize the various aux. information storage buffers.
148 : //
149 : // Resize the storage. Retain the value of the existing pixels.
150 : // New pixels due to resize, if any, are assigned a new Array<T>
151 : // pointer.
152 196 : cfCells_p.resize(nChan_p, nW_p, nPol_p, true);
153 :
154 776 : for (uInt i=0;i<cfCells_p.shape()(0);i++) // nChan_p
155 1190 : for (uInt j=0;j<cfCells_p.shape()(1);j++) // nW_p
156 : {
157 610 : Int k=0; // nPol_p
158 1830 : for(uInt prow=0;prow<muellerElements_p.nelements();prow++)
159 2440 : for(uInt pcol=0;pcol<muellerElements_p(prow).nelements();pcol++)
160 : {
161 1220 : if (cfCells_p(i,j,k).null()) cfCells_p(i,j,k) = new CFCell;
162 1220 : if (cfCells_p(i,j,k)->storage_p.null()) cfCells_p(i,j,k)->storage_p=new Array<TT>;
163 1220 : cfCells_p(i,j,k)->freqValue_p = freqValues(i);
164 1220 : cfCells_p(i,j,k)->freqIncr_p = freqIncr;
165 1220 : cfCells_p(i,j,k)->wValue_p = wValues(j);
166 1220 : cfCells_p(i,j,k)->wIncr_p = wValIncr_p;
167 1220 : cfCells_p(i,j,k)->muellerElement_p = muellerElements_p(prow)(pcol);
168 1220 : cfCells_p(i,j,k)->xSupport_p = 0;
169 1220 : cfCells_p(i,j,k)->ySupport_p = 0;
170 1220 : k++;
171 : }
172 : }
173 196 : }
174 : //
175 : //---------------------------------------------------------------
176 : //
177 : // template <class T> void CFBuffer<T>
178 : // RigidVector<Int, 3> CFBuffer::setParams(const Int& inu, const Int& iw, const Int& muellerElement,
179 : // const TableRecord& miscInfo)
180 : // {
181 : // RigidVector<Int,3> ndx = setParams(inu, iw,
182 : // miscInfo.coordSys, miscInfo.sampling,
183 : // miscInfo.xSupport,miscInfo.ySupport,
184 : // miscInfo.freqValue, miscInfo.wValue,
185 : // muellerElements,
186 : // miscInfo.fileName, miscInfo.conjfFreq,
187 : // miscinfo.conjPoln);
188 : // cfCells_p(ndx(0), ndx(1), ndx(2))->telescopeName = micsInfo.telescopeName;
189 : // }
190 : //
191 : //---------------------------------------------------------------
192 : //
193 : // template <class T> void CFBuffer<T>
194 1020 : RigidVector<Int, 3> CFBuffer::setParams(const Int& inu, const Int& iw, const Int& ipx, const Int& ipy,
195 : const Double& freqValue,
196 : const Double& wValue,
197 : const Int& muellerElement,
198 : CoordinateSystem& cs,
199 : const TableRecord& miscInfo)
200 : {
201 1020 : float sampling; miscInfo.get("Sampling",sampling);
202 1020 : int xSupport, ySupport; miscInfo.get("Xsupport",xSupport);miscInfo.get("Ysupport",ySupport);
203 2040 : String fileName; miscInfo.get("Name",fileName);
204 1020 : double conjFreq; miscInfo.get("ConjFreq", conjFreq);
205 1020 : int conjPoln; miscInfo.get("ConjPoln", conjPoln);
206 2040 : String telescopeName; miscInfo.get("TelescopeName", telescopeName);
207 2040 : String bandName;
208 1020 : if (miscInfo.isDefined("BandName")) miscInfo.get("BandName", bandName);
209 0 : else bandName="UNKNOWABLE";
210 1020 : float diameter; miscInfo.get("Diameter", diameter);
211 : // In the absense of evidence, assume that users are sensible and
212 : // are using AWProjection where it is really need it and not for
213 : // using it as a replacement for rotatially symmetric stuff. So
214 : // by default, the CFs are assumed to be rotationally asymmetric.
215 1020 : bool isRotationallySymmetric=True;
216 :
217 1020 : if (miscInfo.isDefined("OpCode"))
218 1020 : miscInfo.get("OpCode",isRotationallySymmetric);
219 :
220 : RigidVector<Int,3> ndx=setParams(inu, iw, ipx, ipy, freqValue, bandName, wValue, muellerElement, cs,
221 : sampling, xSupport, ySupport, fileName, conjFreq, conjPoln, telescopeName,
222 1020 : diameter);
223 1020 : cfCells_p(ndx(0),ndx(1),ndx(2))->isRotationallySymmetric_p = isRotationallySymmetric;
224 2040 : return ndx;
225 : }
226 1420 : RigidVector<Int, 3> CFBuffer::setParams(const Int& inu, const Int& iw, const Int& /*ipx*/, const Int& /*ipy*/,
227 : const Double& freqValue,
228 : const String& bandName,
229 : const Double& wValue,
230 : const Int& muellerElement,
231 : CoordinateSystem& cs,
232 : Float& sampling,
233 : Int& xSupport, Int& ySupport,
234 : const String& fileName,
235 : const Double& conjFreq,
236 : const Int& conjPoln,
237 : const String& telescopeName,
238 : const Float& diameter)
239 : {
240 1420 : RigidVector<Int,3> ndx=getIndex(freqValue, wValue, muellerElement);
241 1420 : ndx(0)=inu; ndx(1)=iw;//ndx(2) = muellerElements_p(ipx)(ipy);
242 1420 : cfCells_p(ndx(0),ndx(1),ndx(2))->sampling_p = sampling;
243 1420 : cfCells_p(ndx(0),ndx(1),ndx(2))->xSupport_p = xSupport;
244 1420 : cfCells_p(ndx(0),ndx(1),ndx(2))->ySupport_p = ySupport;
245 1420 : cfCells_p(ndx(0),ndx(1),ndx(2))->muellerElement_p = muellerElement;
246 1420 : cfCells_p(ndx(0),ndx(1),ndx(2))->conjFreq_p = conjFreq;
247 1420 : cfCells_p(ndx(0),ndx(1),ndx(2))->conjPoln_p = conjPoln;
248 1420 : if (fileName != "")
249 1020 : cfCells_p(ndx(0),ndx(1),ndx(2))->fileName_p = fileName;
250 1420 : cfCells_p(ndx(0),ndx(1),ndx(2))->telescopeName_p = telescopeName;//String("EVLA");
251 1420 : cfCells_p(ndx(0),ndx(1),ndx(2))->diameter_p = diameter;//25.0;
252 1420 : if (bandName != "") cfCells_p(ndx(0),ndx(1),ndx(2))->bandName_p = bandName;
253 :
254 1420 : Int index=cs.findCoordinate(Coordinate::SPECTRAL);
255 2840 : SpectralCoordinate spCS = cs.spectralCoordinate(index);
256 2840 : Vector<Double> val; val.resize(1);val(0)=freqValues_p(ndx(0));
257 1420 : spCS.setReferenceValue(val);
258 : //val.resize();
259 : //val = spCS.increment();
260 : // val = cfCells_p(ndx(0),ndx(1),ndx(2))->freqIncr_p;
261 : //spCS.setIncrement(val);
262 1420 : cs.replaceCoordinate(spCS,index);
263 :
264 : // cfCells_p(ndx(0),ndx(1),ndx(2))->freqValue_p =
265 : // cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL)).referenceValue()(0);
266 : // cfCells_p(ndx(0),ndx(1),ndx(2))->freqIncr_p =
267 : // cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL)).increment()(0);
268 :
269 1420 : cfCells_p(ndx(0),ndx(1),ndx(2))->coordSys_p = cs;
270 2840 : return ndx;
271 : }
272 : //
273 : //---------------------------------------------------------------
274 : //
275 170 : void CFBuffer::setPA(Float& pa)
276 : {
277 170 : RigidVector<Int,3> ndx(-1);
278 170 : Int nF=cfCells_p.shape()(0), nW=cfCells_p.shape()(1), nM=cfCells_p.shape()(2);
279 680 : for (Int i=0; i<nF; i++)
280 1020 : for(Int j=0; j<nW; j++)
281 1530 : for(Int k=0; k<nM; k++)
282 1020 : cfCells_p(i,j,k)->pa_p=Quantity(pa,"deg");
283 170 : }
284 :
285 0 : void CFBuffer::setParams(Int& /*nx*/, Int& /*ny*/, CoordinateSystem& /*cs*/, Float& /*sampling*/,
286 : Int& /*xSupport*/, Int& /*ySupport*/, const Double& /*freqValue*/,
287 : const Double& /*wValue*/, const Int& /*muellerElement*/,
288 : const String& /*fileName*/)
289 : {
290 : /*
291 : RigidVector<Int,3> ndx=setParams(cs, sampling, xSupport, ySupport,
292 : freqValue, wValue, muellerElement);
293 : // Adding degenerate axis here to make life easier when using these
294 : // buffers to write them as Images later (with a 4D co-ordinate
295 : // system - the degenerate axis become the Freq. and Polarization
296 : // axis).
297 :
298 : // cfCells_p(ndx(0),ndx(1),ndx(2))->storage_p->resize(IPosition(4,nx,ny,1,1));
299 : cfCells_p(ndx(0),ndx(1),ndx(2))->shape_p=IPosition(4,nx,ny,1,1);
300 : */
301 0 : }
302 : //
303 : //---------------------------------------------------------------
304 : //
305 : // template<class T> Array<T>& CFBuffer<T>::
306 0 : CFCell& CFBuffer::getCFCell(const Int& i, const Int& j, const Int& k)
307 : {
308 0 : return *cfCells_p(i,j,k); // Nfreq x Nw x Nmueller
309 : }
310 : //
311 : //---------------------------------------------------------------
312 : //
313 : // template <class T> Array<T>& CFBuffer<T>
314 0 : CFCell& CFBuffer::getCFCell(const Double& freqVal,
315 : const Double& wValue,
316 : const Int& muellerElement) // muellerElement: (i,j) of the Mueller Matrix
317 : {
318 0 : RigidVector<Int,3> ndx=getIndex(freqVal, wValue, muellerElement);
319 0 : return getCFCell(ndx(0), ndx(1),ndx(2));
320 : }
321 : //
322 : //---------------------------------------------------------------
323 : //
324 : // template<class T> Array<T>& CFBuffer<T>::
325 7292352 : CountedPtr<CFCell>& CFBuffer::getCFCellPtr(const Int& i/*FreqNdx*/, const Int& j/*WNdx*/, const Int& k/*MuellerNdx*/)
326 : {
327 : // IPosition shp=cfCells_p.shape();
328 : // if (cfHitsStats.shape().product()==0)
329 : // {
330 : // cfHitsStats.resize(shp);
331 : // cfHitsStats = 0;
332 : // }
333 : // try
334 : // {
335 : // AlwaysAssert(((i<shp[0]) && (j<shp[1]) && (k<shp[2])) , AipsError);
336 : // }
337 : // catch (AipsError x)
338 : // {
339 : // cerr << "#### " << i << " " << j << " " << k << " " << shp << endl;
340 : // throw(x);
341 : // }
342 :
343 : // cfHitsStats(i,j,k)++;
344 : // // cerr << "CFBuffer: Cell: " << i << " " << j << " " << k << " " << cfCells_p(i,j,k)->xSupport_p << endl;
345 7292352 : return cfCells_p(i,j,k); // Nfreq x Nw x Nmueller
346 : }
347 : //
348 : //---------------------------------------------------------------
349 : //
350 : // template <class T> Array<T>& CFBuffer<T>
351 1068 : CountedPtr<CFCell>& CFBuffer::getCFCellPtr(const Double& freqVal,
352 : const Double& wValue,
353 : const Int & muellerElement)
354 : {
355 1068 : RigidVector<Int,3> ndx=getIndex(freqVal, wValue, muellerElement);
356 1068 : return cfCells_p(ndx(0), ndx(1),ndx(2));
357 : }
358 : //
359 : //---------------------------------------------------------------
360 : //
361 : // template<class T> Vector<Int>& CFBuffer<T>
362 : // This should be made more efficient with binary search.
363 2700 : RigidVector<Int,3> CFBuffer::getIndex(const Double& freqValue,
364 : const Double& wValue,
365 : const Int& muellerElement)
366 : {
367 2700 : RigidVector<Int,3> ndx(-1);
368 2700 : Int nF=cfCells_p.shape()(0), nW=cfCells_p.shape()(1), nM=cfCells_p.shape()(2);
369 : Int i,j,k;
370 : //UNUSED: Int di;
371 :
372 : // Double dfMin=0;di=0;
373 : // for (i=0;i<nF;i++)
374 : // {
375 : // Double df=fabs(cfCells_p(i,0,0)->freqValue_p - freqValue);
376 : // if (df < dfMin) {dfMin=df; di=i;}
377 : // }
378 : // i=di;
379 5052 : for (i=0; i<nF; i++)
380 : // if ((fabs(cfCells_p(i,0,0)->freqValue_p - freqValue) < cfCells_p(i,0,0)->freqIncr_p))
381 5052 : if (cfCells_p(i,0,0)->freqValue_p == freqValue) break;
382 5100 : for(j=0; j<nW; j++) if (cfCells_p(i,j,0)->wValue_p == wValue) break;
383 4050 : for (k=0; k<nM; k++) if (cfCells_p(i,j,k)->muellerElement_p == muellerElement) break;
384 :
385 2700 : ndx(0)=i;ndx(1)=j;ndx(2)=k;
386 : // for (Int i=0;i<nF;i++)
387 : // for(Int j=0;j<nW;j++)
388 : // for(Int k=0;k<nM;k++)
389 : // {
390 : // if (
391 : // (fabs(cfCells_p(i,j,k)->freqValue_p - freqValue) < cfCells_p(i,j,k)->freqIncr_p) &&
392 : // (cfCells_p(i,j,k)->wValue_p == wValue) &&
393 : // (cfCells_p(i,j,k)->muellerElement_p == muellerElement)
394 : // )
395 : // {ndx(0)=i;ndx(1)=j;ndx(2)=k;break;}
396 : // }
397 :
398 : // cerr << "@#$%#@ " << freqValue
399 : // << " " << cfCells_p(ndx(0),ndx(1),ndx(2))->freqValue_p
400 : // << " " << cfCells_p(ndx(0),ndx(1),ndx(2))->freqIncr_p
401 : // << " " << ndx
402 : // << endl;
403 : // for(uInt i=0;i<nF;i++)
404 : // if (freqValue==cfCells_pfreqValues_p(i)) {ndx(0)=i;break;}
405 :
406 : // for(uInt i=0;i<wValues_p.nelements();i++)
407 : // if (wValue==wValues_p(i)) {ndx(1)=i;break;}
408 :
409 : // ndx(2)=muellerElements_p(muellerElement(0),muellerElement(1));
410 :
411 2700 : return ndx;
412 : }
413 : //
414 : //---------------------------------------------------------------
415 : //
416 212 : void CFBuffer::getParams(CoordinateSystem& cs, Float& sampling,
417 : Int& xSupport, Int& ySupport,
418 : String& bandName,
419 : const Double& freqVal, const Double& wValue,
420 : const Int& muellerElement)
421 : {
422 212 : RigidVector<Int,3> ndx=getIndex(freqVal, wValue, muellerElement);
423 212 : getParams(cs, sampling, xSupport, ySupport, bandName, ndx(0), ndx(1), ndx(2));
424 212 : }
425 : //
426 : //---------------------------------------------------------------
427 : //
428 0 : Double CFBuffer::nearest(Bool& found, const Double& val,
429 : const Vector<Double>& valList,
430 : const Double& /*incr*/)
431 : {
432 0 : Int n=valList.nelements();
433 0 : if (n==1) {found=true;return valList[0];}
434 0 : Int where = binarySearch(found, valList, val, n);
435 0 : if (found) return valList(where);
436 0 : else return -1.0;
437 : }
438 : //
439 : //---------------------------------------------------------------
440 : //
441 0 : Int CFBuffer::nearestNdx(const Double& val,
442 : const Vector<Double>& /*valList*/,
443 : const Double& incr)
444 : {
445 0 : return SynthesisUtils::nint(incr*val);
446 :
447 : // Int n=valList.nelements();
448 : // if (n==1) return 0;
449 : // else return -1;
450 :
451 : // Int where = binarySearch(found, valList, val, n);
452 : // if (found) return valList(where);
453 : // else return -1.0;
454 : }
455 : //
456 : //---------------------------------------------------------------
457 : //
458 : // template <class T> void CFBuffer<T>
459 0 : void CFBuffer::show(const char *Mesg,ostream &os)
460 : {
461 0 : LogIO log_l(LogOrigin("CFBuffer","show[R&D]"));
462 :
463 0 : if (Mesg != NULL) os << Mesg << endl;
464 0 : os << "---------------------------------------------------------" << endl
465 0 : << "Shapes: " << cfCells_p.shape() << endl;
466 0 : for (Int i=0;i<cfCells_p.shape()(0);i++)
467 0 : for (Int j=0;j<cfCells_p.shape()(1);j++)
468 0 : for (Int k=0;k<cfCells_p.shape()(2);k++)
469 : {
470 0 : os << "CFCell["<<i<<","<<j<<","<<k<<"]:" << endl;
471 0 : cfCells_p(i,j,k)->show(Mesg, os);
472 0 : os << "Pointing offset: " << pointingOffset_p << endl;
473 : }
474 0 : os << "---------------------------------------------------------" << endl;
475 0 : }
476 :
477 38 : void CFBuffer::makePersistent(const char *dir, const char *cfName)
478 : {
479 144 : for (Int i=0;i<cfCells_p.shape()(0);i++)
480 242 : for (Int j=0;j<cfCells_p.shape()(1);j++)
481 408 : for (Int k=0;k<cfCells_p.shape()(2);k++)
482 : {
483 272 : ostringstream name;
484 272 : name << String(cfName) << "_CF_" << i << "_" << j << "_" << k << ".im";
485 : //cerr << "CFB Name : " << name.str() << endl;
486 : // const char *formedName;
487 : // if (cfName != "" ) formedName = name.str().c_str();
488 : // else formedName = cfName;
489 : // cerr << "Formed name = " << formedName << endl;
490 272 : cfCells_p(i,j,k)->makePersistent(dir, name.str().c_str());
491 : }
492 38 : }
493 :
494 988 : Int CFBuffer::nearestFreqNdx(const Double& freqVal)
495 : {
496 : Int index;
497 988 : SynthesisUtils::nearestValue(freqValues_p, freqVal, index);
498 988 : return index;
499 : // // The algorithm below has a N*log(N) cost.
500 : // Vector<Double> diff = fabs(freqValues_p - freqVal);
501 : // Bool dummy;
502 : // Sort sorter(diff.getStorage(dummy), sizeof(Double));
503 : // sorter.sortKey((uInt)0,TpDouble);
504 : // Int nch=freqValues_p.nelements();
505 : // Vector<uInt> sortIndx;
506 : // sorter.sort(sortIndx, nch);
507 :
508 : // return sortIndx(0);
509 :
510 : // Int ndx=min(freqValues_p.nelements()-1,max(0,SynthesisUtils::nint((freqVal-freqValues_p[0])/freqValIncr_p)));
511 : // return ndx;
512 : }
513 :
514 :
515 170 : void CFBuffer::primeTheCache()
516 : {
517 : //
518 : // In each CFCell of this CFBuffer, cache things that might be
519 : // required in tight loops and can be expensive to extract
520 : // otherwise.
521 : //
522 340 : IPosition cfCShape=cfCells_p.shape();
523 680 : for (Int i=0; i < cfCShape(0); i++)
524 1020 : for (Int j=0; j < cfCShape(1); j++)
525 1530 : for (Int k=0; k < cfCShape(2); k++)
526 1020 : getCFCellPtr(i,j,k)->initCache();
527 170 : }
528 :
529 170 : void CFBuffer::initMaps(const VisBuffer2&, // vb,
530 : const Matrix<Double>& freqSelection, const Double& imRefFreq)
531 : {
532 340 : Vector<Double> spwList=freqSelection.column(0);
533 170 : Int maxSpw=(Int)(max(spwList));
534 170 : freqNdxMap_p.resize(maxSpw+1);
535 170 : conjFreqNdxMap_p.resize(maxSpw+1);
536 :
537 664 : for (Int i=0;i<(Int)spwList.nelements(); i++)
538 : {
539 494 : Int spw=(Int)freqSelection(i,0);
540 494 : Double fmin=freqSelection(i,1), fmax=freqSelection(i,2), finc=freqSelection(i,3);
541 494 : Int nchan = (Int)((fmax-fmin)/finc + 1);
542 494 : freqNdxMap_p[spw].resize(nchan,True);
543 494 : conjFreqNdxMap_p[spw].resize(nchan,True);
544 988 : for (Int c=0;c<nchan;c++)
545 : {
546 494 : Double freq=fmin+c*finc;
547 494 : Double conjFreq=sqrt(2*imRefFreq*imRefFreq - freq*freq);
548 494 : freqNdxMap_p[spw][c]=nearestFreqNdx(freq);
549 494 : conjFreqNdxMap_p[spw][c]=nearestFreqNdx(conjFreq);
550 : }
551 : }
552 :
553 :
554 : // cerr << "CFBuffer::initMaps: "
555 : // << freqSelection << endl
556 : // << freqValues_p << endl
557 : // << freqNdxMap_p << endl
558 : // << conjFreqNdxMap_p << endl;
559 :
560 170 : }
561 :
562 144 : void CFBuffer::initPolMaps(PolMapType& polMap, PolMapType& conjPolMap)
563 : {
564 144 : muellerElementsIndex_p.assign(polMap);
565 144 : conjMuellerElementsIndex_p.assign(conjPolMap);
566 144 : }
567 :
568 0 : void CFBuffer::getFreqNdxMaps(Vector<Vector<Int> >& freqNdx, Vector<Vector<Int> >& conjFreqNdx)
569 : {
570 : Int nspw;
571 0 : nspw=freqNdxMap_p.nelements();
572 0 : freqNdx.resize(nspw);
573 0 : for (Int s=0;s<nspw;s++)
574 0 : freqNdx[s].assign(freqNdxMap_p[s]);
575 :
576 0 : nspw=conjFreqNdxMap_p.nelements();
577 0 : conjFreqNdx.resize(nspw);
578 0 : for (Int s=0;s<nspw;s++)
579 0 : conjFreqNdx[s].assign(conjFreqNdxMap_p[s]);
580 0 : }
581 :
582 0 : void CFBuffer::ASSIGNVVofI(Int** &target,Vector<Vector<Int> >& source, Bool& doAlloc)
583 : {
584 : // cerr << "Source: " << target << " " << source << endl;
585 :
586 : Int nx,ny;
587 0 : Bool b=(doAlloc||(target==NULL));
588 0 : nx=source.nelements();
589 0 : if (b) target=(int **)malloc(nx*sizeof(int*));
590 0 : for (int i=0;i<nx;i++)
591 : {
592 0 : ny = source[i].nelements();
593 0 : if (b) (target[i]) = (int *)malloc(ny*sizeof(int));
594 0 : for (int j=0;j<ny;j++)
595 0 : (target)[i][j]=source[i][j];
596 : }
597 :
598 : // cerr << "Target: ";
599 : // for (int ii=0;ii<source.nelements();ii++)
600 : // {
601 : // Int ny=source[ii].nelements();
602 : // for(Int jj=0;jj<ny;jj++)
603 : // cerr << target[ii][jj] << " ";
604 : // cerr << endl;
605 : // }
606 0 : }
607 :
608 0 : void CFBuffer::getAsStruct(CFBStruct& st)
609 : {
610 0 : Vector<Int> shp=cfCells_p.shape().asVector();
611 :
612 0 : Bool doAlloc=(st.CFBStorage == NULL);
613 :
614 0 : st.shape[0]=shp[0];
615 0 : st.shape[1]=shp[1];
616 0 : st.shape[2]=shp[2];
617 0 : st.fIncr = freqValIncr_p; st.wIncr = wValIncr_p;
618 :
619 :
620 : Bool dummy;
621 0 : Int n=cfCells_p.nelements();
622 0 : if (doAlloc) st.CFBStorage = (CFCStruct *)malloc(n*sizeof(CFCStruct));
623 0 : CountedPtr<CFCell> *cfstore=cfCells_p.getStorage(dummy);
624 0 : for (int i=0;i<n;i++)
625 : {
626 : // st.CFBStorage[i]=(CFCStruct *)malloc(sizeof(CFCStruct));
627 0 : (cfstore[i]).operator->()->getAsStruct(st.CFBStorage[i]);
628 : }
629 : // st.CFBStorage[i] = (cfstore[i]).operator->()->getStorage()->getStorage(dummy);
630 :
631 : // if (doAlloc) st.pointingOffset=(Double *)malloc(pointingOffset_p.nelements()*sizeof(Double));
632 : // for (uInt i=0;i<pointingOffset_p.nelements();i++) st.pointingOffset[i]=pointingOffset_p[i];
633 :
634 0 : if (doAlloc) st.freqValues=(Double *)malloc(freqValues_p.nelements()*sizeof(Double));
635 0 : for (uInt i=0;i<freqValues_p.nelements();i++) st.freqValues[i]=freqValues_p[i];
636 :
637 0 : if (doAlloc) st.wValues=(Double *)malloc(wValues_p.nelements()*sizeof(Double));
638 0 : for (uInt i=0;i<wValues_p.nelements();i++) st.wValues[i]=wValues_p[i];
639 :
640 0 : ASSIGNVVofI(st.muellerElements, muellerElements_p, doAlloc);
641 0 : ASSIGNVVofI(st.muellerElementsIndex, muellerElementsIndex_p,doAlloc);
642 0 : ASSIGNVVofI(st.conjMuellerElements, conjMuellerElements_p,doAlloc);
643 0 : ASSIGNVVofI(st.conjMuellerElementsIndex, conjMuellerElementsIndex_p,doAlloc);
644 :
645 0 : st.nMueller = muellerElements_p.nelements();
646 :
647 0 : Vector<Vector<Int> > freqNdx, conjFreqNdx;
648 0 : getFreqNdxMaps(freqNdx, conjFreqNdx);
649 0 : ASSIGNVVofI(st.freqNdxMap, freqNdx, doAlloc);
650 0 : ASSIGNVVofI(st.conjFreqNdxMap, conjFreqNdx, doAlloc);
651 0 : }
652 :
653 : //
654 : //----------------------------------------------------------------------
655 : //
656 0 : void CFBuffer::fill(const Int& /*nx*/, const Int& /*ny*/,
657 : const Vector<Double>& freqValues,
658 : const Vector<Double>& wValues,
659 : const PolMapType& muellerElements)
660 : {
661 :
662 0 : LogIO log_l(LogOrigin("CFBuffer", "fillBuffer[R&D]"));
663 0 : for (uInt imx=0;imx<muellerElements.nelements();imx++) // Loop over all MuellerElements
664 0 : for (uInt imy=0;imy<muellerElements(imx).nelements();imy++)
665 : {
666 0 : for (uInt inu=0;inu<freqValues.nelements();inu++) // All freq. channels
667 : {
668 0 : for (uInt iw=0;iw<wValues.nelements();iw++) // All w-planes
669 : {
670 : log_l << " CF("
671 0 : << "M:"<<muellerElements(imx)(imy)
672 : << ",C:" << inu
673 0 : << ",W:" << iw << "): ";
674 0 : Vector<Double> ftRef(2);
675 :
676 : // ftRef(0)=cfWtBuf.shape()(0)/2.0;
677 : // ftRef(1)=cfWtBuf.shape()(1)/2.0;
678 0 : CoordinateSystem ftCoords;//=cs_l;
679 : // SynthesisUtils::makeFTCoordSys(cs_l, cfWtBuf.shape()(0), ftRef, ftCoords);
680 :
681 : // cfb.setParams(inu,iw,imx,imy,//muellerElements(imx)(imy),
682 : // ftCoords, sampling, xSupport, ySupport,
683 : // freqValues(inu), wValues(iw), muellerElements(imx)(imy));
684 : // cfb.getCFCellPtr(freqValues(inu), wValues(iw),
685 : // muellerElements(imx)(imy))->pa_p=Quantity(vbPA,"rad");
686 : //
687 : // Now tha the CFs have been computed, cache its
688 : // paramters in CFCell for quick access in tight
689 : // loops (in the re-sampler, e.g.).
690 : //
691 0 : (getCFCellPtr(freqValues(inu), wValues(iw),
692 0 : muellerElements(imx)(imy)))->initCache();
693 : }
694 : }
695 : }
696 0 : }
697 : //
698 : //----------------------------------------------------------------------
699 : //
700 191 : int CFBuffer::getMaxCFSize()
701 : {
702 191 : if (maxCFSize_p < 0)
703 : {
704 210 : IPosition shp(cfCells_p.shape());
705 412 : for (Int i=0;i < shp[0]; i++)
706 644 : for (Int j=0; j < shp[1]; j++)
707 1011 : for (Int k=0; k < shp[2]; k++)
708 : {
709 674 : int cfNX = SynthesisUtils::getCFShape(getCFCacheDir(), cfCells_p(i,j,k)->fileName_p)[0];
710 674 : if (cfNX > maxCFSize_p) maxCFSize_p=cfNX;
711 : }
712 : }
713 191 : return maxCFSize_p;
714 : }
715 : //
716 : //----------------------------------------------------------------------
717 : //
718 33789 : bool CFBuffer::finitePointingOffsets()
719 : {
720 67416 : return ((fabs(pointingOffset_p(0)(0))>0) ||
721 67416 : (fabs(pointingOffset_p(0)(1))>0));
722 :
723 : }
724 :
725 : } // end casa namespace
726 :
727 :
728 :
|
