Line data Source code
1 : //# MultiTermFTNew.cc: Implementation of MultiTermFTNew class
2 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
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 :
28 : #include <msvis/MSVis/VisBuffer2.h>
29 : #include <casacore/images/Images/ImageInterface.h>
30 : #include <casacore/images/Images/PagedImage.h>
31 : #include <casacore/casa/Containers/Block.h>
32 : #include <casacore/casa/Containers/Record.h>
33 : #include <casacore/casa/Arrays/ArrayLogical.h>
34 : #include <casacore/casa/Arrays/ArrayMath.h>
35 : #include <casacore/casa/Arrays/Array.h>
36 : #include <casacore/casa/Arrays/MaskedArray.h>
37 : #include <casacore/casa/Arrays/Vector.h>
38 : #include <casacore/casa/Arrays/Matrix.h>
39 : #include <casacore/casa/Arrays/Cube.h>
40 : #include <casacore/casa/BasicSL/String.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 : #include <casacore/casa/Exceptions/Error.h>
43 : #include <casacore/casa/OS/Timer.h>
44 : #include <sstream>
45 :
46 : #include <synthesis/TransformMachines/StokesImageUtil.h>
47 : #include <synthesis/TransformMachines2/VisModelData.h>
48 : #include <casacore/lattices/Lattices/LatticeLocker.h>
49 : #include <casacore/images/Images/ImageInterface.h>
50 : #include <casacore/images/Images/SubImage.h>
51 :
52 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
53 :
54 : #include <synthesis/TransformMachines2/MultiTermFTNew.h>
55 : #include <synthesis/TransformMachines2/MosaicFTNew.h>
56 : #include <synthesis/TransformMachines2/Utils.h>
57 :
58 : // This is the list of FTMachine types supported by MultiTermFTNew
59 : //#include <synthesis/TransformMachines/GridFT.h>
60 : //#include <synthesis/TransformMachines/AWProjectFT.h>
61 : //#include <synthesis/TransformMachines/AWProjectWBFT.h>
62 :
63 : #include<synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
64 :
65 : using namespace casacore;
66 : namespace casa { //# NAMESPACE CASA - BEGIN
67 : namespace refim { // namespace for refactor
68 :
69 : using namespace casacore;
70 : using namespace casa;
71 : using namespace casacore;
72 : using namespace casa::refim;
73 : using namespace casacore;
74 : using namespace casa::vi;
75 :
76 : #define PSOURCE false
77 : #define psource (IPosition(4,1536,1536,0,0))
78 :
79 : //----------------------------------------------------------------------
80 : //-------------------- constructors and descructors ----------------------
81 : //----------------------------------------------------------------------
82 238 : MultiTermFTNew::MultiTermFTNew(CountedPtr<FTMachine>&subftm, Int nterms, Bool forward)
83 238 : :FTMachine(), nterms_p(nterms),
84 238 : reffreq_p(0.0), imweights_p(Matrix<Float>(0,0)), machineName_p("MultiTermFTNew")
85 : // donePSF_p(false)
86 : {
87 :
88 238 : this->setBasePrivates(*subftm);
89 238 : canComputeResiduals_p = subftm->canComputeResiduals();
90 :
91 238 : if( forward ) psfnterms_p = nterms_p;
92 119 : else psfnterms_p = 2*nterms_p-1;
93 :
94 238 : subftms_p.resize(psfnterms_p);
95 824 : for(uInt termindex=0;termindex<psfnterms_p;termindex++)
96 : {
97 : // cout << "Creating new FTM of type : " << subftm->name() << endl;
98 586 : if( termindex==0 ){
99 238 : subftms_p[termindex] = subftm;
100 : }
101 : else {
102 348 : subftms_p[termindex] = getNewFTM(subftm);
103 348 : if((subftms_p[termindex]->name())=="MosaicFTNew"){
104 69 : (static_cast<MosaicFTNew *>(subftms_p[termindex].get()))->setConvFunc( (static_cast<MosaicFTNew * >(subftm.get()))->getConvFunc());
105 :
106 : }
107 : }
108 :
109 586 : subftms_p[termindex]->setMiscInfo(termindex);
110 : }
111 :
112 : // printFTTypes();
113 238 : }
114 :
115 : //----------------------------------------------------------------------
116 : // Construct from the input state record
117 0 : MultiTermFTNew::MultiTermFTNew(const RecordInterface& stateRec)
118 0 : : FTMachine()
119 : {
120 0 : String error;
121 0 : if (!fromRecord(error, stateRec)) {
122 0 : throw (AipsError("Failed to create gridder: " + error));
123 : };
124 0 : }
125 :
126 : //----------------------------------------------------------------------
127 : // Copy constructor
128 6 : MultiTermFTNew::MultiTermFTNew(const MultiTermFTNew& other) : FTMachine(), machineName_p("MultiTermFTNew")
129 : {
130 6 : operator=(other);
131 6 : }
132 :
133 6 : MultiTermFTNew& MultiTermFTNew::operator=(const MultiTermFTNew& other)
134 : {
135 :
136 6 : if(this!=&other)
137 : {
138 6 : FTMachine::operator=(other);
139 :
140 : // Copy local privates
141 6 : machineName_p = other.machineName_p;
142 6 : nterms_p = other.nterms_p;
143 6 : psfnterms_p = other.psfnterms_p;
144 6 : reffreq_p = other.reffreq_p;
145 : // donePSF_p = other.donePSF_p;
146 :
147 : // Make the list of subftms
148 6 : subftms_p.resize(other.subftms_p.nelements());
149 21 : for (uInt termindex=0;termindex<other.subftms_p.nelements();termindex++)
150 : {
151 15 : subftms_p[termindex] = getNewFTM( (other.subftms_p[termindex]) );
152 15 : subftms_p[termindex]->setMiscInfo(termindex);
153 : }
154 : // subftms_p[termindex] = getNewFTM( &(*(other.subftms_p[termindex])) );
155 :
156 : // Just checking....
157 6 : AlwaysAssert( subftms_p.nelements()>0 , AipsError );
158 :
159 : // Check if the sub ftm type can calculate its own residuals....
160 6 : canComputeResiduals_p = subftms_p[0]->canComputeResiduals();
161 :
162 : }
163 :
164 : // cout << "Checking FTtypes at the end of MTFT operator= for " << ( (other.subftms_p.nelements() > nterms_p)?String("grid"):String("degrid") ) << endl;
165 : // printFTTypes();
166 :
167 6 : return *this;
168 :
169 : }
170 :
171 363 : CountedPtr<FTMachine> MultiTermFTNew::getNewFTM(const CountedPtr<FTMachine>& ftm)
172 : {
173 :
174 363 : return ftm->cloneFTM();
175 :
176 : /*
177 : if(ftm->name()=="GridFT")
178 : {
179 : return new GridFT(static_cast<const GridFT&>(*ftm));
180 : }
181 : else if(ftm->name()=="AWProjectWBFT")
182 : { return new AWProjectWBFT(static_cast<const AWProjectWBFT&>(*ftm)); }
183 : else
184 : {throw(AipsError("FTMachine "+ftm->name()+" is not supported with MS-MFS")); }
185 :
186 : return NULL;
187 : */
188 :
189 : }
190 :
191 6 : FTMachine* MultiTermFTNew::cloneFTM()
192 : {
193 6 : return new MultiTermFTNew(*this);
194 : }
195 :
196 :
197 : //----------------------------------------------------------------------
198 488 : MultiTermFTNew::~MultiTermFTNew()
199 : {
200 488 : }
201 0 : void MultiTermFTNew::setMovingSource(const String& sourcename, const String& ephemtable){
202 0 : for (uInt k=0; k < subftms_p.nelements(); ++k)
203 0 : (subftms_p[k])->setMovingSource(sourcename, ephemtable);
204 0 : }
205 0 : void MultiTermFTNew::setMovingSource(const MDirection& mdir){
206 0 : for (uInt k=0; k < subftms_p.nelements(); ++k)
207 0 : (subftms_p[k])->setMovingSource(mdir);
208 0 : }
209 0 : void MultiTermFTNew::setLocation(const MPosition& mloc){
210 0 : for (uInt k=0; k < subftms_p.nelements(); ++k)
211 0 : (subftms_p[k])->setLocation(mloc);
212 0 : }
213 : //---------------------------------------------------------------------------------------------------
214 : //------------ Multi-Term Specific Functions --------------------------------------------------------
215 : //---------------------------------------------------------------------------------------------------
216 :
217 : // Multiply the imaging weights by Taylor functions - in place
218 : // This function MUST be called in ascending Taylor-term order
219 : // NOTE : Add checks to ensure this.
220 93495 : Bool MultiTermFTNew::modifyVisWeights(VisBuffer2& vb,uInt thisterm)
221 : {
222 : {
223 : ///There is no setImagingWeight in vb2 !
224 93495 : Matrix<Float>& imwgt=const_cast<Matrix<Float>& >(vb.imagingWeight());
225 : ////remove the above line when using setImagingWeight
226 :
227 93495 : if(imweights_p.shape() != vb.imagingWeight().shape())
228 113 : imweights_p.resize(vb.imagingWeight().shape());
229 93495 : imweights_p = vb.imagingWeight();
230 :
231 93495 : Float freq=0.0,mulfactor=1.0;
232 93495 : Vector<Double> selfreqlist(vb.getFrequencies(0));
233 :
234 27652446 : for (rownr_t row=0; row<vb.nRows(); row++)
235 170453710 : for (Int chn=0; chn<vb.nChannels(); chn++)
236 : {
237 142894759 : freq = selfreqlist(IPosition(1,chn));
238 142894759 : mulfactor = ((freq-reffreq_p)/reffreq_p);
239 142894759 : (imwgt)(chn,row) *= pow( mulfactor,(Int)thisterm);
240 : // sumwt_p += (vb.imagingWeight())(chn,row);
241 : }
242 93495 : vb.setImagingWeight(imwgt);
243 93495 : }
244 : /* // For debugging.
245 : else
246 : {
247 : for (Int row=0; row<vb.nRow(); row++)
248 : for (Int chn=0; chn<vb.nChannel(); chn++)
249 : {
250 : sumwt_p += (vb.imagingWeight())(chn,row);
251 : }
252 : }
253 : */
254 93495 : return true;
255 : }
256 :
257 0 : void MultiTermFTNew::initMaps(const VisBuffer2& vb){
258 :
259 0 : for (uInt k=0; k < subftms_p.nelements(); ++k){
260 0 : (subftms_p[k])->setFrameValidity( freqFrameValid_p);
261 0 : (subftms_p[k])->freqInterpMethod_p=this->freqInterpMethod_p;
262 0 : (subftms_p[k])->initMaps(vb);
263 : }
264 0 : }
265 : // Reset the imaging weights back to their original values
266 : // to be called just after "put"
267 93495 : void MultiTermFTNew::restoreImagingWeights(VisBuffer2 &vb)
268 : {
269 93495 : AlwaysAssert( imweights_p.shape() == vb.imagingWeight().shape() ,AipsError);
270 : ///There is no setImagingWeight in vb2 !
271 : // Matrix<Float>& imwgt=const_cast<Matrix<Float>& >(vb.imagingWeight());
272 : // ////remove the above line when using setImagingWeight
273 : // imwgt = imweights_p;
274 93495 : vb.setImagingWeight(imweights_p);
275 93495 : }
276 :
277 :
278 : // Multiply the model visibilities by the Taylor functions - in place.
279 22881 : Bool MultiTermFTNew::modifyModelVis(VisBuffer2& vb, uInt thisterm)
280 : {
281 22881 : Float freq=0.0,mulfactor=1.0;
282 22881 : Vector<Double> selfreqlist(vb.getFrequencies(0));
283 :
284 : // DComplex modcount=0.0;
285 :
286 22881 : Cube<Complex> modelVisCube(vb.visCubeModel().shape());
287 22881 : modelVisCube=vb.visCubeModel();
288 68787 : for (uInt pol=0; pol< uInt((vb.visCubeModel()).shape()[0]); pol++)
289 199660 : for (uInt chn=0; chn< uInt(vb.nChannels()); chn++)
290 48953648 : for (uInt row=0; row< uInt(vb.nRows()); row++)
291 : {
292 : // modcount += ( vb.modelVisCube())(pol,chn,row);
293 48799894 : freq = selfreqlist(IPosition(1,chn));
294 48799894 : mulfactor = ((freq-reffreq_p)/reffreq_p);
295 48799894 : modelVisCube(pol,chn,row) *= pow(mulfactor, (Int) thisterm);
296 : }
297 22881 : vb.setVisCubeModel(modelVisCube);
298 : // cout << "field : " << vb.fieldId() << " spw : "
299 : // << vb.spectralWindow() << " --- predicted model before taylor wt mult :"
300 : // << thisterm << " sumvis : " << modcount << endl;
301 :
302 22881 : return true;
303 22881 : }
304 :
305 :
306 : //---------------------------------------------------------------------------------------------------
307 : //---------------------- Prediction and De-gridding -----------------------------------
308 : //---------------------------------------------------------------------------------------------------
309 :
310 : // void MultiTermFTNew::initializeToVis(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec,PtrBlock<SubImage<Float> *> & modelImageVec, PtrBlock<SubImage<Float> *>& weightImageVec, PtrBlock<SubImage<Float> *>& fluxScaleVec,Block<Matrix<Float> >& weightsVec, const VisBuffer& vb)
311 :
312 79 : void MultiTermFTNew::initializeToVisNew(const VisBuffer2& vb,
313 : CountedPtr<SIImageStore> imstore)
314 : {
315 :
316 : // Convert Stokes planes to correlation planes..
317 235 : for(uInt taylor=0;taylor<nterms_p;taylor++)
318 : {
319 :
320 156 : if(!(imstore->forwardGrid(taylor)).get())
321 0 : throw(AipsError("MultiTermFTNew::InitializeToVisNew error imagestore has no valid grid initialized for taylor term "+String::toString(taylor)));
322 156 : stokesToCorrelation( *(imstore->model(taylor)) , *(imstore->forwardGrid(taylor) ) );
323 :
324 156 : if(vb.polarizationFrame()==MSIter::Linear) {
325 0 : StokesImageUtil::changeCStokesRep( *(imstore->forwardGrid(taylor) ), StokesImageUtil::LINEAR);
326 : } else {
327 156 : StokesImageUtil::changeCStokesRep( *(imstore->forwardGrid(taylor) ) , StokesImageUtil::CIRCULAR);
328 : }
329 : }
330 :
331 79 : reffreq_p = imstore->getReferenceFrequency();
332 :
333 235 : for(uInt taylor=0;taylor<nterms_p;taylor++)
334 : {
335 156 : subftms_p[taylor]->initializeToVis(*(imstore->forwardGrid(taylor)),vb);
336 : }
337 :
338 79 : }// end of initializeToVis
339 :
340 :
341 :
342 23841 : void MultiTermFTNew::get(VisBuffer2& vb, Int row)
343 : {
344 :
345 : // De-grid the model for the zeroth order Taylor term
346 23841 : subftms_p[0]->get(vb,row);
347 : // Save the model visibilities in a local cube
348 23841 : modviscube_p.assign( vb.visCubeModel() );
349 :
350 46722 : for(uInt tix=1;tix<nterms_p;tix++) // Only nterms.... not 2nterms-1
351 : {
352 : // Reset the model visibilities to zero
353 22881 : vb.setVisCubeModel(Complex(0.0,0.0));
354 : // De-grid the model onto the modelviscube (other Taylor terms)
355 22881 : subftms_p[tix]->get(vb,row);
356 : // Multiply visibilities by taylor-weights
357 22881 : modifyModelVis(vb,tix);
358 : // Accumulate model visibilities across Taylor terms
359 22881 : modviscube_p += vb.visCubeModel();
360 : }
361 : // Set the vb.modelviscube to what has been accumulated
362 23841 : vb.setVisCubeModel(modviscube_p);
363 23841 : }
364 :
365 79 : void MultiTermFTNew::finalizeToVis()
366 : {
367 79 : AlwaysAssert(subftms_p.nelements() >= nterms_p , AipsError);
368 235 : for(uInt taylor=0;taylor<nterms_p;taylor++) subftms_p[taylor]->finalizeToVis();
369 79 : }
370 :
371 : //---------------------------------------------------------------------------------------------------
372 : //---------------------- Calculate Residual Visibilities -------------------------------
373 : //---------------------------------------------------------------------------------------------------
374 0 : void MultiTermFTNew::ComputeResiduals(VisBuffer2 &vb, Bool useCorrected)
375 : {
376 :
377 0 : if(subftms_p[0]->canComputeResiduals()) subftms_p[0]->ComputeResiduals(vb,useCorrected);
378 0 : else throw(AipsError("MultiTerm::ComputeResiduals : subftm of MultiTermFTNew cannot compute its own residuals !"));
379 :
380 0 : }
381 :
382 : //---------------------------------------------------------------------------------------------------
383 : //---------------------- Gridding --------------------------------------------------------------
384 : //---------------------------------------------------------------------------------------------------
385 :
386 : /// void MultiTermFTNew::initializeToSky(Block<CountedPtr<ImageInterface<Complex> > > & compImageV
387 : //ec, Block<Matrix<Float> >& weightsVec, const VisBuffer& vb, const Bool dopsf)
388 324 : Long MultiTermFTNew::estimateRAM(const CountedPtr<SIImageStore>& imstor){
389 324 : Long mem=0;
390 1178 : for(uInt k=0; k < subftms_p.nelements(); ++k){
391 :
392 854 : mem+=subftms_p[k]->estimateRAM(imstor);
393 : }
394 :
395 324 : return mem;
396 : }
397 :
398 289 : void MultiTermFTNew::initializeToSkyNew(const Bool dopsf,
399 : const VisBuffer2& vb,
400 : CountedPtr<SIImageStore> imstore)
401 : {
402 :
403 : // If PSF is already done, don't ask again !
404 : // AlwaysAssert( !(donePSF_p && dopsf) , AipsError );
405 :
406 : // The PSF needs to be the first thing made (because of weight images)
407 : // AlwaysAssert( !(dopsf==false && donePSF_p==false) , AipsError);
408 :
409 : // if(donePSF_p==true)
410 289 : if(dopsf==false)
411 : {
412 181 : if( subftms_p.nelements() != nterms_p )
413 : {
414 111 : subftms_p.resize( nterms_p ,true);
415 : // cout << "MTFT::initializeToSky : resizing to " << nterms_p << " terms" << endl;
416 : }
417 : }
418 :
419 : // Make the relevant float grid.
420 : // This is needed mainly for facetting (to set facet shapes), but is harmless for non-facetting.
421 289 : if( dopsf ) { imstore->psf(0); } else { imstore->residual(0); }
422 :
423 289 : reffreq_p = imstore->getReferenceFrequency();
424 :
425 289 : Matrix<Float> sumWeight;
426 964 : for(uInt taylor=0;taylor< (dopsf ? psfnterms_p : nterms_p);taylor++)
427 : {
428 :
429 675 : if(! (imstore->backwardGrid(taylor)).get())
430 0 : throw(AipsError("MultiTermFTNew::InitializeToSkyNew error imagestore has no valid grid initialized for taylor term "+String::toString(taylor)));
431 675 : subftms_p[taylor]->initializeToSky(*(imstore->backwardGrid(taylor) ), sumWeight,vb);
432 : }
433 :
434 289 : }// end of initializeToSky
435 :
436 289 : void MultiTermFTNew::initBriggsWeightor(vi::VisibilityIterator2& vi){
437 :
438 :
439 964 : for (uInt k=0; k < subftms_p.nelements(); ++k){
440 675 : subftms_p[k]->initBriggsWeightor(vi);
441 : }
442 :
443 289 : }
444 :
445 72877 : void MultiTermFTNew::put(VisBuffer2& vb, Int row, Bool dopsf, refim::FTMachine::Type type)
446 : {
447 :
448 72877 : subftms_p[0]->put(vb,row,dopsf,type);
449 72877 : if (!dryRun())
450 : {
451 72877 : Int gridnterms=nterms_p;
452 72877 : if(dopsf==true) // && donePSF_p==false)
453 : {
454 25808 : gridnterms=2*nterms_p-1;
455 : }
456 :
457 : //cerr << " Calling put for " << gridnterms << " terms, nelements : " << subftms_p.nelements() << " and dopsf " << dopsf << " reffreq " << reffreq_p << endl;
458 :
459 166372 : for(Int tix=1;tix<gridnterms;tix++)
460 : {
461 93495 : modifyVisWeights(vb,tix);
462 93495 : subftms_p[tix]->put(vb,row,dopsf,type);
463 93495 : restoreImagingWeights(vb);
464 : }
465 : }
466 :
467 72877 : }// end of put
468 :
469 : //----------------------------------------------------------------------
470 :
471 : // void MultiTermFTNew::finalizeToSky(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec, PtrBlock<SubImage<Float> *> & resImageVec, PtrBlock<SubImage<Float> *>& weightImageVec, PtrBlock<SubImage<Float> *>& fluxScaleVec, Bool dopsf, Block<Matrix<Float> >& weightsVec, const VisBuffer& /*vb*/)
472 :
473 289 : void MultiTermFTNew::finalizeToSkyNew(Bool dopsf,
474 : const VisBuffer2& /*vb*/,
475 : CountedPtr<SIImageStore> imstore )
476 : {
477 :
478 : // Collect images and weights from all FTMs
479 964 : for(uInt taylor=0;taylor< (dopsf ? psfnterms_p : nterms_p) ;taylor++)
480 : {
481 675 : Matrix<Float> sumWeights;
482 675 : subftms_p[taylor]->finalizeToSky();
483 675 : shared_ptr<ImageInterface<Float> > theim=dopsf ? imstore->psf(taylor) : imstore->residual(taylor);
484 675 : { LatticeLocker lock1 (*theim, FileLocker::Write);
485 675 : correlationToStokes( subftms_p[taylor]->getImage(sumWeights, false) , *theim, dopsf);
486 675 : }
487 675 : if( subftms_p[taylor]->useWeightImage() && dopsf ) {
488 118 : LatticeLocker lock1 (*(imstore->weight(taylor)), FileLocker::Write);
489 118 : subftms_p[taylor]->getWeightImage(*(imstore->weight(taylor)), sumWeights);
490 118 : }
491 :
492 : // Take sumWeights from corrToStokes here....
493 1350 : Matrix<Float> sumWeightStokes( (imstore->sumwt())->shape()[2], (imstore->sumwt())->shape()[3] );
494 675 : StokesImageUtil::ToStokesSumWt( sumWeightStokes, sumWeights );
495 :
496 675 : AlwaysAssert( ( (imstore->sumwt(taylor))->shape()[2] == sumWeightStokes.shape()[0] ) &&
497 : ((imstore->sumwt(taylor))->shape()[3] == sumWeightStokes.shape()[1] ) , AipsError );
498 675 : LatticeLocker lock1 (*(imstore->sumwt(taylor)), FileLocker::Write);
499 675 : (imstore->sumwt(taylor))->put( sumWeightStokes.reform((imstore->sumwt(taylor))->shape()) );
500 :
501 : // cout << "taylor : " << taylor << " sumwt : " << sumWeights << endl;
502 :
503 675 : }// end for taylor
504 :
505 : // if( dopsf ) donePSF_p = true;
506 :
507 289 : }//end of finalizeToSkyNew
508 :
509 :
510 : //---------------------------------------------------------------------------------------------------
511 : //----------------------------- Obtain Images -----------------------------------------------------
512 : //---------------------------------------------------------------------------------------------------
513 : //----------------------------------------------------------------------
514 0 : void MultiTermFTNew::makeImage(refim::FTMachine::Type type, VisibilityIterator2& vi,
515 : ImageInterface<Complex>& theImage, Matrix<Float>& weight)
516 : {
517 : // cout << "MTFT :: makeImage for taylor 0 only "<< endl;
518 0 : subftms_p[0]->makeImage(type, vi, theImage, weight);
519 0 : }
520 0 : void MultiTermFTNew::makeMTImages(refim::FTMachine::Type type,
521 : vi::VisibilityIterator2& vi,
522 : casacore::Vector<casacore::CountedPtr<casacore::ImageInterface<casacore::Complex> > >& theImage,
523 : casacore::Vector<casacore::CountedPtr<casacore::Matrix<casacore::Float> > >& weight){
524 0 : Int ntaylor= (type== refim::FTMachine::PSF) ? psfnterms_p : nterms_p;
525 :
526 :
527 :
528 0 : vi::VisBuffer2* vb=vi.getVisBuffer();
529 :
530 : // Initialize put (i.e. transform to Sky) for this model
531 0 : vi.origin();
532 0 : for(Int taylor=0;taylor< ntaylor ; ++taylor) {
533 0 : if(vb->polarizationFrame()==MSIter::Linear) {
534 0 : StokesImageUtil::changeCStokesRep(*(theImage[taylor]), StokesImageUtil::LINEAR);
535 : }
536 : else {
537 0 : StokesImageUtil::changeCStokesRep(*(theImage[taylor]), StokesImageUtil::CIRCULAR);
538 : }
539 0 : subftms_p[taylor]->initializeToSky(*(theImage[taylor]), *(weight[taylor]),*vb);
540 :
541 : }
542 : {
543 0 : Vector<Double> refpix = (theImage[0]->coordinates().spectralCoordinate()).referencePixel();
544 0 : (theImage[0]->coordinates().spectralCoordinate()).toWorld( reffreq_p, refpix[0] );
545 0 : }
546 0 : Bool useCorrected= !(MSColumns(vi.ms()).correctedData().isNull());
547 0 : if((type==FTMachine::CORRECTED) && (!useCorrected))
548 0 : type=FTMachine::OBSERVED;
549 :
550 :
551 : // Loop over the visibilities, putting VisBuffers
552 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
553 0 : for (vi.origin(); vi.more(); vi.next()) {
554 :
555 0 : switch(type) {
556 0 : case FTMachine::RESIDUAL:
557 0 : vb->setVisCube(vb->visCubeCorrected());
558 0 : vb->setVisCube(vb->visCube()-vb->visCubeModel());
559 0 : put(*vb, -1, false);
560 0 : break;
561 0 : case FTMachine::CORRECTED:
562 0 : put(*vb, -1, false, FTMachine::CORRECTED);
563 0 : break;
564 0 : case FTMachine::PSF:
565 0 : vb->setVisCube(Complex(1.0,0.0));
566 0 : put(*vb, -1, true, FTMachine::PSF);
567 0 : break;
568 0 : case FTMachine::OBSERVED:
569 0 : put(*vb, -1, false, FTMachine::OBSERVED);
570 0 : break;
571 0 : default:
572 0 : throw(AipsError("Cannot make multiterm image of type requested"));
573 : break;
574 : }
575 : }
576 : }
577 : ///
578 0 : for(Int taylor=0;taylor< ntaylor ; ++taylor) {
579 0 : subftms_p[taylor]->finalizeToSky();
580 0 : subftms_p[taylor]->getImage(*(weight[taylor]), false);
581 :
582 : }
583 0 : }
584 : //---------------------------------------------------------------------------------------------------
585 : //------------------------ To / From Records ---------------------------------------------------------
586 : //---------------------------------------------------------------------------------------------------
587 2 : Bool MultiTermFTNew::toRecord(String& error, RecordInterface& outRec, Bool withImage, const String diskimage)
588 : {
589 : // cout << "MTFTNew :: toRecord for " << subftms_p.nelements() << " subftms" << endl;
590 2 : Bool retval = true;
591 2 : outRec.define("name", this->name());
592 2 : outRec.define("nterms",nterms_p);
593 2 : outRec.define("reffreq",reffreq_p);
594 2 : outRec.define("machinename",machineName_p);
595 2 : outRec.define("psfnterms",psfnterms_p);
596 : // outRec.define("donePSF_p",donePSF_p);
597 :
598 2 : outRec.define("numfts", (Int)subftms_p.nelements() ); // Since the forward and reverse ones are different.
599 :
600 6 : for(uInt tix=0;tix<subftms_p.nelements();tix++)
601 : {
602 4 : Record subFTContainer;
603 4 : String elimage="";
604 4 : if(diskimage != ""){
605 4 : elimage=diskimage+String("_term_")+String::toString(tix);
606 : }
607 4 : subftms_p[tix]->toRecord(error, subFTContainer,withImage, elimage);
608 4 : outRec.defineRecord("subftm_"+String::toString(tix),subFTContainer);
609 4 : }
610 :
611 2 : return retval;
612 : }
613 :
614 : //---------------------------------------------------------------------------------------------------
615 0 : Bool MultiTermFTNew::fromRecord(String& error, const RecordInterface& inRec)
616 : {
617 : //cout << "MTFTNew :: fromRecord "<< endl;
618 0 : Bool retval = true;
619 :
620 0 : inRec.get("nterms",nterms_p);
621 0 : inRec.get("reffreq",reffreq_p);
622 0 : inRec.get("machinename",machineName_p);
623 0 : inRec.get("psfnterms",psfnterms_p);
624 : // inRec.get("donePSF_p",donePSF_p);
625 :
626 0 : Int nftms=1;
627 0 : inRec.get("numfts",nftms);
628 :
629 0 : subftms_p.resize(nftms);
630 : //cerr << "number of ft " << nftms << endl;
631 0 : for(Int tix=0;tix<nftms;tix++)
632 : {
633 0 : Record subFTMRec=inRec.asRecord("subftm_"+String::toString(tix));
634 0 : subftms_p[tix]=VisModelData::NEW_FT(subFTMRec);
635 0 : if(subftms_p[tix].null())
636 0 : throw(AipsError("Could not recover from record term "+String::toString(tix)+" ftmachine"));
637 0 : retval = (retval || subftms_p[tix]->fromRecord(error, subFTMRec));
638 0 : if(!retval) throw(AipsError("Could not recover term "+String::toString(tix)+" ftmachine; \n Error being "+error));
639 0 : }
640 :
641 :
642 0 : return retval;
643 : }
644 : //---------------------------------------------------------------------------------------------------
645 :
646 0 : Bool MultiTermFTNew::storeAsImg(String fileName, ImageInterface<Float>& theImg)
647 : {
648 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), fileName);
649 0 : LatticeExpr<Float> le(theImg);
650 0 : tmp.copyData(le);
651 0 : return true;
652 0 : }
653 : //
654 : // Set the supplied CFCache for all internal FTMs if they do use CFCache mechanism
655 : //
656 0 : void MultiTermFTNew::setCFCache(casacore::CountedPtr<CFCache>& cfc, const casacore::Bool resetCFC)
657 : {
658 0 : for (unsigned int i=0;i<subftms_p.nelements();i++)
659 : {
660 0 : if (subftms_p[i]->isUsingCFCache())
661 0 : subftms_p[i]->setCFCache(cfc,resetCFC);
662 : }
663 0 : }
664 :
665 :
666 : //---------------------------------------------------------------------------------------------------
667 : //---------------------------------------------------------------------------------------------------
668 : //---------------------------------------------------------------------------------------------------
669 : } // end namespace refim
670 : } //# NAMESPACE CASA - END
671 :
|
