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