Line data Source code
1 : //# FTMachine.cc: Implementation of FTMachine class
2 : //# Copyright (C) 1997-2014
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 : #include <cmath>
28 : #include <msvis/MSVis/VisibilityIterator.h>
29 : #include <casacore/casa/Quanta/Quantum.h>
30 : #include <casacore/casa/Quanta/UnitMap.h>
31 : #include <casacore/casa/Quanta/UnitVal.h>
32 : #include <casacore/measures/Measures/Stokes.h>
33 : #include <casacore/casa/Quanta/Euler.h>
34 : #include <casacore/casa/Quanta/RotMatrix.h>
35 : #include <casacore/measures/Measures/MFrequency.h>
36 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
37 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
39 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
40 : #include <casacore/coordinates/Coordinates/Projection.h>
41 : #include <casacore/ms/MeasurementSets/MSColumns.h>
42 : #include <casacore/casa/BasicSL/Constants.h>
43 : #include <synthesis/TransformMachines/FTMachine.h>
44 : #include <casacore/scimath/Mathematics/RigidVector.h>
45 : #include <msvis/MSVis/StokesVector.h>
46 : #include <synthesis/TransformMachines/StokesImageUtil.h>
47 : #include <synthesis/TransformMachines/VisModelData.h>
48 : #include <synthesis/TransformMachines/Utils.h>
49 : #include <msvis/MSVis/VisBuffer.h>
50 : #include <msvis/MSVis/VisSet.h>
51 : #include <casacore/images/Images/ImageInterface.h>
52 : #include <casacore/images/Images/PagedImage.h>
53 : #include <casacore/images/Images/ImageUtilities.h>
54 : #include <casacore/casa/Containers/Block.h>
55 : #include <casacore/casa/Containers/Record.h>
56 : #include <casacore/casa/Arrays/ArrayIter.h>
57 : #include <casacore/casa/Arrays/ArrayLogical.h>
58 : #include <casacore/casa/Arrays/ArrayMath.h>
59 : #include <casacore/casa/Arrays/MatrixMath.h>
60 : #include <casacore/casa/Arrays/MaskedArray.h>
61 : #include <casacore/casa/Arrays/Array.h>
62 : #include <casacore/casa/Arrays/Vector.h>
63 : #include <casacore/casa/Arrays/Matrix.h>
64 : #include <casacore/casa/Arrays/MatrixIter.h>
65 : #include <casacore/casa/BasicSL/String.h>
66 : #include <casacore/casa/Utilities/Assert.h>
67 : #include <casacore/casa/Utilities/BinarySearch.h>
68 : #include <casacore/casa/Exceptions/Error.h>
69 : #include <casacore/scimath/Mathematics/NNGridder.h>
70 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
71 : #include <casacore/measures/Measures/UVWMachine.h>
72 : #include <synthesis/TransformMachines/CFStore2.h>
73 :
74 : #include <casacore/casa/System/ProgressMeter.h>
75 :
76 : #include <casacore/casa/OS/Timer.h>
77 : #include <sstream>
78 : #include <iostream>
79 :
80 : using namespace casacore;
81 : namespace casa { //# NAMESPACE CASA - BEGIN
82 :
83 130 : FTMachine::FTMachine() : isDryRun(false), image(0), uvwMachine_p(0),
84 130 : tangentSpecified_p(false), fixMovingSource_p(false),
85 65 : distance_p(0.0), lastFieldId_p(-1),lastMSId_p(-1),
86 65 : useDoubleGrid_p(false),
87 65 : freqFrameValid_p(false),
88 65 : freqInterpMethod_p(InterpolateArray1D<Double,Complex>::nearestNeighbour),
89 65 : pointingDirCol_p("DIRECTION"),
90 65 : cfStokes_p(), cfCache_p(), cfs_p(), cfwts_p(), cfs2_p(), cfwts2_p(),
91 130 : canComputeResiduals_p(false), toVis_p(true), numthreads_p(-1),pbLimit_p(0.05),
92 195 : sj_p(0),cmplxImage_p(), phaseCenterTime_p(-1.0)
93 : {
94 65 : spectralCoord_p=SpectralCoordinate();
95 65 : isIOnly=false;
96 65 : spwChanSelFlag_p=0;
97 65 : polInUse_p=0;
98 65 : pop_p = new PolOuterProduct;
99 : // cerr << "Called FTMachine()" << endl;
100 65 : }
101 :
102 0 : FTMachine::FTMachine(CountedPtr<CFCache>& cfcache,CountedPtr<ConvolutionFunction>& cf):
103 0 : isDryRun(false), image(0), uvwMachine_p(0),
104 0 : tangentSpecified_p(false), fixMovingSource_p(false),
105 0 : distance_p(0.0), lastFieldId_p(-1),lastMSId_p(-1),
106 0 : useDoubleGrid_p(false),
107 0 : freqFrameValid_p(false),
108 0 : freqInterpMethod_p(InterpolateArray1D<Double,Complex>::nearestNeighbour),
109 0 : pointingDirCol_p("DIRECTION"),
110 0 : cfStokes_p(), cfCache_p(cfcache), cfs_p(), cfwts_p(), cfs2_p(), cfwts2_p(),
111 0 : convFuncCtor_p(cf),canComputeResiduals_p(false), toVis_p(true), numthreads_p(-1),
112 0 : pbLimit_p(0.05),sj_p(0),cmplxImage_p( ), phaseCenterTime_p(-1.0)
113 : {
114 0 : spectralCoord_p=SpectralCoordinate();
115 0 : isIOnly=false;
116 0 : spwChanSelFlag_p=0;
117 0 : polInUse_p=0;
118 0 : pop_p = new PolOuterProduct;
119 : //cerr << "Called FTMachine(CPT<CFCi>),...)" << endl;
120 0 : }
121 :
122 163 : LogIO& FTMachine::logIO() {return logIO_p;};
123 :
124 : //----------------------------------------------------------------------
125 28 : FTMachine& FTMachine::operator=(const FTMachine& other)
126 : {
127 28 : if(this!=&other) {
128 28 : image=other.image;
129 : //generic selection stuff and state
130 28 : nAntenna_p=other.nAntenna_p;
131 28 : distance_p=other.distance_p;
132 28 : lastFieldId_p=other.lastFieldId_p;
133 28 : lastMSId_p=other.lastMSId_p;
134 :
135 28 : tangentSpecified_p=other.tangentSpecified_p;
136 28 : mTangent_p=other.mTangent_p;
137 28 : mImage_p=other.mImage_p;
138 28 : mFrame_p=other.mFrame_p;
139 :
140 28 : nx=other.nx;
141 28 : ny=other.ny;
142 28 : npol=other.npol;
143 28 : nchan=other.nchan;
144 28 : nvischan=other.nvischan;
145 28 : nvispol=other.nvispol;
146 28 : mLocation_p=other.mLocation_p;
147 28 : if(uvwMachine_p)
148 0 : delete uvwMachine_p;
149 28 : if(other.uvwMachine_p)
150 0 : uvwMachine_p=new UVWMachine(*other.uvwMachine_p);
151 : else
152 28 : uvwMachine_p=0;
153 28 : doUVWRotation_p=other.doUVWRotation_p;
154 : //Spectral and pol stuff
155 28 : freqInterpMethod_p=other.freqInterpMethod_p;
156 28 : spwChanSelFlag_p.resize();
157 28 : spwChanSelFlag_p=other.spwChanSelFlag_p;
158 28 : freqFrameValid_p=other.freqFrameValid_p;
159 28 : selectedSpw_p.resize();
160 28 : selectedSpw_p=other.selectedSpw_p;
161 28 : imageFreq_p.resize();
162 28 : imageFreq_p=other.imageFreq_p;
163 28 : lsrFreq_p.resize();
164 28 : lsrFreq_p=other.lsrFreq_p;
165 28 : interpVisFreq_p.resize();
166 28 : interpVisFreq_p=other.interpVisFreq_p;
167 28 : multiChanMap_p=other.multiChanMap_p;
168 28 : chanMap.resize();
169 28 : chanMap=other.chanMap;
170 28 : polMap.resize();
171 28 : polMap=other.polMap;
172 28 : nVisChan_p.resize();
173 28 : nVisChan_p=other.nVisChan_p;
174 28 : spectralCoord_p=other.spectralCoord_p;
175 28 : doConversion_p.resize();
176 28 : doConversion_p=other.doConversion_p;
177 28 : pointingDirCol_p=other.pointingDirCol_p;
178 : //moving source stuff
179 28 : movingDir_p=other.movingDir_p;
180 28 : fixMovingSource_p=other.fixMovingSource_p;
181 28 : firstMovingDir_p=other.firstMovingDir_p;
182 :
183 : //Double precision gridding for those FTMachines that can do
184 28 : useDoubleGrid_p=other.useDoubleGrid_p;
185 28 : cfStokes_p = other.cfStokes_p;
186 28 : polInUse_p = other.polInUse_p;
187 28 : cfs_p = other.cfs_p;
188 28 : cfwts_p = other.cfwts_p;
189 28 : cfs2_p = other.cfs2_p;
190 28 : cfwts2_p = other.cfwts2_p;
191 28 : canComputeResiduals_p = other.canComputeResiduals_p;
192 :
193 28 : pop_p = other.pop_p;
194 28 : toVis_p = other.toVis_p;
195 28 : spwFreqSel_p.resize();
196 28 : spwFreqSel_p = other.spwFreqSel_p;
197 28 : expandedSpwFreqSel_p = other.expandedSpwFreqSel_p;
198 28 : expandedSpwConjFreqSel_p = other.expandedSpwConjFreqSel_p;
199 28 : cmplxImage_p=other.cmplxImage_p;
200 28 : numthreads_p=other.numthreads_p;
201 28 : pbLimit_p=other.pbLimit_p;
202 28 : convFuncCtor_p = other.convFuncCtor_p;
203 28 : sj_p.resize();
204 28 : sj_p=other.sj_p;
205 28 : isDryRun=other.isDryRun;
206 28 : phaseCenterTime_p=other.phaseCenterTime_p;
207 : };
208 28 : return *this;
209 : };
210 :
211 :
212 0 : FTMachine* FTMachine::cloneFTM(){
213 0 : Record rec;
214 0 : String err;
215 0 : if(!(this->toRecord(err, rec)))
216 0 : throw(AipsError("Error in cloning FTMachine"));
217 0 : return VisModelData::NEW_FT(rec);
218 0 : }
219 :
220 : //=================
221 :
222 : /* template <typename T> void FTMachine::getGrid(Array<T>& thegrid){
223 : thegrid.resize();
224 : if(whatType<Array<T>>()==TpArrayComplex)
225 : thegrid.assign(griddedData);
226 : else if((whatType<Array<T>>()==TpArrayDComplex))
227 : thegrid.assign(griddedData2);
228 : else if(((whatType<Array<T>>()==TpArrayFloat))){
229 : thegrid.resize(griddedData.shape());
230 : thegrid=real(griddedData);
231 : }
232 : else if(((whatType<Array<T>>()==TpArrayDouble))){
233 : thegrid.resize(griddedData2.shape());
234 : thegrid=real(griddedData2);
235 : }
236 :
237 :
238 : }
239 : */
240 :
241 : //----------------------------------------------------------------------
242 303 : Bool FTMachine::changed(const VisBuffer&) {
243 303 : return false;
244 : }
245 :
246 : //----------------------------------------------------------------------
247 0 : FTMachine::FTMachine(const FTMachine& other)
248 : {
249 0 : operator=(other);
250 0 : }
251 :
252 0 : Bool FTMachine::doublePrecGrid(){
253 0 : return useDoubleGrid_p;
254 : }
255 :
256 : //----------------------------------------------------------------------
257 63 : void FTMachine::initPolInfo(const VisBuffer& vb)
258 : {
259 : //
260 : // Need to figure out where to compute the following arrays/ints
261 : // in the re-factored code.
262 : // ----------------------------------------------------------------
263 : {
264 63 : polInUse_p = 0;
265 63 : uInt N=0;
266 191 : for(uInt i=0;i<polMap.nelements();i++) if (polMap(i) > -1) polInUse_p++;
267 63 : cfStokes_p.resize(polInUse_p);
268 191 : for(uInt i=0;i<polMap.nelements();i++)
269 128 : if (polMap(i) > -1) {cfStokes_p(N) = vb.corrType()(i);N++;}
270 : }
271 63 : }
272 : //----------------------------------------------------------------------
273 63 : void FTMachine::initMaps(const VisBuffer& vb) {
274 :
275 63 : logIO() << LogOrigin("FTMachine", "initMaps") << LogIO::NORMAL;
276 :
277 63 : AlwaysAssert(image, AipsError);
278 :
279 : // Set the frame for the UVWMachine
280 63 : mFrame_p=MeasFrame(MEpoch(Quantity(vb.time()(0), "s"), vb.msColumns().timeMeas()(0).getRef()), mLocation_p);
281 :
282 : // First get the CoordinateSystem for the image and then find
283 : // the DirectionCoordinate
284 63 : CoordinateSystem coords=image->coordinates();
285 63 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
286 63 : AlwaysAssert(directionIndex>=0, AipsError);
287 : DirectionCoordinate
288 63 : directionCoord=coords.directionCoordinate(directionIndex);
289 :
290 : // get the first position of moving source
291 63 : if(fixMovingSource_p){
292 :
293 : //First convert to HA-DEC or AZEL for parallax correction
294 0 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
295 0 : MDirection tmphadec=MDirection::Convert(movingDir_p, outref1)();
296 0 : MDirection::Ref outref(directionCoord.directionType(), mFrame_p);
297 0 : firstMovingDir_p=MDirection::Convert(tmphadec, outref)();
298 :
299 0 : }
300 :
301 :
302 : // Now we need MDirection of the image phase center. This is
303 : // what we define it to be. So we define it to be the
304 : // center pixel. So we have to do the conversion here.
305 : // This is independent of padding since we just want to know
306 : // what the world coordinates are for the phase center
307 : // pixel
308 : {
309 63 : Vector<Double> pixelPhaseCenter(2);
310 63 : pixelPhaseCenter(0) = Double( image->shape()(0) / 2 );
311 63 : pixelPhaseCenter(1) = Double( image->shape()(1) / 2 );
312 63 : directionCoord.toWorld(mImage_p, pixelPhaseCenter);
313 :
314 63 : }
315 :
316 :
317 : // Decide if uvwrotation is not necessary, if phasecenter and
318 : // image center are with in one pixel distance; Save some
319 : // computation time especially for spectral cubes.
320 : {
321 126 : Vector<Double> equal= (mImage_p.getAngle()-
322 189 : vb.phaseCenter().getAngle()).getValue();
323 189 : if((abs(equal(0)) < abs(directionCoord.increment()(0)))
324 189 : && (abs(equal(1)) < abs(directionCoord.increment()(1)))){
325 31 : doUVWRotation_p=false;
326 : }
327 : else{
328 32 : doUVWRotation_p=true;
329 : }
330 63 : }
331 : // Get the object distance in meters
332 63 : Record info(image->miscInfo());
333 63 : if(info.isDefined("distance")) {
334 0 : info.get("distance", distance_p);
335 0 : if(abs(distance_p)>0.0) {
336 0 : logIO() << "Distance to object is set to " << distance_p/1000.0
337 0 : << "km: applying focus correction" << LogIO::POST;
338 : }
339 : }
340 :
341 : // Set up the UVWMachine.
342 63 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
343 :
344 :
345 63 : String observatory=vb.msColumns().observation().telescopeName()(0);
346 126 : if(observatory.contains("ATCA") || observatory.contains("DRAO")
347 126 : || observatory.contains("WSRT")){
348 0 : uvwMachine_p=new UVWMachine(mImage_p, vb.phaseCenter(phaseCenterTime_p), mFrame_p,
349 0 : true, false);
350 : }
351 : else{
352 126 : uvwMachine_p=new UVWMachine(mImage_p, vb.phaseCenter(phaseCenterTime_p), mFrame_p,
353 63 : false, tangentSpecified_p);
354 : }
355 63 : AlwaysAssert(uvwMachine_p, AipsError);
356 :
357 63 : lastFieldId_p=-1;
358 63 : lastMSId_p=vb.msId();
359 :
360 : // Set up maps
361 63 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
362 63 : AlwaysAssert(spectralIndex>-1, AipsError);
363 63 : spectralCoord_p=coords.spectralCoordinate(spectralIndex);
364 :
365 : //Store the image/grid channels freq values
366 : {
367 63 : Int chanNumbre=image->shape()(3);
368 63 : Vector<Double> pixindex(chanNumbre);
369 63 : imageFreq_p.resize(chanNumbre);
370 63 : Vector<Double> tempStorFreq(chanNumbre);
371 63 : indgen(pixindex);
372 : // pixindex=pixindex+1.0;
373 126 : for (Int ll=0; ll< chanNumbre; ++ll){
374 63 : if( !spectralCoord_p.toWorld(tempStorFreq(ll), pixindex(ll))){
375 0 : logIO() << "Cannot get imageFreq " << LogIO::EXCEPTION;
376 :
377 : }
378 : }
379 63 : convertArray(imageFreq_p,tempStorFreq);
380 63 : }
381 : //Destroy any conversion layer Freq coord if freqframe is not valid
382 63 : if(!freqFrameValid_p){
383 49 : MFrequency::Types imageFreqType=spectralCoord_p.frequencySystem();
384 49 : spectralCoord_p.setFrequencySystem(imageFreqType);
385 49 : spectralCoord_p.setReferenceConversion(imageFreqType,
386 98 : MEpoch(Quantity(vb.time()(0), "s")),
387 49 : mLocation_p,
388 49 : mImage_p);
389 : }
390 :
391 : // Channel map: do this properly by looking up the frequencies
392 : // If a visibility channel does not map onto an image
393 : // pixel then we set the corresponding chanMap to -1.
394 : // This means that put and get must always check for this
395 : // value (see e.g. GridFT)
396 :
397 63 : nvischan = vb.frequency().nelements();
398 63 : interpVisFreq_p.resize();
399 63 : interpVisFreq_p=vb.frequency();
400 63 : if(selectedSpw_p.nelements() < 1){
401 33 : Vector<Int> myspw(1);
402 33 : myspw[0]=vb.spectralWindow();
403 33 : setSpw(myspw, freqFrameValid_p);
404 33 : }
405 :
406 63 : matchAllSpwChans(vb);
407 :
408 63 : chanMap.resize();
409 :
410 63 : chanMap=multiChanMap_p[vb.spectralWindow()];
411 63 : if(chanMap.nelements() == 0)
412 0 : chanMap=Vector<Int>(vb.frequency().nelements(), -1);
413 :
414 : {
415 : //logIO() << LogIO::DEBUGGING << "Channel Map: " << chanMap << LogIO::POST;
416 : }
417 : // Should never get here
418 63 : if(max(chanMap)>=nchan||min(chanMap)<-3) {
419 0 : logIO() << "Illegal Channel Map: " << chanMap << LogIO::EXCEPTION;
420 : }
421 :
422 : // Polarization map
423 63 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
424 63 : AlwaysAssert(stokesIndex>-1, AipsError);
425 63 : StokesCoordinate stokesCoord=coords.stokesCoordinate(stokesIndex);
426 :
427 63 : Vector<Int> visPolMap(vb.corrType());
428 63 : nvispol=visPolMap.nelements();
429 63 : AlwaysAssert(nvispol>0, AipsError);
430 63 : polMap.resize(nvispol);
431 63 : polMap=-1;
432 63 : isIOnly=false;
433 63 : Int pol=0;
434 63 : Bool found=false;
435 : // First we try matching Stokes in the visibilities to
436 : // Stokes in the image that we are gridding into.
437 191 : for (pol=0;pol<nvispol;pol++) {
438 128 : Int p=0;
439 128 : if(stokesCoord.toPixel(p, Stokes::type(visPolMap(pol)))) {
440 0 : AlwaysAssert(p<npol, AipsError);
441 0 : polMap(pol)=p;
442 0 : found=true;
443 : }
444 : }
445 : // If this fails then perhaps we were looking to grid I
446 : // directly. If so then we need to check that the parallel
447 : // hands are present in the visibilities.
448 63 : if(!found) {
449 63 : Int p=0;
450 63 : if(stokesCoord.toPixel(p, Stokes::I)) {
451 63 : polMap=-1;
452 63 : if(vb.polFrame()==VisibilityIterator::Linear) {
453 62 : p=0;
454 186 : for (pol=0;pol<nvispol;pol++) {
455 124 : if(Stokes::type(visPolMap(pol))==Stokes::XX)
456 62 : {polMap(pol)=0;p++;found=true;};
457 124 : if(Stokes::type(visPolMap(pol))==Stokes::YY)
458 62 : {polMap(pol)=0;p++;found=true;};
459 : }
460 : }
461 : else {
462 1 : p=0;
463 5 : for (pol=0;pol<nvispol;pol++) {
464 4 : if(Stokes::type(visPolMap(pol))==Stokes::LL)
465 1 : {polMap(pol)=0;p++;found=true;};
466 4 : if(Stokes::type(visPolMap(pol))==Stokes::RR)
467 1 : {polMap(pol)=0;p++;found=true;};
468 : }
469 : }
470 63 : if(!found) {
471 : logIO() << "Cannot find polarization map: visibility polarizations = "
472 0 : << visPolMap << LogIO::EXCEPTION;
473 : }
474 : else {
475 63 : isIOnly=true;
476 : //logIO() << LogIO::DEBUGGING << "Transforming I only" << LogIO::POST;
477 : }
478 : };
479 : }
480 : //logIO() << LogIO::DEBUGGING << "Polarization map = "<< polMap
481 : // << LogIO::POST;
482 :
483 63 : initPolInfo(vb);
484 63 : pop_p->initCFMaps(visPolMap, polMap);
485 63 : }
486 :
487 64 : FTMachine::~FTMachine()
488 : {
489 64 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
490 64 : }
491 :
492 463 : Bool FTMachine::interpolateFrequencyTogrid(const VisBuffer& vb,
493 : const Matrix<Float>& wt,
494 : Cube<Complex>& data,
495 : Cube<Int>& flags,
496 : Matrix<Float>& weight,
497 : FTMachine::Type type){
498 463 : Cube<Complex> origdata;
499 463 : Cube<Bool> modflagCube;
500 463 : Vector<Double> visFreq(vb.frequency().nelements());
501 463 : if(doConversion_p[vb.spectralWindow()]){
502 0 : visFreq.resize(lsrFreq_p.shape());
503 0 : convertArray(visFreq, lsrFreq_p);
504 : }
505 : else{
506 463 : convertArray(visFreq, vb.frequency());
507 463 : lsrFreq_p.resize();
508 463 : lsrFreq_p=vb.frequency();
509 : }
510 463 : if(type==FTMachine::MODEL){
511 0 : origdata.reference(vb.modelVisCube());
512 : }
513 463 : else if(type==FTMachine::CORRECTED){
514 0 : origdata.reference(vb.correctedVisCube());
515 : }
516 463 : else if(type==FTMachine::OBSERVED){
517 200 : origdata.reference(vb.visCube());
518 : }
519 263 : else if(type==FTMachine::PSF){
520 : // make sure its a size 0 data ...psf
521 : //so avoid reading any data from disk
522 263 : origdata.resize();
523 :
524 : }
525 : else{
526 0 : throw(AipsError("Don't know which column is being regridded"));
527 : }
528 463 : if((imageFreq_p.nelements()==1) || (freqInterpMethod_p== InterpolateArray1D<Double, Complex>::nearestNeighbour) || (vb.nChannel()==1)){
529 463 : data.reference(origdata);
530 : // do something here for apply flag based on spw chan sels
531 : // e.g.
532 : // setSpecFlag(vb, chansels_p) -> newflag cube
533 463 : setSpectralFlag(vb,modflagCube);
534 : //flags.resize(vb.flagCube().shape());
535 463 : flags.resize(modflagCube.shape());
536 463 : flags=0;
537 : //flags(vb.flagCube())=true;
538 463 : flags(modflagCube)=true;
539 463 : weight.reference(wt);
540 463 : interpVisFreq_p.resize();
541 463 : interpVisFreq_p=lsrFreq_p;
542 :
543 463 : return false;
544 : }
545 :
546 0 : Cube<Bool>flag;
547 :
548 : //okay at this stage we have at least 2 channels
549 0 : Double width=fabs(imageFreq_p[1]-imageFreq_p[0])/fabs(visFreq[1]-visFreq[0]);
550 : //if width is smaller than number of points needed for interpolation ...do it directly
551 : //
552 : // If image chan width is more than twice the data chan width, make a new list of
553 : // data frequencies on which to interpolate. This new list is sync'd with the starting image chan
554 : // and have the same width as the data chans.
555 0 : if(((width >2.0) && (freqInterpMethod_p==InterpolateArray1D<Double, Complex>::linear)) ||
556 0 : ((width >4.0) && (freqInterpMethod_p !=InterpolateArray1D<Double, Complex>::linear))){
557 0 : Double minVF=min(visFreq);
558 0 : Double maxVF=max(visFreq);
559 0 : Double minIF=min(imageFreq_p);
560 0 : Double maxIF=max(imageFreq_p);
561 0 : if( ((minIF-fabs(imageFreq_p[1]-imageFreq_p[0])/2.0) > maxVF) ||
562 0 : ((maxIF+fabs(imageFreq_p[1]-imageFreq_p[0])/2.0) < minVF)){
563 : //This function should not have been called with image
564 : //being out of bound of data...but still
565 0 : interpVisFreq_p.resize(imageFreq_p.nelements());
566 0 : interpVisFreq_p=imageFreq_p;
567 0 : chanMap.resize(interpVisFreq_p.nelements());
568 0 : chanMap.set(-1);
569 : }
570 : else{ // Make a new list of frequencies.
571 : Bool found;
572 0 : uInt where=0;
573 0 : Double interpwidth=visFreq[1]-visFreq[0];
574 0 : if(minIF < minVF){ // Need to find the first image-channel with data in it
575 0 : where=binarySearchBrackets(found, imageFreq_p, minVF, imageFreq_p.nelements());
576 0 : if(where != imageFreq_p.nelements()){
577 0 : minIF=imageFreq_p[where];
578 : }
579 : }
580 :
581 0 : if(maxIF > maxVF){
582 0 : where=binarySearchBrackets(found, imageFreq_p, maxVF, imageFreq_p.nelements());
583 0 : if(where!= imageFreq_p.nelements()){
584 0 : maxIF=imageFreq_p[where];
585 : }
586 :
587 : }
588 :
589 : // This new list of frequencies starts at the first image channel minus half image channel.
590 : // It ends at the last image channel plus half image channel.
591 0 : Int ninterpchan=(Int)ceil((maxIF-minIF+fabs(imageFreq_p[1]-imageFreq_p[0]))/fabs(interpwidth));
592 0 : chanMap.resize(ninterpchan);
593 0 : chanMap.set(-1);
594 0 : interpVisFreq_p.resize(ninterpchan);
595 0 : interpVisFreq_p[0]=(interpwidth > 0) ? minIF : maxIF;
596 0 : interpVisFreq_p[0] =(interpwidth >0) ? (interpVisFreq_p[0]-fabs(imageFreq_p[1]-imageFreq_p[0])/2.0):
597 0 : (interpVisFreq_p[0]+fabs(imageFreq_p[1]-imageFreq_p[0])/2.0);
598 0 : for (Int k=1; k < ninterpchan; ++k){
599 0 : interpVisFreq_p[k] = interpVisFreq_p[k-1]+ interpwidth;
600 : }
601 :
602 0 : for (Int k=0; k < ninterpchan; ++k){
603 : ///chanmap with width
604 0 : Double nearestchanval = interpVisFreq_p[k]- (imageFreq_p[1]-imageFreq_p[0])/2.0;
605 0 : where=binarySearchBrackets(found, imageFreq_p, nearestchanval, imageFreq_p.nelements());
606 0 : if(where != imageFreq_p.nelements())
607 0 : chanMap[k]=where;
608 : }
609 :
610 : }// By now, we have a new list of frequencies, synchronized with image channels, but with data chan widths.
611 0 : }// end of ' if (we have to make new frequencies) '
612 : else{
613 : // Interpolate directly onto output image frequencies.
614 0 : interpVisFreq_p.resize(imageFreq_p.nelements());
615 0 : convertArray(interpVisFreq_p, imageFreq_p);
616 0 : chanMap.resize(interpVisFreq_p.nelements());
617 0 : indgen(chanMap);
618 : }
619 :
620 : // Read flags from the vb.
621 0 : setSpectralFlag(vb,modflagCube);
622 :
623 0 : if(type != FTMachine::PSF){ // Interpolating the data
624 : //Need to get new interpolate functions that interpolate explicitly on the 2nd axis
625 : //2 swap of axes needed
626 0 : Cube<Complex> flipdata;
627 0 : Cube<Bool> flipflag;
628 :
629 : // Interpolate the data.
630 : // Input flags are from the previous step ( setSpectralFlag ).
631 : // Output flags contain info about channels that could not be interpolated
632 : // (for example, linear interp with only one data point)
633 0 : swapyz(flipflag,modflagCube);
634 0 : swapyz(flipdata,origdata);
635 : InterpolateArray1D<Double,Complex>::
636 0 : interpolate(data,flag,interpVisFreq_p,visFreq,flipdata,flipflag,freqInterpMethod_p);
637 0 : flipdata.resize();
638 0 : swapyz(flipdata,data);
639 0 : data.resize();
640 0 : data.reference(flipdata);
641 0 : flipflag.resize();
642 0 : swapyz(flipflag,flag);
643 0 : flag.resize();
644 0 : flag.reference(flipflag);
645 : // Note : 'flag' will get augmented with the flags coming out of weight interpolation
646 0 : }
647 : else
648 : { // get the flag array to the correct shape.
649 : // This will get filled at the end of weight-interpolation.
650 0 : flag.resize(vb.nCorr(), interpVisFreq_p.nelements(), vb.nRow());
651 0 : flag.set(false);
652 : }
653 : // Now, interpolate the weights also.
654 : // (1) Read in the flags from the vb ( setSpectralFlags -> modflagCube )
655 : // (2) Collapse the flags along the polarization dimension to match shape of weight.
656 0 : Matrix<Bool> chanflag(wt.shape());
657 0 : AlwaysAssert( chanflag.shape()[0]==modflagCube.shape()[1], AipsError);
658 0 : AlwaysAssert( chanflag.shape()[1]==modflagCube.shape()[2], AipsError);
659 0 : chanflag=false;
660 0 : for(uInt pol=0;pol<modflagCube.shape()[0];pol++)
661 0 : chanflag = chanflag | modflagCube.yzPlane(pol);
662 :
663 : // (3) Interpolate the weights.
664 : // Input flags are the collapsed vb flags : 'chanflag'
665 : // Output flags are in tempoutputflag
666 : // - contains info about channels that couldn't be interpolated.
667 0 : Matrix<Float> flipweight;
668 0 : flipweight=transpose(wt);
669 0 : Matrix<Bool> flipchanflag;
670 0 : flipchanflag=transpose(chanflag);
671 0 : Matrix<Bool> tempoutputflag;
672 : InterpolateArray1D<Double,Float>::
673 0 : interpolate(weight,tempoutputflag, interpVisFreq_p, visFreq,flipweight,flipchanflag,freqInterpMethod_p);
674 0 : flipweight.resize();
675 0 : flipweight=transpose(weight);
676 0 : weight.resize();
677 0 : weight.reference(flipweight);
678 0 : flipchanflag.resize();
679 0 : flipchanflag=transpose(tempoutputflag);
680 0 : tempoutputflag.resize();
681 0 : tempoutputflag.reference(flipchanflag);
682 :
683 : // (4) Now, fill these flags back into the flag cube
684 : // so that they get USED while gridding the PSF (and data)
685 : // Taking the OR of the flags that came out of data-interpolation
686 : // and weight-interpolation, in case they're different.
687 : // Expanding flags across polarization. This will destroy any
688 : // pol-dependent flags for imaging, but msvis::VisImagingWeight
689 : // uses the OR of flags across polarization anyway
690 : // so we don't lose anything.
691 :
692 0 : AlwaysAssert( tempoutputflag.shape()[0]==flag.shape()[1], AipsError);
693 0 : AlwaysAssert( tempoutputflag.shape()[1]==flag.shape()[2], AipsError);
694 0 : for(uInt pol=0;pol<flag.shape()[0];pol++)
695 0 : flag.yzPlane(pol) = tempoutputflag | flag.yzPlane(pol);
696 :
697 : // Fill the output array of image-channel flags.
698 0 : flags.resize(flag.shape());
699 0 : flags=0;
700 0 : flags(flag)=true;
701 :
702 0 : return true;
703 463 : }
704 :
705 358 : void FTMachine::getInterpolateArrays(const VisBuffer& vb,
706 : Cube<Complex>& data, Cube<Int>& flags){
707 :
708 :
709 358 : if((imageFreq_p.nelements()==1) || (freqInterpMethod_p== InterpolateArray1D<Double, Complex>::nearestNeighbour)|| (vb.nChannel()==1)){
710 358 : Cube<Bool> modflagCube;
711 358 : setSpectralFlag(vb,modflagCube);
712 358 : data.reference(vb.modelVisCube());
713 : //flags.resize(vb.flagCube().shape());
714 358 : flags.resize(modflagCube.shape());
715 358 : flags=0;
716 : //flags(vb.flagCube())=true;
717 358 : flags(modflagCube)=true;
718 358 : interpVisFreq_p.resize();
719 358 : interpVisFreq_p=vb.frequency();
720 358 : return;
721 358 : }
722 :
723 0 : data.resize(vb.nCorr(), imageFreq_p.nelements(), vb.nRow());
724 0 : flags.resize(vb.nCorr(), imageFreq_p.nelements(), vb.nRow());
725 0 : data.set(Complex(0.0,0.0));
726 0 : flags.set(0);
727 : //no need to degrid channels that does map over this vb
728 0 : Int maxchan=max(chanMap);
729 0 : for (uInt k =0 ; k < chanMap.nelements() ; ++k){
730 0 : if(chanMap(k)==-1)
731 0 : chanMap(k)=maxchan;
732 : }
733 0 : Int minchan=min(chanMap);
734 0 : if(minchan==maxchan)
735 0 : minchan=-1;
736 :
737 :
738 0 : for(Int k = 0; k < minchan; ++k)
739 0 : flags.xzPlane(k).set(1);
740 :
741 0 : for(uInt k = maxchan + 1; k < imageFreq_p.nelements(); ++k)
742 0 : flags.xzPlane(k).set(1);
743 0 : interpVisFreq_p.resize(imageFreq_p.nelements());
744 0 : convertArray(interpVisFreq_p, imageFreq_p);
745 0 : chanMap.resize(imageFreq_p.nelements());
746 0 : indgen(chanMap);
747 : }
748 :
749 358 : Bool FTMachine::interpolateFrequencyFromgrid(VisBuffer& vb,
750 : Cube<Complex>& data,
751 : FTMachine::Type type){
752 :
753 : Cube<Complex> *origdata;
754 358 : Vector<Double> visFreq(vb.frequency().nelements());
755 :
756 358 : if(doConversion_p[vb.spectralWindow()]){
757 0 : convertArray(visFreq, lsrFreq_p);
758 : }
759 : else{
760 358 : convertArray(visFreq, vb.frequency());
761 : }
762 :
763 358 : if(type==FTMachine::MODEL){
764 358 : origdata=&(vb.modelVisCube());
765 : }
766 0 : else if(type==FTMachine::CORRECTED){
767 0 : origdata=&(vb.correctedVisCube());
768 : }
769 : else{
770 0 : origdata=&(vb.visCube());
771 : }
772 :
773 : //
774 : // If visibility data (vb) has only one channel, or the image cube
775 : // has only one channel, resort to nearestNeighbour interpolation.
776 : // Honour user selection of nearestNeighbour.
777 : //
778 358 : if((imageFreq_p.nelements()==1) ||
779 358 : (vb.nChannel()==1) ||
780 0 : (freqInterpMethod_p== InterpolateArray1D<Double, Complex>::nearestNeighbour)){
781 358 : origdata->reference(data);
782 358 : return false;
783 : }
784 :
785 : //Need to get new interpolate functions that interpolate explicitly on the 2nd axis
786 : //2 swap of axes needed
787 0 : Cube<Complex> flipgrid;
788 0 : flipgrid.resize();
789 0 : swapyz(flipgrid,data);
790 :
791 0 : Cube<Complex> flipdata((origdata->shape())(0),(origdata->shape())(2),
792 0 : (origdata->shape())(1)) ;
793 0 : flipdata.set(Complex(0.0));
794 : InterpolateArray1D<Double,Complex>::
795 0 : interpolate(flipdata,visFreq, imageFreq_p, flipgrid,freqInterpMethod_p);
796 :
797 :
798 :
799 0 : Cube<Bool> copyOfFlag;
800 : //cerr << "spw " << vb.spectralWindow() << " chanMap " << multiChanMap_p[vb.spectralWindow()] << endl;
801 0 : Vector<Int> mychanmap;
802 0 : mychanmap=multiChanMap_p[vb.spectralWindow()];
803 0 : copyOfFlag.assign(vb.flagCube());
804 0 : for (uInt k=0; k< mychanmap.nelements(); ++ k)
805 0 : if(mychanmap(k) ==-1)
806 0 : copyOfFlag.xzPlane(k).set(true);
807 :
808 0 : swapyz(vb.modelVisCube(), copyOfFlag, flipdata);
809 : //swapyz(vb.modelVisCube(), flipdata);
810 :
811 0 : return true;
812 358 : }
813 0 : void FTMachine::girarUVW(Matrix<Double>& uvw, Vector<Double>& dphase,
814 : const VisBuffer& vb)
815 : {
816 :
817 :
818 :
819 : //the uvw rotation is done for common tangent reprojection or if the
820 : //image center is different from the phasecenter
821 : // UVrotation is false only if field never changes
822 :
823 :
824 0 : if((vb.fieldId()!=lastFieldId_p) || (vb.msId()!=lastMSId_p))
825 0 : doUVWRotation_p=true;
826 0 : if(doUVWRotation_p || fixMovingSource_p){
827 :
828 0 : mFrame_p.epoch() != 0 ?
829 0 : mFrame_p.resetEpoch(MEpoch(Quantity(vb.time()(0), "s"))):
830 0 : mFrame_p.set(mLocation_p, MEpoch(Quantity(vb.time()(0), "s"), vb.msColumns().timeMeas()(0).getRef()));
831 : MDirection::Types outType;
832 0 : MDirection::getType(outType, mImage_p.getRefString());
833 0 : MDirection phasecenter=MDirection::Convert(vb.phaseCenter(phaseCenterTime_p), MDirection::Ref(outType, mFrame_p))();
834 :
835 :
836 0 : if(fixMovingSource_p){
837 :
838 : //First convert to HA-DEC or AZEL for parallax correction
839 0 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
840 0 : MDirection tmphadec=MDirection::Convert(movingDir_p, outref1)();
841 0 : MDirection::Ref outref(mImage_p.getRef().getType(), mFrame_p);
842 0 : MDirection sourcenow=MDirection::Convert(tmphadec, outref)();
843 : //cerr << "Rotating to fixed moving source " << MVDirection(phasecenter.getAngle()-firstMovingDir_p.getAngle()+sourcenow.getAngle()) << endl;
844 0 : phasecenter.set(MVDirection(phasecenter.getAngle()+firstMovingDir_p.getAngle()-sourcenow.getAngle()));
845 :
846 0 : }
847 :
848 :
849 : // Set up the UVWMachine only if the field id has changed. If
850 : // the tangent plane is specified then we need a UVWMachine that
851 : // will reproject to that plane iso the image plane
852 0 : if((vb.fieldId()!=lastFieldId_p) || (vb.msId()!=lastMSId_p) || fixMovingSource_p) {
853 :
854 0 : String observatory=vb.msColumns().observation().telescopeName()(0);
855 0 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
856 0 : if(observatory.contains("ATCA") || observatory.contains("WSRT")){
857 : //Tangent specified is being wrongly used...it should be for a
858 : //Use the safest way for now.
859 0 : uvwMachine_p=new UVWMachine(phasecenter, vb.phaseCenter(phaseCenterTime_p), mFrame_p,
860 0 : true, false);
861 0 : phaseShifter_p=new UVWMachine(mImage_p, phasecenter, mFrame_p,
862 0 : true, false);
863 : }
864 : else{
865 0 : uvwMachine_p=new UVWMachine(phasecenter, vb.phaseCenter(phaseCenterTime_p), mFrame_p,
866 0 : false, false);
867 0 : phaseShifter_p=new UVWMachine(mImage_p, phasecenter, mFrame_p,
868 0 : false, false);
869 : }
870 0 : }
871 :
872 0 : lastFieldId_p=vb.fieldId();
873 0 : lastMSId_p=vb.msId();
874 :
875 :
876 0 : AlwaysAssert(uvwMachine_p, AipsError);
877 :
878 : // Always force a recalculation
879 0 : uvwMachine_p->reCalculate();
880 0 : phaseShifter_p->reCalculate();
881 :
882 : // Now do the conversions
883 0 : uInt nrows=dphase.nelements();
884 0 : Vector<Double> thisRow(3);
885 0 : thisRow=0.0;
886 : //CoordinateSystem csys=image->coordinates();
887 : //DirectionCoordinate dc=csys.directionCoordinate(0);
888 : //Vector<Double> thePix(2);
889 : //dc.toPixel(thePix, phasecenter);
890 : //cerr << "field id " << vb.fieldId() << " the Pix " << thePix << endl;
891 : //Vector<Float> scale(2);
892 : //scale(0)=dc.increment()(0);
893 : //scale(1)=dc.increment()(1);
894 0 : for (uInt irow=0; irow<nrows;++irow) {
895 0 : thisRow.assign(uvw.column(irow));
896 : //cerr << " uvw " << thisRow ;
897 : // This is for frame change
898 0 : uvwMachine_p->convertUVW(dphase(irow), thisRow);
899 : // This is for correlator phase center change
900 0 : MVPosition rotphase=phaseShifter_p->rotationPhase() ;
901 : //cerr << " rotPhase " << rotphase << " oldphase "<< rotphase*(uvw.column(irow)) << " newphase " << (rotphase)*thisRow ;
902 : // cerr << " phase " << dphase(irow) << " new uvw " << uvw.column(irow);
903 : //dphase(irow)+= (thePix(0)-nx/2.0)*thisRow(0)*scale(0)+(thePix(1)-ny/2.0)*thisRow(1)*scale(1);
904 : //Double pixphase=(thePix(0)-nx/2.0)*uvw.column(irow)(0)*scale(0)+(thePix(1)-ny/2.0)*uvw.column(irow)(1)*scale(1);
905 : //Double pixphase2=(thePix(0)-nx/2.0)*thisRow(0)*scale(0)+(thePix(1)-ny/2.0)*thisRow(1)*scale(1);
906 : //cerr << " pixphase " << pixphase << " pixphase2 " << pixphase2<< endl;
907 : //dphase(irow)=pixphase;
908 : //dphase(irow)+= rotphase(0)*thisRow(0)+rotphase(1)*thisRow(1);
909 0 : RotMatrix rotMat=phaseShifter_p->rotationUVW();
910 : //cerr << "rot 0 " << rotMat(0,0) << " " << rotMat(1,0) << " " << rotMat(2,0) << endl;
911 : //cerr << "rot 1 " << rotMat(0,1) << " " << rotMat(1,1) << " " << rotMat(2,1) <<
912 0 : uvw.column(irow)(0)=thisRow(0)*rotMat(0,0)+thisRow(1)*rotMat(1,0);
913 0 : uvw.column(irow)(1)=thisRow(0)*rotMat(0,1)+thisRow(1)*rotMat(1,1);
914 :
915 0 : uvw.column(irow)(2)=thisRow(0)*rotMat(0,2)+thisRow(1)*rotMat(1,2)+thisRow(2)*rotMat(2,2);
916 : //cerr << "w term " << thisRow(2) << " aft " << uvw.column(irow)(2) << endl;
917 0 : dphase(irow)+= rotphase(0)*uvw.column(irow)(0)+rotphase(1)*uvw.column(irow)(1);
918 0 : }
919 :
920 :
921 0 : }
922 0 : }
923 :
924 503 : void FTMachine::rotateUVW(Matrix<Double>& uvw, Vector<Double>& dphase,
925 : const VisBuffer& vb)
926 : {
927 :
928 :
929 :
930 : //the uvw rotation is done for common tangent reprojection or if the
931 : //image center is different from the phasecenter
932 : // UVrotation is false only if field never changes
933 503 : if((vb.fieldId()!=lastFieldId_p) || (vb.msId()!=lastMSId_p))
934 104 : doUVWRotation_p=true;
935 503 : if(doUVWRotation_p || tangentSpecified_p || fixMovingSource_p){
936 503 : ok();
937 :
938 503 : mFrame_p.epoch() != 0 ?
939 1006 : mFrame_p.resetEpoch(MEpoch(Quantity(vb.time()(0), "s"))):
940 503 : mFrame_p.set(mLocation_p, MEpoch(Quantity(vb.time()(0), "s"), vb.msColumns().timeMeas()(0).getRef()));
941 :
942 503 : MDirection phasecenter=mImage_p;
943 503 : if(fixMovingSource_p){
944 :
945 : //First convert to HA-DEC or AZEL for parallax correction
946 0 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
947 0 : MDirection tmphadec=MDirection::Convert(movingDir_p, outref1)();
948 0 : MDirection::Ref outref(mImage_p.getRef().getType(), mFrame_p);
949 0 : MDirection sourcenow=MDirection::Convert(tmphadec, outref)();
950 : //cerr << "Rotating to fixed moving source " << MVDirection(phasecenter.getAngle()-firstMovingDir_p.getAngle()+sourcenow.getAngle()) << endl;
951 0 : phasecenter.set(MVDirection(phasecenter.getAngle()+firstMovingDir_p.getAngle()-sourcenow.getAngle()));
952 :
953 0 : }
954 :
955 :
956 : // Set up the UVWMachine only if the field id has changed. If
957 : // the tangent plane is specified then we need a UVWMachine that
958 : // will reproject to that plane iso the image plane
959 503 : if((vb.fieldId()!=lastFieldId_p) || (vb.msId()!=lastMSId_p) || fixMovingSource_p) {
960 :
961 104 : String observatory=vb.msColumns().observation().telescopeName()(0);
962 104 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
963 104 : if(observatory.contains("ATCA") || observatory.contains("WSRT")){
964 : //Tangent specified is being wrongly used...it should be for a
965 : //Use the safest way for now.
966 0 : uvwMachine_p=new UVWMachine(phasecenter, vb.phaseCenter(phaseCenterTime_p), mFrame_p,
967 0 : true, false);
968 : }
969 : else{
970 208 : uvwMachine_p=new UVWMachine(phasecenter, vb.phaseCenter(phaseCenterTime_p), mFrame_p,
971 104 : false,tangentSpecified_p);
972 : }
973 104 : }
974 :
975 503 : lastFieldId_p=vb.fieldId();
976 503 : lastMSId_p=vb.msId();
977 :
978 :
979 503 : AlwaysAssert(uvwMachine_p, AipsError);
980 :
981 : // Always force a recalculation
982 503 : uvwMachine_p->reCalculate();
983 :
984 : // Now do the conversions
985 503 : uInt nrows=dphase.nelements();
986 503 : Vector<Double> thisRow(3);
987 503 : thisRow=0.0;
988 : uInt irow;
989 : //#pragma omp parallel default(shared) private(irow,thisRow)
990 : {
991 : //#pragma omp for
992 360958 : for (irow=0; irow<nrows;++irow) {
993 360455 : thisRow.reference(uvw.column(irow));
994 360455 : convUVW(dphase(irow), thisRow);
995 : }
996 :
997 : }//end pragma
998 503 : }
999 503 : }
1000 :
1001 360455 : void FTMachine::convUVW(Double& dphase, Vector<Double>& thisrow){
1002 : //for (uInt i=0;i<3;i++) thisRow(i)=uvw(i,row);
1003 360455 : uvwMachine_p->convertUVW(dphase, thisrow);
1004 : //for (uint i=0;i<3;i++) uvw(i,row)=thisRow(i);
1005 :
1006 360455 : }
1007 :
1008 :
1009 0 : void FTMachine::locateuvw(const Double*& uvw, const Double*& dphase,
1010 : const Double*& freq, const Int& nvchan,
1011 : const Double*& scale, const Double*& offset, const Int& sampling, Int*& loc, Int*& off, Complex*& phasor, const Int& row, const bool& doW){
1012 :
1013 0 : Int rowoff=row*nvchan;
1014 : Double phase;
1015 : Double pos;
1016 0 : Int nel= doW ? 3 : 2;
1017 0 : for (Int f=0; f<nvchan; ++f){
1018 0 : for (Int k=0; k <2; ++k){
1019 0 : pos=(scale[k])*uvw[3*row+k]*(freq[f])/C::c+((offset[k])+1.0);
1020 0 : loc[(rowoff+f)*nel+k]=std::lround(pos);
1021 0 : off[(rowoff+f)*nel+k]=std::lround((Double(loc[(rowoff+f)*nel+k])-pos)*Double(sampling));
1022 : //off[(rowoff+f)*2+k]=(loc[(rowoff+f)*2+k]-pos(k))*sampling;
1023 : }
1024 0 : phase=-Double(2.0)*C::pi*dphase[row]*(freq[f])/C::c;
1025 0 : phasor[rowoff+f]=Complex(cos(phase), sin(phase));
1026 :
1027 : ///This is for W-Projection
1028 0 : if(doW){
1029 0 : pos=sqrt(fabs(scale[2]*uvw[3*row+2]*(freq[f])/C::c))+offset[2]+1.0;
1030 0 : loc[(rowoff+f)*nel+2]=std::lround(pos);
1031 0 : off[(rowoff+f)*nel+2]=0;
1032 : }
1033 : }
1034 :
1035 :
1036 :
1037 :
1038 0 : }
1039 :
1040 14 : void FTMachine::setnumthreads(Int num){
1041 14 : numthreads_p=num;
1042 14 : }
1043 0 : Int FTMachine::getnumthreads(){
1044 0 : return numthreads_p;
1045 : }
1046 :
1047 : //
1048 : // Refocus the array on a point at finite distance
1049 : //
1050 503 : void FTMachine::refocus(Matrix<Double>& uvw, const Vector<Int>& ant1,
1051 : const Vector<Int>& ant2,
1052 : Vector<Double>& dphase, const VisBuffer& vb)
1053 : {
1054 :
1055 503 : ok();
1056 :
1057 503 : if(abs(distance_p)>0.0) {
1058 :
1059 0 : nAntenna_p=max(vb.antenna2())+1;
1060 :
1061 : // Positions of antennas
1062 0 : Matrix<Double> antPos(3,nAntenna_p);
1063 0 : antPos=0.0;
1064 0 : Vector<Int> nAntPos(nAntenna_p);
1065 0 : nAntPos=0;
1066 :
1067 0 : uInt aref = min(ant1);
1068 :
1069 : // Now find the antenna locations: for this we just reference to a common
1070 : // point. We ignore the time variation within this buffer.
1071 0 : uInt nrows=dphase.nelements();
1072 0 : for (uInt row=0;row<nrows;row++) {
1073 0 : uInt a1=ant1(row);
1074 0 : uInt a2=ant2(row);
1075 0 : for(uInt dim=0;dim<3;dim++) {
1076 0 : antPos(dim, a1)+=uvw(dim, row);
1077 0 : antPos(dim, a2)-=uvw(dim, row);
1078 : }
1079 0 : nAntPos(a1)+=1;
1080 0 : nAntPos(a2)+=1;
1081 : }
1082 :
1083 : // Now remove the reference location
1084 0 : Vector<Double> center(3);
1085 0 : for(uInt dim=0;dim<3;dim++) {
1086 0 : center(dim) = antPos(dim,aref)/nAntPos(aref);
1087 : }
1088 :
1089 : // Now normalize
1090 0 : for (uInt ant=0; ant<nAntenna_p; ant++) {
1091 0 : if(nAntPos(ant)>0) {
1092 0 : for(uInt dim=0;dim<3;dim++) {
1093 0 : antPos(dim,ant)/=nAntPos(ant);
1094 0 : antPos(dim,ant)-=center(dim);
1095 : }
1096 : }
1097 : }
1098 :
1099 : // Now calculate the offset needed to focus the array,
1100 : // including the w term. This must have the correct asymptotic
1101 : // form so that at infinity no net change occurs
1102 0 : for (uInt row=0;row<nrows;row++) {
1103 0 : uInt a1=ant1(row);
1104 0 : uInt a2=ant2(row);
1105 :
1106 0 : Double d1=distance_p*distance_p-2*distance_p*antPos(2,a1);
1107 0 : Double d2=distance_p*distance_p-2*distance_p*antPos(2,a2);
1108 0 : for(uInt dim=0;dim<3;dim++) {
1109 0 : d1+=antPos(dim,a1)*antPos(dim,a1);
1110 0 : d2+=antPos(dim,a2)*antPos(dim,a2);
1111 : }
1112 0 : d1=sqrt(d1);
1113 0 : d2=sqrt(d2);
1114 0 : for(uInt dim=0;dim<2;dim++) {
1115 0 : dphase(row)-=(antPos(dim,a1)*antPos(dim,a1)-antPos(dim,a2)*antPos(dim,a2))/(2*distance_p);
1116 : }
1117 0 : uvw(0,row)=distance_p*(antPos(0,a1)/d1-antPos(0,a2)/d2);
1118 0 : uvw(1,row)=distance_p*(antPos(1,a1)/d1-antPos(1,a2)/d2);
1119 0 : uvw(2,row)=distance_p*(antPos(2,a1)/d1-antPos(2,a2)/d2)+dphase(row);
1120 : }
1121 0 : }
1122 503 : }
1123 :
1124 0 : void FTMachine::ok() {
1125 0 : AlwaysAssert(image, AipsError);
1126 0 : AlwaysAssert(uvwMachine_p, AipsError);
1127 0 : }
1128 :
1129 0 : Bool FTMachine::toRecord(String& error, RecordInterface& outRecord,
1130 : Bool withImage, const String diskimage) {
1131 : // Save the FTMachine to a Record
1132 : //
1133 0 : outRecord.define("name", this->name());
1134 0 : if(withImage){
1135 0 : CoordinateSystem cs=image->coordinates();
1136 0 : DirectionCoordinate dircoord=cs.directionCoordinate(cs.findCoordinate(Coordinate::DIRECTION));
1137 0 : dircoord.setReferenceValue(mImage_p.getAngle().getValue());
1138 0 : if(diskimage != ""){
1139 : try{
1140 0 : PagedImage<Complex> imCopy(TiledShape(toVis_p ? griddedData.shape(): image->shape()), image->coordinates(), diskimage);
1141 0 : toVis_p ? imCopy.put(griddedData) : imCopy.copyData(*image);
1142 0 : ImageUtilities::copyMiscellaneous(imCopy, *image);
1143 0 : Vector<Double> pixcen(2);
1144 0 : pixcen(0)=Double(imCopy.shape()(0)/2); pixcen(1)=Double(imCopy.shape()(1)/2);
1145 0 : dircoord.setReferencePixel(pixcen);
1146 0 : cs.replaceCoordinate(dircoord, cs.findCoordinate(Coordinate::DIRECTION));
1147 0 : imCopy.setCoordinateInfo(cs);
1148 :
1149 0 : }
1150 0 : catch(...){
1151 0 : throw(AipsError(String("Failed to save model image "+diskimage+String(" to disk"))));
1152 0 : }
1153 0 : outRecord.define("diskimage", diskimage);
1154 :
1155 : }
1156 : else{
1157 0 : Record imrec;
1158 0 : Vector<Double> pixcen(2);
1159 0 : pixcen(0)=Double(griddedData.shape()(0)/2); pixcen(1)=Double(griddedData.shape()(1)/2);
1160 0 : dircoord.setReferencePixel(pixcen);
1161 0 : cs.replaceCoordinate(dircoord, cs.findCoordinate(Coordinate::DIRECTION));
1162 0 : TempImage<Complex> imCopy(griddedData.shape(), cs);
1163 0 : imCopy.put(griddedData) ;
1164 0 : ImageUtilities::copyMiscellaneous(imCopy, *image);
1165 0 : if(imCopy.toRecord(error, imrec))
1166 0 : outRecord.defineRecord("image", imrec);
1167 0 : }
1168 0 : }
1169 0 : outRecord.define("nantenna", nAntenna_p);
1170 0 : outRecord.define("distance", distance_p);
1171 0 : outRecord.define("lastfieldid", lastFieldId_p);
1172 0 : outRecord.define("lastmsid", lastMSId_p);
1173 0 : outRecord.define("tangentspecified", tangentSpecified_p);
1174 0 : saveMeasure(outRecord, String("mtangent_rec"), error, mTangent_p);
1175 0 : saveMeasure(outRecord, "mimage_rec", error, mImage_p);
1176 : //mFrame_p not necessary to save as it is generated from mLocation_p
1177 0 : outRecord.define("nx", nx);
1178 0 : outRecord.define("ny", ny);
1179 0 : outRecord.define("npol", npol);
1180 0 : outRecord.define("nchan", nchan);
1181 0 : outRecord.define("nvischan", nvischan);
1182 0 : outRecord.define("nvispol", nvispol);
1183 0 : saveMeasure(outRecord, "mlocation_rec", error, mLocation_p);
1184 : //no need to save uvwMachine_p
1185 0 : outRecord.define("douvwrotation", doUVWRotation_p);
1186 0 : outRecord.define("freqinterpmethod", static_cast<Int>(freqInterpMethod_p));
1187 0 : outRecord.define("spwchanselflag", spwChanSelFlag_p);
1188 0 : outRecord.define("freqframevalid", freqFrameValid_p);
1189 0 : outRecord.define("selectedspw", selectedSpw_p);
1190 0 : outRecord.define("imagefreq", imageFreq_p);
1191 0 : outRecord.define("lsrfreq", lsrFreq_p);
1192 0 : outRecord.define("interpvisfreq", interpVisFreq_p);
1193 0 : Record multichmaprec;
1194 0 : for (uInt k=0; k < multiChanMap_p.nelements(); ++k)
1195 0 : multichmaprec.define(k, multiChanMap_p[k]);
1196 0 : outRecord.defineRecord("multichanmaprec", multichmaprec);
1197 0 : outRecord.define("chanmap", chanMap);
1198 0 : outRecord.define("polmap", polMap);
1199 0 : outRecord.define("nvischanmulti", nVisChan_p);
1200 0 : spectralCoord_p.save(outRecord, "spectralcoord");
1201 0 : outRecord.define("doconversion", doConversion_p);
1202 0 : outRecord.define("pointingdircol", pointingDirCol_p);
1203 0 : saveMeasure(outRecord, "movingdir_rec", error, movingDir_p);
1204 0 : outRecord.define("fixmovingsource", fixMovingSource_p);
1205 0 : saveMeasure(outRecord, "firstmovingdir_rec", error, firstMovingDir_p);
1206 0 : outRecord.define("usedoublegrid", useDoubleGrid_p);
1207 0 : outRecord.define("cfstokes", cfStokes_p);
1208 0 : outRecord.define("polinuse", polInUse_p);
1209 0 : outRecord.define("tovis", toVis_p);
1210 0 : outRecord.define("sumweight", sumWeight);
1211 0 : outRecord.define("numthreads", numthreads_p);
1212 0 : outRecord.define("phasecentertime", phaseCenterTime_p);
1213 : //Need to serialized sj_p...the user has to set the sj_p after recovering from record
1214 0 : return true;
1215 0 : };
1216 :
1217 0 : Bool FTMachine::saveMeasure(RecordInterface& rec, const String& name, String& err, const Measure& meas){
1218 0 : Record tmprec;
1219 0 : MeasureHolder mh(meas);
1220 0 : if(mh.toRecord(err, tmprec)){
1221 0 : rec.defineRecord(name, tmprec);
1222 0 : return true;
1223 : }
1224 0 : return false;
1225 0 : }
1226 :
1227 0 : Bool FTMachine::fromRecord(String& error,
1228 : const RecordInterface& inRecord
1229 : ) {
1230 : // Restore an FTMachine from a Record
1231 : //
1232 0 : uvwMachine_p=0; //when null it is reconstructed from mImage_p and mFrame_p
1233 : //mFrame_p is not necessary as it is generated in initMaps from mLocation_p
1234 0 : inRecord.get("nx", nx);
1235 0 : inRecord.get("ny", ny);
1236 0 : inRecord.get("npol", npol);
1237 0 : inRecord.get("nchan", nchan);
1238 0 : inRecord.get("nvischan", nvischan);
1239 0 : inRecord.get("nvispol", nvispol);
1240 0 : cmplxImage_p=NULL;
1241 0 : inRecord.get("tovis", toVis_p);
1242 0 : if(inRecord.isDefined("image")){
1243 0 : cmplxImage_p=new TempImage<Complex>();
1244 0 : image=&(*cmplxImage_p);
1245 :
1246 0 : const Record rec=inRecord.asRecord("image");
1247 0 : if(!cmplxImage_p->fromRecord(error, rec))
1248 0 : return false;
1249 :
1250 0 : }
1251 0 : else if(inRecord.isDefined("diskimage")){
1252 0 : String theDiskImage;
1253 0 : inRecord.get("diskimage", theDiskImage);
1254 : try{
1255 0 : File eldir(theDiskImage);
1256 0 : if(! eldir.exists()){
1257 0 : String subPathname[30];
1258 0 : String sep = "/";
1259 0 : uInt nw = split(theDiskImage, subPathname, 20, sep);
1260 :
1261 0 : String theposs=(subPathname[nw-1]);
1262 0 : Bool isExistant=File(theposs).exists();
1263 0 : if(isExistant)
1264 0 : theDiskImage=theposs;
1265 0 : for (uInt i=nw-2 ; i>0; --i){
1266 0 : theposs=subPathname[i]+"/"+theposs;
1267 0 : File newEldir(theposs);
1268 0 : if(newEldir.exists()){
1269 0 : isExistant=true;
1270 0 : theDiskImage=theposs;
1271 : }
1272 0 : }
1273 0 : if(!isExistant)
1274 0 : throw(AipsError("Could not locate mage"));
1275 0 : }
1276 0 : cmplxImage_p=new PagedImage<Complex> (theDiskImage);
1277 0 : image=&(*cmplxImage_p);
1278 0 : }
1279 0 : catch(...){
1280 0 : throw(AipsError(String("Failure to load ")+theDiskImage+String(" image from disk")));
1281 0 : }
1282 0 : }
1283 0 : if(toVis_p && !cmplxImage_p.null()) {
1284 0 : griddedData.resize(image->shape());
1285 0 : griddedData=image->get();
1286 : }
1287 0 : else if(!toVis_p){
1288 0 : IPosition gridShape(4, nx, ny, npol, nchan);
1289 0 : griddedData.resize(gridShape);
1290 0 : griddedData=Complex(0.0);
1291 0 : }
1292 :
1293 0 : nAntenna_p=inRecord.asuInt("nantenna");
1294 0 : distance_p=inRecord.asDouble("distance");
1295 0 : lastFieldId_p=inRecord.asInt("lastfieldid");
1296 0 : lastMSId_p=inRecord.asInt("lastmsid");
1297 0 : inRecord.get("tangentspecified", tangentSpecified_p);
1298 0 : { const Record rec=inRecord.asRecord("mtangent_rec");
1299 0 : MeasureHolder mh;
1300 0 : if(!mh.fromRecord(error, rec))
1301 0 : return false;
1302 0 : mTangent_p=mh.asMDirection();
1303 0 : }
1304 0 : { const Record rec=inRecord.asRecord("mimage_rec");
1305 0 : MeasureHolder mh;
1306 0 : if(!mh.fromRecord(error, rec))
1307 0 : return false;
1308 0 : mImage_p=mh.asMDirection();
1309 0 : }
1310 :
1311 :
1312 0 : { const Record rec=inRecord.asRecord("mlocation_rec");
1313 0 : MeasureHolder mh;
1314 0 : if(!mh.fromRecord(error, rec))
1315 0 : return false;
1316 0 : mLocation_p=mh.asMPosition();
1317 0 : }
1318 0 : inRecord.get("douvwrotation", doUVWRotation_p);
1319 :
1320 : //inRecord.get("spwchanselflag", spwChanSelFlag_p);
1321 : //We won't respect the chanselflag as the vister may have different selections
1322 0 : spwChanSelFlag_p.resize();
1323 0 : inRecord.get("freqframevalid", freqFrameValid_p);
1324 0 : inRecord.get("selectedspw", selectedSpw_p);
1325 0 : inRecord.get("imagefreq", imageFreq_p);
1326 0 : inRecord.get("lsrfreq", lsrFreq_p);
1327 0 : inRecord.get("interpvisfreq", interpVisFreq_p);
1328 0 : const Record multichmaprec=inRecord.asRecord("multichanmaprec");
1329 0 : multiChanMap_p.resize(multichmaprec.nfields(), true, false);
1330 0 : for (uInt k=0; k < multichmaprec.nfields(); ++k)
1331 0 : multichmaprec.get(k, multiChanMap_p[k]);
1332 0 : inRecord.get("chanmap", chanMap);
1333 0 : inRecord.get("polmap", polMap);
1334 0 : inRecord.get("nvischanmulti", nVisChan_p);
1335 0 : SpectralCoordinate *tmpSpec=SpectralCoordinate::restore(inRecord, "spectralcoord");
1336 0 : if(tmpSpec){
1337 0 : spectralCoord_p=*tmpSpec;
1338 0 : delete tmpSpec;
1339 : }
1340 0 : inRecord.get("doconversion", doConversion_p);
1341 0 : inRecord.get("pointingdircol", pointingDirCol_p);
1342 0 : { const Record rec=inRecord.asRecord("movingdir_rec");
1343 0 : MeasureHolder mh;
1344 0 : if(!mh.fromRecord(error, rec))
1345 0 : return false;
1346 0 : movingDir_p=mh.asMDirection();
1347 0 : }
1348 0 : inRecord.get("fixmovingsource", fixMovingSource_p);
1349 0 : { const Record rec=inRecord.asRecord("firstmovingdir_rec");
1350 0 : MeasureHolder mh;
1351 0 : if(!mh.fromRecord(error, rec))
1352 0 : return false;
1353 0 : firstMovingDir_p=mh.asMDirection();
1354 0 : }
1355 0 : inRecord.get("usedoublegrid", useDoubleGrid_p);
1356 0 : inRecord.get("cfstokes", cfStokes_p);
1357 0 : inRecord.get("polinuse", polInUse_p);
1358 :
1359 :
1360 0 : inRecord.get("sumweight", sumWeight);
1361 0 : if(toVis_p){
1362 0 : freqInterpMethod_p=InterpolateArray1D<Double, Complex>::nearestNeighbour;
1363 : }
1364 : else{
1365 : Int tmpInt;
1366 0 : inRecord.get("freqinterpmethod", tmpInt);
1367 0 : freqInterpMethod_p=static_cast<InterpolateArray1D<Double, Complex >::InterpolationMethod>(tmpInt);
1368 : }
1369 0 : inRecord.get("numthreads", numthreads_p);
1370 0 : inRecord.get("phasecentertime", phaseCenterTime_p);
1371 0 : return true;
1372 0 : };
1373 :
1374 : // Make a plain straightforward honest-to-FSM image. This returns
1375 : // a complex image, without conversion to Stokes. The representation
1376 : // is that required for the visibilities.
1377 : //----------------------------------------------------------------------
1378 0 : void FTMachine::makeImage(FTMachine::Type type,
1379 : ROVisibilityIterator& vi,
1380 : ImageInterface<Complex>& theImage,
1381 : Matrix<Float>& weight) {
1382 :
1383 :
1384 0 : logIO() << LogOrigin("FTMachine", "makeImage0") << LogIO::NORMAL;
1385 :
1386 : // Loop over all visibilities and pixels
1387 0 : VisBuffer vb(vi);
1388 :
1389 : // Initialize put (i.e. transform to Sky) for this model
1390 0 : vi.origin();
1391 :
1392 0 : if(vb.polFrame()==MSIter::Linear) {
1393 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
1394 : }
1395 : else {
1396 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
1397 : }
1398 :
1399 0 : initializeToSky(theImage,weight,vb);
1400 0 : Bool useCorrected= !(vi.msColumns().correctedData().isNull());
1401 0 : if((type==FTMachine::CORRECTED) && (!useCorrected))
1402 0 : type=FTMachine::OBSERVED;
1403 0 : Bool normalize=true;
1404 0 : if(type==FTMachine::COVERAGE)
1405 0 : normalize=false;
1406 :
1407 : // Loop over the visibilities, putting VisBuffers
1408 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1409 0 : for (vi.origin(); vi.more(); vi++) {
1410 :
1411 0 : switch(type) {
1412 0 : case FTMachine::RESIDUAL:
1413 0 : vb.visCube()=vb.correctedVisCube();
1414 0 : vb.visCube()-=vb.modelVisCube();
1415 0 : put(vb, -1, false);
1416 0 : break;
1417 0 : case FTMachine::MODEL:
1418 0 : put(vb, -1, false, FTMachine::MODEL);
1419 0 : break;
1420 0 : case FTMachine::CORRECTED:
1421 0 : put(vb, -1, false, FTMachine::CORRECTED);
1422 0 : break;
1423 0 : case FTMachine::PSF:
1424 0 : vb.visCube()=Complex(1.0,0.0);
1425 0 : put(vb, -1, true, FTMachine::PSF);
1426 0 : break;
1427 0 : case FTMachine::COVERAGE:
1428 0 : vb.visCube()=Complex(1.0);
1429 0 : put(vb, -1, true, FTMachine::COVERAGE);
1430 0 : break;
1431 0 : case FTMachine::OBSERVED:
1432 : default:
1433 0 : put(vb, -1, false, FTMachine::OBSERVED);
1434 0 : break;
1435 : }
1436 : }
1437 : }
1438 0 : finalizeToSky();
1439 : // Normalize by dividing out weights, etc.
1440 0 : getImage(weight, normalize);
1441 0 : }
1442 :
1443 :
1444 : // Make a plain straightforward honest-to-God image. This returns
1445 : // a complex image, without conversion to Stokes. The representation
1446 : // is that required for the visibilities. This version always works
1447 : // but has unnecessary operations for synthesis.
1448 : //----------------------------------------------------------------------
1449 0 : void FTMachine::makeImage(FTMachine::Type type,
1450 : VisSet& vs,
1451 : ImageInterface<Complex>& theImage,
1452 : Matrix<Float>& weight) {
1453 :
1454 0 : logIO() << LogOrigin("FTMachine", "makeImage") << LogIO::NORMAL;
1455 :
1456 : // If we want to calculate the PSF, we'll have to fill in the MODEL_DATA
1457 : // column
1458 0 : if(type==FTMachine::PSF) {
1459 :
1460 0 : VisIter& vi(vs.iter());
1461 :
1462 : // Loop over all visibilities and pixels
1463 0 : VisBuffer vb(vi);
1464 :
1465 : // Initialize put (i.e. transform to Sky) for this model
1466 0 : vi.origin();
1467 :
1468 0 : logIO() << "Calculating MODEL_DATA column from point source model" << LogIO::POST;
1469 0 : TempImage<Float> pointImage(theImage.shape(), theImage.coordinates());
1470 0 : IPosition start(4, theImage.shape()(0)/2, theImage.shape()(1)/2, 0, 0);
1471 0 : IPosition shape(4, 1, 1, 1, theImage.shape()(3));
1472 0 : Array<Float> line(shape);
1473 0 : pointImage.set(0.0);
1474 0 : line=1.0;
1475 0 : pointImage.putSlice(line, start);
1476 0 : TempImage<Complex> cPointImage(theImage.shape(), theImage.coordinates());
1477 0 : StokesImageUtil::From(cPointImage, pointImage);
1478 0 : if(vb.polFrame()==MSIter::Linear) {
1479 0 : StokesImageUtil::changeCStokesRep(cPointImage, StokesImageUtil::LINEAR);
1480 : }
1481 : else {
1482 0 : StokesImageUtil::changeCStokesRep(cPointImage, StokesImageUtil::CIRCULAR);
1483 : }
1484 0 : initializeToVis(cPointImage, vb);
1485 : // Loop over all visibilities
1486 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1487 0 : for (vi.origin(); vi.more(); vi++) {
1488 0 : get(vb, -1);
1489 0 : vi.setVis(vb.modelVisCube(),VisibilityIterator::Model);
1490 : }
1491 : }
1492 0 : finalizeToVis();
1493 0 : }
1494 :
1495 0 : ROVisIter& vi(vs.iter());
1496 :
1497 : // Loop over all visibilities and pixels
1498 0 : VisBuffer vb(vi);
1499 :
1500 : // Initialize put (i.e. transform to Sky) for this model
1501 0 : vi.origin();
1502 :
1503 : // Initialize put (i.e. transform to Sky) for this model
1504 0 : vi.origin();
1505 :
1506 0 : if(vb.polFrame()==MSIter::Linear) {
1507 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
1508 : }
1509 : else {
1510 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
1511 : }
1512 :
1513 0 : initializeToSky(theImage,weight,vb);
1514 0 : Bool useCorrected= !(vi.msColumns().correctedData().isNull());
1515 0 : if((type==FTMachine::CORRECTED) && (!useCorrected))
1516 0 : type=FTMachine::OBSERVED;
1517 :
1518 :
1519 : // Loop over the visibilities, putting VisBuffers
1520 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1521 0 : for (vi.origin(); vi.more(); vi++) {
1522 :
1523 0 : switch(type) {
1524 0 : case FTMachine::RESIDUAL:
1525 0 : vb.visCube()=vb.correctedVisCube();
1526 0 : vb.visCube()-=vb.modelVisCube();
1527 0 : put(vb, -1, false);
1528 0 : break;
1529 0 : case FTMachine::PSF:
1530 : case FTMachine::MODEL:
1531 : //vb.visCube()=vb.modelVisCube();
1532 0 : put(vb, -1, false, FTMachine::MODEL);
1533 0 : break;
1534 0 : case FTMachine::CORRECTED:
1535 : //vb.visCube()=vb.correctedVisCube();
1536 0 : put(vb, -1, false, FTMachine::CORRECTED);
1537 0 : break;
1538 0 : case FTMachine::COVERAGE:
1539 0 : vb.visCube()=Complex(1.0);
1540 0 : put(vb, -1, true);
1541 0 : break;
1542 0 : case FTMachine::OBSERVED:
1543 : default:
1544 0 : put(vb, -1, false);
1545 0 : break;
1546 : }
1547 : }
1548 : }
1549 0 : finalizeToSky();
1550 : // Normalize by dividing out weights, etc.
1551 0 : getImage(weight, true);
1552 :
1553 0 : }
1554 :
1555 :
1556 47 : Bool FTMachine::setSpw(Vector<Int>& spws, Bool validFrame){
1557 :
1558 47 : freqFrameValid_p=validFrame;
1559 47 : if(spws.nelements() >= 1){
1560 33 : selectedSpw_p.resize();
1561 33 : selectedSpw_p=spws;
1562 33 : multiChanMap_p.resize(max(spws)+1);
1563 33 : return true;
1564 : }
1565 :
1566 14 : return false;
1567 : }
1568 :
1569 63 : Bool FTMachine::matchAllSpwChans(const VisBuffer& vb){
1570 :
1571 63 : vb.allSelectedSpectralWindows(selectedSpw_p, nVisChan_p);
1572 :
1573 :
1574 63 : doConversion_p.resize(max(selectedSpw_p)+1);
1575 63 : doConversion_p.set(false);
1576 :
1577 63 : multiChanMap_p.resize(max(selectedSpw_p)+1, true);
1578 :
1579 63 : Bool anymatchChan=false;
1580 63 : Bool anyTopo=false;
1581 126 : for (uInt k=0; k < selectedSpw_p.nelements(); ++k){
1582 63 : Bool matchthis=matchChannel(selectedSpw_p[k], vb);
1583 63 : anymatchChan= (anymatchChan || matchthis);
1584 63 : anyTopo=anyTopo || ((MFrequency::castType(vb.msColumns().spectralWindow().measFreqRef()(selectedSpw_p[k]))==MFrequency::TOPO) && freqFrameValid_p);
1585 : }
1586 :
1587 : // if TOPO and valid frame things may match later but not now thus we'll go
1588 : // through the data
1589 : // hoping the user made the right choice
1590 63 : if (!anymatchChan && !anyTopo){
1591 : logIO() << "No overlap in frequency between image channels and selected data found for this FTMachine \n"
1592 : << " Check your data selection and image parameters if you end up with a blank image"
1593 0 : << LogIO::WARN << LogIO::POST;
1594 :
1595 : }
1596 :
1597 63 : return true;
1598 :
1599 : }
1600 :
1601 63 : Bool FTMachine::matchChannel(const Int& spw,
1602 : const VisBuffer& vb){
1603 :
1604 63 : if(nVisChan_p[spw] < 0)
1605 : logIO() << " Spectral window " << spw
1606 0 : << " does not seem to have been selected" << LogIO::EXCEPTION;
1607 63 : nvischan = nVisChan_p[spw];
1608 63 : chanMap.resize(nvischan);
1609 63 : chanMap.set(-1);
1610 63 : Vector<Double> lsrFreq(0);
1611 63 : Bool condoo=false;
1612 :
1613 : //cerr << "doConve " << spw << " " << doConversion_p[spw] << " freqframeval " << freqFrameValid_p << endl;
1614 :
1615 63 : vb.lsrFrequency(spw, lsrFreq, condoo, !freqFrameValid_p);
1616 63 : doConversion_p[spw]=condoo;
1617 :
1618 63 : if(lsrFreq.nelements() ==0){
1619 0 : return false;
1620 : }
1621 63 : lsrFreq_p.resize(lsrFreq.nelements());
1622 63 : lsrFreq_p=lsrFreq;
1623 63 : Vector<Double> c(1);
1624 63 : c=0.0;
1625 63 : Vector<Double> f(1);
1626 63 : Int nFound=0;
1627 :
1628 :
1629 : //cout.precision(10);
1630 241 : for (Int chan=0;chan<nvischan;chan++) {
1631 178 : f(0)=lsrFreq[chan];
1632 178 : if(spectralCoord_p.toPixel(c, f)) {
1633 178 : Int pixel=Int(floor(c(0)+0.5)); // round to chan freq at chan center
1634 : //cerr << "spw " << spw << " f " << f(0) << " pixel "<< c(0) << " " << pixel << endl;
1635 : /////////////
1636 : //c(0)=pixel;
1637 : //spectralCoord_p.toWorld(f, c);
1638 : //cerr << "f1 " << f(0) << " pixel "<< c(0) << " " << pixel << endl;
1639 : ////////////////
1640 178 : if(pixel>-1&&pixel<nchan) {
1641 79 : chanMap(chan)=pixel;
1642 79 : nFound++;
1643 79 : if(nvischan>1&&(chan==0||chan==nvischan-1)) {
1644 : /*logIO() << LogIO::DEBUGGING
1645 : << "Selected visibility channel : " << chan+1
1646 : << " has frequency "
1647 : << MFrequency(Quantity(f(0), "Hz")).get("GHz").getValue()
1648 : << " GHz and maps to image pixel " << pixel+1 << LogIO::POST;
1649 : */
1650 : }
1651 : }
1652 : else{
1653 :
1654 :
1655 99 : if(nvischan > 1){
1656 99 : Double fwidth=lsrFreq[1]-lsrFreq[0];
1657 99 : Double limit=0;
1658 99 : Double where=c(0)*fabs(spectralCoord_p.increment()(0));
1659 99 : if( freqInterpMethod_p==InterpolateArray1D<Double,Complex>::linear)
1660 0 : limit=1;
1661 99 : else if( freqInterpMethod_p==InterpolateArray1D<Double,Complex>::cubic || freqInterpMethod_p==InterpolateArray1D<Double,Complex>::spline)
1662 0 : limit=2;
1663 99 : if(((pixel<0) && (where >= (0-limit*fabs(fwidth)))) )
1664 0 : chanMap(chan)=-2;
1665 99 : if((pixel>=nchan) ) {
1666 99 : where=f(0);
1667 : Double fend;
1668 99 : spectralCoord_p.toWorld(fend, Double(nchan));
1669 99 : if( ( (fwidth >0) &&where < (fend+limit*fwidth)) || ( (fwidth <0) &&where > (fend+limit*fwidth)) )
1670 0 : chanMap(chan)=-2;
1671 : }
1672 : }
1673 : }
1674 : }
1675 : }
1676 :
1677 63 : multiChanMap_p[spw].resize();
1678 63 : multiChanMap_p[spw]=chanMap;
1679 :
1680 63 : if(nFound==0) {
1681 : /*
1682 : logIO() << "Visibility channels in spw " << spw+1
1683 : << " of ms " << vb.msId() << " is not being used "
1684 : << LogIO::WARN << LogIO::POST;
1685 : */
1686 0 : return false;
1687 : }
1688 :
1689 :
1690 :
1691 :
1692 63 : return true;
1693 :
1694 63 : }
1695 :
1696 :
1697 :
1698 :
1699 :
1700 821 : void FTMachine::gridOk(Int convSupport){
1701 :
1702 821 : if (nx <= 2*convSupport) {
1703 0 : logIO_p
1704 : << "number of pixels "
1705 : << nx << " on x axis is smaller than or equals to the gridding support "
1706 : << 2*convSupport << " Please use a larger value "
1707 0 : << LogIO::EXCEPTION;
1708 : }
1709 :
1710 821 : if (ny <= 2*convSupport) {
1711 0 : logIO_p
1712 : << "number of pixels "
1713 : << ny << " on y axis is smaller than or equals to the gridding support "
1714 : << 2*convSupport << " Please use a larger value "
1715 0 : << LogIO::EXCEPTION;
1716 : }
1717 :
1718 821 : }
1719 :
1720 0 : void FTMachine::setLocation(const MPosition& loc){
1721 :
1722 0 : mLocation_p=loc;
1723 :
1724 0 : }
1725 :
1726 0 : MPosition& FTMachine::getLocation(){
1727 :
1728 0 : return mLocation_p;
1729 : }
1730 :
1731 :
1732 0 : void FTMachine::setMovingSource(const String& sourcename){
1733 :
1734 0 : fixMovingSource_p=true;
1735 0 : movingDir_p=MDirection(Quantity(0.0,"deg"), Quantity(90.0, "deg"));
1736 0 : movingDir_p.setRefString(sourcename);
1737 :
1738 0 : }
1739 :
1740 0 : void FTMachine::setMovingSource(const MDirection& mdir){
1741 :
1742 0 : fixMovingSource_p=true;
1743 0 : movingDir_p=mdir;
1744 :
1745 0 : }
1746 :
1747 14 : void FTMachine::setFreqInterpolation(const String& method){
1748 :
1749 14 : String meth=method;
1750 14 : meth.downcase();
1751 14 : if(meth.contains("linear")){
1752 0 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::linear;
1753 : }
1754 14 : else if(meth.contains("splin")){
1755 0 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::spline;
1756 : }
1757 14 : else if(meth.contains("cub")){
1758 0 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::cubic;
1759 : }
1760 : else{
1761 14 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::nearestNeighbour;
1762 : }
1763 :
1764 14 : }
1765 :
1766 :
1767 : // helper function to swap the y and z axes of a Cube
1768 0 : void FTMachine::swapyz(Cube<Complex>& out, const Cube<Complex>& in)
1769 : {
1770 0 : IPosition inShape=in.shape();
1771 0 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
1772 : //resize breaks references...so out better have the right shape
1773 : //if references is not to be broken
1774 0 : if(out.nelements()==0)
1775 0 : out.resize(nxx,nyy,nzz);
1776 : Bool deleteIn,deleteOut;
1777 0 : const Complex* pin = in.getStorage(deleteIn);
1778 0 : Complex* pout = out.getStorage(deleteOut);
1779 0 : uInt i=0, zOffset=0;
1780 0 : for (uInt iz=0; iz<nzz; ++iz, zOffset+=nxx) {
1781 0 : Int yOffset=zOffset;
1782 0 : for (uInt iy=0; iy<nyy; ++iy, yOffset+=nxx*nzz) {
1783 0 : for (uInt ix=0; ix<nxx; ++ix){
1784 0 : pout[i++] = pin[ix+yOffset];
1785 : }
1786 : }
1787 : }
1788 0 : out.putStorage(pout,deleteOut);
1789 0 : in.freeStorage(pin,deleteIn);
1790 0 : }
1791 0 : void FTMachine::swapyz(Cube<Complex>& out, const Cube<Bool>& outFlag, const Cube<Complex>& in)
1792 : {
1793 0 : IPosition inShape=in.shape();
1794 0 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
1795 : //resize breaks references...so out better have the right shape
1796 : //if references is not to be broken
1797 0 : if(out.nelements()==0)
1798 0 : out.resize(nxx,nyy,nzz);
1799 : Bool deleteIn,deleteOut, delFlag;
1800 0 : const Complex* pin = in.getStorage(deleteIn);
1801 0 : const Bool* poutflag= outFlag.getStorage(delFlag);
1802 0 : Complex* pout = out.getStorage(deleteOut);
1803 0 : uInt i=0, zOffset=0;
1804 0 : for (uInt iz=0; iz<nzz; ++iz, zOffset+=nxx) {
1805 0 : Int yOffset=zOffset;
1806 0 : for (uInt iy=0; iy<nyy; ++iy, yOffset+=nxx*nzz) {
1807 0 : for (uInt ix=0; ix<nxx; ++ix){
1808 0 : if(!poutflag[i] )
1809 0 : pout[i] = pin[ix+yOffset];
1810 0 : ++i;
1811 : }
1812 : }
1813 : }
1814 0 : out.putStorage(pout,deleteOut);
1815 0 : in.freeStorage(pin,deleteIn);
1816 0 : outFlag.freeStorage(poutflag, delFlag);
1817 0 : }
1818 : // helper function to swap the y and z axes of a Cube
1819 0 : void FTMachine::swapyz(Cube<Bool>& out, const Cube<Bool>& in)
1820 : {
1821 0 : IPosition inShape=in.shape();
1822 0 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
1823 0 : if(out.nelements()==0)
1824 0 : out.resize(nxx,nyy,nzz);
1825 : Bool deleteIn,deleteOut;
1826 0 : const Bool* pin = in.getStorage(deleteIn);
1827 0 : Bool* pout = out.getStorage(deleteOut);
1828 0 : uInt i=0, zOffset=0;
1829 0 : for (uInt iz=0; iz<nzz; iz++, zOffset+=nxx) {
1830 0 : Int yOffset=zOffset;
1831 0 : for (uInt iy=0; iy<nyy; iy++, yOffset+=nxx*nzz) {
1832 0 : for (uInt ix=0; ix<nxx; ix++) pout[i++] = pin[ix+yOffset];
1833 : }
1834 : }
1835 0 : out.putStorage(pout,deleteOut);
1836 0 : in.freeStorage(pin,deleteIn);
1837 0 : }
1838 :
1839 0 : void FTMachine::setPointingDirColumn(const String& column){
1840 0 : pointingDirCol_p=column;
1841 0 : pointingDirCol_p.upcase();
1842 0 : if( (pointingDirCol_p != "DIRECTION") &&(pointingDirCol_p != "TARGET") && (pointingDirCol_p != "ENCODER") && (pointingDirCol_p != "POINTING_OFFSET") && (pointingDirCol_p != "SOURCE_OFFSET")){
1843 :
1844 : //basically at this stage you don't know what you're doing...so you get the default
1845 :
1846 0 : pointingDirCol_p="DIRECTION";
1847 :
1848 : }
1849 0 : }
1850 :
1851 0 : String FTMachine::getPointingDirColumnInUse(){
1852 :
1853 0 : return pointingDirCol_p;
1854 :
1855 : }
1856 :
1857 14 : void FTMachine::setSpwChanSelection(const Cube<Int>& spwchansels) {
1858 14 : spwChanSelFlag_p.resize();
1859 14 : spwChanSelFlag_p=spwchansels;
1860 14 : }
1861 :
1862 14 : void FTMachine::setSpwFreqSelection(const Matrix<Double>& spwFreqs)
1863 : {
1864 14 : spwFreqSel_p.assign(spwFreqs);
1865 14 : SynthesisUtils::expandFreqSelection(spwFreqs,expandedSpwFreqSel_p, expandedSpwConjFreqSel_p);
1866 : //cerr << expandedSpwFreqSel_p << endl;
1867 14 : }
1868 :
1869 821 : void FTMachine::setSpectralFlag(const VisBuffer& vb, Cube<Bool>& modflagcube){
1870 :
1871 821 : modflagcube.resize(vb.flagCube().shape());
1872 821 : modflagcube=vb.flagCube();
1873 821 : uInt nchan = vb.nChannel();
1874 821 : uInt msid = vb.msId();
1875 821 : uInt selspw = vb.spectralWindow();
1876 821 : Bool spwFlagIsSet=( (spwChanSelFlag_p.nelements() > 0) && (spwChanSelFlag_p.shape()(1) > selspw) &&
1877 821 : (spwChanSelFlag_p.shape()(0) > msid) &&
1878 0 : (spwChanSelFlag_p.shape()(2) >=nchan));
1879 : //cerr << "spwFlagIsSet " << spwFlagIsSet << endl;
1880 23291 : for (uInt i=0;i<nchan;i++) {
1881 : //Flag those channels that did not get selected...
1882 : //respect the flags from vb if selected or
1883 : //if spwChanSelFlag is wrong shape
1884 22470 : if ((spwFlagIsSet) && (spwChanSelFlag_p(msid,selspw,i)!=1)) {
1885 0 : modflagcube.xzPlane(i).set(true);
1886 : }
1887 : }
1888 821 : }
1889 :
1890 : //-----------------------------------------------------------------------------------------------------------------
1891 : //------------ Vectorized versions of initializeToVis, initializeToSky, finalizeToSky
1892 : //------------ that are called from CubeSkyEquation.
1893 : //------------ They call getImage,getWeightImage, which are implemented in all FTMs
1894 : //------------ Also, Correlation / Stokes conversions and gS/ggS normalizations.
1895 :
1896 : // Vectorized InitializeToVis
1897 0 : void FTMachine::initializeToVis(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec,
1898 : PtrBlock<SubImage<Float> *> & modelImageVec,
1899 : PtrBlock<SubImage<Float> *>& weightImageVec,
1900 : PtrBlock<SubImage<Float> *>& fluxScaleVec,
1901 : Block<Matrix<Float> >& weightsVec,
1902 : const VisBuffer& vb)
1903 : {
1904 0 : AlwaysAssert(compImageVec.nelements()==1, AipsError);
1905 0 : AlwaysAssert(modelImageVec.nelements()==1, AipsError);
1906 0 : AlwaysAssert(weightImageVec.nelements()==1, AipsError);
1907 0 : AlwaysAssert(fluxScaleVec.nelements()==1, AipsError);
1908 0 : AlwaysAssert(weightsVec.nelements()==1, AipsError);
1909 0 : Matrix<Float> tempWts;
1910 :
1911 : // Convert from Stokes planes to Correlation planes
1912 0 : stokesToCorrelation(*(modelImageVec[0]), *(compImageVec[0]));
1913 0 : if(sj_p.nelements() >0 ){
1914 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
1915 0 : (sj_p(k))->apply(*(compImageVec[0]), *(compImageVec[0]), vb, 0, true);
1916 : }
1917 : }
1918 : // Call initializeToVis
1919 0 : initializeToVis(*(compImageVec[0]), vb); // Pure virtual
1920 :
1921 0 : };
1922 :
1923 : // Vectorized finalizeToVis is not implemented because it does nothing and is never called.
1924 :
1925 : // Vectorized InitializeToSky
1926 0 : void FTMachine::initializeToSky(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec,
1927 : Block<Matrix<Float> >& weightsVec,
1928 : const VisBuffer& vb,
1929 : const Bool /*dopsf*/)
1930 :
1931 : {
1932 0 : AlwaysAssert(compImageVec.nelements()==1, AipsError);
1933 0 : AlwaysAssert(weightsVec.nelements()==1, AipsError);
1934 :
1935 0 : initializeToSky(*(compImageVec[0]) , weightsVec[0] , vb);
1936 0 : };
1937 :
1938 : // Vectorized finalizeToSky
1939 0 : void FTMachine::finalizeToSky(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec,
1940 : PtrBlock<SubImage<Float> *> & resImageVec,
1941 : PtrBlock<SubImage<Float> *>& weightImageVec,
1942 : PtrBlock<SubImage<Float> *>& fluxScaleVec,
1943 : Bool dopsf,
1944 : Block<Matrix<Float> >& weightsVec, const VisBuffer& vb)
1945 : {
1946 : // Call default finalizeToSky
1947 0 : finalizeToSky(); // Pure virtual
1948 :
1949 : // Check vector lengths.
1950 0 : AlwaysAssert(compImageVec.nelements()==1, AipsError);
1951 0 : AlwaysAssert(resImageVec.nelements()==1, AipsError);
1952 0 : AlwaysAssert(weightImageVec.nelements()==1, AipsError);
1953 0 : AlwaysAssert(fluxScaleVec.nelements()==1, AipsError);
1954 0 : AlwaysAssert(weightsVec.nelements()==1, AipsError);
1955 :
1956 : // Get the gridded image
1957 0 : (*(compImageVec[0])).copyData(getImage(weightsVec[0],false)); // getImage is Pure virtual
1958 0 : Bool doSky=(sj_p.nelements()>0) && (!dopsf);
1959 0 : TempImage<Float> work;
1960 : // Convert from correlation (complex) to Stokes (float)
1961 0 : if (doSky){
1962 0 : work=TempImage<Float>((compImageVec[0])->shape(), (compImageVec[0])->coordinates());
1963 0 : work.set(0);
1964 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
1965 0 : (sj_p(k))->apply(*(compImageVec[0]), *(compImageVec[0]), vb, 0, false);
1966 : }
1967 0 : correlationToStokes(*(compImageVec[0]), work , dopsf);
1968 0 : LatticeExpr<Float> le((*(resImageVec[0]))+work);
1969 0 : (resImageVec[0])->copyData(le);
1970 0 : getWeightImage(work , weightsVec[0]);
1971 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
1972 0 : (sj_p(k))->applySquare(work, work, vb, 0);
1973 : }
1974 0 : (weightImageVec[0])->copyData((LatticeExpr<Float>)((*(weightImageVec[0]))+work)) ;
1975 0 : }
1976 : else{
1977 0 : correlationToStokes(*(compImageVec[0]), *(resImageVec[0]), dopsf);
1978 0 : getWeightImage((*(weightImageVec[0])), weightsVec[0]); // Pure virtual
1979 : }
1980 :
1981 : // ForSB //
1982 : // For FTM, normalizeImage should get called with normtype=0 (only sumwt normalization)
1983 : // For AWP, call getWeightImage, and then normalizeImage with normtype=2 (sumwt and pb)
1984 : // Currently, both call it with normtype=2 because for GridFT, "pb" is filled with ones.
1985 :
1986 : // Normalize the output image by the sum of wts and sensitivity image
1987 : // cerr << "FTM::finalizeToSky: Weights = " << weightsVec[0] << endl;
1988 : // storeImg(String("wt.im"),*(weightImageVec[0]));
1989 : // storeImg(String("stokes.im"),*(resImageVec[0]));
1990 :
1991 0 : normalizeImage( *(resImageVec[0]) , weightsVec[0], *(weightImageVec[0]) , dopsf,
1992 : // (Float)1e-03,
1993 0 : (Float)pbLimit_p,
1994 : doSky? (Int)6 : (Int)0); // Normalize by sum-of-wts.
1995 : // (Int)2); // Normalize by (sum-of-wts*avgPB)
1996 :
1997 : // storeImg(String("stokes1.im"),*(resImageVec[0]));
1998 :
1999 0 : return;
2000 0 : };
2001 :
2002 0 : void FTMachine::setSkyJones(Vector<CountedPtr<SkyJones> >& sj){
2003 0 : sj_p.resize();
2004 0 : sj_p=sj;
2005 : // cout << "FTM : Set Sky Jones of length " << sj_p.nelements() << " and types ";
2006 0 : for( uInt i=0;i<sj_p.nelements();i++) cout << sj_p[i]->type() << " ";
2007 0 : cout << endl;
2008 0 : }
2009 : // Convert complex correlation planes to float Stokes planes
2010 0 : void FTMachine::correlationToStokes(ImageInterface<Complex>& compImage,
2011 : ImageInterface<Float>& resImage,
2012 : const Bool dopsf)
2013 : {
2014 : // Convert correlation image to IQUV format
2015 0 : AlwaysAssert(compImage.shape()[0]==resImage.shape()[0], AipsError);
2016 0 : AlwaysAssert(compImage.shape()[1]==resImage.shape()[1], AipsError);
2017 0 : AlwaysAssert(compImage.shape()[3]==resImage.shape()[3], AipsError);
2018 :
2019 0 : if(dopsf)
2020 : {
2021 : // For the PSF, choose only those stokes planes that have a valid PSF
2022 0 : StokesImageUtil::ToStokesPSF(resImage,compImage);
2023 : }
2024 : else
2025 : {
2026 0 : StokesImageUtil::To(resImage,compImage);
2027 : }
2028 0 : };
2029 :
2030 : // Convert float Stokes planes to complex correlation planes
2031 0 : void FTMachine::stokesToCorrelation(ImageInterface<Float>& modelImage,
2032 : ImageInterface<Complex>& compImage)
2033 : {
2034 : /*
2035 : StokesCoordinate stcomp=compImage.coordinates().stokesCoordinate(compImage.coordinates().findCoordinate(Coordinate::STOKES));
2036 : StokesCoordinate stfloat = modelImage.coordinates().stokesCoordinate(modelImage.coordinates().findCoordinate(Coordinate::STOKES));
2037 :
2038 : cout << "Stokes types : complex : " << stcomp.stokes() << " float : " << stfloat.stokes() << endl;
2039 : cout << "Shapes : complex : " << compImage.shape() << " float : " << modelImage.shape() << endl;
2040 : */
2041 :
2042 : //Pol axis need not be same
2043 0 : AlwaysAssert(modelImage.shape()[0]==compImage.shape()[0], AipsError);
2044 0 : AlwaysAssert(modelImage.shape()[1]==compImage.shape()[1], AipsError);
2045 0 : AlwaysAssert(modelImage.shape()[3]==compImage.shape()[3], AipsError);
2046 : // Convert from Stokes to Complex
2047 0 : StokesImageUtil::From(compImage, modelImage);
2048 0 : };
2049 :
2050 : //------------------------------------------------------------------------------------------------------------------
2051 :
2052 0 : void FTMachine::normalizeImage(ImageInterface<Float>& skyImage,
2053 : Matrix<Float>& sumOfWts,
2054 : ImageInterface<Float>& sensitivityImage,
2055 : Bool dopsf, Float pblimit, Int normtype)
2056 : {
2057 :
2058 : //Normalize the sky Image
2059 0 : Int nXX=(skyImage).shape()(0);
2060 0 : Int nYY=(skyImage).shape()(1);
2061 0 : Int npola= (skyImage).shape()(2);
2062 0 : Int nchana= (skyImage).shape()(3);
2063 :
2064 0 : IPosition pcentre(4,nXX/2,nYY/2,0,0);
2065 : // IPosition psource(4,nXX/2+22,nYY/2,0,0);
2066 :
2067 : // storeImg(String("norm_resimage.im") , skyImage);
2068 : // storeImg(String("norm_sensitivity.im"), sensitivityImage);
2069 :
2070 : ///// cout << "FTM::norm : pblimit : " << pblimit << endl;
2071 :
2072 : // Note : This is needed because initial prediction has no info about sumwt.
2073 : // Not a clean solution. // ForSB -- if you see a better way, go for it.
2074 0 : if(sumOfWts.shape() != IPosition(2,npola,nchana))
2075 : {
2076 0 : cout << "Empty SumWt.... resizing and filling with 1.0" << endl;
2077 0 : sumOfWts.resize(IPosition(2,npola,nchana));
2078 0 : sumOfWts=1.0;
2079 : }
2080 :
2081 : // if(dopsf) cout << "*** FTM::normalizeImage : Image Center : " << skyImage.getAt(pcentre) << " and weightImage : " << sensitivityImage.getAt(pcentre) << " SumWt : " << sumOfWts[0,0];
2082 : // else cout << "*** FTM::normalizeImage : Source Loc : " << skyImage.getAt(psource) << " and weightImage : " << sensitivityImage.getAt(psource) << " SumWt : " << sumOfWts[0,0];
2083 :
2084 :
2085 :
2086 0 : IPosition blc(4,nXX, nYY, npola, nchana);
2087 0 : IPosition trc(4, nXX, nYY, npola, nchana);
2088 0 : blc(0)=0; blc(1)=0; trc(0)=nXX-1; trc(1)=nYY-1;
2089 : //max weights per plane
2090 0 : for (Int pol=0; pol < npola; ++pol){
2091 0 : for (Int chan=0; chan < nchana ; ++chan){
2092 :
2093 0 : blc(2)=pol; trc(2)=pol;
2094 0 : blc(3)=chan; trc(3)=chan;
2095 0 : Slicer sl(blc, trc, Slicer::endIsLast);
2096 0 : SubImage<Float> subSkyImage(skyImage, sl, false);
2097 0 : SubImage<Float> subSensitivityImage(sensitivityImage, sl, false);
2098 0 : SubImage<Float> subOutput(skyImage, sl, true);
2099 0 : Float sumWt = sumOfWts(pol,chan);
2100 0 : if(sumWt > 0.0){
2101 0 : switch(normtype)
2102 : {
2103 0 : case 0: // only sum Of Weights - FTM only (ForSB)
2104 0 : subOutput.copyData( (LatticeExpr<Float>)
2105 0 : ((dopsf?1.0:-1.0)*subSkyImage/(sumWt)));
2106 0 : break;
2107 :
2108 0 : case 1: // only sensitivityImage Ic/avgPB (ForSB)
2109 0 : subOutput.copyData( (LatticeExpr<Float>)
2110 0 : (iif(subSensitivityImage > (pblimit),
2111 0 : (subSkyImage/(subSensitivityImage)),
2112 : (subSkyImage))));
2113 : // 0.0)));
2114 0 : break;
2115 :
2116 0 : case 2: // sum of Weights and sensitivityImage IGridded/(SoW*avgPB) and PSF --> Id (ForSB)
2117 0 : subOutput.copyData( (LatticeExpr<Float>)
2118 0 : (iif(subSensitivityImage > (pblimit),
2119 0 : ((dopsf?1.0:-1.0)*subSkyImage/(subSensitivityImage*sumWt)),
2120 : //((dopsf?1.0:-1.0)*subSkyImage))));
2121 : 0.0)));
2122 0 : break;
2123 :
2124 0 : case 3: // MULTIPLY by the sensitivityImage avgPB
2125 0 : subOutput.copyData( (LatticeExpr<Float>) (subSkyImage * subSensitivityImage) );
2126 0 : break;
2127 :
2128 0 : case 4: // DIVIDE by sqrt of sensitivityImage
2129 0 : subOutput.copyData( (LatticeExpr<Float>)
2130 0 : (iif((subSensitivityImage) > (pblimit),
2131 0 : (subSkyImage/(sqrt(subSensitivityImage))),
2132 : (subSkyImage))));
2133 : //0.0)));
2134 0 : break;
2135 :
2136 0 : case 5: // MULTIPLY by sqrt of sensitivityImage
2137 0 : subOutput.copyData( (LatticeExpr<Float>)
2138 0 : (iif((subSensitivityImage) > (pblimit),
2139 0 : (subSkyImage * (sqrt(subSensitivityImage))),
2140 : (subSkyImage))));
2141 :
2142 0 : break;
2143 :
2144 0 : case 6: // divide by non normalized sensitivity image
2145 : {
2146 0 : Float elpblimit=max( subSensitivityImage).getFloat() * pblimit;
2147 0 : subOutput.copyData( (LatticeExpr<Float>)
2148 0 : (iif(subSensitivityImage > (elpblimit),
2149 0 : ((dopsf?1.0:-1.0)*subSkyImage/(subSensitivityImage)),
2150 : 0.0)));
2151 : }
2152 0 : break;
2153 0 : default:
2154 0 : throw(AipsError("Unrecognized norm-type in FTM::normalizeImage : "+String::toString(normtype)));
2155 : }
2156 : }
2157 : else{
2158 0 : subOutput.set(0.0);
2159 : }
2160 0 : }
2161 : }
2162 :
2163 : //if(dopsf) cout << " Normalized (" << normtype << ") Image Center : " << skyImage.getAt(pcentre) << endl;
2164 : // else cout << " Normalized (" << normtype << ") Source Loc : " << skyImage.getAt(psource) << endl;
2165 :
2166 0 : }
2167 :
2168 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2169 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2170 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2171 : ///// For use with the new framework
2172 : ///// (Sorry about these copies, but need to keep old system working)
2173 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2174 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2175 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2176 :
2177 : // Vectorized InitializeToVis
2178 0 : void FTMachine::initializeToVisNew(const VisBuffer& vb,
2179 : CountedPtr<SIImageStore> imstore)
2180 : {
2181 0 : AlwaysAssert(imstore->getNTaylorTerms(false)==1, AipsError);
2182 :
2183 0 : Matrix<Float> tempWts;
2184 :
2185 : // Convert from Stokes planes to Correlation planes
2186 0 : stokesToCorrelation(*(imstore->model()), *(imstore->forwardGrid()));
2187 :
2188 0 : if(vb.polFrame()==MSIter::Linear) {
2189 0 : StokesImageUtil::changeCStokesRep(*(imstore->forwardGrid()),
2190 : StokesImageUtil::LINEAR);
2191 : }
2192 : else {
2193 0 : StokesImageUtil::changeCStokesRep(*(imstore->forwardGrid()),
2194 : StokesImageUtil::CIRCULAR);
2195 : }
2196 :
2197 : //------------------------------------------------------------------------------------
2198 : // Image Mosaic only : Multiply the input model with the Primary Beam
2199 0 : if(sj_p.nelements() >0 ){
2200 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2201 0 : (sj_p(k))->apply(*(imstore->forwardGrid()), *(imstore->forwardGrid()), vb, 0, true);
2202 : }
2203 : }
2204 : //------------------------------------------------------------------------------------
2205 :
2206 : // Call initializeToVis
2207 0 : initializeToVis(*(imstore->forwardGrid()), vb); // Pure virtual
2208 :
2209 0 : };
2210 :
2211 : // Vectorized finalizeToVis is not implemented because it does nothing and is never called.
2212 :
2213 : // Vectorized InitializeToSky
2214 0 : void FTMachine::initializeToSkyNew(const Bool dopsf,
2215 : const VisBuffer& vb,
2216 : CountedPtr<SIImageStore> imstore)
2217 :
2218 : {
2219 0 : AlwaysAssert(imstore->getNTaylorTerms(false)==1, AipsError);
2220 :
2221 : // Make the relevant float grid.
2222 : // This is needed mainly for facetting (to set facet shapes), but is harmless for non-facetting.
2223 0 : if( dopsf ) { imstore->psf(); } else { imstore->residual(); }
2224 :
2225 : // Initialize the complex grid (i.e. tell FTMachine what array to use internally)
2226 0 : Matrix<Float> sumWeight;
2227 0 : initializeToSky(*(imstore->backwardGrid()) , sumWeight , vb);
2228 :
2229 0 : };
2230 :
2231 : // Vectorized finalizeToSky
2232 0 : void FTMachine::finalizeToSkyNew(Bool dopsf,
2233 : const VisBuffer& vb,
2234 : CountedPtr<SIImageStore> imstore )
2235 : {
2236 : // Check vector lengths.
2237 0 : AlwaysAssert( imstore->getNTaylorTerms(false)==1, AipsError);
2238 :
2239 0 : Matrix<Float> sumWeights;
2240 0 : finalizeToSky();
2241 :
2242 : //------------------------------------------------------------------------------------
2243 : // Straightforward case. No extra primary beams. No image mosaic
2244 0 : if(sj_p.nelements() == 0 )
2245 : {
2246 0 : correlationToStokes( getImage(sumWeights, false) , ( dopsf ? *(imstore->psf()) : *(imstore->residual()) ), dopsf);
2247 :
2248 0 : if( useWeightImage() && dopsf ) {
2249 0 : getWeightImage( *(imstore->weight()) , sumWeights);
2250 :
2251 : //cerr << "weightim val " << (imstore->weight())->getAt(IPosition(4, nx/2, ny/2, 0, 0)) << " nx, ny " << ny <<", " << ny << " sumWeights " << sumWeights << endl;
2252 :
2253 : // Fill weight image only once, during PSF generation. Remember.... it is normalized only once
2254 : // during PSF generation.
2255 : }
2256 :
2257 : // Take sumWeights from corrToStokes here....
2258 0 : Matrix<Float> sumWeightStokes( (imstore->sumwt())->shape()[2], (imstore->sumwt())->shape()[3] );
2259 0 : StokesImageUtil::ToStokesSumWt( sumWeightStokes, sumWeights );
2260 :
2261 0 : AlwaysAssert( ( (imstore->sumwt())->shape()[2] == sumWeightStokes.shape()[0] ) &&
2262 : ((imstore->sumwt())->shape()[3] == sumWeightStokes.shape()[1] ) , AipsError );
2263 :
2264 0 : (imstore->sumwt())->put( sumWeightStokes.reform((imstore->sumwt())->shape()) );
2265 :
2266 :
2267 0 : }
2268 : //------------------------------------------------------------------------------------
2269 : // Image Mosaic only : Multiply the residual, and weight image by the PB.
2270 : else
2271 : {
2272 :
2273 : // Take the FT of the gridded values. Writes into backwardGrid().
2274 0 : getImage(sumWeights, false);
2275 :
2276 : // Multiply complex image grid by PB.
2277 0 : if( !dopsf )
2278 : {
2279 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2280 0 : (sj_p(k))->apply(*(imstore->backwardGrid()), *(imstore->backwardGrid()), vb, 0, true);
2281 : }
2282 : }
2283 :
2284 : // Convert from correlation to Stokes onto a new temporary grid.
2285 0 : SubImage<Float> targetImage( ( dopsf ? *(imstore->psf()) : *(imstore->residual()) ) , true);
2286 0 : TempImage<Float> temp( targetImage.shape(), targetImage.coordinates() );
2287 0 : correlationToStokes( *(imstore->backwardGrid()), temp, false);
2288 :
2289 : // Add the temporary Stokes image to the residual or PSF, whichever is being made.
2290 0 : LatticeExpr<Float> addToRes( targetImage + temp );
2291 0 : targetImage.copyData(addToRes);
2292 :
2293 : // Now, do the same with the weight image and sumwt ( only on the first pass )
2294 0 : if( dopsf )
2295 : {
2296 0 : SubImage<Float> weightImage( *(imstore->weight()) , true);
2297 0 : TempImage<Float> temp(weightImage.shape(), weightImage.coordinates());
2298 0 : getWeightImage(temp, sumWeights);
2299 :
2300 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2301 0 : (sj_p(k))->applySquare(temp,temp, vb, -1);
2302 : }
2303 :
2304 0 : LatticeExpr<Float> addToWgt( weightImage + temp );
2305 0 : weightImage.copyData(addToWgt);
2306 :
2307 :
2308 0 : AlwaysAssert( ( (imstore->sumwt())->shape()[2] == sumWeights.shape()[0] ) &&
2309 : ((imstore->sumwt())->shape()[3] == sumWeights.shape()[1] ) , AipsError );
2310 :
2311 0 : SubImage<Float> sumwtImage( *(imstore->sumwt()) , true);
2312 0 : TempImage<Float> temp2(sumwtImage.shape(), sumwtImage.coordinates());
2313 0 : temp2.put( sumWeights.reform(sumwtImage.shape()) );
2314 0 : LatticeExpr<Float> addToWgt2( sumwtImage + temp2 );
2315 0 : sumwtImage.copyData(addToWgt2);
2316 :
2317 :
2318 : //cout << "In finalizeGridCoreMos : sumwt : " << sumwtImage.get() << endl;
2319 :
2320 0 : }
2321 :
2322 0 : }///image mosaic only
2323 : //------------------------------------------------------------------------------------
2324 :
2325 :
2326 :
2327 0 : return;
2328 0 : };
2329 :
2330 0 : Bool FTMachine::changedSkyJonesLogic(const VisBuffer& vb, Bool& firstRow, Bool& internalRow)
2331 : {
2332 0 : firstRow=false;
2333 0 : internalRow=false;
2334 :
2335 0 : if( sj_p.nelements()==0 )
2336 0 : {throw(AipsError("Internal Error : Checking changedSkyJones, but it is not yet set."));}
2337 :
2338 0 : CountedPtr<SkyJones> ej = sj_p[0];
2339 0 : if(ej.null())
2340 0 : return false;
2341 0 : if(ej->changed(vb,0))
2342 0 : firstRow=true;
2343 0 : Int row2temp=0;
2344 0 : if(ej->changedBuffer(vb,0,row2temp)) {
2345 0 : internalRow=true;
2346 : }
2347 0 : return (firstRow || internalRow) ;
2348 0 : }
2349 :
2350 :
2351 0 : void FTMachine::setCFCache(CountedPtr<CFCache>& /*cfc*/, const Bool /*loadCFC*/)
2352 : {
2353 0 : throw(AipsError("FTMachine::setCFCache() directly called!"));
2354 : }
2355 : /*
2356 : /// Move to individual FTMs............ make it pure virtual.
2357 : Bool FTMachine::useWeightImage()
2358 : {
2359 : if( name() == "GridFT" || name() == "WProjectFT" )
2360 : { return false; }
2361 : else
2362 : { return true; }
2363 : }
2364 : */
2365 :
2366 : } //# NAMESPACE CASA - END
2367 :
|
