Line data Source code
1 : //# SIImageStoreMultiTerm.cc: Implementation of Imager.h
2 : //# Copyright (C) 1997-2008
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program 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 General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@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 <casacore/casa/Exceptions/Error.h>
29 : #include <iostream>
30 : #include <sstream>
31 :
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 :
36 : #include <casacore/casa/Logging.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : #include <casacore/casa/Logging/LogMessage.h>
39 : #include <casacore/casa/Logging/LogSink.h>
40 : #include <casacore/casa/Logging/LogMessage.h>
41 :
42 : #include <casacore/casa/OS/DirectoryIterator.h>
43 : #include <casacore/casa/OS/File.h>
44 : #include <casacore/casa/OS/Path.h>
45 :
46 : #include <casacore/casa/OS/HostInfo.h>
47 : #include <casacore/images/Images/TempImage.h>
48 : #include <casacore/images/Images/PagedImage.h>
49 : #include <casacore/ms/MeasurementSets/MSHistoryHandler.h>
50 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
51 : #include <synthesis/TransformMachines/StokesImageUtil.h>
52 : #include <casacore/images/Images/TempImage.h>
53 : #include <casacore/images/Images/SubImage.h>
54 : #include <casacore/images/Regions/ImageRegion.h>
55 :
56 : #include <synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
57 :
58 : #include <casacore/casa/Arrays/MatrixMath.h>
59 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
60 :
61 : #include <sys/types.h>
62 : #include <unistd.h>
63 : using namespace std;
64 :
65 : using namespace casacore;
66 : namespace casa { //# NAMESPACE CASA - BEGIN
67 :
68 : //////////////////////////////////////////////////////////////////////////////////////////////////////
69 : //////////////////////////////////////////////////////////////////////////////////////////////////////
70 :
71 0 : SIImageStoreMultiTerm::SIImageStoreMultiTerm():SIImageStore()
72 : {
73 0 : itsNTerms=0;
74 :
75 0 : itsPsfs.resize(0);
76 0 : itsModels.resize(0);
77 0 : itsResiduals.resize(0);
78 0 : itsWeights.resize(0);
79 0 : itsImages.resize(0);
80 0 : itsSumWts.resize(0);
81 0 : itsImagePBcors.resize(0);
82 0 : itsPBs.resize(0);
83 :
84 0 : itsForwardGrids.resize(0);
85 0 : itsBackwardGrids.resize(0);
86 :
87 0 : itsUseWeight=false;
88 :
89 0 : init();
90 :
91 0 : validate();
92 :
93 0 : }
94 :
95 : // Used from SynthesisNormalizer::makeImageStore()
96 0 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(const String &imagename,
97 : const CoordinateSystem &imcoordsys,
98 : const IPosition &imshape,
99 : const String &objectname,
100 : const Record &miscinfo,
101 : const Int /*nfacets*/,
102 : const Bool /*overwrite*/,
103 : uInt ntaylorterms,
104 0 : Bool useweightimage)
105 : {
106 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","Open new Images",WHERE) );
107 :
108 0 : itsNTerms = ntaylorterms;
109 :
110 0 : itsPsfs.resize(2 * itsNTerms - 1);
111 0 : itsModels.resize(itsNTerms);
112 0 : itsResiduals.resize(itsNTerms);
113 0 : itsWeights.resize(2 * itsNTerms - 1);
114 0 : itsImages.resize(itsNTerms);
115 0 : itsSumWts.resize(2 * itsNTerms - 1);
116 0 : itsPBs.resize(itsNTerms);
117 0 : itsImagePBcors.resize(itsNTerms);
118 :
119 0 : itsMask.reset( );
120 0 : itsGridWt.reset( );
121 :
122 0 : itsForwardGrids.resize(itsNTerms);
123 0 : itsBackwardGrids.resize(2 * itsNTerms - 1);
124 :
125 : // itsNFacets = nfacets; // So that sumwt shape happens properly, via checkValidity
126 : // itsFacetId = -1;
127 0 : itsNFacets=1;
128 0 : itsFacetId=0;
129 0 : itsNChanChunks = 1;
130 0 : itsChanId = 0;
131 0 : itsNPolChunks = 1;
132 0 : itsPolId = 0;
133 :
134 0 : itsImageName = imagename;
135 0 : itsCoordSys = imcoordsys;
136 0 : itsImageShape = imshape;
137 0 : itsObjectName = objectname;
138 0 : itsMiscInfo = miscinfo;
139 :
140 0 : itsUseWeight = useweightimage;
141 :
142 0 : init();
143 :
144 0 : validate();
145 :
146 0 : }
147 :
148 : // Used from SynthesisNormalizer::makeImageStore()
149 0 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(const String &imagename, uInt ntaylorterms,
150 0 : const Bool ignorefacets, const Bool ignoresumwt)
151 : {
152 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","Open existing Images",WHERE) );
153 :
154 0 : itsNTerms = ntaylorterms;
155 :
156 0 : auto fltPtrReset = [](Block<shared_ptr<ImageInterface<Float> > >&a) {for(uInt i=0; i < a.nelements(); ++i) a[i].reset(); };
157 0 : itsPsfs.resize(2 * itsNTerms - 1); fltPtrReset(itsPsfs);
158 0 : itsModels.resize(itsNTerms); fltPtrReset(itsModels);
159 0 : itsResiduals.resize(itsNTerms); fltPtrReset(itsResiduals);
160 0 : itsWeights.resize(2 * itsNTerms - 1); fltPtrReset(itsWeights);
161 0 : itsImages.resize(itsNTerms); fltPtrReset(itsImages);
162 0 : itsPBs.resize(itsNTerms); fltPtrReset(itsPBs);
163 0 : itsSumWts.resize(2 * itsNTerms - 1); fltPtrReset(itsSumWts);
164 0 : itsMask.reset( );
165 0 : itsGridWt.reset( );
166 0 : itsImagePBcors.resize(itsNTerms); fltPtrReset(itsImagePBcors);
167 :
168 :
169 :
170 0 : itsMiscInfo=Record();
171 0 : auto cmplxPtrReset = [](Block<shared_ptr<ImageInterface<Complex> > >&a) {for(uInt i=0; i < a.nelements(); ++i) a[i].reset(); };
172 0 : itsForwardGrids.resize(itsNTerms); cmplxPtrReset(itsForwardGrids);
173 0 : itsBackwardGrids.resize(2 * itsNTerms - 1); cmplxPtrReset(itsBackwardGrids);
174 :
175 0 : itsImageName = imagename;
176 :
177 0 : itsNFacets=1;
178 0 : itsFacetId=0;
179 0 : itsNChanChunks = 1;
180 0 : itsChanId = 0;
181 0 : itsNPolChunks = 1;
182 0 : itsPolId = 0;
183 :
184 0 : Bool exists=true;
185 0 : Bool sumwtexists=true;
186 0 : Bool modelexists=true;
187 0 : for(uInt tix=0;tix<2*itsNTerms-1;tix++)
188 : {
189 0 : if( tix<itsNTerms ) {
190 0 : exists &= ( doesImageExist( itsImageName+String(".residual.tt")+String::toString(tix) ) ||
191 0 : doesImageExist( itsImageName+String(".psf.tt")+String::toString(tix) ) );
192 0 : modelexists &= ( doesImageExist( itsImageName+String(".model.tt")+String::toString(tix) ) ||
193 0 : doesImageExist( itsImageName+String(".model.tt")+String::toString(tix) ) );
194 0 : sumwtexists &= ( doesImageExist( itsImageName+String(".sumwt.tt")+String::toString(tix) ) );
195 : }else {
196 0 : exists &= ( doesImageExist( itsImageName+String(".psf.tt")+String::toString(tix) ) );
197 0 : sumwtexists &= ( doesImageExist( itsImageName+String(".sumwt.tt")+String::toString(tix) ) );
198 : }
199 : }
200 :
201 : // The PSF or Residual images must exist. ( or the gridwt image)
202 : // All this is just for the shape and coordinate system.
203 0 : if( exists || modelexists || doesImageExist(itsImageName+String(".gridwt")) )
204 : {
205 0 : std::shared_ptr<ImageInterface<Float> > imptr;
206 0 : if( doesImageExist(itsImageName+String(".psf.tt0")) )
207 : {
208 0 : buildImage( imptr, (itsImageName+String(".psf.tt0")) );
209 :
210 : //cout << "Opening PSF image to read csys" << endl;
211 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".psf.tt0")) );
212 : }
213 0 : else if( doesImageExist(itsImageName+String(".residual.tt0")) )
214 : {
215 0 : buildImage( imptr, (itsImageName+String(".residual.tt0")) );
216 : //cout << "Opening Residual image to read csys" << endl;
217 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".residual.tt0")) );
218 : }
219 0 : else if( doesImageExist(itsImageName+String(".model.tt0")) )
220 : {
221 0 : buildImage( imptr, (itsImageName+String(".model.tt0")) );
222 : //cout << "Opening Model image to read csys" << endl;
223 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".model.tt0")) );
224 : }
225 : else
226 : {
227 : // How can this be right ?
228 : //cout << "Opening Sumwt image to read csys" << endl;
229 0 : buildImage( imptr, (itsImageName+String(".gridwt")) );
230 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".gridwt")) );
231 : }
232 :
233 0 : itsObjectName=imptr->imageInfo().objectName();
234 0 : itsImageShape = imptr->shape();
235 0 : itsCoordSys = imptr->coordinates();
236 0 : itsMiscInfo=imptr->miscInfo();
237 :
238 : }
239 : else
240 : {
241 0 : throw( AipsError( "Multi-term PSF, Residual or Model Images do not exist. Please create one of them." ) );
242 : }
243 :
244 0 : if( doesImageExist(itsImageName+String(".residual.tt0")) ||
245 0 : doesImageExist(itsImageName+String(".psf.tt0")) )
246 : {
247 0 : if( sumwtexists )
248 : {
249 0 : std::shared_ptr<ImageInterface<Float> > imptr;
250 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt.tt0")) );
251 0 : buildImage( imptr, (itsImageName+String(".sumwt.tt0")) );
252 0 : itsNFacets = imptr->shape()[0];
253 0 : itsFacetId = 0;
254 0 : itsUseWeight = getUseWeightImage( *imptr );
255 : /////redo this here as psf may have different coordinates
256 0 : itsCoordSys = imptr->coordinates();
257 0 : itsMiscInfo=imptr->miscInfo();
258 0 : if( itsUseWeight && ! doesImageExist(itsImageName+String(".weight.tt0")) )
259 : {
260 0 : throw(AipsError("Internal error : MultiTerm Sumwt has a useweightimage=true but the weight image does not exist."));
261 : }
262 : }
263 : else
264 : {
265 0 : if(!ignoresumwt)
266 0 : {throw( AipsError( "Multi-term SumWt does not exist. Please create PSFs or Residuals." ) );}
267 : else
268 : {
269 0 : os << "SumWt.ttx do not exist. Proceeding only with PSFs" << LogIO::POST;
270 0 : std::shared_ptr<ImageInterface<Float> > imptr;
271 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt.tt0")) );
272 0 : if( doesImageExist(itsImageName+String(".residual.tt0")) )
273 0 : {buildImage( imptr, (itsImageName+String(".residual.tt0")) );}
274 : else
275 0 : {buildImage( imptr, (itsImageName+String(".psf.tt0")) );}
276 :
277 0 : itsNFacets = 1;
278 0 : itsFacetId = 0;
279 0 : itsUseWeight = False;
280 0 : itsCoordSys = imptr->coordinates();
281 0 : itsMiscInfo=imptr->miscInfo();
282 : }
283 : }
284 : }// if psf0 or res0 exist
285 :
286 0 : if( ignorefacets==true ) itsNFacets=1;
287 :
288 0 : init();
289 0 : validate();
290 :
291 0 : }
292 :
293 : /*
294 : /////////////Constructor with pointers already created else where but taken over here
295 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(Block<std::shared_ptr<ImageInterface<Float> > > modelims,
296 : Block<std::shared_ptr<ImageInterface<Float> > >residims,
297 : Block<std::shared_ptr<ImageInterface<Float> > >psfims,
298 : Block<std::shared_ptr<ImageInterface<Float> > >weightims,
299 : Block<std::shared_ptr<ImageInterface<Float> > >restoredims,
300 : Block<std::shared_ptr<ImageInterface<Float> > >sumwtims,
301 : std::shared_ptr<ImageInterface<Float> > newmask,
302 : std::shared_ptr<ImageInterface<Float> > newalpha,
303 : std::shared_ptr<ImageInterface<Float> > newbeta)
304 : {
305 :
306 : itsPsfs=psfims;
307 : itsModels=modelims;
308 : itsResiduals=residims;
309 : itsWeights=weightims;
310 : itsImages=restoredims;
311 : itsSumWts=sumwtims;
312 : itsMask = newmask;
313 : itsAlpha = newalpha;
314 : itsBeta = newbeta;
315 :
316 : itsNTerms = itsResiduals.nelements();
317 : itsMiscInfo=Record();
318 :
319 : AlwaysAssert( itsPsfs.nelements() == 2*itsNTerms-1 , AipsError );
320 : AlwaysAssert( itsPsfs.nelements()>0 && itsPsfs[0] , AipsError );
321 : AlwaysAssert( itsSumWts.nelements()>0 && itsSumWts[0] , AipsError );
322 :
323 : itsForwardGrids.resize( itsNTerms );
324 : itsBackwardGrids.resize( 2 * itsNTerms - 1 );
325 :
326 : itsImageShape=psfims[0]->shape();
327 : itsCoordSys = psfims[0]->coordinates();
328 : itsMiscInfo = psfims[0]->miscInfo();
329 :
330 : itsNFacets=sumwtims[0]->shape()[0];
331 : itsUseWeight=getUseWeightImage( *(sumwtims[0]) );
332 :
333 : itsImageName = String("reference"); // This is what the access functions use to guard against allocs...
334 :
335 : init();
336 : validate();
337 :
338 : }
339 : */
340 :
341 : //////////////////////////////////////////////////////////////////////////////////////////////////////
342 : /////////////Constructor with pointers already created else where but taken over here
343 : // used from getSubImageStore(), for example when creating the facets list
344 : // this would be safer if it was refactored as a copy constructor of the generic stuff +
345 : // initialization of the facet related parameters
346 0 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(const Block<std::shared_ptr<ImageInterface<Float> > > &modelims,
347 : const Block<std::shared_ptr<ImageInterface<Float> > > &residims,
348 : const Block<std::shared_ptr<ImageInterface<Float> > > &psfims,
349 : const Block<std::shared_ptr<ImageInterface<Float> > > &weightims,
350 : const Block<std::shared_ptr<ImageInterface<Float> > > &restoredims,
351 : const Block<std::shared_ptr<ImageInterface<Float> > > &sumwtims,
352 : const Block<std::shared_ptr<ImageInterface<Float> > > &pbims,
353 : const Block<std::shared_ptr<ImageInterface<Float> > > &restoredpbcorims,
354 : const std::shared_ptr<ImageInterface<Float> > &newmask,
355 : const std::shared_ptr<ImageInterface<Float> > &newalpha,
356 : const std::shared_ptr<ImageInterface<Float> > &newbeta,
357 : const std::shared_ptr<ImageInterface<Float> > &newalphaerror,
358 : const std::shared_ptr<ImageInterface<Float> > &newalphapbcor,
359 : const std::shared_ptr<ImageInterface<Float> > &newbetapbcor,
360 : const CoordinateSystem& csys,
361 : const IPosition &imshape,
362 : const String &imagename,
363 : const String &objectname,
364 : const Record &miscinfo,
365 : const Int facet, const Int nfacets,
366 : const Int chan, const Int nchanchunks,
367 0 : const Int pol, const Int npolchunks)
368 : {
369 0 : itsPsfs=psfims;
370 0 : itsModels=modelims;
371 0 : itsResiduals=residims;
372 0 : itsWeights=weightims;
373 0 : itsImages=restoredims;
374 0 : itsSumWts=sumwtims;
375 0 : itsMask = newmask;
376 0 : itsAlpha = newalpha;
377 0 : itsBeta = newbeta;
378 0 : itsAlphaError = newalphaerror;
379 0 : itsAlphaPBcor = newalphapbcor;
380 0 : itsBetaPBcor = newbetapbcor;
381 :
382 0 : itsPBs=pbims;
383 0 : itsImagePBcors=restoredpbcorims;
384 :
385 0 : itsNTerms = itsResiduals.nelements();
386 0 : itsMiscInfo=Record();
387 :
388 0 : AlwaysAssert( itsPsfs.nelements() == 2*itsNTerms-1 , AipsError );
389 : // AlwaysAssert( itsPsfs.nelements()>0 && itsPsfs[0] , AipsError );
390 : // AlwaysAssert( itsSumWts.nelements()>0 && itsSumWts[0] , AipsError );
391 :
392 0 : itsForwardGrids.resize( itsNTerms );
393 0 : itsBackwardGrids.resize( 2 * itsNTerms - 1 );
394 :
395 0 : itsObjectName = objectname;
396 0 : itsMiscInfo = miscinfo;
397 :
398 0 : itsNFacets = nfacets;
399 0 : itsFacetId = facet;
400 0 : itsNChanChunks = nchanchunks;
401 0 : itsChanId = chan;
402 0 : itsNPolChunks = npolchunks;
403 0 : itsPolId = pol;
404 :
405 0 : itsParentImageShape = imshape;
406 0 : itsImageShape = imshape;
407 : ///// itsImageShape = IPosition(4,0,0,0,0);
408 :
409 0 : itsCoordSys = csys; // Hopefully this doesn't change for a reference image
410 0 : itsImageName = imagename;
411 :
412 :
413 : //-----------------------
414 0 : init(); // Connect parent pointers to the images.
415 : //-----------------------
416 :
417 : // Set these to null, to be set later upon first access.
418 : // Setting to null will hopefully set all elements of each array, to NULL.
419 0 : itsPsfs=std::shared_ptr<ImageInterface<Float> >();
420 0 : itsModels=std::shared_ptr<ImageInterface<Float> >();
421 0 : itsResiduals=std::shared_ptr<ImageInterface<Float> >();
422 0 : itsWeights=std::shared_ptr<ImageInterface<Float> >();
423 0 : itsImages=std::shared_ptr<ImageInterface<Float> >();
424 0 : itsSumWts=std::shared_ptr<ImageInterface<Float> >();
425 0 : itsPBs=std::shared_ptr<ImageInterface<Float> >();
426 :
427 0 : itsMask.reset( );
428 :
429 0 : validate();
430 :
431 0 : }
432 :
433 : //////////////////////////////////////////////////////////////////////////////////////////////////////
434 :
435 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
436 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
437 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
438 :
439 0 : uInt SIImageStoreMultiTerm::getNTaylorTerms(Bool dopsf)
440 : {
441 0 : return dopsf ? (2*itsNTerms-1) : itsNTerms;
442 : }
443 :
444 : // Check if images that are asked-for are ready and all have the same shape.
445 : /*
446 : Bool SIImageStoreMultiTerm::checkValidity(const Bool ipsf, const Bool iresidual,
447 : const Bool iweight, const Bool imodel, const Bool irestored,
448 : const Bool imask,const Bool isumwt,
449 : const Bool ialpha, const Bool ibeta)
450 : {
451 :
452 : // cout << "In MT::checkValidity imask is " << imask << endl;
453 :
454 : Bool valid = true;
455 :
456 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
457 : {
458 :
459 : if(ipsf==true)
460 : { psf(tix);
461 : valid = valid & ( itsPsfs[tix] && itsPsfs[tix]->shape()==itsImageShape ); }
462 : if(iweight==true)
463 : { weight(tix);
464 : valid = valid & ( itsWeights[tix] && itsWeights[tix]->shape()==itsImageShape ); }
465 :
466 : if(isumwt==true) {
467 : IPosition useShape(itsImageShape);
468 : useShape[0]=itsNFacets; useShape[1]=itsNFacets;
469 : sumwt(tix);
470 : valid = valid & ( itsSumWts[tix] && itsSumWts[tix]->shape()==useShape );
471 : }
472 :
473 : if( tix< itsNTerms )
474 : {
475 : if(iresidual==true)
476 : { residual(tix);
477 : valid = valid & ( itsResiduals[tix] && itsResiduals[tix]->shape()==itsImageShape ); }
478 : if(imodel==true)
479 : { model(tix);
480 : valid = valid & ( itsModels[tix] && itsModels[tix]->shape()==itsImageShape); }
481 : if(irestored==true)
482 : { image(tix);
483 : valid = valid & ( itsImages[tix] && itsImages[tix]->shape()==itsImageShape); }
484 : }
485 : }
486 :
487 : if(imask==true)
488 : { mask(); valid = valid & ( itsMask && itsMask->shape()==itsImageShape);
489 : // cout << " Mask null ? " << (bool) itsMask << endl;
490 : }
491 : if(ialpha==true)
492 : { alpha(); valid = valid & ( itsAlpha && itsAlpha->shape()==itsImageShape ); }
493 : if(ibeta==true)
494 : { beta(); valid = valid & ( itsBeta && itsBeta->shape()==itsImageShape ); }
495 :
496 : return valid;
497 :
498 : }
499 : */
500 :
501 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
502 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
503 : //////////////////////////////////////////////////////////////////////////////////////////////////////
504 : //////////////////////////////////////////////////////////////////////////////////////////////////////
505 :
506 0 : SIImageStoreMultiTerm::~SIImageStoreMultiTerm()
507 : {
508 0 : }
509 :
510 : //////////////////////////////////////////////////////////////////////////////////////////////////////
511 : //////////////////////////////////////////////////////////////////////////////////////////////////////
512 :
513 0 : Bool SIImageStoreMultiTerm::releaseLocks()
514 : {
515 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","releaseLocks",WHERE) );
516 :
517 0 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
518 : {
519 0 : if( itsPsfs[tix] ) releaseImage( itsPsfs[tix] );
520 0 : if( itsWeights[tix] ) releaseImage( itsWeights[tix] );
521 0 : if( itsSumWts[tix] ) releaseImage( itsSumWts[tix] );
522 0 : if( tix < itsNTerms )
523 : {
524 0 : if( itsModels[tix] ) releaseImage( itsModels[tix] );
525 0 : if( itsResiduals[tix] ) releaseImage( itsResiduals[tix] );
526 0 : if( itsImages[tix] ) releaseImage( itsImages[tix] );
527 0 : if( itsPBs[tix] ) releaseImage( itsPBs[tix] );
528 0 : if( itsImagePBcors[tix] ) releaseImage( itsImagePBcors[tix] );
529 : }
530 : }
531 0 : if( itsMask ) releaseImage( itsMask );
532 0 : if( itsAlpha ) releaseImage( itsAlpha );
533 0 : if( itsBeta ) releaseImage( itsBeta );
534 0 : if( itsAlphaError ) releaseImage( itsAlphaError );
535 0 : if( itsAlphaPBcor ) releaseImage( itsAlphaPBcor );
536 0 : if( itsBetaPBcor ) releaseImage( itsBetaPBcor );
537 0 : if( itsGridWt ) releaseImage( itsGridWt );
538 :
539 0 : return true; // do something more intelligent here.
540 : }
541 :
542 0 : Bool SIImageStoreMultiTerm::releaseComplexGrids()
543 : {
544 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","releaseComplexGrids",WHERE) );
545 :
546 0 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
547 : {
548 0 : if( itsBackwardGrids[tix] ) releaseImage( itsBackwardGrids[tix] );
549 0 : if( tix < itsNTerms )
550 : {
551 0 : if( itsForwardGrids[tix] ) releaseImage( itsForwardGrids[tix] );
552 : }
553 : }
554 0 : return True; // do something more intelligent here.
555 : }
556 :
557 :
558 0 : Double SIImageStoreMultiTerm::getReferenceFrequency()
559 : {
560 : Double theRefFreq;
561 :
562 0 : Vector<Double> refpix = (itsCoordSys.spectralCoordinate()).referencePixel();
563 0 : AlwaysAssert( refpix.nelements()>0, AipsError );
564 0 : (itsCoordSys.spectralCoordinate()).toWorld( theRefFreq, refpix[0] );
565 : // cout << "Reading ref freq as : " << theRefFreq << endl;
566 0 : return theRefFreq;
567 : }
568 :
569 : //////////////////////////////////////////////////////////////////////////////////////////////////////
570 : //////////////////////////////////////////////////////////////////////////////////////////////////////
571 0 : Vector<String> SIImageStoreMultiTerm::getModelImageName()
572 : {
573 0 : Vector<String> mods(itsNTerms);
574 0 : for(uInt tix=0;tix<itsNTerms;tix++)
575 0 : {mods[tix]=itsImageName + imageExts(MODEL)+".tt"+String::toString(tix);}
576 0 : return mods;
577 : }
578 :
579 0 : void SIImageStoreMultiTerm::setModelImage( const Vector<String> &modelnames )
580 : {
581 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","setModelImage",WHERE) );
582 :
583 0 : if( modelnames.nelements() > itsNTerms )
584 0 : { throw(AipsError("We currently cannot support more than nterms images as the starting model. "));
585 : }
586 :
587 0 : if( modelnames.nelements() > 0 && modelnames.nelements() <= itsNTerms )
588 : {
589 0 : for(uInt tix=0;tix<modelnames.nelements();tix++)
590 : {
591 0 : setModelImageOne( modelnames[tix], tix );
592 : }
593 : }
594 :
595 0 : }
596 :
597 : /*
598 : void SIImageStoreMultiTerm::setModelImage( String modelname )
599 : {
600 : LogIO os( LogOrigin("SIImageStoreMultiTerm","setModelImage",WHERE) );
601 :
602 : for(uInt tix=0;tix<itsNTerms;tix++)
603 : {
604 :
605 : Directory immodel( modelname+String(".model.tt")+String::toString(tix) );
606 : if( !immodel.exists() )
607 : {
608 : os << "Starting model image does not exist for term : " << tix << LogIO::POST;
609 : }
610 : else
611 : {
612 : std::shared_ptr<PagedImage<Float> > newmodel( new PagedImage<Float>( modelname+String(".model.tt")+String::toString(tix) ) );
613 : // Check shapes, coordsys with those of other images. If different, try to re-grid here.
614 :
615 : if( newmodel->shape() != model(tix)->shape() )
616 : {
617 : os << "Regridding input model to target coordinate system for term " << tix << LogIO::POST;
618 : regridToModelImage( *newmodel , tix);
619 : // For now, throw an exception.
620 : //throw( AipsError( "Input model image "+modelname+".model.tt"+String::toString(tix)+" is not the same shape as that defined for output in "+ itsImageName + ".model" ) );
621 : }
622 : else
623 : {
624 : os << "Setting " << modelname << " as model for term " << tix << LogIO::POST;
625 : // Then, add its contents to itsModel.
626 : //itsModel->put( itsModel->get() + model->get() );
627 : /////( model(tix) )->put( newmodel->get() );
628 : model(tix)->copyData( LatticeExpr<Float> (*newmodel) );
629 : }
630 : }
631 : }//nterms
632 : }
633 : */
634 :
635 : //////////////////////////////////////////////////////////////////////////////////////////////////////
636 : //////////////////////////////////////////////////////////////////////////////////////////////////////
637 :
638 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::psf(uInt term)
639 : {
640 0 : AlwaysAssert( itsPsfs.nelements() > term, AipsError );
641 0 : accessImage( itsPsfs[term], itsParentPsfs[term], imageExts(PSF)+".tt"+String::toString(term) );
642 0 : return itsPsfs[term];
643 : }
644 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::residual(uInt term)
645 : {
646 0 : accessImage( itsResiduals[term], itsParentResiduals[term], imageExts(RESIDUAL)+".tt"+String::toString(term) );
647 : // Record mi = itsResiduals[term]->miscInfo(); ostringstream oss;mi.print(oss);cout<<"MiscInfo(res) " << term << " : " << oss.str() << endl;
648 :
649 0 : return itsResiduals[term];
650 : }
651 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::weight(uInt term)
652 : {
653 0 : accessImage( itsWeights[term], itsParentWeights[term], imageExts(WEIGHT)+".tt"+String::toString(term) );
654 0 : return itsWeights[term];
655 : }
656 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::sumwt(uInt term)
657 : {
658 0 : accessImage( itsSumWts[term], itsParentSumWts[term], imageExts(SUMWT)+".tt"+String::toString(term) );
659 :
660 :
661 0 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
662 0 : {itsUseWeight = getUseWeightImage( *itsParentSumWts[0] );}
663 0 : setUseWeightImage( *(itsSumWts[term]) , itsUseWeight); // Sets a flag in the SumWt image.
664 :
665 0 : return itsSumWts[term];
666 : }
667 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::model(uInt term)
668 : {
669 :
670 0 : accessImage( itsModels[term], itsParentModels[term], imageExts(MODEL)+".tt"+String::toString(term) );
671 :
672 0 : itsModels[term]->setUnits("Jy/pixel");
673 0 : return itsModels[term];
674 : }
675 :
676 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::image(uInt term)
677 : {
678 :
679 0 : accessImage( itsImages[term], itsParentImages[term], imageExts(IMAGE)+".tt"+String::toString(term) );
680 0 : itsImages[term]->setUnits("Jy/beam");
681 0 : return itsImages[term];
682 : }
683 :
684 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::pb(uInt term)
685 : {
686 :
687 0 : accessImage( itsPBs[term], itsParentPBs[term], imageExts(PB)+".tt"+String::toString(term) );
688 0 : return itsPBs[term];
689 : }
690 :
691 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::imagepbcor(uInt term)
692 : {
693 :
694 0 : accessImage( itsImagePBcors[term], itsParentImagePBcors[term], imageExts(IMAGE)+".tt"+String::toString(term)+ ".pbcor" );
695 0 : itsImagePBcors[term]->setUnits("Jy/beam");
696 0 : return itsImagePBcors[term];
697 : }
698 :
699 0 : std::shared_ptr<ImageInterface<Complex> > SIImageStoreMultiTerm::forwardGrid(uInt term){
700 0 : if( itsForwardGrids[term] )// && (itsForwardGrids[term]->shape() == itsImageShape))
701 0 : return itsForwardGrids[term];
702 0 : Vector<Int> whichStokes(0);
703 0 : IPosition cimageShape;
704 0 : cimageShape=itsImageShape;
705 0 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
706 0 : whichStokes, itsDataPolRep);
707 0 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
708 0 : cimageShape(2)=whichStokes.nelements();
709 0 : itsForwardGrids[term].reset(new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()));
710 0 : return itsForwardGrids[term];
711 : }
712 0 : std::shared_ptr<ImageInterface<Complex> > SIImageStoreMultiTerm::backwardGrid(uInt term){
713 :
714 0 : if( itsBackwardGrids[term] && (itsBackwardGrids[term]->shape() == itsImageShape))
715 0 : return itsBackwardGrids[term];
716 : // cout << "MT : Making backward grid of shape : " << itsImageShape << endl;
717 0 : Vector<Int> whichStokes(0);
718 0 : IPosition cimageShape;
719 0 : cimageShape=itsImageShape;
720 0 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
721 0 : whichStokes, itsDataPolRep);
722 0 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
723 0 : cimageShape(2)=whichStokes.nelements();
724 0 : itsBackwardGrids[term].reset(new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()));
725 0 : return itsBackwardGrids[term];
726 : }
727 :
728 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::alpha()
729 : {
730 0 : if( itsAlpha && itsAlpha->shape() == itsImageShape ) { return itsAlpha; }
731 : // checkRef( itsAlpha , "alpha" );
732 0 : itsAlpha = openImage( itsImageName+String(".alpha"), false );
733 : // itsAlpha->setUnits("Alpha");
734 0 : return itsAlpha;
735 : }
736 :
737 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::beta()
738 : {
739 0 : if( itsBeta && itsBeta->shape() == itsImageShape ) { return itsBeta; }
740 : // checkRef( itsBeta , "beta" );
741 0 : itsBeta = openImage( itsImageName+String(".beta"), false );
742 : // itsBeta->setUnits("Beta");
743 0 : return itsBeta;
744 : }
745 :
746 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::alphaerror()
747 : {
748 0 : if( itsAlphaError && itsAlphaError->shape() == itsImageShape ) { return itsAlphaError; }
749 : // checkRef( itsAlpha , "alpha" );
750 0 : itsAlphaError = openImage( itsImageName+String(".alpha.error"), false );
751 : // itsAlpha->setUnits("Alpha");
752 0 : return itsAlphaError;
753 : }
754 :
755 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::alphapbcor()
756 : {
757 0 : if( itsAlphaPBcor && itsAlphaPBcor->shape() == itsImageShape ) { return itsAlphaPBcor; }
758 : // checkRef( itsAlphaPBcor , "alpha" );
759 0 : itsAlphaPBcor = openImage( itsImageName+String(".alpha.pbcor"), False );
760 : // itsAlphaPBcor->setUnits("Alpha");
761 0 : return itsAlphaPBcor;
762 : }
763 :
764 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::betapbcor()
765 : {
766 0 : if( itsBetaPBcor && itsBetaPBcor->shape() == itsImageShape ) { return itsBetaPBcor; }
767 : // checkRef( itsBetaPBcor , "beta" );
768 0 : itsBetaPBcor = openImage( itsImageName+String(".beta.pbcor"), False );
769 : // itsBetaPBcor->setUnits("Beta");
770 0 : return itsBetaPBcor;
771 : }
772 :
773 :
774 :
775 : // TODO : Move to an image-wrapper class ? Same function exists in SynthesisDeconvolver.
776 0 : Bool SIImageStoreMultiTerm::doesImageExist(String imagename)
777 : {
778 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","doesImageExist",WHERE) );
779 0 : Directory image( imagename );
780 0 : return image.exists();
781 : }
782 :
783 :
784 0 : void SIImageStoreMultiTerm::resetImages( Bool resetpsf, Bool resetresidual, Bool resetweight )
785 : {
786 0 : for(uInt tix=0;tix<2*itsNTerms-1;tix++)
787 : {
788 0 : if( resetpsf ) psf(tix)->set(0.0);
789 :
790 0 : if( tix < itsNTerms ) {
791 0 : if( resetresidual ) {
792 : //removeMask( residual(tix) );
793 0 : residual(tix)->set(0.0);
794 : }
795 : }
796 0 : if( resetweight && itsWeights[tix] ) weight(tix)->set(0.0);
797 0 : if( resetweight ) sumwt(tix)->set(0.0);
798 : }//nterms
799 0 : }
800 :
801 0 : void SIImageStoreMultiTerm::addImages( std::shared_ptr<SIImageStore> imagestoadd,
802 : Bool addpsf, Bool addresidual, Bool addweight, Bool adddensity)
803 : {
804 0 : for(uInt tix=0;tix<2*itsNTerms-1;tix++)
805 : {
806 :
807 0 : if(addpsf)
808 : {
809 0 : LatticeExpr<Float> adderPsf( *(psf(tix)) + *(imagestoadd->psf(tix)) );
810 0 : psf(tix)->copyData(adderPsf);
811 : }
812 0 : if(addweight)
813 : {
814 :
815 0 : if(getUseWeightImage( *(imagestoadd->sumwt(tix)) ) ) // Access and add weight only if it is needed.
816 : {
817 0 : LatticeExpr<Float> adderWeight( *(weight(tix)) + *(imagestoadd->weight(tix)) );
818 0 : weight(tix)->copyData(adderWeight);
819 : }
820 :
821 0 : LatticeExpr<Float> adderSumWt( *(sumwt(tix)) + *(imagestoadd->sumwt(tix)) );
822 0 : sumwt(tix)->copyData(adderSumWt);
823 0 : setUseWeightImage( *sumwt(tix), getUseWeightImage(*(imagestoadd->sumwt(tix)) ) );
824 : }
825 :
826 0 : if(tix < itsNTerms && addresidual)
827 : {
828 0 : LatticeExpr<Float> adderRes( *(residual(tix)) + *(imagestoadd->residual(tix)) );
829 0 : residual(tix)->copyData(adderRes);
830 : }
831 :
832 0 : if( tix==0 && adddensity )
833 : {
834 0 : LatticeExpr<Float> adderDensity( *(gridwt()) + *(imagestoadd->gridwt()) );
835 0 : gridwt()->copyData(adderDensity);
836 : }
837 :
838 : }
839 0 : }
840 :
841 0 : void SIImageStoreMultiTerm::dividePSFByWeight(const Float /*pblimit*/)
842 : {
843 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","dividePSFByWeight",WHERE) );
844 :
845 : //// for(uInt tix=0;tix<2*itsNTerms-1;tix++)
846 0 : for(Int tix=2*itsNTerms-1-1;tix>-1;tix--) // AAH go backwards so that zeroth term is normalized last..... sigh sigh sigh.
847 : {
848 :
849 : //cout << "npsfs : " << itsPsfs.nelements() << " and tix : " << tix << endl;
850 :
851 0 : normPSF(tix);
852 :
853 0 : if ( itsUseWeight) {
854 0 : divideImageByWeightVal( *weight(tix) );
855 : }
856 :
857 : }
858 0 : }
859 :
860 0 : void SIImageStoreMultiTerm::normalizePrimaryBeam(const Float pblimit)
861 : {
862 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","normalizePrimaryBeam",WHERE) );
863 0 : if ( itsUseWeight) {
864 :
865 0 : makePBFromWeight(pblimit);
866 : }
867 0 : else { makePBImage(pblimit); }
868 : // calcSensitivity();
869 0 : }
870 :
871 :
872 : // Make another for the PSF too.
873 0 : void SIImageStoreMultiTerm::divideResidualByWeight(Float pblimit, const String normtype)
874 : {
875 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","divideResidualByWeight",WHERE) );
876 :
877 0 : if( itsUseWeight )
878 : {
879 0 : itsPBScaleFactor = getPbMax();
880 : }
881 :
882 0 : for(uInt tix=0;tix<itsNTerms;tix++)
883 : {
884 :
885 0 : divideImageByWeightVal( *residual(tix) );
886 :
887 : // if(doesImageExist(itsImageName+String(".weight.tt0")) )
888 0 : if( itsUseWeight )
889 : {
890 0 : LatticeExpr<Float> ratio;
891 0 : Float scalepb = fabs(pblimit) * itsPBScaleFactor * itsPBScaleFactor ;
892 0 : if( normtype=="flatnoise"){
893 0 : LatticeExpr<Float> deno = LatticeExpr<Float> ( sqrt( abs(*(weight(0)) ) ) * itsPBScaleFactor );
894 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".residual.tt")+String::toString(tix) ;
895 0 : os << " by [ sqrt(weightimage) * " << itsPBScaleFactor ;
896 0 : os << " ] to get flat noise with unit pb peak."<< LogIO::POST;
897 0 : LatticeExpr<Float> mask( iif( (deno) > scalepb , 1.0, 0.0 ) );
898 0 : LatticeExpr<Float> maskinv( iif( (deno) > scalepb , 0.0, 1.0 ) );
899 0 : ratio=LatticeExpr<Float> ( ( (*(residual(tix))) * mask ) / ( deno + maskinv ) );
900 : }
901 0 : else if(normtype=="pbsquare"){
902 0 : Float deno = itsPBScaleFactor*itsPBScaleFactor ;
903 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".residual.tt")+String::toString(tix) ;
904 0 : os << deno ;
905 0 : os << " ] to get optimal sig/noise with unit pb peak."<< LogIO::POST;
906 :
907 0 : ratio=LatticeExpr<Float> ( ( *(residual(tix)) ) / ( deno ) );
908 :
909 :
910 :
911 : }
912 0 : else if( normtype=="flatsky") {
913 0 : LatticeExpr<Float> deno( *(weight(0)) );
914 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".residual.tt")+String::toString(tix) ;
915 0 : os << " by [ weight ] to get flat sky"<< LogIO::POST;
916 0 : LatticeExpr<Float> mask( iif( (deno) > scalepb , 1.0, 0.0 ) );
917 0 : LatticeExpr<Float> maskinv( iif( (deno) > scalepb , 0.0, 1.0 ) );
918 0 : ratio=LatticeExpr<Float> ( ( (*(residual(tix))) * mask ) / ( deno + maskinv ) );
919 : }
920 : else{
921 0 : throw(AipsError("Don't know how to proceed with normtype "+normtype));
922 : }
923 :
924 :
925 :
926 0 : residual(tix)->copyData(ratio);
927 : }
928 :
929 0 : if( (residual(tix)->getDefaultMask()=="") && hasPB() && pblimit >= 0.0 )
930 0 : {copyMask(pb(),residual(tix));}
931 :
932 0 : if( pblimit <0.0 && (residual(tix)->getDefaultMask()).matches("mask0") ) removeMask( residual(tix) );
933 :
934 : }
935 0 : }
936 :
937 0 : void SIImageStoreMultiTerm::divideModelByWeight(Float pblimit, const String normtype)
938 : {
939 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","divideModelByWeight",WHERE) );
940 :
941 0 : if( itsUseWeight // only when needed
942 0 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
943 : {
944 :
945 0 : if( normtype=="flatsky") {
946 0 : Array<Float> arrmod;
947 0 : model(0)->get( arrmod, true );
948 :
949 0 : os << "Model is already flat sky with peak flux : " << max(arrmod);
950 0 : os << ". No need to divide before prediction" << LogIO::POST;
951 :
952 0 : return;
953 : }
954 :
955 0 : itsPBScaleFactor = getPbMax();
956 :
957 0 : for(uInt tix=0;tix<itsNTerms;tix++)
958 0 : { LatticeExpr<Float> deno;
959 0 : if(normtype=="flatnoise"){
960 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".model")+String::toString(tix);
961 0 : os << " by [ sqrt(weight) / " << itsPBScaleFactor ;
962 0 : os <<" ] to get to flat sky model before prediction" << LogIO::POST;
963 :
964 0 : deno = LatticeExpr<Float> ( sqrt( abs(*(weight(0))) ) / itsPBScaleFactor );
965 : }
966 0 : else if(normtype=="pbsquare"){
967 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".model")+String::toString(tix);
968 0 : os << " by [ (weight) / " << itsPBScaleFactor*itsPBScaleFactor ;
969 0 : os <<" ] to get an optimal sig/noise model before prediction" << LogIO::POST;
970 :
971 0 : deno = LatticeExpr<Float> ( abs(*(weight(0))) / (itsPBScaleFactor*itsPBScaleFactor) );
972 : }
973 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
974 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
975 0 : LatticeExpr<Float> ratio( ( (*(model(tix))) * mask ) / ( deno + maskinv ) );
976 :
977 0 : itsModels[tix]->copyData(ratio);
978 : }
979 : }
980 : }
981 :
982 :
983 0 : void SIImageStoreMultiTerm::multiplyModelByWeight(Float pblimit, const String normtype)
984 : {
985 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","multiplyModelByWeight",WHERE) );
986 :
987 :
988 0 : if( itsUseWeight // only when needed
989 0 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
990 : {
991 :
992 0 : if( normtype=="flatsky") {
993 0 : os << "Model is already flat sky. No need to multiply back after prediction" << LogIO::POST;
994 0 : return;
995 : }
996 :
997 0 : itsPBScaleFactor = getPbMax();
998 :
999 0 : for(uInt tix=0;tix<itsNTerms;tix++)
1000 : {
1001 0 : LatticeExpr<Float> deno;
1002 0 : if( normtype=="flatnoise") {
1003 0 : os << LogIO::NORMAL1 << "Multiplying " << itsImageName+String(".model")+String::toString(tix);
1004 0 : os << " by [ sqrt(weight) / " << itsPBScaleFactor;
1005 0 : os << " ] to take model back to flat noise with unit pb peak." << LogIO::POST;
1006 :
1007 0 : deno=LatticeExpr<Float> ( sqrt( abs(*(weight(0)) ) ) / itsPBScaleFactor );
1008 : }
1009 0 : else if ( normtype=="pbsquare"){
1010 0 : os << LogIO::NORMAL1 << "Multiplying " << itsImageName+String(".model")+String::toString(tix);
1011 0 : os << " by [ weight / " << itsPBScaleFactor*itsPBScaleFactor;
1012 0 : os << " ] to take model back to optima sig/noise with unit pb peak." << LogIO::POST;
1013 :
1014 0 : deno=LatticeExpr<Float> ( abs(*(weight(0)) ) / (itsPBScaleFactor*itsPBScaleFactor) );
1015 : }
1016 : else{
1017 0 : throw(AipsError("No idea of what to do for "+normtype));
1018 : }
1019 :
1020 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1021 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1022 0 : LatticeExpr<Float> ratio( ( (*(model(tix))) * mask ) * ( deno + maskinv ) );
1023 :
1024 0 : itsModels[tix]->copyData(ratio);
1025 : }
1026 : }
1027 : }
1028 :
1029 :
1030 0 : void SIImageStoreMultiTerm::restore(GaussianBeam& rbeam, String& usebeam, uInt /*term*/, Float psfcutoff)
1031 : {
1032 :
1033 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","restore",WHERE) );
1034 :
1035 0 : for(uInt tix=0; tix<itsNTerms; tix++)
1036 : {
1037 0 : SIImageStore::restore(rbeam, usebeam, tix, psfcutoff);
1038 : }
1039 :
1040 0 : calculateAlphaBeta("image");
1041 :
1042 0 : }// restore
1043 :
1044 0 : void SIImageStoreMultiTerm::calculateAlphaBeta(String imtype)
1045 : {
1046 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","calculateAlphaBeta",WHERE) );
1047 : try
1048 : {
1049 :
1050 : // Check if this is Stokes I.
1051 0 : Bool isOnlyStokesI=True;
1052 :
1053 0 : Vector<Int> stokes ( (itsParentCoordSys.stokesCoordinate()).stokes() );
1054 0 : AlwaysAssert( stokes.nelements()>0 , AipsError);
1055 0 : if( stokes.nelements()>1 || stokes[0]!=1 ) isOnlyStokesI=False;
1056 :
1057 0 : if( ! isOnlyStokesI )
1058 0 : { os << "Alpha and Beta images will not be computed for images that contain anything other than stokes I in them." << LogIO::POST; }
1059 :
1060 : //Calculate alpha, beta
1061 0 : if( itsNTerms > 1 && isOnlyStokesI)
1062 : {
1063 : // Calculate alpha and beta
1064 0 : LatticeExprNode leMaxRes = max( *( residual(0) ) );
1065 0 : Float maxres = leMaxRes.getFloat();
1066 0 : Float specthreshold = maxres/10.0; //////////// do something better here.....
1067 :
1068 0 : os << "Calculating spectral parameters for Intensity > peakresidual/10 = " << specthreshold << " Jy/beam" << LogIO::POST;
1069 :
1070 : /////////////////////////////////////////////////////////
1071 : /////// Calculate alpha
1072 0 : if(imtype=="image")
1073 : {
1074 0 : LatticeExpr<Float> mask1(iif(((*(image(0))))>(specthreshold),1.0,0.0));
1075 0 : LatticeExpr<Float> mask0(iif(((*(image(0))))>(specthreshold),0.0,1.0));
1076 0 : LatticeExpr<Float> alphacalc( (((*(image(1))))*mask1)/(((*(image(0))))+(mask0)) );
1077 0 : alpha()->copyData(alphacalc);
1078 :
1079 0 : ImageInfo ii = image(0)->imageInfo();
1080 : // Set the restoring beam for alpha
1081 0 : alpha()->setImageInfo(ii);
1082 : //alpha()->table().unmarkForDelete();
1083 :
1084 : // Make a mask for the alpha image
1085 0 : LatticeExpr<Bool> lemask(iif(((*(image(0))) > specthreshold) , True, False));
1086 0 : removeMask( alpha() );
1087 0 : createMask( lemask, (alpha()) );
1088 :
1089 : /////// Calculate alpha error
1090 0 : alphaerror()->set(0.0);
1091 :
1092 0 : LatticeExpr<Float> alphacalcerror( abs(alphacalc) * sqrt( ( (*residual(0)*mask1)/(*image(0)+mask0) )*( (*residual(0)*mask1)/(*image(0)+mask0) ) + ( (*residual(1)*mask1)/(*image(1)+mask0) )*( (*residual(1)*mask1)/(*image(1)+mask0) ) ) );
1093 0 : alphaerror()->copyData(alphacalcerror);
1094 0 : alphaerror()->setImageInfo(ii);
1095 0 : removeMask(alphaerror());
1096 0 : createMask(lemask, alphaerror());
1097 : // alphaerror()->table().unmarkForDelete();
1098 0 : os << "Written Spectral Index Error Image : " << alphaerror()->name() << LogIO::POST;
1099 :
1100 0 : if(itsNTerms>2) // calculate beta too.
1101 : {
1102 0 : beta()->set(0.0);
1103 0 : LatticeExpr<Float> betacalc( (*image(2)*mask1)/((*image(0))+(mask0))-0.5*(*alpha())*((*alpha())-1.0) );
1104 0 : beta()->copyData(betacalc);
1105 0 : beta()->setImageInfo(ii);
1106 : //imbeta.setUnits(Unit("Spectral Curvature"));
1107 0 : removeMask(beta());
1108 0 : createMask(lemask, beta());
1109 0 : os << "Written Spectral Curvature Image : " << beta()->name() << LogIO::POST;
1110 : }
1111 :
1112 : }
1113 : else
1114 : {
1115 0 : LatticeExpr<Float> mask1(iif(((*(imagepbcor(0))))>(specthreshold),1.0,0.0));
1116 0 : LatticeExpr<Float> mask0(iif(((*(imagepbcor(0))))>(specthreshold),0.0,1.0));
1117 0 : LatticeExpr<Float> alphacalc( (((*(imagepbcor(1))))*mask1)/(((*(imagepbcor(0))))+(mask0)) );
1118 0 : alphapbcor()->copyData(alphacalc);
1119 :
1120 0 : ImageInfo ii = image(0)->imageInfo();
1121 : // Set the restoring beam for alpha
1122 0 : alphapbcor()->setImageInfo(ii);
1123 : //alpha()->table().unmarkForDelete();
1124 :
1125 : // Make a mask for the alpha image
1126 0 : LatticeExpr<Bool> lemask(iif(((*(imagepbcor(0))) > specthreshold) , True, False));
1127 0 : removeMask( alphapbcor() );
1128 0 : createMask( lemask, (alphapbcor()) );
1129 :
1130 : /////////////////////////////////////////////////////////
1131 0 : if(itsNTerms>2) // calculate beta too.
1132 : {
1133 0 : betapbcor()->set(0.0);
1134 0 : LatticeExpr<Float> betacalc( (*imagepbcor(2)*mask1)/((*imagepbcor(0))+(mask0))-0.5*(*alphapbcor())*((*alphapbcor())-1.0) );
1135 0 : betapbcor()->copyData(betacalc);
1136 0 : betapbcor()->setImageInfo(ii);
1137 : //imbeta.setUnits(Unit("Spectral Curvature"));
1138 0 : removeMask(betapbcor());
1139 0 : createMask(lemask, betapbcor());
1140 0 : os << "Written Spectral Curvature Image : " << betapbcor()->name() << LogIO::POST;
1141 : }
1142 :
1143 : }// pbcor
1144 :
1145 : }//if nterms>1
1146 : }
1147 0 : catch(AipsError &x)
1148 : {
1149 0 : throw( AipsError("Error in computing Alpha (and Beta) : " + x.getMesg() ) );
1150 : }
1151 :
1152 0 : }// calculateAlphaBeta
1153 :
1154 0 : void SIImageStoreMultiTerm::pbcor()
1155 : {
1156 :
1157 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","pbcorPlane",WHERE) );
1158 :
1159 : /// Temp Code to prevent this approximate PBCOR from happening for EVLA data
1160 : if(1)
1161 : {
1162 0 : String telescope = itsCoordSys.obsInfo().telescope();
1163 0 : if ( telescope != "ALMA" )
1164 : {
1165 0 : os << LogIO::WARN << "Wideband (multi-term) PB correction is not yet available via tclean in the 4.7 release. Please use the widebandpbcor task instead. "<< LogIO::POST;
1166 0 : return;
1167 : }
1168 : else
1169 : {
1170 0 : os << LogIO::WARN << "Wideband (multi-term) PB Correction is currently only an approximation. It assumes no PB frequency dependence. This code has been added for the 4.7 release to support the current ALMA pipeline, which does not apply corrections for the frequency dependence of the primary beam across small fractional bandwidths. Please look at the help for the 'pbcor' parameter and use the widebandpbcor task if needed. " <<LogIO::POST;
1171 : }
1172 :
1173 :
1174 : }
1175 :
1176 :
1177 : // message saying that it's only stokes I for now...
1178 :
1179 0 : for(uInt tix=0; tix<itsNTerms; tix++)
1180 : {
1181 0 : SIImageStore::pbcor(tix);
1182 : }
1183 :
1184 0 : calculateAlphaBeta("pbcor");
1185 :
1186 : }
1187 :
1188 0 : void SIImageStoreMultiTerm::calcSensitivity()
1189 : {
1190 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","calcSensitivity",WHERE) );
1191 :
1192 : // Construct Hessian.
1193 :
1194 0 : Matrix<Float> hess(IPosition(2,itsNTerms,itsNTerms));
1195 0 : for(uInt tay1=0;tay1<itsNTerms;tay1++)
1196 0 : for(uInt tay2=0;tay2<itsNTerms;tay2++)
1197 : {
1198 : //uInt taymin = (tay1<=tay2)? tay1 : tay2;
1199 : //uInt taymax = (tay1>=tay2)? tay1 : tay2;
1200 : //uInt ind = (taymax*(taymax+1)/2)+taymin;
1201 :
1202 0 : uInt ind = tay1+tay2;
1203 0 : AlwaysAssert( ind < 2*itsNTerms-1, AipsError );
1204 :
1205 0 : Array<Float> lsumwt;
1206 0 : sumwt( ind )->get( lsumwt, false );
1207 : // cout << "lsumwt shape : " << lsumwt.shape() << endl;
1208 0 : AlwaysAssert( lsumwt.shape().nelements()==4, AipsError );
1209 0 : AlwaysAssert( lsumwt.shape()[0]>0, AipsError );
1210 :
1211 : // hess(tay1,tay2) = lsumwt(IPosition(1,0)); //Always pick sumwt from first facet only.
1212 0 : hess(tay1,tay2) = lsumwt(IPosition(4,0,0,0,0)); //Always pick sumwt from first facet only.
1213 : }
1214 :
1215 0 : os << "Multi-Term Hessian Matrix : " << hess << LogIO::POST;
1216 :
1217 : // Invert Hessian.
1218 : try
1219 : {
1220 0 : Float deter=0.0;
1221 0 : Matrix<Float> invhess;
1222 0 : invertSymPosDef(invhess,deter,hess);
1223 0 : os << "Multi-Term Covariance Matrix : " << invhess << LogIO::POST;
1224 :
1225 : // Just print the sqrt of the diagonal elements.
1226 :
1227 0 : for(uInt tix=0;tix<itsNTerms;tix++)
1228 : {
1229 0 : os << "[" << itsImageName << "][Taylor"<< tix << "] Theoretical sensitivity (Jy/bm):" ;
1230 0 : if( invhess(tix,tix) > 0.0 ) { os << sqrt(invhess(tix,tix)) << LogIO::POST; }
1231 0 : else { os << "none" << LogIO::POST; }
1232 : }
1233 : }
1234 0 : catch(AipsError &x)
1235 : {
1236 0 : os << LogIO::WARN << "Cannot invert Hessian Matrix : " << x.getMesg() << " || Calculating approximate sensitivity " << LogIO::POST;
1237 :
1238 : // Approximate : 1/h^2
1239 0 : for(uInt tix=0;tix<itsNTerms;tix++)
1240 : {
1241 0 : Array<Float> lsumwt;
1242 0 : AlwaysAssert( 2*tix < 2*itsNTerms-1, AipsError );
1243 0 : sumwt(2*tix)->get( lsumwt , false );
1244 :
1245 0 : IPosition shp( lsumwt.shape() );
1246 : //cout << "Sumwt shape : " << shp << " : " << lsumwt << endl;
1247 : //AlwaysAssert( shp.nelements()==4 , AipsError );
1248 :
1249 0 : os << "[" << itsImageName << "][Taylor"<< tix << "] Approx Theoretical sensitivity (Jy/bm):" ;
1250 :
1251 0 : IPosition it(4,0,0,0,0);
1252 0 : for ( it[0]=0; it[0]<shp[0]; it[0]++)
1253 0 : for ( it[1]=0; it[1]<shp[1]; it[1]++)
1254 0 : for ( it[2]=0; it[2]<shp[2]; it[2]++)
1255 0 : for ( it[3]=0; it[3]<shp[3]; it[3]++)
1256 : {
1257 0 : if( shp[0]>1 ){os << "f"<< it[0]+(it[1]*shp[0]) << ":" ;}
1258 0 : if( lsumwt( it ) > 1e-07 )
1259 : {
1260 0 : os << 1.0/(sqrt(lsumwt(it))) << " " ;
1261 : }
1262 : else
1263 : {
1264 0 : os << "none" << " ";
1265 : }
1266 : }
1267 :
1268 0 : os << LogIO::POST;
1269 : }
1270 : }
1271 :
1272 0 : calcFractionalBandwidth();
1273 0 : }
1274 :
1275 0 : Double SIImageStoreMultiTerm::calcFractionalBandwidth()
1276 : {
1277 :
1278 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","calcFractionalBandwidth",WHERE) );
1279 :
1280 0 : Double fbw=1.0;
1281 :
1282 0 : for(uInt i=0; i<itsCoordSys.nCoordinates(); i++)
1283 : {
1284 0 : if( itsCoordSys.type(i) == Coordinate::SPECTRAL )
1285 : {
1286 0 : SpectralCoordinate speccoord(itsCoordSys.spectralCoordinate(i));
1287 0 : Double startfreq=0.0,startpixel=-0.5;
1288 0 : Double endfreq=0.0,endpixel=+0.5;
1289 0 : speccoord.toWorld(startfreq,startpixel);
1290 0 : speccoord.toWorld(endfreq,endpixel);
1291 0 : Double midfreq = (endfreq+startfreq)/2.0;
1292 0 : fbw = ((endfreq - startfreq)/midfreq) * 100.0;
1293 0 : os << "MFS frequency range : " << startfreq/1.0e+9 << " GHz -> " << endfreq/1.0e+9 << "GHz.";
1294 0 : os << "Fractional Bandwidth : " << fbw << " %.";
1295 0 : os << "Reference Frequency for Taylor Expansion : "<< getReferenceFrequency()/1.0e+9 << "GHz." << LogIO::POST;
1296 : }
1297 : }
1298 0 : return fbw;
1299 : }
1300 :
1301 :
1302 : // Check for non-zero model (this is different from getting model flux, for derived SIIMMT)
1303 0 : Bool SIImageStoreMultiTerm::isModelEmpty()
1304 : {
1305 : /// There MUST be a more efficient way to do this !!!!! I hope.
1306 : /// Maybe save this info and change it anytime model() is accessed....
1307 0 : Bool emptymodel=true;
1308 0 : for(uInt tix=0;tix<itsNTerms;tix++)
1309 : {
1310 : //if( fabs( getModelFlux(tix) ) > 1e-08 ) emptymodel=false;
1311 0 : if( doesImageExist(itsImageName+String(".model.tt")+String::toString(tix)) &&
1312 0 : fabs( getModelFlux(tix) ) > 1e-08 ) emptymodel=false;
1313 : }
1314 0 : return emptymodel;
1315 : }
1316 :
1317 0 : void SIImageStoreMultiTerm::printImageStats()
1318 : {
1319 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","printImageStats",WHERE) );
1320 : // FIXME minresmask needs to be initialized here, or else compiler complains
1321 0 : Float minresmask = 0, maxresmask, minres, maxres;
1322 0 : ArrayLattice<Bool> pixelmask(residual()->getMask());
1323 :
1324 : // findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
1325 :
1326 0 : if (hasMask())
1327 : {
1328 : // findMinMaxLattice(*residual(), *mask() , maxres,maxresmask, minres, minresmask);
1329 0 : findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
1330 : }
1331 : else
1332 : {
1333 0 : LatticeExpr<Float> reswithpixmask(iif(pixelmask, *residual(), 0));
1334 : //LatticeExprNode pres( max( *residual() ) );
1335 0 : LatticeExprNode pres( max( reswithpixmask ) );
1336 0 : maxres = pres.getFloat();
1337 : //LatticeExprNode pres2( min( *residual() ) );
1338 0 : LatticeExprNode pres2( min( reswithpixmask ) );
1339 0 : minres = pres2.getFloat();
1340 : }
1341 :
1342 0 : os << "[" << itsImageName << "]" ;
1343 0 : os << " Peak residual (max,min) " ;
1344 0 : if( minresmask!=0.0 || maxresmask!=0.0 )
1345 0 : { os << "within mask : (" << maxresmask << "," << minresmask << ") "; }
1346 0 : os << "over full image : (" << maxres << "," << minres << ")" << LogIO::POST;
1347 :
1348 0 : os << "[" << itsImageName << "] Total Model Flux : " ;
1349 0 : for(uInt tix=0;tix<itsNTerms;tix++)
1350 0 : {os << getModelFlux(tix) << "(tt" << String::toString(tix) << ")"; }
1351 0 : os<<LogIO::POST;
1352 :
1353 0 : }
1354 :
1355 0 : std::shared_ptr<SIImageStore> SIImageStoreMultiTerm::getSubImageStore(const Int facet, const Int nfacets,
1356 : const Int chan, const Int nchanchunks,
1357 : const Int pol, const Int npolchunks)
1358 : {
1359 : std::shared_ptr<SIImageStore> multiTermStore =
1360 0 : std::make_shared<SIImageStoreMultiTerm>(itsModels, itsResiduals, itsPsfs, itsWeights, itsImages, itsSumWts, itsPBs, itsImagePBcors, itsMask, itsAlpha, itsBeta, itsAlphaError,itsAlphaPBcor, itsBetaPBcor, itsCoordSys, itsParentImageShape, itsImageName, itsObjectName, itsMiscInfo, facet, nfacets, chan, nchanchunks, pol, npolchunks);
1361 0 : return multiTermStore;
1362 : }
1363 :
1364 :
1365 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1366 :
1367 : //
1368 : //-------------------------------------------------------------
1369 : // Initialize the internals of the object. Perhaps other such
1370 : // initializations in the constructors can be moved here too.
1371 : //
1372 0 : void SIImageStoreMultiTerm::init()
1373 : {
1374 0 : imageExts.resize(MAX_IMAGE_IDS);
1375 :
1376 0 : imageExts(MASK)=".mask";
1377 0 : imageExts(PSF)=".psf";
1378 0 : imageExts(MODEL)=".model";
1379 0 : imageExts(RESIDUAL)=".residual";
1380 0 : imageExts(WEIGHT)=".weight";
1381 0 : imageExts(IMAGE)=".image";
1382 0 : imageExts(SUMWT)=".sumwt";
1383 0 : imageExts(GRIDWT)=".gridwt";
1384 0 : imageExts(PB)=".pb";
1385 0 : imageExts(FORWARDGRID)=".forward";
1386 0 : imageExts(BACKWARDGRID)=".backward";
1387 : // imageExts(IMAGEPBCOR)=".image.pbcor";
1388 :
1389 0 : itsParentPsfs.resize(itsPsfs.nelements());
1390 0 : itsParentModels.resize(itsModels.nelements());
1391 0 : itsParentResiduals.resize(itsResiduals.nelements());
1392 0 : itsParentWeights.resize(itsWeights.nelements());
1393 0 : itsParentImages.resize(itsImages.nelements());
1394 0 : itsParentSumWts.resize(itsSumWts.nelements());
1395 0 : itsParentImagePBcors.resize(itsImagePBcors.nelements());
1396 0 : itsParentPBs.resize(itsPBs.nelements());
1397 :
1398 0 : itsParentPsfs = itsPsfs;
1399 0 : itsParentModels=itsModels;
1400 0 : itsParentResiduals=itsResiduals;
1401 0 : itsParentWeights=itsWeights;
1402 0 : itsParentImages=itsImages;
1403 0 : itsParentSumWts=itsSumWts;
1404 0 : itsParentImagePBcors=itsImagePBcors;
1405 0 : itsParentPBs=itsPBs;
1406 :
1407 0 : itsParentMask=itsMask;
1408 :
1409 0 : itsParentImageShape = itsImageShape;
1410 0 : itsParentCoordSys = itsCoordSys;
1411 0 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 ) { itsImageShape=IPosition(4,0,0,0,0); }
1412 :
1413 0 : itsOpened=0;
1414 :
1415 0 : }
1416 :
1417 : /*
1418 : Bool SIImageStoreMultiTerm::getUseWeightImage()
1419 : { if( itsParentSumWts.nelements()==0 || ! itsParentSumWts[0] )
1420 : {return false;}
1421 : else
1422 : {
1423 : Bool ret = SIImageStore::getUseWeightImage( *(itsParentSumWts[0]) );
1424 : return ret;
1425 : }
1426 : }
1427 : */
1428 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1429 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1430 :
1431 : } //# NAMESPACE CASA - END
1432 :
|
