Line data Source code
1 : // -*- C++ -*-
2 : //# AWProjectWBFT.cc: Implementation of AWProjectWBFT 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 :
29 : #include <synthesis/TransformMachines/AWProjectWBFT.h>
30 : #include <synthesis/TransformMachines/AWVisResampler.h>
31 : #include <synthesis/TransformMachines/StokesImageUtil.h>
32 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
33 : #include <casacore/scimath/Mathematics/FFTServer.h>
34 : #include <casacore/scimath/Mathematics/Convolver.h>
35 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
36 : #include <casacore/images/Images/ImageInterface.h>
37 : #include <casacore/images/Images/PagedImage.h>
38 : #include <msvis/MSVis/VisBuffer.h>
39 : #include <casacore/casa/Arrays/Vector.h>
40 : #include <casacore/casa/Arrays/Slice.h>
41 : #include <casacore/casa/Arrays/ArrayMath.h>
42 : #include <casacore/casa/Arrays/Array.h>
43 : #include <casacore/casa/OS/HostInfo.h>
44 : #include <sstream>
45 : #include <casacore/casa/Utilities/CompositeNumber.h>
46 :
47 : //#define CONVSIZE (1024*4)
48 : //#define OVERSAMPLING 10
49 : #define USETABLES 0 // If equal to 1, use tabulated exp() and
50 : // complex exp() functions.
51 : #define MAXPOINTINGERROR 250.0 // Max. pointing error in arcsec used to
52 : // determine the resolution of the
53 : // tabulated exp() function.
54 : using namespace casacore;
55 : namespace casa { //# NAMESPACE CASA - BEGIN
56 : //
57 : //---------------------------------------------------------------
58 : //
59 :
60 0 : AWProjectWBFT::AWProjectWBFT(Int nWPlanes, Long icachesize,
61 : CountedPtr<CFCache>& cfcache,
62 : CountedPtr<ConvolutionFunction>& cf,
63 : CountedPtr<VisibilityResamplerBase>& visResampler,
64 : Bool applyPointingOffset,
65 : Bool doPBCorr,
66 : Int itilesize,
67 : Float paSteps,
68 : Float pbLimit,
69 : Bool usezero,
70 : Bool conjBeams,
71 0 : Bool doublePrecGrid)
72 : : AWProjectFT(nWPlanes,icachesize,cfcache,cf,visResampler,applyPointingOffset,doPBCorr, itilesize,pbLimit,usezero,conjBeams,doublePrecGrid),
73 0 : avgPBReady_p(false),resetPBs_p(true),wtImageFTDone_p(false),fieldIds_p(0),rotatedCFWts_p(),visResamplerWt_p(),oneTimeMessage_p(false)
74 : {
75 : (void)paSteps;
76 : //
77 : // Set the function pointer for FORTRAN call for GCF services.
78 : // This is a pure C function pointer which FORTRAN can call.
79 : // The function itself uses GCF object services.
80 : //
81 : // convSize=0;
82 0 : paChangeDetector.reset();
83 0 : pbLimit_p=pbLimit;
84 0 : if (applyPointingOffset) doPointing=1; else doPointing=0;
85 0 : maxConvSupport=-1;
86 : //
87 : // Set up the Conv. Func. disk cache manager object.
88 : //
89 : //convSampling=OVERSAMPLING;
90 : // convSize=CONVSIZE;
91 : //use memory size defined in aipsrc if exists
92 0 : Long hostRAM = (HostInfo::memoryTotal(true)*1024); // In bytes
93 0 : hostRAM = hostRAM/(sizeof(Float)*2); // In complex pixels
94 0 : if (cachesize > hostRAM) cachesize=hostRAM;
95 :
96 : // visResampler_p->init(useDoubleGrid_p);
97 0 : lastPAUsedForWtImg = MAGICPAVALUE;
98 0 : }
99 : //
100 : //---------------------------------------------------------------
101 : //
102 0 : AWProjectWBFT::AWProjectWBFT(const RecordInterface& stateRec)
103 0 : : AWProjectFT(stateRec), oneTimeMessage_p(false)
104 : {
105 0 : LogIO log_l(LogOrigin("AWProjectWBFT", "AWProjectWBFT[R&D]"));
106 : // Construct from the input state record
107 :
108 : // if (!fromRecord(error, stateRec))
109 0 : if (!fromRecord(stateRec))
110 0 : log_l << "Failed to create " << name() << " object." << LogIO::EXCEPTION;
111 :
112 0 : maxConvSupport=-1;
113 : //convSampling=OVERSAMPLING;
114 : //convSize=CONVSIZE;
115 0 : visResampler_p->init(useDoubleGrid_p);
116 0 : }
117 : //
118 : //---------------------------------------------------------------
119 : //
120 0 : AWProjectWBFT& AWProjectWBFT::operator=(const AWProjectWBFT& other)
121 : {
122 0 : if(this!=&other)
123 : {
124 0 : AWProjectFT::operator=(other);
125 0 : padding_p = other.padding_p;
126 0 : nWPlanes_p = other.nWPlanes_p;
127 0 : imageCache = other.imageCache;
128 0 : cachesize = other.cachesize;
129 0 : tilesize = other.tilesize;
130 0 : gridder = other.gridder;
131 0 : isTiled = other.isTiled;
132 0 : lattice = other.lattice;
133 0 : arrayLattice = other.arrayLattice;
134 0 : maxAbsData = other.maxAbsData;
135 0 : centerLoc = other.centerLoc;
136 0 : offsetLoc = other.offsetLoc;
137 0 : pointingToImage = other.pointingToImage;
138 0 : usezero_p = other.usezero_p;
139 0 : doPBCorrection = other.doPBCorrection;
140 0 : maxConvSupport = other.maxConvSupport;
141 0 : avgPBReady_p = other.avgPBReady_p;
142 0 : resetPBs_p = other.resetPBs_p;
143 0 : wtImageFTDone_p = other.wtImageFTDone_p;
144 0 : rotatedCFWts_p = other.rotatedCFWts_p;
145 0 : visResamplerWt_p=other.visResamplerWt_p; // Copy the counted pointer
146 0 : oneTimeMessage_p = other.oneTimeMessage_p;
147 : // *visResamplerWt_p = *other.visResamplerWt_p; // Call the appropriate operator=()
148 :
149 : // visResampler_p=other.visResampler_p->clone();
150 : };
151 0 : return *this;
152 : };
153 : //
154 : //----------------------------------------------------------------------
155 : //
156 0 : Int AWProjectWBFT::findPointingOffsets(const VisBuffer& vb,
157 : Array<Float> &l_off,
158 : Array<Float> &m_off,
159 : Bool Evaluate)
160 : {
161 0 : LogIO log_l(LogOrigin("AWProjectWBFT","findPointingOffsets[R&D]"));
162 0 : Int NAnt=0;
163 : //
164 : // This will return 0 if EPJ Table is not given. Otherwise will
165 : // return the number of antennas it detected (from the EPJ table)
166 : // and the offsets in l_off and m_off.
167 : //
168 0 : NAnt = AWProjectFT::findPointingOffsets(vb,l_off,m_off,Evaluate);
169 : // NAnt = l_off.shape()(2);
170 :
171 : // Resize the offset arrays if no pointing table was given.
172 : //
173 0 : if (NAnt <=0 )
174 : {
175 0 : NAnt=vb.msColumns().antenna().nrow();
176 0 : l_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
177 0 : m_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
178 0 : l_off = m_off = 0.0;
179 : }
180 : //
181 : // Add field offsets to the pointing errors.
182 : //
183 :
184 : // Float dayHr=2*3.141592653589793116;
185 : // MVDirection ref(directionCoord.referenceValue()(0),
186 : // directionCoord.referenceValue()(1)),
187 : // vbDir(vb.direction1()(0).getAngle().getValue()(0),
188 : // vb.direction1()(0).getAngle().getValue()(1));
189 :
190 : // if (0)
191 : // {
192 : // l_off = l_off - (Float)(directionCoord.referenceValue()(0) -
193 : // dayHr-vb.direction1()(0).getAngle().getValue()(0));
194 : // m_off = m_off - (Float)(directionCoord.referenceValue()(1) -
195 : // vb.direction1()(0).getAngle().getValue()(1));
196 : // }
197 :
198 : //
199 : // Convert the direction from image co-ordinate system and the VB
200 : // to MDirection. Then convert the MDirection to Quantity so that
201 : // arithematic operation (subtraction) can be done. Then use the
202 : // subtracted Quantity to construct another MDirection and use
203 : // *it's* getAngle().getValue() to extract the difference in the
204 : // sky direction (from image co-ordinate system) and the phase
205 : // center direction of the VB in radians!
206 : //
207 : // If only VisBuffer, and DirectionCoordinate class could return
208 : // MDirection, and MDirection class had operator-(), the code
209 : // below could look more readable as:
210 : // MDirection diff=vb.mDirection()-directionCoord.mDirection();
211 : //
212 0 : CoordinateSystem coords(image->coordinates());
213 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
214 0 : AlwaysAssert(directionIndex>=0, AipsError);
215 0 : DirectionCoordinate directionCoord=coords.directionCoordinate(directionIndex);
216 0 : MDirection vbMDir(vb.direction1()(0)),skyMDir, diff;
217 0 : directionCoord.toWorld(skyMDir, directionCoord.referencePixel());
218 :
219 0 : diff = MDirection(skyMDir.getAngle()-vbMDir.getAngle());
220 0 : l_off = l_off - (Float)diff.getAngle().getValue()(0);
221 0 : m_off = m_off - (Float)diff.getAngle().getValue()(1);
222 :
223 : static int firstPass=0;
224 : // static int fieldID=-1;
225 0 : static Vector<Float> offsets0,offsets1;
226 0 : if (firstPass==0)
227 : {
228 0 : offsets0.resize(NAnt);
229 0 : offsets1.resize(NAnt);
230 : // MLCG mlcg((Int)(vb.time()(0)));
231 : // Normal nrand(&mlcg,0.0,10.0);
232 : // for(Int i=0;i<NAnt;i++) offsets0(i) = (Float)(nrand());
233 : // for(Int i=0;i<NAnt;i++) offsets1(i) = (Float)(nrand());
234 0 : offsets0 = offsets1 = 0.0;
235 : }
236 0 : for(Int i=0;i<NAnt;i++)
237 : {
238 0 : l_off(IPosition(3,0,0,i)) = l_off(IPosition(3,0,0,i)) + offsets0(i)/2.062642e+05;
239 0 : m_off(IPosition(3,0,0,i)) = m_off(IPosition(3,0,0,i)) + offsets1(i)/2.062642e+05;
240 : }
241 :
242 : //m_off=l_off=0.0;
243 : // if (fieldID != vb.fieldId())
244 : // {
245 : // fieldID = vb.fieldId();
246 : // cout << l_off*2.062642e5 << endl;
247 : // cout << m_off*2.062642e5 << endl;
248 : // }
249 0 : if (firstPass==0)
250 : {
251 : // cout << (Float)(directionCoord.referenceValue()(0)) << " "
252 : // << (Float)(directionCoord.referenceValue()(1)) << endl;
253 : // cout << vb.direction1()(0).getAngle().getValue()(0) << " "
254 : // << vb.direction1()(0).getAngle().getValue()(1) << endl;
255 : // cout << l_off << endl;
256 : // cout << m_off << endl;
257 : }
258 0 : firstPass++;
259 0 : return NAnt;
260 0 : }
261 : //
262 : //---------------------------------------------------------------
263 : //
264 0 : void AWProjectWBFT::normalizeAvgPB()
265 : {
266 0 : LogIO log_l(LogOrigin("AWProjectWBFT","normalizeAvgPB[R&D]"));
267 : // We accumulated normalized PBs. So don't normalize the average
268 : // PB.
269 0 : pbNormalized_p = false;
270 :
271 0 : }
272 : //
273 : //---------------------------------------------------------------
274 : // For wide-band case, this just tells the user that the sensitivity
275 : // image will be computed during the first gridding cycle.
276 : //
277 : // Weights image (sum of convolution functions) is accumulated
278 : // during gridding. At the end of the gridding cycle, the weight
279 : // image is FT'ed (using ftWeightImage()) and properly normalized
280 : // and converted to sensitivity image (using
281 : // makeSensitivityImage()). These methods are triggred in the first
282 : // call to getImage(). In subsequent calls to getImage() these
283 : // calls are NoOps.
284 : //
285 0 : void AWProjectWBFT::makeSensitivityImage(const VisBuffer&,
286 : const ImageInterface<Complex>& /*imageTemplate*/,
287 : ImageInterface<Float>& /*sensitivityImage*/)
288 : {
289 0 : if (oneTimeMessage_p == false)
290 : {
291 0 : LogIO log_l(LogOrigin("AWProjectWBFT", "makeSensitivityImage[R&D]"));
292 0 : log_l << "Setting up for weights accumulation ";
293 0 : if (sensitivityPatternQualifierStr_p != "") log_l << "for term " << sensitivityPatternQualifier_p << " ";
294 : log_l << "during gridding to compute sensitivity pattern."
295 : << endl
296 : << "Consequently, the first gridding cycle will be slower than the subsequent ones."
297 0 : << LogIO::WARN;
298 0 : oneTimeMessage_p=true;
299 0 : }
300 0 : }
301 : //
302 : //---------------------------------------------------------------
303 : //
304 0 : void AWProjectWBFT::ftWeightImage(Lattice<Complex>& wtImage,
305 : const Matrix<Float>& sumWt,
306 : const Bool& doFFTNorm)
307 : {
308 0 : LogIO log_l(LogOrigin("AWProjectWBFT", "ftWeightImage[R&D]"));
309 0 : if (wtImageFTDone_p) return;
310 :
311 : //cerr << "From ftWeightImage" << endl;
312 : // not used
313 : // Bool doSumWtNorm=true;
314 : // if (sumWt.shape().nelements()==0) doSumWtNorm=false;
315 0 : if ((sumWt.shape().nelements() < 2) ||
316 0 : (sumWt.shape()(0) != wtImage.shape()(2)) ||
317 0 : (sumWt.shape()(1) != wtImage.shape()(3)))
318 0 : log_l << "Sum of weights per poln and chan is required" << LogIO::EXCEPTION;
319 : //
320 : // Use only the amplitude of the gridded weights. LatticeExpr
321 : // classes while useful, appear to be too strict about types (the
322 : // code below will not compile if the abs(wtImage) is not
323 : // converted back to a complex type).
324 :
325 : // LatticeExpr<Complex> le(abs(wtImage)*Complex(1,0));
326 : // wtImage.copyData(le);
327 :
328 : // {
329 : // String name("wtimg.im");
330 : // storeArrayAsImage(name,griddedWeights.coordinates(),wtImage.get());
331 : // }
332 0 : LatticeFFT::cfft2d(wtImage,false);
333 : // {
334 : // String name("ftwtimg.im");
335 : // storeArrayAsImage(name,griddedWeights.coordinates(),wtImage.get());
336 : // }
337 0 : wtImageFTDone_p=true;
338 :
339 0 : Int sizeX=wtImage.shape()(0), sizeY=wtImage.shape()(1);
340 :
341 0 : Array<Complex> wtBuf; wtImage.get(wtBuf,false);
342 0 : ArrayLattice<Complex> wtLat(wtBuf,true);
343 : //
344 : // Copy one 2D plane at a time, normalizing by the sum of weights
345 : // and possibly 2D FFT.
346 : //
347 : // Set up Lattice iterators on wtImage and sensitivityImage
348 : //
349 0 : IPosition axisPath(4, 0, 1, 2, 3);
350 0 : IPosition cursorShape(4, sizeX, sizeY, 1, 1);
351 :
352 0 : LatticeStepper wtImStepper(wtImage.shape(), cursorShape, axisPath);
353 0 : LatticeIterator<Complex> wtImIter(wtImage, wtImStepper);
354 : //
355 : // Iterate over channel and polarization axis
356 : //
357 0 : if (!doFFTNorm) sizeX=sizeY=1;
358 : //
359 : // Normalize each frequency and polarization plane of the complex
360 : // sensitivity pattern
361 : //
362 : // getSumOfCFWeights() returns the sum-of-weight appropriate for
363 : // computing sensitivity patterns (which could be different from
364 : // the data-sum-of-weights.
365 : //
366 :
367 : // USEFUL DEBUG MESSAGE
368 : //cerr << "SumCFWt: " << getSumOfCFWeights() << " " << max(wtBuf) << " " << sensitivityPatternQualifier_p << endl;
369 0 : for(wtImIter.reset(); !wtImIter.atEnd(); wtImIter++)
370 : {
371 0 : Int pol_l=wtImIter.position()(2), chan_l=wtImIter.position()(3);
372 0 : Double sumwt_l=1.0;;
373 : // Lets write some mildly obfuscated code ~[8-)
374 : //if ((sensitivityPatternQualifier_p == -1) && (doSumWtNorm))
375 : // sumwt_l = ((sumwt_l = getSumOfCFWeights()(pol_l,chan_l))==0)?1.0:sumwt_l;
376 :
377 0 : sumwt_l = getSumOfCFWeights()(pol_l,chan_l);
378 :
379 0 : wtImIter.rwCursor() = (wtImIter.rwCursor()
380 0 : *Float(sizeX)*Float(sizeY)
381 0 : /sumwt_l
382 0 : );
383 :
384 : //////////////////// wtImIter.rwCursor() = sqrt( fabs(wtImIter.rwCursor()) );
385 :
386 : // Double maxval = fabs( max( wtImIter.rwCursor() ) );
387 : // sumwt_l = getSumOfCFWeights()(pol_l,chan_l);
388 : // weightRatio_p = maxval * Float(sizeX)*Float(sizeY) / sumwt_l;
389 : }
390 :
391 0 : }
392 : //
393 : //---------------------------------------------------------------
394 : //
395 0 : void AWProjectWBFT::makeSensitivitySqImage(Lattice<Complex>& wtImage,
396 : ImageInterface<Complex>& sensitivitySqImage,
397 : const Matrix<Float>& sumWt,
398 : const Bool& doFFTNorm)
399 : {
400 : // if (avgPBReady_p) return;
401 0 : LogIO log_l(LogOrigin("AWProjectWBFT", "makeSensitivitySqImage[R&D]"));
402 :
403 : // avgPBReady_p=true;
404 :
405 0 : ftWeightImage(wtImage, sumWt, doFFTNorm);
406 :
407 0 : sensitivitySqImage.resize(griddedWeights.shape());
408 0 : sensitivitySqImage.setCoordinateInfo(griddedWeights.coordinates());
409 :
410 0 : Int sizeX=wtImage.shape()(0), sizeY=wtImage.shape()(1);
411 0 : Array<Complex> senSqBuf; sensitivitySqImage.get(senSqBuf,false);
412 0 : ArrayLattice<Complex> senSqLat(senSqBuf, true);
413 :
414 0 : Array<Complex> wtBuf; wtImage.get(wtBuf,false);
415 0 : ArrayLattice<Complex> wtLat(wtBuf,true);
416 : //
417 : // Copy one 2D plane at a time, normalizing by the sum of weights
418 : // and possibly 2D FFT.
419 : //
420 : // Set up Lattice iterators on wtImage and sensitivityImage
421 : //
422 0 : IPosition axisPath(4, 0, 1, 2, 3);
423 0 : IPosition cursorShape(4, sizeX, sizeY, 1, 1);
424 :
425 0 : LatticeStepper wtImStepper(wtImage.shape(), cursorShape, axisPath);
426 0 : LatticeIterator<Complex> wtImIter(wtImage, wtImStepper);
427 :
428 0 : LatticeStepper senSqImStepper(senSqLat.shape(), cursorShape, axisPath);
429 0 : LatticeIterator<Complex> senSqImIter(senSqLat, senSqImStepper);
430 : //
431 : // Iterate over channel and polarization axis
432 : //
433 : //
434 : // The following code is averaging RR and LL planes and writing
435 : // the result to back to both planes. This needs to be
436 : // generalized for full-pol case.
437 : //
438 0 : IPosition start0(4,0,0,0,0), start1(4,0,0,1,0), length(4,sizeX,sizeY,1,1);
439 0 : Slicer slicePol0(start0,length), slicePol1(start1,length);
440 0 : Array<Complex> polPlane0C, polPlane1C, polPlaneSq0C, polPlaneSq1C;
441 :
442 0 : senSqLat.getSlice(polPlaneSq0C, slicePol0);
443 0 : senSqLat.getSlice(polPlaneSq1C,slicePol1);
444 0 : wtLat.getSlice(polPlane0C, slicePol0);
445 0 : wtLat.getSlice(polPlane1C, slicePol1);
446 :
447 0 : polPlaneSq0C = polPlane0C*polPlane1C;
448 0 : polPlaneSq1C = polPlaneSq0C;
449 :
450 0 : pbNormalized_p=false;
451 :
452 0 : resetPBs_p=false;
453 : // String name("avgPBSq.im");
454 : // storeImg(name,sensitivitySqImage);
455 0 : }
456 : //
457 : //---------------------------------------------------------------
458 : //
459 0 : void AWProjectWBFT::makeSensitivityImage(Lattice<Complex>& wtImage,
460 : ImageInterface<Float>& sensitivityImage,
461 : const Matrix<Float>& sumWt,
462 : const Bool& doFFTNorm)
463 : {
464 0 : if (avgPBReady_p) return;
465 0 : LogIO log_l(LogOrigin("AWProjectWBFT", "makeSensitivityImage[R&D]"));
466 :
467 :
468 : // Matrix<Float> cfWts(sumWt.shape());
469 : // convertArray(cfWts,sumCFWeight);
470 : // ftWeightImage(wtImage, cfWts, doFFTNorm);
471 : // {
472 : // String name("wtimg.im");
473 : // storeArrayAsImage(name,griddedWeights.coordinates(),wtImage.get());
474 : // }
475 0 : ftWeightImage(wtImage, sumWt, doFFTNorm);
476 : // {
477 : // String name("ftwtimg.im");
478 : // storeArrayAsImage(name,griddedWeights.coordinates(),wtImage.get());
479 : // }
480 : // if (tt_pp != "")
481 : // { // uncommented by UUU
482 : // String name("ftwtimg"+sensitivityPatternQualifierStr_p+".im.conj");
483 : // storeArrayAsImage(name,griddedWeights.coordinates(),wtImage.get());
484 : // }
485 : // else
486 : // { // uncommented by UUU
487 : // String name("ftwtimg"+sensitivityPatternQualifierStr_p+".im");
488 : // storeArrayAsImage(name,griddedWeights.coordinates(),wtImage.get());
489 : // }
490 :
491 0 : sensitivityImage.resize(griddedWeights.shape());
492 0 : sensitivityImage.setCoordinateInfo(griddedWeights.coordinates());
493 :
494 0 : Int sizeX=wtImage.shape()(0), sizeY=wtImage.shape()(1);
495 0 : Array<Float> senBuf; sensitivityImage.get(senBuf,false);
496 0 : ArrayLattice<Float> senLat(senBuf, true);
497 :
498 0 : Array<Complex> wtBuf; wtImage.get(wtBuf,false);
499 0 : ArrayLattice<Complex> wtLat(wtBuf,true);
500 : //
501 : // Set up Lattice iteratos on wtImage and sensitivityImage
502 : //
503 0 : IPosition axisPath(4, 0, 1, 2, 3);
504 0 : IPosition cursorShape(4, sizeX, sizeY, 1, 1);
505 :
506 0 : LatticeStepper senImStepper(senLat.shape(), cursorShape, axisPath);
507 0 : LatticeIterator<Float> senImIter(senLat, senImStepper);
508 : //
509 : // The following code is averaging RR and LL planes and writing
510 : // the result to back to both planes. This needs to be
511 : // generalized for full-pol case.
512 : //
513 0 : IPosition start0(4,0,0,0,0), start1(4,0,0,1,0), length(4,sizeX,sizeY,1,1);
514 0 : Slicer slicePol0(start0,length), slicePol1(start1,length);
515 0 : Array<Float> polPlane0F, polPlane1F;
516 0 : Array<Complex> polPlane0C, polPlane1C;
517 :
518 : // Use StokesImageUtil to convert from Complex sensitivity pattern
519 : // to real value image planes.
520 : // StokesImageUtil::To(sensitivityImage, griddedWeights);
521 :
522 : //
523 : // Copy the real part of the average of all planes of the wtImage to
524 : // all the planes of the sensitivity image (senImage).
525 : //
526 0 : LatticeStepper wtImStepper(wtImage.shape(), cursorShape, axisPath);
527 0 : LatticeIterator<Complex> wtImIter(wtImage, wtImStepper);
528 :
529 : // Matrix<Complex> tmp(senImIter.rwMatrixCursor().shape());
530 0 : Matrix<Complex> tmp(senImIter.rwMatrixCursor().shape());
531 0 : tmp = 1.0;
532 0 : Int n=0;
533 0 : for(wtImIter.reset(); !wtImIter.atEnd(); wtImIter++)
534 : {
535 : // tmp += real(wtImIter.rwMatrixCursor().nonDegenerate());
536 0 : tmp = (wtImIter.rwMatrixCursor().nonDegenerate());
537 0 : n++;
538 : }
539 : // tmp = tmp/n;
540 : //UUU// tmp = fabs(tmp); // fabs(Array<Complex>&) returns a complex array
541 : // tmp = (sqrt(fabs(tmp))); // fabs(Array<Complex>&) returns a complex array
542 0 : for(senImIter.reset(); !senImIter.atEnd(); senImIter++)
543 0 : senImIter.rwMatrixCursor() = real(tmp);
544 :
545 0 : if (tt_pp == "")
546 0 : cfCache_p->flush(sensitivityImage,sensitivityPatternQualifierStr_p);
547 :
548 0 : pbNormalized_p=false;
549 0 : resetPBs_p=false;
550 :
551 0 : avgPBReady_p=true;
552 0 : }
553 : //
554 : //---------------------------------------------------------------
555 : //
556 : // Convert the gridded visibilities to image. This implements the
557 : // basic 2D transform (vC-Z Theorem) using the FFT.
558 : //
559 : // This specialization of getImage() exists since the sensitivity
560 : // pattern in this FTMachine is computed as the FT of the gridded
561 : // weights (convolution functions). The gridded weights are
562 : // available along with the gridded data at the end of the gridding
563 : // cycle. This method first converts the gridded weights to
564 : // sensitivity image and then calls AWProjectFT::getImage(), which
565 : // in-turn calls normalizeImage() with the sensitivityImage and the
566 : // sum of weights.
567 : //
568 0 : ImageInterface<Complex>& AWProjectWBFT::getImage(Matrix<Float>& weights,
569 : Bool fftNormalization)
570 : {
571 0 : AlwaysAssert(image, AipsError);
572 0 : LogIO log_l(LogOrigin("AWProjectWBFT", "getImage[R&D]"));
573 :
574 0 : weights.resize(sumWeight.shape());
575 0 : convertArray(weights, sumWeight);//I suppose this converts a
576 : //Matrix<Double> (sumWeights) to
577 : //Matrix<Float> (weights). Why
578 : //is this conversion required?
579 : //--SB (Dec. 2010)
580 : //
581 : // Compute the sensitivity image as FT of the griddedWegiths.
582 : // Return the result in avgPB_p. Cache it in the cfCache_p
583 : // object. And raise the avgPBReady_p flag so that this becomes a
584 : // null call the next time around.
585 : // // REMOVE THIS CODE
586 : // {
587 : // tt_pp="_conj_";
588 : // makeSensitivityImage(griddedConjWeights, *avgPB_p, weights, true);
589 : // tt_pp="";
590 : // wtImageFTDone_p = avgPBReady_p=false;
591 : // }
592 : // // REMOVE THIS CODE
593 0 : makeSensitivityImage(griddedWeights, *avgPB_p, weights, true);
594 :
595 : // if (avgPBSq_p.null()) avgPBSq_p = new TempImage<Complex>();
596 : // makeSensitivitySqImage(griddedWeights, *avgPBSq_p, weights, true);
597 :
598 : //
599 : // This calls the overloadable method normalizeImage() which
600 : // normalizes the raw image by the sensitivty pattern (avgPB_p).
601 : //
602 0 : AWProjectFT::getImage(weights,fftNormalization);
603 : // if (makingPSF)
604 : // {
605 : // String name("psf.im");
606 : // image->put(griddedData);
607 : // storeImg(name,*image);
608 : // }
609 :
610 0 : return *image;
611 0 : }
612 :
613 :
614 : // #define NEED_UNDERSCORES
615 : // #if defined(NEED_UNDERSCORES)
616 : // #define gpbmos gpbmos_
617 : // #define dpbmos dpbmos_
618 : // #define dpbmosgrad dpbmosgrad_
619 : // #define dpbwgrad dpbwgrad_
620 : // #endif
621 :
622 : // extern "C" {
623 : // void gpbmos(Double *uvw,
624 : // Double *dphase,
625 : // const Complex *values,
626 : // Int *nvispol,
627 : // Int *nvischan,
628 : // Int *dopsf,
629 : // const Int *flag,
630 : // const Int *rflag,
631 : // const Float *weight,
632 : // Int *nrow,
633 : // Int *rownum,
634 : // Double *scale,
635 : // Double *offset,
636 : // Complex *grid,
637 : // Int *nx,
638 : // Int *ny,
639 : // Int *npol,
640 : // Int *nchan,
641 : // const Double *freq,
642 : // const Double *c,
643 : // Int *support,
644 : // Int *convsize,
645 : // Int *convwtsize,
646 : // Int *sampling,
647 : // Int *wconvsize,
648 : // Complex *convfunc,
649 : // Int *chanmap,
650 : // Int *polmap,
651 : // Int *polused,
652 : // Double *sumwt,
653 : // Int *ant1,
654 : // Int *ant2,
655 : // Int *nant,
656 : // Int *scanno,
657 : // Double *sigma,
658 : // Float *raoff,
659 : // Float *decoff,
660 : // Double *area,
661 : // Int *doGrad,
662 : // Int *doPointingCorrection,
663 : // Int *nPA,
664 : // Int *paIndex,
665 : // Int *CFMap,
666 : // Int *ConjCFMap,
667 : // Double *currentCFPA, Double *actualPA,
668 : // Int *avgPBReady_p,
669 : // Complex *avgPB, Double *cfRefFreq_p,
670 : // Complex *convWeights,
671 : // Int *convWtSupport);
672 : // void dpbmos(Double *uvw,
673 : // Double *dphase,
674 : // Complex *values,
675 : // Int *nvispol,
676 : // Int *nvischan,
677 : // const Int *flag,
678 : // const Int *rflag,
679 : // Int *nrow,
680 : // Int *rownum,
681 : // Double *scale,
682 : // Double *offset,
683 : // const Complex *grid,
684 : // Int *nx,
685 : // Int *ny,
686 : // Int *npol,
687 : // Int *nchan,
688 : // const Double *freq,
689 : // const Double *c,
690 : // Int *support,
691 : // Int *convsize,
692 : // Int *sampling,
693 : // Int *wconvsize,
694 : // Complex *convfunc,
695 : // Int *chanmap,
696 : // Int *polmap,
697 : // Int *polused,
698 : // Int *ant1,
699 : // Int *ant2,
700 : // Int *nant,
701 : // Int *scanno,
702 : // Double *sigma,
703 : // Float *raoff, Float *decoff,
704 : // Double *area,
705 : // Int *dograd,
706 : // Int *doPointingCorrection,
707 : // Int *nPA,
708 : // Int *paIndex,
709 : // Int *CFMap,
710 : // Int *ConjCFMap,
711 : // Double *currentCFPA, Double *actualPA, Double *cfRefFreq_p);
712 : // void dpbmosgrad(Double *uvw,
713 : // Double *dphase,
714 : // Complex *values,
715 : // Int *nvispol,
716 : // Int *nvischan,
717 : // Complex *gazvalues,
718 : // Complex *gelvalues,
719 : // Int *doconj,
720 : // const Int *flag,
721 : // const Int *rflag,
722 : // Int *nrow,
723 : // Int *rownum,
724 : // Double *scale,
725 : // Double *offset,
726 : // const Complex *grid,
727 : // Int *nx,
728 : // Int *ny,
729 : // Int *npol,
730 : // Int *nchan,
731 : // const Double *freq,
732 : // const Double *c,
733 : // Int *support,
734 : // Int *convsize,
735 : // Int *sampling,
736 : // Int *wconvsize,
737 : // Complex *convfunc,
738 : // Int *chanmap,
739 : // Int *polmap,
740 : // Int *polused,
741 : // Int *ant1,
742 : // Int *ant2,
743 : // Int *nant,
744 : // Int *scanno,
745 : // Double *sigma,
746 : // Float *raoff, Float *decoff,
747 : // Double *area,
748 : // Int *dograd,
749 : // Int *doPointingCorrection,
750 : // Int *nPA,
751 : // Int *paIndex,
752 : // Int *CFMap,
753 : // Int *ConjCFMap,
754 : // Double *currentCFPA, Double *actualPA,Double *cfRefFreq_p);
755 : // }
756 : // //
757 : // //----------------------------------------------------------------------
758 : // //
759 : // void AWProjectWBFT::runFortranGet(Matrix<Double>& uvw,Vector<Double>& dphase,
760 : // Cube<Complex>& visdata,
761 : // IPosition& s,
762 : // Int& Conj,
763 : // Cube<Int>& flags,Vector<Int>& rowFlags,
764 : // Int& rownr,Vector<Double>& actualOffset,
765 : // Array<Complex>* dataPtr,
766 : // Int& aNx, Int& aNy, Int& npol, Int& nchan,
767 : // VisBuffer& vb,Int& Nant_p, Int& scanNo,
768 : // Double& sigma,
769 : // Array<Float>& l_off,
770 : // Array<Float>& m_off,
771 : // Double area,
772 : // Int& doGrad,
773 : // Int paIndex)
774 : // {
775 : // (void)Conj;
776 : // LogIO log_l(LogOrigin("AWProjectWBFT","runFortranGet[R&D]"));
777 : // enum whichGetStorage {RAOFF,DECOFF,UVW,DPHASE,VISDATA,GRADVISAZ,GRADVISEL,
778 : // FLAGS,ROWFLAGS,UVSCALE,ACTUALOFFSET,DATAPTR,VBFREQ,
779 : // CONVSUPPORT,CONVFUNC,CHANMAP,POLMAP,VBANT1,VBANT2,CONJCFMAP,CFMAP};
780 : // Vector<Bool> deleteThem(21);
781 :
782 : // Double *uvw_p, *dphase_p, *actualOffset_p, *vb_freq_p, *uvScale_p;
783 : // Complex *visdata_p, *dataPtr_p, *f_convFunc_p;
784 : // Int *flags_p, *rowFlags_p, *chanMap_p, *polMap_p, *convSupport_p, *vb_ant1_p, *vb_ant2_p,
785 : // *ConjCFMap_p, *CFMap_p;
786 : // Float *l_off_p, *m_off_p;
787 : // Double actualPA;
788 :
789 : // Vector<Int> ConjCFMap, CFMap;
790 : // Int N;
791 : // actualPA = getVBPA(vb);
792 :
793 : // N=polMap.nelements();
794 : // CFMap = polMap; ConjCFMap = polMap;
795 : // for(Int i=0;i<N;i++) CFMap[i] = polMap[N-i-1];
796 :
797 : // Array<Complex> rotatedConvFunc;
798 :
799 : // // Rotate the convolution function using Image rotation and
800 : // // disable rotation in the gridder
801 : // SynthesisUtils::rotateComplexArray(log_l, *cfs_p.data,/*convFunc_p*/ cfs_p.coordSys,
802 : // rotatedConvFunc,(currentCFPA-actualPA),"LINEAR");
803 : // actualPA = currentCFPA;
804 :
805 : // // SynthesisUtils::rotateComplexArray(log_l, convFunc_p, cfs_p.coordSys,
806 : // // rotatedConvFunc,0.0);
807 :
808 : // ConjCFMap = polMap;
809 : // makeCFPolMap(vb,cfStokes_p,CFMap);
810 : // makeConjPolMap(vb,CFMap,ConjCFMap);
811 :
812 :
813 : // ConjCFMap_p = ConjCFMap.getStorage(deleteThem(CONJCFMAP));
814 : // CFMap_p = CFMap.getStorage(deleteThem(CFMAP));
815 :
816 : // uvw_p = uvw.getStorage(deleteThem(UVW));
817 : // dphase_p = dphase.getStorage(deleteThem(DPHASE));
818 : // visdata_p = visdata.getStorage(deleteThem(VISDATA));
819 : // flags_p = flags.getStorage(deleteThem(FLAGS));
820 : // rowFlags_p = rowFlags.getStorage(deleteThem(ROWFLAGS));
821 : // uvScale_p = uvScale.getStorage(deleteThem(UVSCALE));
822 : // actualOffset_p = actualOffset.getStorage(deleteThem(ACTUALOFFSET));
823 : // dataPtr_p = dataPtr->getStorage(deleteThem(DATAPTR));
824 : // vb_freq_p = vb.frequency().getStorage(deleteThem(VBFREQ));
825 : // convSupport_p = cfs_p.xSupport.getStorage(deleteThem(CONVSUPPORT));
826 : // // f_convFunc_p = convFunc_p.getStorage(deleteThem(CONVFUNC));
827 : // f_convFunc_p = rotatedConvFunc.getStorage(deleteThem(CONVFUNC));
828 : // chanMap_p = chanMap.getStorage(deleteThem(CHANMAP));
829 : // polMap_p = polMap.getStorage(deleteThem(POLMAP));
830 : // vb_ant1_p = vb.antenna1().getStorage(deleteThem(VBANT1));
831 : // vb_ant2_p = vb.antenna2().getStorage(deleteThem(VBANT2));
832 : // l_off_p = l_off.getStorage(deleteThem(RAOFF));
833 : // m_off_p = m_off.getStorage(deleteThem(DECOFF));
834 :
835 : // // static int ttt=0;
836 : // // if (ttt==0) cout << l_off(IPosition(3,0,0,0)) << " " << m_off(IPosition(3,0,0,0)) << endl;
837 : // // ttt++;
838 :
839 : // Int npa=1,actualConvSize;
840 : // Int paIndex_Fortran = paIndex;
841 : // // actualConvSize = convFunc_p.shape()(0);
842 : // actualConvSize = cfs_p.data->shape()(0);
843 :
844 : // // IPosition shp=convSupport.shape();
845 : // Int alwaysDoPointing=1;
846 : // alwaysDoPointing=doPointing;
847 : // dpbmos(uvw_p,
848 : // dphase_p,
849 : // visdata_p,
850 : // &s.asVector()(0),
851 : // &s.asVector()(1),
852 : // flags_p,
853 : // rowFlags_p,
854 : // &s.asVector()(2),
855 : // &rownr,
856 : // uvScale_p,
857 : // actualOffset_p,
858 : // dataPtr_p,
859 : // &aNx,
860 : // &aNy,
861 : // &npol,
862 : // &nchan,
863 : // vb_freq_p,
864 : // &C::c,
865 : // convSupport_p,
866 : // &actualConvSize,
867 : // &convSampling,
868 : // &wConvSize,
869 : // f_convFunc_p,
870 : // chanMap_p,
871 : // polMap_p,
872 : // &polInUse_p,
873 : // vb_ant1_p,
874 : // vb_ant2_p,
875 : // &Nant_p,
876 : // &scanNo,
877 : // &sigma,
878 : // l_off_p, m_off_p,
879 : // &area,
880 : // &doGrad,
881 : // &alwaysDoPointing,
882 : // &npa,
883 : // &paIndex_Fortran,
884 : // CFMap_p,
885 : // ConjCFMap_p,
886 : // ¤tCFPA
887 : // ,&actualPA, &cfRefFreq_p
888 : // );
889 :
890 : // ConjCFMap.freeStorage((const Int *&)ConjCFMap_p,deleteThem(CONJCFMAP));
891 : // CFMap.freeStorage((const Int *&)CFMap_p,deleteThem(CFMAP));
892 :
893 : // l_off.freeStorage((const Float*&)l_off_p,deleteThem(RAOFF));
894 : // m_off.freeStorage((const Float*&)m_off_p,deleteThem(DECOFF));
895 : // uvw.freeStorage((const Double*&)uvw_p,deleteThem(UVW));
896 : // dphase.freeStorage((const Double*&)dphase_p,deleteThem(DPHASE));
897 : // visdata.putStorage(visdata_p,deleteThem(VISDATA));
898 : // flags.freeStorage((const Int*&) flags_p,deleteThem(FLAGS));
899 : // rowFlags.freeStorage((const Int *&)rowFlags_p,deleteThem(ROWFLAGS));
900 : // actualOffset.freeStorage((const Double*&)actualOffset_p,deleteThem(ACTUALOFFSET));
901 : // dataPtr->freeStorage((const Complex *&)dataPtr_p,deleteThem(DATAPTR));
902 : // uvScale.freeStorage((const Double*&) uvScale_p,deleteThem(UVSCALE));
903 : // vb.frequency().freeStorage((const Double*&)vb_freq_p,deleteThem(VBFREQ));
904 : // cfs_p.xSupport.freeStorage((const Int*&)convSupport_p,deleteThem(CONVSUPPORT));
905 : // // convFunc_p.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVFUNC));
906 : // chanMap.freeStorage((const Int*&)chanMap_p,deleteThem(CHANMAP));
907 : // polMap.freeStorage((const Int*&) polMap_p,deleteThem(POLMAP));
908 : // vb.antenna1().freeStorage((const Int*&) vb_ant1_p,deleteThem(VBANT1));
909 : // vb.antenna2().freeStorage((const Int*&) vb_ant2_p,deleteThem(VBANT2));
910 : // }
911 : // //
912 : // //----------------------------------------------------------------------
913 : // //
914 : // void AWProjectWBFT::runFortranGetGrad(Matrix<Double>& uvw,Vector<Double>& dphase,
915 : // Cube<Complex>& visdata,
916 : // IPosition& s,
917 : // Cube<Complex>& gradVisAzData,
918 : // Cube<Complex>& gradVisElData,
919 : // Int& Conj,
920 : // Cube<Int>& flags,Vector<Int>& rowFlags,
921 : // Int& rownr,Vector<Double>& actualOffset,
922 : // Array<Complex>* dataPtr,
923 : // Int& aNx, Int& aNy, Int& npol, Int& nchan,
924 : // VisBuffer& vb,Int& Nant_p, Int& scanNo,
925 : // Double& sigma,
926 : // Array<Float>& l_off,
927 : // Array<Float>& m_off,
928 : // Double area,
929 : // Int& doGrad,
930 : // Int paIndex)
931 : // {
932 : // LogIO log_l(LogOrigin("AWProjectWBFT","runFortranGetGrad[R&D]"));
933 : // enum whichGetStorage {RAOFF,DECOFF,UVW,DPHASE,VISDATA,GRADVISAZ,GRADVISEL,
934 : // FLAGS,ROWFLAGS,UVSCALE,ACTUALOFFSET,DATAPTR,VBFREQ,
935 : // CONVSUPPORT,CONVFUNC,CHANMAP,POLMAP,VBANT1,VBANT2,CONJCFMAP,CFMAP};
936 : // Vector<Bool> deleteThem(21);
937 :
938 : // Double *uvw_p, *dphase_p, *actualOffset_p, *vb_freq_p, *uvScale_p;
939 : // Complex *visdata_p, *dataPtr_p, *f_convFunc_p;
940 : // Complex *gradVisAzData_p, *gradVisElData_p;
941 : // Int *flags_p, *rowFlags_p, *chanMap_p, *polMap_p, *convSupport_p, *vb_ant1_p, *vb_ant2_p,
942 : // *ConjCFMap_p, *CFMap_p;
943 : // Float *l_off_p, *m_off_p;
944 : // Double actualPA;
945 :
946 : // Vector<Int> ConjCFMap, CFMap;
947 : // actualPA = getVBPA(vb);
948 : // ConjCFMap = polMap;
949 : // makeCFPolMap(vb,cfStokes_p,CFMap);
950 : // makeConjPolMap(vb,CFMap,ConjCFMap);
951 :
952 : // Array<Complex> rotatedConvFunc;
953 : // // Rotate the convolution function using Image rotation and
954 : // // disable rotation in the gridder
955 : // SynthesisUtils::rotateComplexArray(log_l, *(cfs_p.data) /*convFunc_p*/, cfs_p.coordSys,
956 : // rotatedConvFunc,(currentCFPA-actualPA));
957 : // actualPA = currentCFPA;
958 :
959 : // // SynthesisUtils::rotateComplexArray(log_l, convFunc_p, cfs_p.coordSys,
960 : // // rotatedConvFunc,0.0);
961 :
962 : // ConjCFMap_p = ConjCFMap.getStorage(deleteThem(CONJCFMAP));
963 : // CFMap_p = CFMap.getStorage(deleteThem(CFMAP));
964 :
965 : // uvw_p = uvw.getStorage(deleteThem(UVW));
966 : // dphase_p = dphase.getStorage(deleteThem(DPHASE));
967 : // visdata_p = visdata.getStorage(deleteThem(VISDATA));
968 : // gradVisAzData_p = gradVisAzData.getStorage(deleteThem(GRADVISAZ));
969 : // gradVisElData_p = gradVisElData.getStorage(deleteThem(GRADVISEL));
970 : // flags_p = flags.getStorage(deleteThem(FLAGS));
971 : // rowFlags_p = rowFlags.getStorage(deleteThem(ROWFLAGS));
972 : // uvScale_p = uvScale.getStorage(deleteThem(UVSCALE));
973 : // actualOffset_p = actualOffset.getStorage(deleteThem(ACTUALOFFSET));
974 : // dataPtr_p = dataPtr->getStorage(deleteThem(DATAPTR));
975 : // vb_freq_p = vb.frequency().getStorage(deleteThem(VBFREQ));
976 : // convSupport_p = cfs_p.xSupport.getStorage(deleteThem(CONVSUPPORT));
977 : // // f_convFunc_p = convFunc_p.getStorage(deleteThem(CONVFUNC));
978 : // f_convFunc_p = rotatedConvFunc.getStorage(deleteThem(CONVFUNC));
979 : // chanMap_p = chanMap.getStorage(deleteThem(CHANMAP));
980 : // polMap_p = polMap.getStorage(deleteThem(POLMAP));
981 : // vb_ant1_p = vb.antenna1().getStorage(deleteThem(VBANT1));
982 : // vb_ant2_p = vb.antenna2().getStorage(deleteThem(VBANT2));
983 : // l_off_p = l_off.getStorage(deleteThem(RAOFF));
984 : // m_off_p = m_off.getStorage(deleteThem(DECOFF));
985 :
986 : // Int npa=1,actualConvSize;
987 : // Int paIndex_Fortran = paIndex;
988 : // // actualConvSize = convFunc_p.shape()(0);
989 : // actualConvSize = cfs_p.data->shape()(0);
990 :
991 : // // IPosition shp=convSupport.shape();
992 : // Int alwaysDoPointing=1;
993 : // alwaysDoPointing = doPointing;
994 : // dpbmosgrad(uvw_p,
995 : // dphase_p,
996 : // visdata_p,
997 : // &s.asVector()(0),
998 : // &s.asVector()(1),
999 : // gradVisAzData_p,
1000 : // gradVisElData_p,
1001 : // &Conj,
1002 : // flags_p,
1003 : // rowFlags_p,
1004 : // &s.asVector()(2),
1005 : // &rownr,
1006 : // uvScale_p,
1007 : // actualOffset_p,
1008 : // dataPtr_p,
1009 : // &aNx,
1010 : // &aNy,
1011 : // &npol,
1012 : // &nchan,
1013 : // vb_freq_p,
1014 : // &C::c,
1015 : // convSupport_p,
1016 : // &actualConvSize,
1017 : // &convSampling,
1018 : // &wConvSize,
1019 : // f_convFunc_p,
1020 : // chanMap_p,
1021 : // polMap_p,
1022 : // &polInUse_p,
1023 : // vb_ant1_p,
1024 : // vb_ant2_p,
1025 : // &Nant_p,
1026 : // &scanNo,
1027 : // &sigma,
1028 : // l_off_p, m_off_p,
1029 : // &area,
1030 : // &doGrad,
1031 : // &alwaysDoPointing,
1032 : // &npa,
1033 : // &paIndex_Fortran,
1034 : // CFMap_p,
1035 : // ConjCFMap_p,
1036 : // ¤tCFPA
1037 : // ,&actualPA,&cfRefFreq_p
1038 : // );
1039 : // ConjCFMap.freeStorage((const Int *&)ConjCFMap_p,deleteThem(CONJCFMAP));
1040 : // CFMap.freeStorage((const Int *&)CFMap_p,deleteThem(CFMAP));
1041 :
1042 : // l_off.freeStorage((const Float*&)l_off_p,deleteThem(RAOFF));
1043 : // m_off.freeStorage((const Float*&)m_off_p,deleteThem(DECOFF));
1044 : // uvw.freeStorage((const Double*&)uvw_p,deleteThem(UVW));
1045 : // dphase.freeStorage((const Double*&)dphase_p,deleteThem(DPHASE));
1046 : // visdata.putStorage(visdata_p,deleteThem(VISDATA));
1047 : // gradVisAzData.putStorage(gradVisAzData_p,deleteThem(GRADVISAZ));
1048 : // gradVisElData.putStorage(gradVisElData_p,deleteThem(GRADVISEL));
1049 : // flags.freeStorage((const Int*&) flags_p,deleteThem(FLAGS));
1050 : // rowFlags.freeStorage((const Int *&)rowFlags_p,deleteThem(ROWFLAGS));
1051 : // actualOffset.freeStorage((const Double*&)actualOffset_p,deleteThem(ACTUALOFFSET));
1052 : // dataPtr->freeStorage((const Complex *&)dataPtr_p,deleteThem(DATAPTR));
1053 : // uvScale.freeStorage((const Double*&) uvScale_p,deleteThem(UVSCALE));
1054 : // vb.frequency().freeStorage((const Double*&)vb_freq_p,deleteThem(VBFREQ));
1055 : // cfs_p.xSupport.freeStorage((const Int*&)convSupport_p,deleteThem(CONVSUPPORT));
1056 : // // convFunc_p.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVFUNC));
1057 : // chanMap.freeStorage((const Int*&)chanMap_p,deleteThem(CHANMAP));
1058 : // polMap.freeStorage((const Int*&) polMap_p,deleteThem(POLMAP));
1059 : // vb.antenna1().freeStorage((const Int*&) vb_ant1_p,deleteThem(VBANT1));
1060 : // vb.antenna2().freeStorage((const Int*&) vb_ant2_p,deleteThem(VBANT2));
1061 : // }
1062 : // //
1063 : // //----------------------------------------------------------------------
1064 : // //
1065 : // void AWProjectWBFT::runFortranPut(Matrix<Double>& uvw,Vector<Double>& dphase,
1066 : // const Complex& visdata,
1067 : // IPosition& s,
1068 : // Int& Conj,
1069 : // Cube<Int>& flags,Vector<Int>& rowFlags,
1070 : // const Matrix<Float>& weight,
1071 : // Int& rownr,Vector<Double>& actualOffset,
1072 : // Array<Complex>& dataPtr,
1073 : // Int& aNx, Int& aNy, Int& npol, Int& nchan,
1074 : // const VisBuffer& vb,Int& Nant_p, Int& scanNo,
1075 : // Double& sigma,
1076 : // Array<Float>& l_off,
1077 : // Array<Float>& m_off,
1078 : // Matrix<Double>& sumWeight,
1079 : // Double& area,
1080 : // Int& doGrad,
1081 : // Int& doPSF,
1082 : // Int paIndex)
1083 : // {
1084 : // (void)Conj;
1085 : // LogIO log_l(LogOrigin("AWProjectWBFT","runFortranPut[R&D]"));
1086 : // enum whichGetStorage {RAOFF,DECOFF,UVW,DPHASE,VISDATA,GRADVISAZ,GRADVISEL,
1087 : // FLAGS,ROWFLAGS,UVSCALE,ACTUALOFFSET,DATAPTR,VBFREQ,
1088 : // CONVSUPPORT,CONVWTSUPPORT,CONVFUNC,CHANMAP,POLMAP,VBANT1,VBANT2,WEIGHT,
1089 : // SUMWEIGHT,CONJCFMAP,CFMAP,AVGPBPTR,CONVWTS};
1090 : // Vector<Bool> deleteThem(25);
1091 :
1092 : // Double *uvw_p, *dphase_p, *actualOffset_p, *vb_freq_p, *uvScale_p;
1093 : // Complex *dataPtr_p, *f_convFunc_p, *f_convWts_p,*avgPBPtr;
1094 : // Int *flags_p, *rowFlags_p, *chanMap_p, *polMap_p, *convSupport_p, *convWtSupport_p,
1095 : // *vb_ant1_p, *vb_ant2_p,
1096 : // *ConjCFMap_p, *CFMap_p;
1097 : // Float *l_off_p, *m_off_p;
1098 : // Float *weight_p;Double *sumwt_p,*isumwt_p;
1099 : // Double actualPA;
1100 : // const Complex *visdata_p=&visdata;
1101 :
1102 : // Matrix<Double> iSumWt(sumWeight.shape());
1103 : // iSumWt=0.0;
1104 : // Vector<Int> ConjCFMap, CFMap;
1105 : // actualPA = getVBPA(vb);
1106 :
1107 : // ConjCFMap = polMap;
1108 : // makeCFPolMap(vb,cfStokes_p,CFMap);
1109 : // makeConjPolMap(vb,CFMap,ConjCFMap);
1110 :
1111 : // // {
1112 : // // IPosition tt;
1113 : // // cfs_p.coordSys.list(log_l,MDoppler::RADIO,tt,tt);
1114 : // // cfwts_p.coordSys.list(log_l,MDoppler::RADIO,tt,tt);
1115 : // // }
1116 : // Array<Complex> rotatedConvFunc_l, rotatedConvWeights_l;
1117 :
1118 : // // Rotate the convolution function using Image rotation and
1119 : // // disable rotation in the gridder
1120 : // SynthesisUtils::rotateComplexArray(log_l, *(cfs_p.data)/*convFunc_p*/, cfs_p.coordSys,
1121 : // rotatedConvFunc_l,(currentCFPA-actualPA));
1122 : // if (!avgPBReady_p)
1123 : // SynthesisUtils::rotateComplexArray(log_l, *(cfwts_p.data), /*convWeights_p,*/ cfwts_p.coordSys,
1124 : // rotatedConvWeights_l,(currentCFPA-actualPA));
1125 : // // Disable rotation in the gridder
1126 : // actualPA = currentCFPA;
1127 :
1128 : // // SynthesisUtils::rotateComplexArray(log_l, convFunc_p, cfs_p.coordSys,
1129 : // // rotatedConvFunc_l,0.0);
1130 : // // if (!avgPBReady_p)
1131 : // // SynthesisUtils::rotateComplexArray(log_l, convWeights_p, cfwts_p.coordSys,
1132 : // // rotatedConvWeights_l,0.0);
1133 :
1134 :
1135 : // ConjCFMap_p = ConjCFMap.getStorage(deleteThem(CONJCFMAP));
1136 : // CFMap_p = CFMap.getStorage(deleteThem(CFMAP));
1137 :
1138 : // uvw_p = uvw.getStorage(deleteThem(UVW));
1139 : // dphase_p = dphase.getStorage(deleteThem(DPHASE));
1140 : // flags_p = flags.getStorage(deleteThem(FLAGS));
1141 : // rowFlags_p = rowFlags.getStorage(deleteThem(ROWFLAGS));
1142 : // uvScale_p = uvScale.getStorage(deleteThem(UVSCALE));
1143 : // actualOffset_p = actualOffset.getStorage(deleteThem(ACTUALOFFSET));
1144 : // dataPtr_p = dataPtr.getStorage(deleteThem(DATAPTR));
1145 : // vb_freq_p = (Double *)(vb.frequency().getStorage(deleteThem(VBFREQ)));
1146 : // convSupport_p = cfs_p.xSupport.getStorage(deleteThem(CONVSUPPORT));
1147 : // convWtSupport_p = cfwts_p.xSupport.getStorage(deleteThem(CONVWTSUPPORT));
1148 : // f_convFunc_p = rotatedConvFunc_l.getStorage(deleteThem(CONVFUNC));
1149 : // f_convWts_p = rotatedConvWeights_l.getStorage(deleteThem(CONVWTS));
1150 : // chanMap_p = chanMap.getStorage(deleteThem(CHANMAP));
1151 : // polMap_p = polMap.getStorage(deleteThem(POLMAP));
1152 : // vb_ant1_p = (Int *)(vb.antenna1().getStorage(deleteThem(VBANT1)));
1153 : // vb_ant2_p = (Int *)(vb.antenna2().getStorage(deleteThem(VBANT2)));
1154 : // l_off_p = l_off.getStorage(deleteThem(RAOFF));
1155 : // m_off_p = m_off.getStorage(deleteThem(DECOFF));
1156 : // weight_p = (Float *)(weight.getStorage(deleteThem(WEIGHT)));
1157 : // sumwt_p = sumWeight.getStorage(deleteThem(SUMWEIGHT));
1158 : // Bool tmp;
1159 : // isumwt_p = iSumWt.getStorage(tmp);
1160 :
1161 : // // Array<Complex> avgPB_p(griddedWeights.get());
1162 : // Array<Complex> avgAperture;
1163 : // if (!avgPBReady_p)
1164 : // {
1165 : // avgAperture.resize(griddedWeights.shape());
1166 : // avgAperture.set(Complex(0,0));
1167 : // avgPBPtr = avgAperture.getStorage(deleteThem(AVGPBPTR));
1168 : // }
1169 : // else
1170 : // avgPBPtr=NULL;
1171 :
1172 : // Int npa=1,actualConvSize, actualConvWtSize;
1173 : // Int paIndex_Fortran = paIndex;
1174 : // // Int doAvgPB=((avgPBReady_p==false) &&
1175 : // // ((fabs(lastPAUsedForWtImg-actualPA)*57.2956 >= DELTAPA) ||
1176 : // // (lastPAUsedForWtImg == MAGICPAVALUE)));
1177 :
1178 : // Int doAvgPB=computeAvgPB(actualPA, lastPAUsedForWtImg);//(avgPBReady_p==false);
1179 : // // actualConvSize = convFunc_p.shape()(0);
1180 : // actualConvSize = cfs_p.data->shape()(0);
1181 : // // actualConvWtSize = convWeights_p.shape()(0);
1182 : // actualConvWtSize = cfwts_p.data->shape()(0);
1183 :
1184 : // if (fabs(lastPAUsedForWtImg-actualPA)*57.2956 >= DELTAPA) lastPAUsedForWtImg = actualPA;
1185 :
1186 : // // IPosition shp=convSupport.shape();
1187 : // Int alwaysDoPointing=1;
1188 : // alwaysDoPointing = doPointing;
1189 :
1190 : // gpbmos(uvw_p,
1191 : // dphase_p,
1192 : // visdata_p,
1193 : // &s.asVector()(0),
1194 : // &s.asVector()(1),
1195 : // &doPSF,
1196 : // flags_p,
1197 : // rowFlags_p,
1198 : // weight_p,
1199 : // &s.asVector()(2),
1200 : // &rownr,
1201 : // uvScale_p,
1202 : // actualOffset_p,
1203 : // dataPtr_p,
1204 : // &aNx,
1205 : // &aNy,
1206 : // &npol,
1207 : // &nchan,
1208 : // vb_freq_p,
1209 : // &C::c,
1210 : // convSupport_p,
1211 : // &actualConvSize,
1212 : // &actualConvWtSize,
1213 : // &convSampling,
1214 : // &wConvSize,
1215 : // f_convFunc_p,
1216 : // chanMap_p,
1217 : // polMap_p,
1218 : // &polInUse_p,
1219 : // // sumwt_p,
1220 : // isumwt_p,
1221 : // vb_ant1_p,
1222 : // vb_ant2_p,
1223 : // &Nant_p,
1224 : // &scanNo,
1225 : // &sigma,
1226 : // l_off_p, m_off_p,
1227 : // &area,
1228 : // &doGrad,
1229 : // &alwaysDoPointing,
1230 : // &npa,
1231 : // &paIndex_Fortran,
1232 : // CFMap_p,
1233 : // ConjCFMap_p,
1234 : // ¤tCFPA,&actualPA,
1235 : // &doAvgPB,
1236 : // avgPBPtr,&cfRefFreq_p,
1237 : // f_convWts_p,convWtSupport_p
1238 : // );
1239 :
1240 : // ConjCFMap.freeStorage((const Int *&)ConjCFMap_p,deleteThem(CONJCFMAP));
1241 : // CFMap.freeStorage((const Int *&)CFMap_p,deleteThem(CFMAP));
1242 :
1243 : // l_off.freeStorage((const Float*&)l_off_p,deleteThem(RAOFF));
1244 : // m_off.freeStorage((const Float*&)m_off_p,deleteThem(DECOFF));
1245 : // uvw.freeStorage((const Double*&)uvw_p,deleteThem(UVW));
1246 : // dphase.freeStorage((const Double*&)dphase_p,deleteThem(DPHASE));
1247 : // flags.freeStorage((const Int*&) flags_p,deleteThem(FLAGS));
1248 : // rowFlags.freeStorage((const Int *&)rowFlags_p,deleteThem(ROWFLAGS));
1249 : // actualOffset.freeStorage((const Double*&)actualOffset_p,deleteThem(ACTUALOFFSET));
1250 : // dataPtr.freeStorage((const Complex *&)dataPtr_p,deleteThem(DATAPTR));
1251 : // uvScale.freeStorage((const Double*&) uvScale_p,deleteThem(UVSCALE));
1252 : // vb.frequency().freeStorage((const Double*&)vb_freq_p,deleteThem(VBFREQ));
1253 : // cfs_p.xSupport.freeStorage((const Int*&)convSupport_p,deleteThem(CONVSUPPORT));
1254 : // // convFunc_p.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVFUNC));
1255 : // // convWeights_p.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVWTS));
1256 : // chanMap.freeStorage((const Int*&)chanMap_p,deleteThem(CHANMAP));
1257 : // polMap.freeStorage((const Int*&) polMap_p,deleteThem(POLMAP));
1258 : // vb.antenna1().freeStorage((const Int*&) vb_ant1_p,deleteThem(VBANT1));
1259 : // vb.antenna2().freeStorage((const Int*&) vb_ant2_p,deleteThem(VBANT2));
1260 : // weight.freeStorage((const Float*&)weight_p,deleteThem(WEIGHT));
1261 : // sumWeight.putStorage(sumwt_p,deleteThem(SUMWEIGHT));
1262 : // iSumWt.putStorage(isumwt_p,tmp);
1263 : // sumWeight += iSumWt;
1264 :
1265 : // if (!avgPBReady_p)
1266 : // {
1267 : // // Get the griddedWeigths as a referenced array
1268 : // Array<Complex> gwts; Bool removeDegenerateAxis=false;
1269 : // griddedWeights.get(gwts, removeDegenerateAxis);
1270 : // // griddedWeights.put(griddedWeights.get()+avgPB_p);
1271 : // gwts = gwts + avgAperture;
1272 : // // if (!reference)
1273 : // // griddedWeights.put(gwts);
1274 : // }
1275 : // }
1276 : //
1277 : //---------------------------------------------------------------
1278 : //
1279 : // Initialize the FFT to the Sky. Here we have to setup and
1280 : // initialize the grid.
1281 : //
1282 0 : void AWProjectWBFT::initializeToSky(ImageInterface<Complex>& iimage,
1283 : Matrix<Float>& weight,
1284 : const VisBuffer& vb)
1285 : {
1286 0 : LogIO log_l(LogOrigin("AWProjectWBFT","initializeToSky[R&D]"));
1287 0 : AWProjectFT::initializeToSky(iimage,weight,vb);
1288 : // The following code is same as that in the parent class call above.
1289 : // image=&iimage;
1290 :
1291 : // init();
1292 : // initMaps(vb);
1293 : // nx = image->shape()(0);
1294 : // ny = image->shape()(1);
1295 : // npol = image->shape()(2);
1296 : // nchan = image->shape()(3);
1297 :
1298 : // isTiled = false;
1299 :
1300 : // sumWeight=0.0;
1301 : // weight.resize(sumWeight.shape());
1302 : // weight=0.0;
1303 :
1304 : // if(isTiled)
1305 : // {
1306 : // imageCache->flush();
1307 : // image->set(Complex(0.0));
1308 : // //lattice=image;
1309 : // lattice=CountedPtr<Lattice<Complex> > (image, false);
1310 : // }
1311 : // else
1312 : // {
1313 : // IPosition gridShape(4, nx, ny, npol, nchan);
1314 : // griddedData.resize(gridShape);
1315 : // griddedData=Complex(0.0);
1316 : // // if(arrayLattice) delete arrayLattice; arrayLattice=0;
1317 : // arrayLattice = new ArrayLattice<Complex>(griddedData);
1318 : // lattice=arrayLattice;
1319 : // visResampler_p->initializeToSky(griddedData, sumWeight);
1320 : // }
1321 :
1322 : //AlwaysAssert(lattice, AipsError);
1323 0 : if (resetPBs_p)
1324 : {
1325 0 : griddedWeights.resize(iimage.shape());
1326 0 : griddedWeights.setCoordinateInfo(iimage.coordinates());
1327 0 : griddedWeights.set(0.0);
1328 0 : pbPeaks.resize(griddedWeights.shape()(2));
1329 0 : pbPeaks.set(0.0);
1330 :
1331 0 : resetPBs_p=false;
1332 : }
1333 0 : avgPBReady_p = (cfCache_p->loadAvgPB(avgPB_p,sensitivityPatternQualifierStr_p) != CFDefs::NOTCACHED);
1334 : // avgPBReady_p = cfCache_p->avgPBReady(sensitivityPatternQualifierStr_p);
1335 : // Need to grid the weighted Convolution Functions to make the sensitivity pattern.
1336 0 : if (!avgPBReady_p)
1337 : {
1338 : // Make a copy of the re-sampler and set it up.
1339 0 : if (visResamplerWt_p.null()) visResamplerWt_p = visResampler_p->clone();
1340 0 : visResamplerWt_p = visResampler_p;
1341 0 : visResamplerWt_p->setMaps(chanMap, polMap);
1342 0 : Array<Complex> gwts; Bool removeDegenerateAxis=false;
1343 0 : griddedWeights.get(gwts, removeDegenerateAxis);
1344 : // cerr << "initializeToSky for gwts" << endl;
1345 0 : visResamplerWt_p->initializeToSky(gwts, sumCFWeight);
1346 0 : }
1347 0 : }
1348 : //
1349 : //---------------------------------------------------------------
1350 : //
1351 0 : void AWProjectWBFT::finalizeToSky()
1352 : {
1353 0 : LogIO log_l(LogOrigin("AWProjectWBFT", "finalizeToSky[R&D]"));
1354 0 : AWProjectFT::finalizeToSky();
1355 : // The following commented code is the same as in the parent class
1356 : // call above.
1357 :
1358 : // if(isTiled)
1359 : // {
1360 : // AlwaysAssert(image, AipsError);
1361 : // AlwaysAssert(imageCache, AipsError);
1362 : // imageCache->flush();
1363 : // ostringstream o;
1364 : // imageCache->showCacheStatistics(o);
1365 : // log_l << o.str() << LogIO::POST;
1366 : // }
1367 :
1368 : // if(pointingToImage) delete pointingToImage; pointingToImage=0;
1369 :
1370 : // paChangeDetector.reset();
1371 : // cfCache_p->flush();
1372 : // visResampler_p->finalizeToSky(griddedData, sumWeight);
1373 :
1374 0 : if (!avgPBReady_p)
1375 : {
1376 0 : Array<Complex> gwts; Bool removeDegenerateAxis=false;
1377 0 : griddedWeights.get(gwts, removeDegenerateAxis);
1378 0 : visResamplerWt_p->finalizeToSky(gwts, sumCFWeight);
1379 0 : visResamplerWt_p->releaseBuffers();
1380 0 : }
1381 : /*
1382 : cerr << "Gridding run time = "
1383 : << " " << visResampler_p->runTimeG_p
1384 : << " " << visResampler_p->runTimeG1_p
1385 : << " " << visResampler_p->runTimeG2_p
1386 : << " " << visResampler_p->runTimeG3_p
1387 : << " " << visResampler_p->runTimeG4_p
1388 : << " " << visResampler_p->runTimeG5_p
1389 : << " " << visResampler_p->runTimeG6_p
1390 : << " " << visResampler_p->runTimeG7_p
1391 : << " C " << runTime1_p
1392 : << endl;
1393 : */
1394 0 : visResampler_p->runTimeG_p=visResampler_p->runTimeG1_p=visResampler_p->runTimeG2_p=visResampler_p->runTimeG3_p=visResampler_p->runTimeG4_p=visResampler_p->runTimeG5_p=visResampler_p->runTimeG6_p=visResampler_p->runTimeG7_p=0.0;
1395 0 : runTime1_p=0;
1396 0 : }
1397 : //
1398 : //---------------------------------------------------------------
1399 : //
1400 0 : void AWProjectWBFT::resampleCFToGrid(Array<Complex>& gwts,
1401 : VBStore& vbs, const VisBuffer& vb)
1402 : {
1403 : //
1404 : // Grid the weighted convolution function as well
1405 : //
1406 0 : LogIO log_l(LogOrigin("AWProjectFT", "resampleCFToGrid[R&D]"));
1407 : //
1408 : // Now rotate and put the rotated convolution weight function
1409 : // in rotatedCFWts_l object.
1410 : //
1411 :
1412 : // makeWBCFWt(*cfwts2_p, imRefFreq_p);
1413 :
1414 0 : timer_p.mark();
1415 0 : visResamplerWt_p->copy(*visResampler_p);
1416 :
1417 0 : Vector<Double> pointingOffset(convFuncCtor_p->findPointingOffset(*image, vb));
1418 : //cerr << "AWPWB: " << pointingOffset << endl;
1419 0 : visResamplerWt_p->makeVBRow2CFMap(*cfwts2_p,*convFuncCtor_p, vb,
1420 0 : paChangeDetector.getParAngleTolerance(),
1421 0 : chanMap,polMap,pointingOffset);
1422 0 : VBRow2CFBMapType& theMap=visResamplerWt_p->getVBRow2CFBMap();
1423 0 : convFuncCtor_p->prepareConvFunction(vb,theMap);
1424 0 : runTime1_p += timer_p.real();
1425 : //
1426 : // Set the uvw array to zero-sized array and dopsf=true.
1427 : // uvw.nelements()==0 is a hint to the re-sampler to put the
1428 : // gridded weights at the origin of the uv-grid. dopsf=true so
1429 : // that CF*Wts are accumulated (as against CF*Wts*Vis).
1430 : //
1431 : // Receive the sum-of-weights in a dummy array.
1432 0 : Matrix<Double> uvwOrigin;
1433 0 : vbs.uvw_p.reference(uvwOrigin);
1434 0 : Bool dopsf_l=true;
1435 0 : vbs.accumCFs_p=((vbs.uvw_p.nelements() == 0) && dopsf_l);
1436 :
1437 : // Array<Complex> gwts; Bool removeDegenerateAxis=false;
1438 : // wtsGrid.get(gwts, removeDegenerateAxis);
1439 0 : Int nDataChan = vbs.flagCube_p.shape()[1];
1440 :
1441 0 : vbs.startChan_p = 0; vbs.endChan_p = nDataChan;
1442 0 : visResamplerWt_p->DataToGrid(gwts, vbs, sumCFWeight, dopsf_l);
1443 0 : }
1444 : //
1445 : //---------------------------------------------------------------
1446 : //
1447 0 : void AWProjectWBFT::resampleDataToGrid(Array<Complex>& griddedData_l,
1448 : VBStore& vbs, const VisBuffer& vb,
1449 : Bool& dopsf)
1450 : {
1451 0 : AWProjectFT::resampleDataToGrid(griddedData_l,vbs,vb,dopsf);
1452 0 : if (!avgPBReady_p)
1453 : {
1454 : //
1455 : // Get a reference to the pixels of griddedWeights (a
1456 : // TempImage!)
1457 : //
1458 0 : Array<Complex> gwts; Bool removeDegenerateAxis=false;
1459 0 : griddedWeights.get(gwts, removeDegenerateAxis);
1460 0 : resampleCFToGrid(gwts, vbs, vb);
1461 0 : }
1462 0 : };
1463 : //
1464 : //---------------------------------------------------------------
1465 : //
1466 0 : void AWProjectWBFT::resampleDataToGrid(Array<DComplex>& griddedData_l,
1467 : VBStore& vbs, const VisBuffer& vb,
1468 : Bool& dopsf)
1469 : {
1470 0 : AWProjectFT::resampleDataToGrid(griddedData_l,vbs,vb,dopsf);
1471 0 : if (!avgPBReady_p)
1472 : {
1473 : //
1474 : // Get a reference to the pixels of griddedWeights (a
1475 : // TempImage!)
1476 : //
1477 0 : Array<Complex> gwts; Bool removeDegenerateAxis=false;
1478 0 : griddedWeights.get(gwts, removeDegenerateAxis);
1479 0 : resampleCFToGrid(gwts, vbs, vb);
1480 0 : }
1481 0 : };
1482 : // void AWProjectWBFT::resampleGridToData(VBStore& vbs, const VisBuffer& vb) {};
1483 :
1484 0 : void AWProjectWBFT::setCFCache(CountedPtr<CFCache>& cfc, const Bool resetCFC)
1485 : {
1486 0 : if (resetCFC) cfCache_p = cfc;
1487 0 : if (!cfCache_p.null())
1488 : {
1489 0 : cfs2_p = CountedPtr<CFStore2>(&cfCache_p->memCache2_p[0],false);//new CFStore2;
1490 0 : cfwts2_p = CountedPtr<CFStore2>(&cfCache_p->memCacheWt2_p[0],false);//new CFStore2;
1491 :
1492 : // cfCache_p->summarize(cfCache_p->memCache2_p,String("New CFC"));
1493 : // cfCache_p->summarize(cfCache_p->memCacheWt2_p,String(""));
1494 0 : avgPBReady_p=false;
1495 : }
1496 0 : }
1497 :
1498 : } //# NAMESPACE CASA - END
|
