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