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/TransformMachines2/VB2CFBMap.h>
30 : #include <synthesis/TransformMachines2/AWProjectWBFT.h>
31 : #include <synthesis/TransformMachines2/AWVisResampler.h>
32 : #include <synthesis/TransformMachines/StokesImageUtil.h>
33 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
34 : #include <casacore/scimath/Mathematics/FFTServer.h>
35 : #include <casacore/scimath/Mathematics/Convolver.h>
36 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
37 : #include <casacore/images/Images/ImageInterface.h>
38 : #include <casacore/images/Images/PagedImage.h>
39 : #include <msvis/MSVis/VisBuffer2.h>
40 : #include <casacore/casa/Arrays/Vector.h>
41 : #include <casacore/casa/Arrays/Slice.h>
42 : #include <casacore/casa/Arrays/ArrayMath.h>
43 : #include <casacore/casa/Arrays/Array.h>
44 : #include <casacore/casa/OS/HostInfo.h>
45 : #include <sstream>
46 : #include <casacore/casa/Utilities/CompositeNumber.h>
47 :
48 : using namespace casacore;
49 : namespace casa { //# NAMESPACE CASA - BEGIN
50 : namespace refim{
51 : //==================================================================================================
52 : // Various template instantiations
53 : //
54 : template
55 : void AWProjectWBFT::ftWeightImage(Lattice<casacore::Complex>& wtImage,
56 : const Matrix<Float>& sumWt,
57 : const Bool& doFFTNorm);
58 : template
59 : void AWProjectWBFT::makeWBSensitivityImage(Lattice<casacore::Complex>& wtImage,
60 : ImageInterface<Float>& sensitivityImage,
61 : const Matrix<Float>& sumWt,
62 : const Bool& doFFTNorm);
63 : //==================================================================================================
64 : template
65 : void AWProjectWBFT::ftWeightImage(Lattice<casacore::DComplex>& wtImage,
66 : const Matrix<Float>& sumWt,
67 : const Bool& doFFTNorm);
68 : template
69 : void AWProjectWBFT::makeWBSensitivityImage(Lattice<casacore::DComplex>& wtImage,
70 : ImageInterface<Float>& sensitivityImage,
71 : const Matrix<Float>& sumWt,
72 : const Bool& doFFTNorm);
73 : //==================================================================================================
74 : template
75 : void AWProjectWBFT::resampleCFToGrid(Array<casacore::Complex>& gwts,
76 : VBStore& vbs, const VisBuffer2& vb);
77 : template
78 : void AWProjectWBFT::resampleCFToGrid(Array<casacore::DComplex>& gwts,
79 : VBStore& vbs, const VisBuffer2& vb);
80 : //==================================================================================================
81 :
82 : //
83 : //---------------------------------------------------------------
84 : //
85 93 : AWProjectWBFT::AWProjectWBFT(Int nWPlanes, Long icachesize,
86 : CountedPtr<CFCache>& cfcache,
87 : CountedPtr<ConvolutionFunction>& cf,
88 : CountedPtr<VisibilityResamplerBase>& visResampler,
89 : Bool applyPointingOffset,
90 : vector<float> pointingOffsetSigDev,
91 : Bool doPBCorr,
92 : Int itilesize,
93 : Float paSteps,
94 : Float pbLimit,
95 : Bool usezero,
96 : Bool conjBeams,
97 93 : Bool doublePrecGrid)
98 : : AWProjectFT(nWPlanes,icachesize,cfcache,cf,visResampler,applyPointingOffset,pointingOffsetSigDev,doPBCorr, itilesize,pbLimit,usezero,conjBeams,doublePrecGrid),
99 : //avgPBReady_p(false),
100 93 : resetPBs_p(true),wtImageFTDone_p(false),fieldIds_p(0),rotatedCFWts_p(),visResamplerWt_p(nullptr),oneTimeMessage_p(false)
101 : {
102 : (void)paSteps;
103 : //
104 : // Set the function pointer for FORTRAN call for GCF services.
105 : // This is a pure C function pointer which FORTRAN can call.
106 : // The function itself uses GCF object services.
107 : //
108 93 : paChangeDetector.reset();
109 93 : pbLimit_p=pbLimit;
110 93 : if (applyPointingOffset) doPointing=1; else doPointing=0;
111 93 : maxConvSupport=-1;
112 : //
113 : // Set up the Conv. Func. disk cache manager object.
114 : //
115 : //use memory size defined in aipsrc if exists
116 93 : Long hostRAM = (HostInfo::memoryTotal(true)*1024); // In bytes
117 93 : hostRAM = hostRAM/(sizeof(Float)*2); // In complex pixels
118 93 : if (cachesize > hostRAM) cachesize=hostRAM;
119 :
120 93 : lastPAUsedForWtImg = MAGICPAVALUE;
121 93 : }
122 : //
123 : //---------------------------------------------------------------
124 : //
125 0 : AWProjectWBFT::AWProjectWBFT(const RecordInterface& stateRec)
126 0 : : AWProjectFT(stateRec), oneTimeMessage_p(false)
127 : {
128 0 : LogIO log_l(LogOrigin("AWProjectWBFT2", "AWProjectWBFT[R&D]"));
129 : // Construct from the input state record
130 :
131 0 : if (!fromRecord(stateRec))
132 0 : log_l << "Failed to create " << name() << " object." << LogIO::EXCEPTION;
133 :
134 0 : maxConvSupport=-1;
135 0 : visResampler_p->init(useDoubleGrid_p);
136 0 : }
137 : //
138 : //---------------------------------------------------------------
139 : //
140 159 : AWProjectWBFT& AWProjectWBFT::operator=(const AWProjectWBFT& other)
141 : {
142 159 : if(this!=&other)
143 : {
144 159 : AWProjectFT::operator=(other);
145 159 : padding_p = other.padding_p;
146 159 : nWPlanes_p = other.nWPlanes_p;
147 159 : imageCache = other.imageCache;
148 159 : cachesize = other.cachesize;
149 159 : tilesize = other.tilesize;
150 159 : gridder = other.gridder;
151 159 : isTiled = other.isTiled;
152 159 : lattice = other.lattice;
153 159 : maxAbsData = other.maxAbsData;
154 159 : centerLoc = other.centerLoc;
155 159 : offsetLoc = other.offsetLoc;
156 159 : pointingToImage = other.pointingToImage;
157 159 : usezero_p = other.usezero_p;
158 159 : doPBCorrection = other.doPBCorrection;
159 159 : maxConvSupport = other.maxConvSupport;
160 159 : avgPBReady_p= other.avgPBReady_p;
161 159 : resetPBs_p = other.resetPBs_p;
162 159 : wtImageFTDone_p = other.wtImageFTDone_p;
163 159 : rotatedCFWts_p = other.rotatedCFWts_p;
164 159 : visResamplerWt_p=other.visResamplerWt_p; // Copy the counted pointer
165 159 : oneTimeMessage_p = other.oneTimeMessage_p;
166 : };
167 159 : return *this;
168 : };
169 : //
170 : //----------------------------------------------------------------------
171 : //
172 0 : Int AWProjectWBFT::findPointingOffsets(const VisBuffer2& vb,
173 : Array<Float> &l_off,
174 : Array<Float> &m_off,
175 : Bool Evaluate)
176 : {
177 0 : LogIO log_l(LogOrigin("AWProjectWBFT2","findPointingOffsets[R&D]"));
178 0 : Int NAnt=0;
179 : //
180 : // This will return 0 if EPJ Table is not given. Otherwise will
181 : // return the number of antennas it detected (from the EPJ table)
182 : // and the offsets in l_off and m_off.
183 : //
184 0 : NAnt = AWProjectFT::findPointingOffsets(vb,l_off,m_off,Evaluate);
185 : //
186 : // Resize the offset arrays if no pointing table was given.
187 : //
188 0 : if (NAnt <=0 )
189 : {
190 0 : NAnt=vb.subtableColumns().antenna().nrow();
191 0 : l_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
192 0 : m_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
193 0 : l_off = m_off = 0.0;
194 : }
195 : //
196 : // Add field offsets to the pointing errors.
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 0 : offsets0 = offsets1 = 0.0;
231 : }
232 0 : for(Int i=0;i<NAnt;i++)
233 : {
234 0 : l_off(IPosition(3,0,0,i)) = l_off(IPosition(3,0,0,i)) + offsets0(i)/2.062642e+05;
235 0 : m_off(IPosition(3,0,0,i)) = m_off(IPosition(3,0,0,i)) + offsets1(i)/2.062642e+05;
236 : }
237 :
238 0 : if (firstPass==0)
239 : {
240 : // cout << (Float)(directionCoord.referenceValue()(0)) << " "
241 : // << (Float)(directionCoord.referenceValue()(1)) << endl;
242 : // cout << vb.direction1()(0).getAngle().getValue()(0) << " "
243 : // << vb.direction1()(0).getAngle().getValue()(1) << endl;
244 : // cout << l_off << endl;
245 : // cout << m_off << endl;
246 : }
247 0 : firstPass++;
248 0 : return NAnt;
249 0 : }
250 : //
251 : //---------------------------------------------------------------
252 : //
253 : // void AWProjectWBFT::normalizeAvgPB()
254 : // {
255 : // LogIO log_l(LogOrigin("AWProjectWBFT2","normalizeAvgPB[R&D]"));
256 : // // We accumulated normalized PBs. So don't normalize the average
257 : // // PB.
258 : // pbNormalized_p = false;
259 :
260 : // }
261 : //
262 : //---------------------------------------------------------------
263 : // For wide-band case, this just tells the user that the sensitivity
264 : // image will be computed during the first gridding cycle.
265 : //
266 : // Weights image (sum of convolution functions) is accumulated
267 : // during gridding. At the end of the gridding cycle, the weight
268 : // image is FT'ed (using ftWeightImage()) and properly normalized
269 : // and converted to sensitivity image (using
270 : // makeSensitivityImage()). These methods are triggred in the first
271 : // call to getImage(). In subsequent calls to getImage() these
272 : // calls are NoOps.
273 : //
274 59 : void AWProjectWBFT::makeSensitivityImage(const VisBuffer2&,
275 : const ImageInterface<Complex>& /*Complex imageTemplate*/,
276 : ImageInterface<Float>& /*sensitivityImage*/)
277 : {
278 59 : if (avgPBReady_p == false)
279 : {
280 59 : if (oneTimeMessage_p == false)
281 : {
282 118 : LogIO log_l(LogOrigin("AWProjectWBFT2", "makeSensitivityImage(Complex)[R&D]"));
283 59 : log_l << "Setting up for weights accumulation ";
284 59 : if (sensitivityPatternQualifierStr_p != "") log_l << "for term " << sensitivityPatternQualifier_p << " ";
285 : log_l << "during gridding to compute sensitivity pattern."
286 : << endl
287 : << "Consequently, the first gridding cycle will be slower than the subsequent ones."
288 59 : << LogIO::WARN;
289 59 : oneTimeMessage_p=true;
290 59 : }
291 : }
292 59 : }
293 0 : void AWProjectWBFT::makeSensitivityImage(const VisBuffer2&,
294 : const ImageInterface<DComplex>& /* DComplex imageTemplate*/,
295 : ImageInterface<Float>& /*sensitivityImage*/)
296 : {
297 0 : if (avgPBReady_p == false)
298 : {
299 0 : if (oneTimeMessage_p == false)
300 : {
301 0 : LogIO log_l(LogOrigin("AWProjectWBFT2", "makeSensitivityImage(DComplex)[R&D]"));
302 0 : log_l << "Setting up for weights accumulation ";
303 0 : if (sensitivityPatternQualifierStr_p != "") log_l << "for term " << sensitivityPatternQualifier_p << " ";
304 : log_l << "during gridding to compute sensitivity pattern."
305 : << endl
306 : << "Consequently, the first gridding cycle will be slower than the subsequent ones."
307 0 : << LogIO::WARN;
308 0 : oneTimeMessage_p=true;
309 0 : }
310 : }
311 0 : }
312 : //
313 : //---------------------------------------------------------------
314 : //
315 : template <class T>
316 52 : void AWProjectWBFT::ftWeightImage(Lattice<T>& wtImage,
317 : const Matrix<Float>& sumWt,
318 : const Bool& doFFTNorm)
319 : {
320 104 : LogIO log_l(LogOrigin("AWProjectWBFT2", "ftWeightImage[R&D]"));
321 52 : if (wtImageFTDone_p) return;
322 :
323 52 : if ((sumWt.shape().nelements() < 2) ||
324 208 : (sumWt.shape()(0) != wtImage.shape()(2)) ||
325 104 : (sumWt.shape()(1) != wtImage.shape()(3)))
326 0 : log_l << "Sum of weights per poln and chan is required" << LogIO::EXCEPTION;
327 : //
328 : // Use only the amplitude of the gridded weights. LatticeExpr
329 : // classes while useful, appear to be too strict about types (the
330 : // code below will not compile if the abs(wtImage) is not
331 : // converted back to a complex type).
332 :
333 52 : LatticeFFT::cfft2d(wtImage,false);
334 52 : wtImageFTDone_p=true;
335 :
336 52 : Int sizeX=wtImage.shape()(0), sizeY=wtImage.shape()(1);
337 :
338 52 : Array<T> wtBuf; wtImage.get(wtBuf,false);
339 52 : ArrayLattice<T> wtLat(wtBuf,true);
340 : //
341 : // Copy one 2D plane at a time, normalizing by the sum of weights
342 : // and possibly 2D FFT.
343 : //
344 : // Set up Lattice iterators on wtImage and sensitivityImage
345 : //
346 52 : IPosition axisPath(4, 0, 1, 2, 3);
347 52 : IPosition cursorShape(4, sizeX, sizeY, 1, 1);
348 :
349 52 : LatticeStepper wtImStepper(wtImage.shape(), cursorShape, axisPath);
350 52 : LatticeIterator<T> wtImIter(wtImage, wtImStepper);
351 : //
352 : // Iterate over channel and polarization axis
353 : //
354 52 : if (!doFFTNorm) sizeX=sizeY=1;
355 : //
356 : // Normalize each frequency and polarization plane of the complex
357 : // sensitivity pattern
358 : //
359 : // getSumOfCFWeights() returns the sum-of-weight appropriate for
360 : // computing sensitivity patterns (which could be different from
361 : // the data-sum-of-weights.
362 : //
363 :
364 130 : for(wtImIter.reset(); !wtImIter.atEnd(); wtImIter++)
365 : {
366 78 : wtImIter.rwCursor() = (wtImIter.rwCursor()
367 78 : *Float(sizeX)*Float(sizeY));
368 : }
369 52 : }
370 : //
371 : //---------------------------------------------------------------
372 : //
373 : template <class T>
374 169 : void AWProjectWBFT::makeWBSensitivityImage(Lattice<T>& wtImage,
375 : ImageInterface<Float>& sensitivityImage,
376 : const Matrix<Float>& sumWt,
377 : const Bool& doFFTNorm)
378 : {
379 169 : if (avgPBReady_p) return;
380 104 : LogIO log_l(LogOrigin("AWProjectWBFT2", "makeWBSensitivityImage[R&D]"));
381 :
382 52 : ftWeightImage(wtImage, sumWt, doFFTNorm);
383 :
384 52 : if (useDoubleGrid_p)
385 : {
386 52 : sensitivityImage.resize(griddedWeights_D.shape());
387 52 : sensitivityImage.setCoordinateInfo(griddedWeights_D.coordinates());
388 : }
389 : else
390 : {
391 0 : sensitivityImage.resize(griddedWeights.shape());
392 0 : sensitivityImage.setCoordinateInfo(griddedWeights.coordinates());
393 : }
394 : //
395 : // Rest of the code below is to (1) average the poln planes of the
396 : // weight image for each freq. channel, and (2) copy the averaged
397 : // value to the all the poln planes for each freq. channels of the
398 : // sensitivity image.
399 : //
400 : // Set up Lattice iterator on wtImage for averaging all poln
401 : // planes and writing the result the only the first poln plane.
402 : //
403 52 : IPosition axisPath(2, 2, 3); // Step through the Poln (2) and Freq (3) axis.
404 104 : IPosition cursorShapeWt(4, wtImage.shape()(0), wtImage.shape()(1), 1, 1);
405 :
406 52 : LatticeStepper wtImStepper(wtImage.shape(), cursorShapeWt, axisPath);
407 52 : LatticeIterator<T> wtImIter(wtImage, wtImStepper);
408 52 : int nPolPlanes=wtImage.shape()(2);
409 : // First average all poln planes and write the result back to the first poln plane.
410 208 : for(wtImIter.reset(); !wtImIter.atEnd(); /*increment is inside the for-loop below*/)
411 : {
412 78 : Matrix<T> tmpWt(wtImIter.rwMatrixCursor().shape());
413 78 : ArrayLattice<T> tt(tmpWt,true);
414 156 : for(int i=0;i<nPolPlanes;i++)
415 : {
416 78 : Matrix<T> tmp0_ref;tmp0_ref.reference(wtImIter.rwMatrixCursor().nonDegenerate());
417 78 : ArrayLattice<T> p0(tmp0_ref,true); //Make ArrayLattice from Array by refence
418 156 : LatticeExpr<T> le(tt+p0);
419 78 : tt.copyData(le);
420 78 : wtImIter++;
421 : }
422 156 : LatticeExpr<T> le(tt/nPolPlanes);
423 78 : tt.copyData(le);
424 : }
425 :
426 104 : IPosition cursorShape(4, sensitivityImage.shape()(0), sensitivityImage.shape()(1), 1, 1);
427 52 : LatticeStepper senImStepper(sensitivityImage.shape(), cursorShape, axisPath);
428 52 : LatticeIterator<Float> senImIter(sensitivityImage, senImStepper);
429 : //
430 : // Copy the real part of the average of all planes of the wtImage
431 : // to all the planes of the sensitivity image (senImage).
432 : //
433 208 : for(wtImIter.reset(),senImIter.reset(); !wtImIter.atEnd();/*increment is inside the first for-loop below*/)
434 : {
435 : // Extract the first polarization plane from the weight image.
436 78 : Matrix<T> wt_ref;wt_ref.reference(wtImIter.rwMatrixCursor().nonDegenerate());
437 78 : ArrayLattice<T> p0(wt_ref,true);
438 156 : for (int i=0;i<nPolPlanes;i++) wtImIter++; // Skip the rest of the poln planes
439 :
440 156 : for(int i=0;i<sensitivityImage.shape()(2);i++)
441 : {
442 : // // Now replace the polarization planes with their average. REVIEW THIS FOR FM CASE.
443 78 : Matrix<Float> tmp;tmp.reference(senImIter.rwMatrixCursor().nonDegenerate());
444 78 : ArrayLattice<Float> tt(tmp,true); //Make ArrayLattice from Array by refence
445 156 : LatticeExpr<Float> le(abs(p0));
446 78 : tt.copyData(le);
447 :
448 78 : senImIter++;
449 : }
450 : }
451 :
452 52 : if (tt_pp == "")
453 52 : cfCache_p->flush(sensitivityImage,sensitivityPatternQualifierStr_p);
454 :
455 52 : pbNormalized_p=false;
456 52 : resetPBs_p=false;
457 :
458 52 : avgPBReady_p=true;
459 52 : }
460 : //
461 : //---------------------------------------------------------------
462 : //
463 : // Convert the gridded visibilities to image. This implements the
464 : // basic 2D transform (vC-Z Theorem) using the FFT.
465 : //
466 : // This specialization of getImage() exists since the sensitivity
467 : // pattern in this FTMachine is computed as the FT of the gridded
468 : // weights (convolution functions). The gridded weights are
469 : // available along with the gridded data at the end of the gridding
470 : // cycle. This method first converts the gridded weights to
471 : // sensitivity image and then calls AWProjectFT::getImage(), which
472 : // in-turn calls normalizeImage() with the sensitivityImage and the
473 : // sum of weights.
474 : //
475 169 : ImageInterface<Complex>& AWProjectWBFT::getImage(Matrix<Float>& weights,
476 : Bool fftNormalization)
477 : {
478 169 : AlwaysAssert(image, AipsError);
479 338 : LogIO log_l(LogOrigin("AWProjectWBFT2", "getImage[R&D]"));
480 :
481 169 : weights.resize(sumWeight.shape());
482 169 : convertArray(weights, sumWeight);//I suppose this converts a
483 : //Matrix<Double> (sumWeights) to
484 : //Matrix<Float> (weights). Why
485 : //is this conversion required?
486 : //--SB (Dec. 2010)
487 : //
488 : // Compute the sensitivity image as FT of the griddedWeigths.
489 : // Return the result in avgPB_p. Cache it in the cfCache_p
490 : // object. And raise the avgPBReady_p flag so that this becomes a
491 : // null call the next time around.
492 169 : if (useDoubleGrid_p)
493 169 : makeWBSensitivityImage(griddedWeights_D, *avgPB_p, weights, true);
494 : else
495 0 : makeWBSensitivityImage(griddedWeights, *avgPB_p, weights, true);
496 : //
497 : // This calls the overloadable method normalizeImage() which
498 : // normalizes the raw image by the sensitivty pattern (avgPB_p).
499 : //
500 169 : AWProjectFT::getImage(weights,fftNormalization);
501 :
502 169 : return *image;
503 169 : }
504 : //
505 : //---------------------------------------------------------------
506 : //
507 : // Initialize the FFT to the Sky. Here we have to setup and
508 : // initialize the grid.
509 : //
510 176 : void AWProjectWBFT::initializeToSky(ImageInterface<Complex>& iimage,
511 : Matrix<Float>& weight,
512 : const VisBuffer2& vb)
513 : {
514 352 : LogIO log_l(LogOrigin("AWProjectWBFT2","initializeToSky[R&D]"));
515 176 : AWProjectFT::initializeToSky(iimage,weight,vb);
516 :
517 176 : if (resetPBs_p)
518 : {
519 117 : if (useDoubleGrid_p)
520 : {
521 117 : griddedWeights_D.resize(iimage.shape());
522 117 : griddedWeights_D.setCoordinateInfo(iimage.coordinates());
523 117 : griddedWeights_D.set(0.0);
524 117 : pbPeaks.resize(griddedWeights_D.shape()(2));
525 117 : pbPeaks.set(0.0);
526 : }
527 : else
528 : {
529 0 : griddedWeights.resize(iimage.shape());
530 0 : griddedWeights.setCoordinateInfo(iimage.coordinates());
531 0 : griddedWeights.set(0.0);
532 0 : pbPeaks.resize(griddedWeights.shape()(2));
533 0 : pbPeaks.set(0.0);
534 : }
535 :
536 117 : resetPBs_p=false;
537 : }
538 :
539 176 : std::tuple<int, double>cubeinfo(1,-1.0);
540 : double freqofBegChan;
541 176 : spectralCoord_p.toWorld(freqofBegChan, 0.0);
542 :
543 176 : cubeinfo=std::make_tuple(iimage.shape()(3),freqofBegChan);
544 :
545 : ///load AVGPB is quite the memory consumer for cubes as it will load the whole cube in memory a couple of times even.
546 : //cerr << "###Avoiding loading of avgPB " << avgPBReady_p << endl;
547 : ////TESTOO need to oveload this init in HPG
548 176 : if(!avgPBReady_p)
549 87 : avgPBReady_p = (cfCache_p->loadAvgPB(avgPB_p,sensitivityPatternQualifierStr_p, cubeinfo) != CFDefs::NOTCACHED);
550 :
551 176 : if(avgPBReady_p){
552 130 : LatticeExprNode le( max( *avgPB_p ) );
553 130 : Float avgPB_max=le.getFloat();
554 :
555 130 : if(avgPB_max <= 0.0) avgPBReady_p = false;
556 130 : }
557 : // Need to grid the weighted Convolution Functions to make the sensitivity pattern.
558 176 : if (!avgPBReady_p)
559 : {
560 : // Make a copy of the re-sampler and set it up.
561 59 : if (visResamplerWt_p.null()) visResamplerWt_p = visResampler_p->clone();
562 59 : visResamplerWt_p = visResampler_p;
563 59 : visResamplerWt_p->setMaps(chanMap, polMap);
564 59 : if (useDoubleGrid_p)
565 : {
566 59 : Array<DComplex> gwts; Bool removeDegenerateAxis=false;
567 59 : griddedWeights_D.get(gwts, removeDegenerateAxis);
568 59 : visResamplerWt_p->initializeToSky(gwts, sumCFWeight); //A NoOp right now.
569 59 : }
570 : else
571 : {
572 0 : Array<Complex> gwts; Bool removeDegenerateAxis=false;
573 0 : griddedWeights.get(gwts, removeDegenerateAxis);
574 0 : visResamplerWt_p->initializeToSky(gwts, sumCFWeight); //A NoOp right now.
575 0 : }
576 : }
577 176 : }
578 : //
579 : //---------------------------------------------------------------
580 : //
581 169 : void AWProjectWBFT::finalizeToSky()
582 : {
583 338 : LogIO log_l(LogOrigin("AWProjectWBFT2", "finalizeToSky[R&D]"));
584 169 : AWProjectFT::finalizeToSky();
585 :
586 :
587 169 : if(name()=="AWProjectWBFTHPG")
588 0 : return;
589 : ////
590 169 : if(!visResamplerWt_p)
591 82 : return;
592 :
593 : //
594 : // SB: Not sure why the follwing if-stmt (not by SB) needs to be
595 : // in this class. This be simply avoided by overloading this
596 : // method in the AWProjectWBFTHPG class making the NoOp obvious
597 : // for that class, and easily extendable & maintainable.
598 : //
599 :
600 :
601 87 : if (!avgPBReady_p)
602 : {
603 52 : if (useDoubleGrid_p)
604 : {
605 52 : Array<DComplex> gwts; Bool removeDegenerateAxis=false;
606 52 : griddedWeights_D.get(gwts, removeDegenerateAxis);
607 52 : visResamplerWt_p->finalizeToSky(gwts, sumCFWeight);
608 52 : visResamplerWt_p->releaseBuffers();
609 52 : }
610 : else
611 : {
612 0 : Array<Complex> gwts; Bool removeDegenerateAxis=false;
613 0 : griddedWeights.get(gwts, removeDegenerateAxis);
614 0 : visResamplerWt_p->finalizeToSky(gwts, sumCFWeight);
615 0 : visResamplerWt_p->releaseBuffers();
616 0 : }
617 : }
618 :
619 174 : visResampler_p->runTimeG_p
620 261 : =visResampler_p->runTimeG1_p
621 87 : =visResampler_p->runTimeG2_p
622 87 : =visResampler_p->runTimeG3_p
623 87 : =visResampler_p->runTimeG4_p
624 87 : =visResampler_p->runTimeG5_p
625 87 : =visResampler_p->runTimeG6_p
626 174 : =visResampler_p->runTimeG7_p
627 87 : =0.0;
628 87 : runTime1_p=0;
629 169 : }
630 : //
631 : //---------------------------------------------------------------
632 : // Method sets up the inputs for the VisResampler::DataToGrid()
633 : // method for gridding the CFs at the origin of the UVW-plane. This
634 : // includes application of pointing offset correction and any other
635 : // term that might be included in VB2CFBMap::makeVBRow2CFBMap()
636 : // method.
637 : template <class T>
638 7150 : void AWProjectWBFT::resampleCFToGrid(Array<T>& gwts,
639 : VBStore& vbs, const VisBuffer2& vb)
640 : {
641 7150 : visResamplerWt_p->copy(*visResampler_p);
642 :
643 7150 : po_p->fetchPointingOffset(*image, vb, doPointing);
644 14300 : vb2CFBMap_p->makeVBRow2CFBMap(*cfwts2_p,vb,
645 : paChangeDetector.getParAngleTolerance(),
646 7150 : chanMap,polMap,po_p);
647 7150 : convFuncCtor_p->prepareConvFunction(vb,*vb2CFBMap_p);
648 : //
649 : // Set the uvw array to zero-sized array and dopsf=true.
650 : // uvw.nelements()==0 is a hint to the re-sampler to put the
651 : // gridded weights at the origin of the uv-grid. dopsf=true so
652 : // that CF*Wts are accumulated (as against CF*Wts*Vis).
653 : //
654 : // Receive the sum-of-weights in a dummy array.
655 7150 : Matrix<Double> uvwOrigin;
656 7150 : vbs.uvw_p.reference(uvwOrigin);
657 7150 : Bool dopsf_l=true;
658 7150 : vbs.accumCFs_p=((vbs.uvw_p.nelements() == 0) && dopsf_l);
659 7150 : vbs.ftmType_p=casa::refim::FTMachine::WEIGHT;
660 7150 : Int nDataChan = vbs.flagCube_p.shape()[1];
661 :
662 7150 : visResamplerWt_p->setVB2CFMap(vb2CFBMap_p);
663 7150 : vbs.startChan_p = 0; vbs.endChan_p = nDataChan;
664 :
665 7150 : visResamplerWt_p->DataToGrid(gwts, vbs, sumCFWeight, dopsf_l);
666 7150 : }
667 : //
668 : //---------------------------------------------------------------
669 : // Methods to accumulate data on the grid. If necessary, also
670 : // trigger the accumulation of CFs via ::resampleCFToGrid() method.
671 : //
672 0 : void AWProjectWBFT::resampleDataToGrid(Array<Complex>& griddedData_l,
673 : VBStore& vbs, const VisBuffer2& vb,
674 : Bool& dopsf)
675 : {
676 0 : AWProjectFT::resampleDataToGrid(griddedData_l,vbs,vb,dopsf);
677 0 : if (!avgPBReady_p)
678 : {
679 : //
680 : // Get a reference to the pixels of griddedWeights (a
681 : // TempImage!)
682 : //
683 0 : Array<Complex> gwts; Bool removeDegenerateAxis=false;
684 0 : griddedWeights.get(gwts, removeDegenerateAxis);
685 0 : resampleCFToGrid(gwts, vbs, vb);
686 0 : }
687 0 : };
688 : //
689 : //---------------------------------------------------------------
690 : //
691 20933 : void AWProjectWBFT::resampleDataToGrid(Array<DComplex>& griddedData_l,
692 : VBStore& vbs, const VisBuffer2& vb,
693 : Bool& dopsf)
694 : {
695 20933 : AWProjectFT::resampleDataToGrid(griddedData_l,vbs,vb,dopsf);
696 :
697 20933 : if (!avgPBReady_p)
698 : {
699 : //
700 : // Get a reference to the pixels of griddedWeights (a
701 : // TempImage!)
702 : //
703 7150 : Array<DComplex> gwts; Bool removeDegenerateAxis=false;
704 7150 : griddedWeights_D.get(gwts, removeDegenerateAxis);
705 7150 : resampleCFToGrid(gwts, vbs, vb);
706 7150 : }
707 20933 : };
708 :
709 : ////-----------------------------------------------------------------
710 : // SB: Is this functionality not served by the
711 : // ::resampleCFToGrid() method above? Recommend that that be used
712 : // (or modifed as required) to ensure all mathematical terms are
713 : // included.
714 : //
715 0 : void AWProjectWBFT::gridImgWeights(const VisBuffer2& vb)
716 : {
717 0 : findConvFunction(*image, vb);
718 : //cerr << "IN gridWTImage" << endl;
719 0 : if(avgPBReady_p)
720 0 : return;
721 : else
722 0 : avgPB_p=nullptr; ///make sure it is not pointing to anything
723 :
724 : // This does WEIGHTs gridding via the framework route
725 : // (set ftmType and use ::put() for *all* gridding)
726 0 : ftmType_p=casa::refim::FTMachine::WEIGHT;
727 0 : put(vb,-1,false);
728 0 : return;
729 :
730 : // try
731 : // {
732 : // findConvFunction(*image, vb);
733 : // if(avgPBReady_p)
734 : // return;
735 : // }
736 : // catch(AipsError& x)
737 : // {
738 : // LogIO log_l(LogOrigin("AWProjectFT2", "put[R&D]"));
739 : // log_l << x.getMesg() << LogIO::WARN;
740 : // return;
741 : // }
742 : // Matrix<Float> imagingweight;
743 : // getImagingWeight(imagingweight, vb);
744 : // matchChannel(vb);
745 : // Cube<Complex> data;
746 : // //Fortran gridder need the flag as ints
747 : // Cube<Int> flags;
748 : // Matrix<Float> elWeight;
749 :
750 : // interpolateFrequencyTogrid(vb, imagingweight,data, flags , elWeight, FTMachine::PSF);
751 : // Matrix<Double> uvw(negateUV(vb));
752 : // Vector<Double> dphase(vb.nRows());
753 : // dphase=0.0;
754 : // doUVWRotation_p=true;
755 : // VBStore vbs;
756 : // Vector<Int> gridShape = griddedData2.shape().asVector();
757 : // setupVBStore(vbs,vb, elWeight,data,uvw,flags, dphase, True /*dopsf*/,gridShape);
758 : // if (useDoubleGrid_p){
759 : // Array<DComplex> gwts; Bool removeDegenerateAxis=false;
760 : // griddedWeights_D.get(gwts, removeDegenerateAxis);
761 : // resampleCFToGrid(gwts, vbs, vb);
762 :
763 : // }
764 : // else{
765 : // Array<Complex> gwts; Bool removeDegenerateAxis=false;
766 : // griddedWeights.get(gwts, removeDegenerateAxis);
767 : // resampleCFToGrid(gwts, vbs, vb);
768 :
769 : // }
770 : }
771 :
772 93 : void AWProjectWBFT::setCFCache(CountedPtr<CFCache>& cfc, const Bool resetCFC)
773 : {
774 93 : if (resetCFC) cfCache_p = cfc;
775 93 : if (!cfCache_p.null())
776 : {
777 93 : cfs2_p = CountedPtr<CFStore2>(&cfCache_p->memCache2_p[0],false);//new CFStore2;
778 93 : cfwts2_p = CountedPtr<CFStore2>(&cfCache_p->memCacheWt2_p[0],false);//new CFStore2;
779 :
780 93 : avgPBReady_p=false;
781 : }
782 93 : }
783 :
784 : } //# NAMESPACE CASA - END
785 : };
|
