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 : #include "SIImageStoreMultiTerm.h"
64 : using namespace std;
65 :
66 : using namespace casacore;
67 : namespace casa { //# NAMESPACE CASA - BEGIN
68 :
69 : //////////////////////////////////////////////////////////////////////////////////////////////////////
70 : //////////////////////////////////////////////////////////////////////////////////////////////////////
71 :
72 0 : SIImageStoreMultiTerm::SIImageStoreMultiTerm():SIImageStore()
73 : {
74 0 : itsNTerms=0;
75 :
76 0 : itsPsfs.resize(0);
77 0 : itsModels.resize(0);
78 0 : itsResiduals.resize(0);
79 0 : itsWeights.resize(0);
80 0 : itsImages.resize(0);
81 0 : itsSumWts.resize(0);
82 0 : itsImagePBcors.resize(0);
83 0 : itsPBs.resize(0);
84 :
85 0 : itsForwardGrids.resize(0);
86 0 : itsBackwardGrids.resize(0);
87 :
88 0 : itsUseWeight=false;
89 :
90 0 : init();
91 :
92 0 : validate();
93 :
94 0 : }
95 :
96 : // Used from SynthesisNormalizer::makeImageStore()
97 0 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(const String &imagename,
98 : const CoordinateSystem &imcoordsys,
99 : const IPosition &imshape,
100 : const String &objectname,
101 : const Record &miscinfo,
102 : const Int /*nfacets*/,
103 : const Bool /*overwrite*/,
104 : uInt ntaylorterms,
105 0 : Bool useweightimage)
106 : {
107 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","Open new Images",WHERE) );
108 :
109 0 : itsNTerms = ntaylorterms;
110 :
111 0 : itsPsfs.resize(2 * itsNTerms - 1);
112 0 : itsModels.resize(itsNTerms);
113 0 : itsResiduals.resize(itsNTerms);
114 0 : itsWeights.resize(2 * itsNTerms - 1);
115 0 : itsImages.resize(itsNTerms);
116 0 : itsSumWts.resize(2 * itsNTerms - 1);
117 0 : itsPBs.resize(itsNTerms);
118 0 : itsImagePBcors.resize(itsNTerms);
119 :
120 0 : itsMask.reset( );
121 0 : itsGridWt.reset( );
122 :
123 0 : itsForwardGrids.resize(itsNTerms);
124 0 : itsBackwardGrids.resize(2 * itsNTerms - 1);
125 :
126 : // itsNFacets = nfacets; // So that sumwt shape happens properly, via checkValidity
127 : // itsFacetId = -1;
128 0 : itsNFacets=1;
129 0 : itsFacetId=0;
130 0 : itsNChanChunks = 1;
131 0 : itsChanId = 0;
132 0 : itsNPolChunks = 1;
133 0 : itsPolId = 0;
134 :
135 0 : itsImageName = imagename;
136 0 : itsCoordSys = imcoordsys;
137 0 : itsImageShape = imshape;
138 0 : itsObjectName = objectname;
139 0 : itsMiscInfo = miscinfo;
140 :
141 0 : itsUseWeight = useweightimage;
142 :
143 0 : init();
144 :
145 0 : validate();
146 :
147 0 : }
148 :
149 : // Used from SynthesisNormalizer::makeImageStore()
150 0 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(const String &imagename, uInt ntaylorterms,
151 0 : const Bool ignorefacets, const Bool ignoresumwt)
152 : {
153 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","Open existing Images",WHERE) );
154 :
155 0 : itsNTerms = ntaylorterms;
156 :
157 0 : auto fltPtrReset = [](Block<shared_ptr<ImageInterface<Float> > >&a) {for(uInt i=0; i < a.nelements(); ++i) a[i].reset(); };
158 0 : itsPsfs.resize(2 * itsNTerms - 1); fltPtrReset(itsPsfs);
159 0 : itsModels.resize(itsNTerms); fltPtrReset(itsModels);
160 0 : itsResiduals.resize(itsNTerms); fltPtrReset(itsResiduals);
161 0 : itsWeights.resize(2 * itsNTerms - 1); fltPtrReset(itsWeights);
162 0 : itsImages.resize(itsNTerms); fltPtrReset(itsImages);
163 0 : itsPBs.resize(itsNTerms); fltPtrReset(itsPBs);
164 0 : itsSumWts.resize(2 * itsNTerms - 1); fltPtrReset(itsSumWts);
165 0 : itsMask.reset( );
166 0 : itsGridWt.reset( );
167 0 : itsImagePBcors.resize(itsNTerms); fltPtrReset(itsImagePBcors);
168 :
169 :
170 :
171 0 : itsMiscInfo=Record();
172 0 : auto cmplxPtrReset = [](Block<shared_ptr<ImageInterface<Complex> > >&a) {for(uInt i=0; i < a.nelements(); ++i) a[i].reset(); };
173 0 : itsForwardGrids.resize(itsNTerms); cmplxPtrReset(itsForwardGrids);
174 0 : itsBackwardGrids.resize(2 * itsNTerms - 1); cmplxPtrReset(itsBackwardGrids);
175 :
176 0 : itsImageName = imagename;
177 :
178 0 : itsNFacets=1;
179 0 : itsFacetId=0;
180 0 : itsNChanChunks = 1;
181 0 : itsChanId = 0;
182 0 : itsNPolChunks = 1;
183 0 : itsPolId = 0;
184 :
185 0 : Bool exists=true;
186 0 : Bool sumwtexists=true;
187 0 : Bool modelexists=true;
188 0 : for(uInt tix=0;tix<2*itsNTerms-1;tix++)
189 : {
190 0 : if( tix<itsNTerms ) {
191 0 : exists &= ( doesImageExist( itsImageName+String(".residual.tt")+String::toString(tix) ) ||
192 0 : doesImageExist( itsImageName+String(".psf.tt")+String::toString(tix) ) );
193 0 : modelexists &= ( doesImageExist( itsImageName+String(".model.tt")+String::toString(tix) ) ||
194 0 : doesImageExist( itsImageName+String(".model.tt")+String::toString(tix) ) );
195 0 : sumwtexists &= ( doesImageExist( itsImageName+String(".sumwt.tt")+String::toString(tix) ) );
196 : }else {
197 0 : exists &= ( doesImageExist( itsImageName+String(".psf.tt")+String::toString(tix) ) );
198 0 : sumwtexists &= ( doesImageExist( itsImageName+String(".sumwt.tt")+String::toString(tix) ) );
199 : }
200 : }
201 :
202 : // The PSF or Residual images must exist. ( or the gridwt image)
203 : // All this is just for the shape and coordinate system.
204 0 : if( exists || modelexists || doesImageExist(itsImageName+String(".gridwt")) )
205 : {
206 0 : std::shared_ptr<ImageInterface<Float> > imptr;
207 0 : if( doesImageExist(itsImageName+String(".psf.tt0")) )
208 : {
209 0 : buildImage( imptr, (itsImageName+String(".psf.tt0")) );
210 :
211 : //cout << "Opening PSF image to read csys" << endl;
212 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".psf.tt0")) );
213 : }
214 0 : else if( doesImageExist(itsImageName+String(".residual.tt0")) )
215 : {
216 0 : buildImage( imptr, (itsImageName+String(".residual.tt0")) );
217 : //cout << "Opening Residual image to read csys" << endl;
218 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".residual.tt0")) );
219 : }
220 0 : else if( doesImageExist(itsImageName+String(".model.tt0")) )
221 : {
222 0 : buildImage( imptr, (itsImageName+String(".model.tt0")) );
223 : //cout << "Opening Model image to read csys" << endl;
224 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".model.tt0")) );
225 : }
226 : else
227 : {
228 : // How can this be right ?
229 : //cout << "Opening Sumwt image to read csys" << endl;
230 0 : buildImage( imptr, (itsImageName+String(".gridwt")) );
231 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".gridwt")) );
232 : }
233 :
234 0 : itsObjectName=imptr->imageInfo().objectName();
235 0 : itsImageShape = imptr->shape();
236 0 : itsCoordSys = imptr->coordinates();
237 0 : itsMiscInfo=imptr->miscInfo();
238 :
239 0 : }
240 : else
241 : {
242 0 : throw( AipsError( "Multi-term PSF, Residual or Model Images do not exist. Please create one of them." ) );
243 : }
244 :
245 0 : if( doesImageExist(itsImageName+String(".residual.tt0")) ||
246 0 : doesImageExist(itsImageName+String(".psf.tt0")) )
247 : {
248 0 : if( sumwtexists )
249 : {
250 0 : std::shared_ptr<ImageInterface<Float> > imptr;
251 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt.tt0")) );
252 0 : buildImage( imptr, (itsImageName+String(".sumwt.tt0")) );
253 0 : itsNFacets = imptr->shape()[0];
254 0 : itsFacetId = 0;
255 0 : itsUseWeight = getUseWeightImage( *imptr );
256 : /////redo this here as psf may have different coordinates
257 0 : itsCoordSys = imptr->coordinates();
258 0 : itsMiscInfo=imptr->miscInfo();
259 0 : if(!ignoresumwt)
260 : {
261 0 : if( itsUseWeight && ! doesImageExist(itsImageName+String(".weight.tt0")) )
262 : {
263 0 : throw(AipsError("Internal error : MultiTerm Sumwt has a useweightimage=true but the weight image does not exist."));
264 : }
265 : }
266 0 : }
267 : else
268 : {
269 0 : if(!ignoresumwt)
270 0 : {throw( AipsError( "Multi-term SumWt does not exist. Please create PSFs or Residuals." ) );}
271 : else
272 : {
273 0 : os << "SumWt.ttx do not exist. Proceeding only with PSFs" << LogIO::POST;
274 0 : std::shared_ptr<ImageInterface<Float> > imptr;
275 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt.tt0")) );
276 0 : if( doesImageExist(itsImageName+String(".residual.tt0")) )
277 0 : {buildImage( imptr, (itsImageName+String(".residual.tt0")) );}
278 : else
279 0 : {buildImage( imptr, (itsImageName+String(".psf.tt0")) );}
280 :
281 0 : itsNFacets = 1;
282 0 : itsFacetId = 0;
283 0 : itsUseWeight = False;
284 0 : itsCoordSys = imptr->coordinates();
285 0 : itsMiscInfo=imptr->miscInfo();
286 0 : }
287 : }
288 : }// if psf0 or res0 exist
289 :
290 0 : if( ignorefacets==true ) itsNFacets=1;
291 :
292 0 : init();
293 0 : validate();
294 :
295 0 : }
296 :
297 : /*
298 : /////////////Constructor with pointers already created else where but taken over here
299 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(Block<std::shared_ptr<ImageInterface<Float> > > modelims,
300 : Block<std::shared_ptr<ImageInterface<Float> > >residims,
301 : Block<std::shared_ptr<ImageInterface<Float> > >psfims,
302 : Block<std::shared_ptr<ImageInterface<Float> > >weightims,
303 : Block<std::shared_ptr<ImageInterface<Float> > >restoredims,
304 : Block<std::shared_ptr<ImageInterface<Float> > >sumwtims,
305 : std::shared_ptr<ImageInterface<Float> > newmask,
306 : std::shared_ptr<ImageInterface<Float> > newalpha,
307 : std::shared_ptr<ImageInterface<Float> > newbeta)
308 : {
309 :
310 : itsPsfs=psfims;
311 : itsModels=modelims;
312 : itsResiduals=residims;
313 : itsWeights=weightims;
314 : itsImages=restoredims;
315 : itsSumWts=sumwtims;
316 : itsMask = newmask;
317 : itsAlpha = newalpha;
318 : itsBeta = newbeta;
319 :
320 : itsNTerms = itsResiduals.nelements();
321 : itsMiscInfo=Record();
322 :
323 : AlwaysAssert( itsPsfs.nelements() == 2*itsNTerms-1 , AipsError );
324 : AlwaysAssert( itsPsfs.nelements()>0 && itsPsfs[0] , AipsError );
325 : AlwaysAssert( itsSumWts.nelements()>0 && itsSumWts[0] , AipsError );
326 :
327 : itsForwardGrids.resize( itsNTerms );
328 : itsBackwardGrids.resize( 2 * itsNTerms - 1 );
329 :
330 : itsImageShape=psfims[0]->shape();
331 : itsCoordSys = psfims[0]->coordinates();
332 : itsMiscInfo = psfims[0]->miscInfo();
333 :
334 : itsNFacets=sumwtims[0]->shape()[0];
335 : itsUseWeight=getUseWeightImage( *(sumwtims[0]) );
336 :
337 : itsImageName = String("reference"); // This is what the access functions use to guard against allocs...
338 :
339 : init();
340 : validate();
341 :
342 : }
343 : */
344 :
345 : //////////////////////////////////////////////////////////////////////////////////////////////////////
346 : /////////////Constructor with pointers already created else where but taken over here
347 : // used from getSubImageStore(), for example when creating the facets list
348 : // this would be safer if it was refactored as a copy constructor of the generic stuff +
349 : // initialization of the facet related parameters
350 0 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(const Block<std::shared_ptr<ImageInterface<Float> > > &modelims,
351 : const Block<std::shared_ptr<ImageInterface<Float> > > &residims,
352 : const Block<std::shared_ptr<ImageInterface<Float> > > &psfims,
353 : const Block<std::shared_ptr<ImageInterface<Float> > > &weightims,
354 : const Block<std::shared_ptr<ImageInterface<Float> > > &restoredims,
355 : const Block<std::shared_ptr<ImageInterface<Float> > > &sumwtims,
356 : const Block<std::shared_ptr<ImageInterface<Float> > > &pbims,
357 : const Block<std::shared_ptr<ImageInterface<Float> > > &restoredpbcorims,
358 : const std::shared_ptr<ImageInterface<Float> > &newmask,
359 : const std::shared_ptr<ImageInterface<Float> > &newalpha,
360 : const std::shared_ptr<ImageInterface<Float> > &newbeta,
361 : const std::shared_ptr<ImageInterface<Float> > &newalphaerror,
362 : const std::shared_ptr<ImageInterface<Float> > &newalphapbcor,
363 : const std::shared_ptr<ImageInterface<Float> > &newbetapbcor,
364 : const CoordinateSystem& csys,
365 : const IPosition &imshape,
366 : const String &imagename,
367 : const String &objectname,
368 : const Record &miscinfo,
369 : const Int facet, const Int nfacets,
370 : const Int chan, const Int nchanchunks,
371 0 : const Int pol, const Int npolchunks)
372 : {
373 0 : itsPsfs=psfims;
374 0 : itsModels=modelims;
375 0 : itsResiduals=residims;
376 0 : itsWeights=weightims;
377 0 : itsImages=restoredims;
378 0 : itsSumWts=sumwtims;
379 0 : itsMask = newmask;
380 0 : itsAlpha = newalpha;
381 0 : itsBeta = newbeta;
382 0 : itsAlphaError = newalphaerror;
383 0 : itsAlphaPBcor = newalphapbcor;
384 0 : itsBetaPBcor = newbetapbcor;
385 :
386 0 : itsPBs=pbims;
387 0 : itsImagePBcors=restoredpbcorims;
388 :
389 0 : itsNTerms = itsResiduals.nelements();
390 0 : itsMiscInfo=Record();
391 :
392 0 : AlwaysAssert( itsPsfs.nelements() == 2*itsNTerms-1 , AipsError );
393 : // AlwaysAssert( itsPsfs.nelements()>0 && itsPsfs[0] , AipsError );
394 : // AlwaysAssert( itsSumWts.nelements()>0 && itsSumWts[0] , AipsError );
395 :
396 0 : itsForwardGrids.resize( itsNTerms );
397 0 : itsBackwardGrids.resize( 2 * itsNTerms - 1 );
398 :
399 0 : itsObjectName = objectname;
400 0 : itsMiscInfo = miscinfo;
401 :
402 0 : itsNFacets = nfacets;
403 0 : itsFacetId = facet;
404 0 : itsNChanChunks = nchanchunks;
405 0 : itsChanId = chan;
406 0 : itsNPolChunks = npolchunks;
407 0 : itsPolId = pol;
408 :
409 0 : itsParentImageShape = imshape;
410 0 : itsImageShape = imshape;
411 : ///// itsImageShape = IPosition(4,0,0,0,0);
412 :
413 0 : itsCoordSys = csys; // Hopefully this doesn't change for a reference image
414 0 : itsImageName = imagename;
415 :
416 :
417 : //-----------------------
418 0 : init(); // Connect parent pointers to the images.
419 : //-----------------------
420 :
421 : // Set these to null, to be set later upon first access.
422 : // Setting to null will hopefully set all elements of each array, to NULL.
423 0 : itsPsfs=std::shared_ptr<ImageInterface<Float> >();
424 0 : itsModels=std::shared_ptr<ImageInterface<Float> >();
425 0 : itsResiduals=std::shared_ptr<ImageInterface<Float> >();
426 0 : itsWeights=std::shared_ptr<ImageInterface<Float> >();
427 0 : itsImages=std::shared_ptr<ImageInterface<Float> >();
428 0 : itsSumWts=std::shared_ptr<ImageInterface<Float> >();
429 0 : itsPBs=std::shared_ptr<ImageInterface<Float> >();
430 :
431 0 : itsMask.reset( );
432 :
433 0 : validate();
434 :
435 0 : }
436 :
437 : //////////////////////////////////////////////////////////////////////////////////////////////////////
438 :
439 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
440 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
441 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
442 :
443 0 : uInt SIImageStoreMultiTerm::getNTaylorTerms(Bool dopsf)
444 : {
445 0 : return dopsf ? (2*itsNTerms-1) : itsNTerms;
446 : }
447 :
448 : // Check if images that are asked-for are ready and all have the same shape.
449 : /*
450 : Bool SIImageStoreMultiTerm::checkValidity(const Bool ipsf, const Bool iresidual,
451 : const Bool iweight, const Bool imodel, const Bool irestored,
452 : const Bool imask,const Bool isumwt,
453 : const Bool ialpha, const Bool ibeta)
454 : {
455 :
456 : // cout << "In MT::checkValidity imask is " << imask << endl;
457 :
458 : Bool valid = true;
459 :
460 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
461 : {
462 :
463 : if(ipsf==true)
464 : { psf(tix);
465 : valid = valid & ( itsPsfs[tix] && itsPsfs[tix]->shape()==itsImageShape ); }
466 : if(iweight==true)
467 : { weight(tix);
468 : valid = valid & ( itsWeights[tix] && itsWeights[tix]->shape()==itsImageShape ); }
469 :
470 : if(isumwt==true) {
471 : IPosition useShape(itsImageShape);
472 : useShape[0]=itsNFacets; useShape[1]=itsNFacets;
473 : sumwt(tix);
474 : valid = valid & ( itsSumWts[tix] && itsSumWts[tix]->shape()==useShape );
475 : }
476 :
477 : if( tix< itsNTerms )
478 : {
479 : if(iresidual==true)
480 : { residual(tix);
481 : valid = valid & ( itsResiduals[tix] && itsResiduals[tix]->shape()==itsImageShape ); }
482 : if(imodel==true)
483 : { model(tix);
484 : valid = valid & ( itsModels[tix] && itsModels[tix]->shape()==itsImageShape); }
485 : if(irestored==true)
486 : { image(tix);
487 : valid = valid & ( itsImages[tix] && itsImages[tix]->shape()==itsImageShape); }
488 : }
489 : }
490 :
491 : if(imask==true)
492 : { mask(); valid = valid & ( itsMask && itsMask->shape()==itsImageShape);
493 : // cout << " Mask null ? " << (bool) itsMask << endl;
494 : }
495 : if(ialpha==true)
496 : { alpha(); valid = valid & ( itsAlpha && itsAlpha->shape()==itsImageShape ); }
497 : if(ibeta==true)
498 : { beta(); valid = valid & ( itsBeta && itsBeta->shape()==itsImageShape ); }
499 :
500 : return valid;
501 :
502 : }
503 : */
504 :
505 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
506 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
507 : //////////////////////////////////////////////////////////////////////////////////////////////////////
508 : //////////////////////////////////////////////////////////////////////////////////////////////////////
509 :
510 0 : SIImageStoreMultiTerm::~SIImageStoreMultiTerm()
511 : {
512 0 : }
513 :
514 : //////////////////////////////////////////////////////////////////////////////////////////////////////
515 : //////////////////////////////////////////////////////////////////////////////////////////////////////
516 :
517 0 : Bool SIImageStoreMultiTerm::releaseLocks()
518 : {
519 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","releaseLocks",WHERE) );
520 :
521 0 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
522 : {
523 0 : if( itsPsfs[tix] ) releaseImage( itsPsfs[tix] );
524 0 : if( itsWeights[tix] ) releaseImage( itsWeights[tix] );
525 0 : if( itsSumWts[tix] ) releaseImage( itsSumWts[tix] );
526 0 : if( tix < itsNTerms )
527 : {
528 0 : if( itsModels[tix] ) releaseImage( itsModels[tix] );
529 0 : if( itsResiduals[tix] ) releaseImage( itsResiduals[tix] );
530 0 : if( itsImages[tix] ) releaseImage( itsImages[tix] );
531 0 : if( itsPBs[tix] ) releaseImage( itsPBs[tix] );
532 0 : if( itsImagePBcors[tix] ) releaseImage( itsImagePBcors[tix] );
533 : }
534 : }
535 0 : if( itsMask ) releaseImage( itsMask );
536 0 : if( itsAlpha ) releaseImage( itsAlpha );
537 0 : if( itsBeta ) releaseImage( itsBeta );
538 0 : if( itsAlphaError ) releaseImage( itsAlphaError );
539 0 : if( itsAlphaPBcor ) releaseImage( itsAlphaPBcor );
540 0 : if( itsBetaPBcor ) releaseImage( itsBetaPBcor );
541 0 : if( itsGridWt ) releaseImage( itsGridWt );
542 :
543 0 : return true; // do something more intelligent here.
544 0 : }
545 :
546 0 : Bool SIImageStoreMultiTerm::releaseComplexGrids()
547 : {
548 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","releaseComplexGrids",WHERE) );
549 :
550 0 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
551 : {
552 0 : if( itsBackwardGrids[tix] ) releaseImage( itsBackwardGrids[tix] );
553 0 : if( tix < itsNTerms )
554 : {
555 0 : if( itsForwardGrids[tix] ) releaseImage( itsForwardGrids[tix] );
556 : }
557 : }
558 0 : return True; // do something more intelligent here.
559 0 : }
560 :
561 :
562 0 : Double SIImageStoreMultiTerm::getReferenceFrequency()
563 : {
564 : Double theRefFreq;
565 :
566 0 : Vector<Double> refpix = (itsCoordSys.spectralCoordinate()).referencePixel();
567 0 : AlwaysAssert( refpix.nelements()>0, AipsError );
568 0 : (itsCoordSys.spectralCoordinate()).toWorld( theRefFreq, refpix[0] );
569 : // cout << "Reading ref freq as : " << theRefFreq << endl;
570 0 : return theRefFreq;
571 0 : }
572 :
573 : //////////////////////////////////////////////////////////////////////////////////////////////////////
574 : //////////////////////////////////////////////////////////////////////////////////////////////////////
575 0 : Vector<String> SIImageStoreMultiTerm::getModelImageName()
576 : {
577 0 : Vector<String> mods(itsNTerms);
578 0 : for(uInt tix=0;tix<itsNTerms;tix++)
579 0 : {mods[tix]=itsImageName + imageExts(MODEL)+".tt"+String::toString(tix);}
580 0 : return mods;
581 0 : }
582 :
583 0 : void SIImageStoreMultiTerm::setModelImage( const Vector<String> &modelnames )
584 : {
585 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","setModelImage",WHERE) );
586 :
587 0 : if( modelnames.nelements() > itsNTerms )
588 0 : { throw(AipsError("We currently cannot support more than nterms images as the starting model. "));
589 : }
590 :
591 0 : if( modelnames.nelements() > 0 && modelnames.nelements() <= itsNTerms )
592 : {
593 0 : for(uInt tix=0;tix<modelnames.nelements();tix++)
594 : {
595 0 : setModelImageOne( modelnames[tix], tix );
596 : }
597 : }
598 :
599 0 : }
600 :
601 : /*
602 : void SIImageStoreMultiTerm::setModelImage( String modelname )
603 : {
604 : LogIO os( LogOrigin("SIImageStoreMultiTerm","setModelImage",WHERE) );
605 :
606 : for(uInt tix=0;tix<itsNTerms;tix++)
607 : {
608 :
609 : Directory immodel( modelname+String(".model.tt")+String::toString(tix) );
610 : if( !immodel.exists() )
611 : {
612 : os << "Starting model image does not exist for term : " << tix << LogIO::POST;
613 : }
614 : else
615 : {
616 : std::shared_ptr<PagedImage<Float> > newmodel( new PagedImage<Float>( modelname+String(".model.tt")+String::toString(tix) ) );
617 : // Check shapes, coordsys with those of other images. If different, try to re-grid here.
618 :
619 : if( newmodel->shape() != model(tix)->shape() )
620 : {
621 : os << "Regridding input model to target coordinate system for term " << tix << LogIO::POST;
622 : regridToModelImage( *newmodel , tix);
623 : // For now, throw an exception.
624 : //throw( AipsError( "Input model image "+modelname+".model.tt"+String::toString(tix)+" is not the same shape as that defined for output in "+ itsImageName + ".model" ) );
625 : }
626 : else
627 : {
628 : os << "Setting " << modelname << " as model for term " << tix << LogIO::POST;
629 : // Then, add its contents to itsModel.
630 : //itsModel->put( itsModel->get() + model->get() );
631 : /////( model(tix) )->put( newmodel->get() );
632 : model(tix)->copyData( LatticeExpr<Float> (*newmodel) );
633 : }
634 : }
635 : }//nterms
636 : }
637 : */
638 :
639 : //////////////////////////////////////////////////////////////////////////////////////////////////////
640 : //////////////////////////////////////////////////////////////////////////////////////////////////////
641 :
642 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::psf(uInt term)
643 : {
644 0 : AlwaysAssert( itsPsfs.nelements() > term, AipsError );
645 0 : accessImage( itsPsfs[term], itsParentPsfs[term], imageExts(PSF)+".tt"+String::toString(term) );
646 0 : return itsPsfs[term];
647 : }
648 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::residual(uInt term)
649 : {
650 0 : accessImage( itsResiduals[term], itsParentResiduals[term], imageExts(RESIDUAL)+".tt"+String::toString(term) );
651 : // Record mi = itsResiduals[term]->miscInfo(); ostringstream oss;mi.print(oss);cout<<"MiscInfo(res) " << term << " : " << oss.str() << endl;
652 :
653 0 : return itsResiduals[term];
654 : }
655 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::weight(uInt term)
656 : {
657 0 : accessImage( itsWeights[term], itsParentWeights[term], imageExts(WEIGHT)+".tt"+String::toString(term) );
658 0 : return itsWeights[term];
659 : }
660 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::sumwt(uInt term)
661 : {
662 0 : accessImage( itsSumWts[term], itsParentSumWts[term], imageExts(SUMWT)+".tt"+String::toString(term) );
663 :
664 :
665 0 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
666 0 : {itsUseWeight = getUseWeightImage( *itsParentSumWts[0] );}
667 0 : setUseWeightImage( *(itsSumWts[term]) , itsUseWeight); // Sets a flag in the SumWt image.
668 :
669 0 : return itsSumWts[term];
670 : }
671 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::model(uInt term)
672 : {
673 :
674 0 : accessImage( itsModels[term], itsParentModels[term], imageExts(MODEL)+".tt"+String::toString(term) );
675 :
676 0 : itsModels[term]->setUnits("Jy/pixel");
677 0 : return itsModels[term];
678 : }
679 :
680 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::image(uInt term)
681 : {
682 :
683 0 : accessImage( itsImages[term], itsParentImages[term], imageExts(IMAGE)+".tt"+String::toString(term) );
684 0 : itsImages[term]->setUnits("Jy/beam");
685 0 : return itsImages[term];
686 : }
687 :
688 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::pb(uInt term)
689 : {
690 :
691 0 : accessImage( itsPBs[term], itsParentPBs[term], imageExts(PB)+".tt"+String::toString(term) );
692 0 : return itsPBs[term];
693 : }
694 :
695 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::imagepbcor(uInt term)
696 : {
697 :
698 0 : accessImage( itsImagePBcors[term], itsParentImagePBcors[term], imageExts(IMAGE)+".tt"+String::toString(term)+ ".pbcor" );
699 0 : itsImagePBcors[term]->setUnits("Jy/beam");
700 0 : return itsImagePBcors[term];
701 : }
702 :
703 0 : std::shared_ptr<ImageInterface<Complex> > SIImageStoreMultiTerm::forwardGrid(uInt term){
704 0 : if( itsForwardGrids[term] )// && (itsForwardGrids[term]->shape() == itsImageShape))
705 0 : return itsForwardGrids[term];
706 0 : Vector<Int> whichStokes(0);
707 0 : IPosition cimageShape;
708 0 : cimageShape=itsImageShape;
709 0 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
710 0 : whichStokes, itsDataPolRep);
711 0 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
712 0 : cimageShape(2)=whichStokes.nelements();
713 0 : itsForwardGrids[term].reset(new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()));
714 0 : return itsForwardGrids[term];
715 0 : }
716 0 : std::shared_ptr<ImageInterface<Complex> > SIImageStoreMultiTerm::backwardGrid(uInt term){
717 :
718 0 : if( itsBackwardGrids[term] && (itsBackwardGrids[term]->shape() == itsImageShape))
719 0 : return itsBackwardGrids[term];
720 : // cout << "MT : Making backward grid of shape : " << itsImageShape << endl;
721 0 : Vector<Int> whichStokes(0);
722 0 : IPosition cimageShape;
723 0 : cimageShape=itsImageShape;
724 0 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
725 0 : whichStokes, itsDataPolRep);
726 0 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
727 0 : cimageShape(2)=whichStokes.nelements();
728 0 : itsBackwardGrids[term].reset(new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()));
729 0 : return itsBackwardGrids[term];
730 0 : }
731 :
732 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::alpha()
733 : {
734 0 : if( itsAlpha && itsAlpha->shape() == itsImageShape ) { return itsAlpha; }
735 : // checkRef( itsAlpha , "alpha" );
736 0 : itsAlpha = openImage( itsImageName+String(".alpha"), false );
737 : // itsAlpha->setUnits("Alpha");
738 0 : return itsAlpha;
739 : }
740 :
741 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::beta()
742 : {
743 0 : if( itsBeta && itsBeta->shape() == itsImageShape ) { return itsBeta; }
744 : // checkRef( itsBeta , "beta" );
745 0 : itsBeta = openImage( itsImageName+String(".beta"), false );
746 : // itsBeta->setUnits("Beta");
747 0 : return itsBeta;
748 : }
749 :
750 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::alphaerror()
751 : {
752 0 : if( itsAlphaError && itsAlphaError->shape() == itsImageShape ) { return itsAlphaError; }
753 : // checkRef( itsAlpha , "alpha" );
754 0 : itsAlphaError = openImage( itsImageName+String(".alpha.error"), false );
755 : // itsAlpha->setUnits("Alpha");
756 0 : return itsAlphaError;
757 : }
758 :
759 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::alphapbcor()
760 : {
761 0 : if( itsAlphaPBcor && itsAlphaPBcor->shape() == itsImageShape ) { return itsAlphaPBcor; }
762 : // checkRef( itsAlphaPBcor , "alpha" );
763 0 : itsAlphaPBcor = openImage( itsImageName+String(".alpha.pbcor"), False );
764 : // itsAlphaPBcor->setUnits("Alpha");
765 0 : return itsAlphaPBcor;
766 : }
767 :
768 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::betapbcor()
769 : {
770 0 : if( itsBetaPBcor && itsBetaPBcor->shape() == itsImageShape ) { return itsBetaPBcor; }
771 : // checkRef( itsBetaPBcor , "beta" );
772 0 : itsBetaPBcor = openImage( itsImageName+String(".beta.pbcor"), False );
773 : // itsBetaPBcor->setUnits("Beta");
774 0 : return itsBetaPBcor;
775 : }
776 :
777 :
778 :
779 : // TODO : Move to an image-wrapper class ? Same function exists in SynthesisDeconvolver.
780 0 : Bool SIImageStoreMultiTerm::doesImageExist(String imagename)
781 : {
782 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","doesImageExist",WHERE) );
783 0 : Directory image( imagename );
784 0 : return image.exists();
785 0 : }
786 :
787 :
788 0 : void SIImageStoreMultiTerm::resetImages( Bool resetpsf, Bool resetresidual, Bool resetweight )
789 : {
790 0 : for(uInt tix=0;tix<2*itsNTerms-1;tix++)
791 : {
792 0 : if( resetpsf ) psf(tix)->set(0.0);
793 :
794 0 : if( tix < itsNTerms ) {
795 0 : if( resetresidual ) {
796 : //removeMask( residual(tix) );
797 0 : residual(tix)->set(0.0);
798 : }
799 : }
800 0 : if( resetweight && itsWeights[tix] ) weight(tix)->set(0.0);
801 0 : if( resetweight ) sumwt(tix)->set(0.0);
802 : }//nterms
803 0 : }
804 :
805 0 : void SIImageStoreMultiTerm::addImages( std::shared_ptr<SIImageStore> imagestoadd,
806 : Bool addpsf, Bool addresidual, Bool addweight, Bool adddensity)
807 : {
808 0 : for(uInt tix=0;tix<2*itsNTerms-1;tix++)
809 : {
810 :
811 0 : if(addpsf)
812 : {
813 0 : LatticeExpr<Float> adderPsf( *(psf(tix)) + *(imagestoadd->psf(tix)) );
814 0 : psf(tix)->copyData(adderPsf);
815 0 : }
816 0 : if(addweight)
817 : {
818 :
819 0 : if(getUseWeightImage( *(imagestoadd->sumwt(tix)) ) ) // Access and add weight only if it is needed.
820 : {
821 0 : LatticeExpr<Float> adderWeight( *(weight(tix)) + *(imagestoadd->weight(tix)) );
822 0 : weight(tix)->copyData(adderWeight);
823 0 : }
824 :
825 0 : LatticeExpr<Float> adderSumWt( *(sumwt(tix)) + *(imagestoadd->sumwt(tix)) );
826 0 : sumwt(tix)->copyData(adderSumWt);
827 0 : setUseWeightImage( *sumwt(tix), getUseWeightImage(*(imagestoadd->sumwt(tix)) ) );
828 0 : }
829 :
830 0 : if(tix < itsNTerms && addresidual)
831 : {
832 0 : LatticeExpr<Float> adderRes( *(residual(tix)) + *(imagestoadd->residual(tix)) );
833 0 : residual(tix)->copyData(adderRes);
834 0 : }
835 :
836 0 : if( tix==0 && adddensity )
837 : {
838 0 : LatticeExpr<Float> adderDensity( *(gridwt()) + *(imagestoadd->gridwt()) );
839 0 : gridwt()->copyData(adderDensity);
840 0 : }
841 :
842 : }
843 0 : }
844 :
845 0 : void SIImageStoreMultiTerm::dividePSFByWeight(const Float /*pblimit*/)
846 : {
847 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","dividePSFByWeight",WHERE) );
848 :
849 : //// for(uInt tix=0;tix<2*itsNTerms-1;tix++)
850 0 : for(Int tix=2*itsNTerms-1-1;tix>-1;tix--) // AAH go backwards so that zeroth term is normalized last..... sigh sigh sigh.
851 : {
852 :
853 : //cout << "npsfs : " << itsPsfs.nelements() << " and tix : " << tix << endl;
854 :
855 0 : normPSF(tix);
856 :
857 : //if ( itsUseWeight) {
858 : // divideImageByWeightVal( *weight(tix) );
859 : //}
860 :
861 : }
862 0 : }
863 :
864 0 : void SIImageStoreMultiTerm::divideWeightBySumWt()
865 : {
866 0 : LogIO os(LogOrigin("SIImageStoreMultiTerm", "divideWeightByWeight", WHERE));
867 :
868 : //// for(uInt tix=0;tix<2*itsNTerms-1;tix++)
869 0 : for (Int tix = 2 * itsNTerms - 1 - 1;tix > -1;tix--) // AAH go backwards so that zeroth term is normalized last..... sigh sigh sigh.
870 : {
871 :
872 0 : if ( itsUseWeight) {
873 0 : divideImageByWeightVal( *weight(tix) );
874 : }
875 :
876 : }
877 0 : }
878 0 : void SIImageStoreMultiTerm::normalizePrimaryBeam(const Float pblimit)
879 : {
880 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","normalizePrimaryBeam",WHERE) );
881 0 : if ( itsUseWeight) {
882 :
883 0 : makePBFromWeight(pblimit);
884 : }
885 0 : else { makePBImage(pblimit); }
886 : // calcSensitivity();
887 0 : }
888 :
889 :
890 : // Make another for the PSF too.
891 0 : void SIImageStoreMultiTerm::divideResidualByWeight(Float pblimit, const String normtype)
892 : {
893 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","divideResidualByWeight",WHERE) );
894 :
895 0 : if( itsUseWeight )
896 : {
897 0 : itsPBScaleFactor = getPbMax();
898 : }
899 :
900 0 : for(uInt tix=0;tix<itsNTerms;tix++)
901 : {
902 :
903 0 : divideImageByWeightVal( *residual(tix) );
904 :
905 : // if(doesImageExist(itsImageName+String(".weight.tt0")) )
906 0 : if( itsUseWeight )
907 : {
908 0 : LatticeExpr<Float> ratio;
909 0 : Float scalepb = fabs(pblimit) * itsPBScaleFactor * itsPBScaleFactor ;
910 0 : if( normtype=="flatnoise"){
911 0 : LatticeExpr<Float> deno = LatticeExpr<Float> ( sqrt( abs(*(weight(0)) ) ) * itsPBScaleFactor );
912 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".residual.tt")+String::toString(tix) ;
913 0 : os << " by [ sqrt(weightimage) * " << itsPBScaleFactor ;
914 0 : os << " ] to get flat noise with unit pb peak."<< LogIO::POST;
915 0 : LatticeExpr<Float> mask( iif( (deno) > scalepb , 1.0, 0.0 ) );
916 0 : LatticeExpr<Float> maskinv( iif( (deno) > scalepb , 0.0, 1.0 ) );
917 0 : ratio=LatticeExpr<Float> ( ( (*(residual(tix))) * mask ) / ( deno + maskinv ) );
918 0 : }
919 0 : else if(normtype=="pbsquare"){
920 0 : Float deno = itsPBScaleFactor*itsPBScaleFactor ;
921 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".residual.tt")+String::toString(tix) ;
922 0 : os << deno ;
923 0 : os << " ] to get optimal sig/noise with unit pb peak."<< LogIO::POST;
924 :
925 0 : ratio=LatticeExpr<Float> ( ( *(residual(tix)) ) / ( deno ) );
926 :
927 :
928 :
929 : }
930 0 : else if( normtype=="flatsky") {
931 0 : LatticeExpr<Float> deno( *(weight(0)) );
932 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".residual.tt")+String::toString(tix) ;
933 0 : os << " by [ weight ] to get flat sky"<< LogIO::POST;
934 0 : LatticeExpr<Float> mask( iif( (deno) > scalepb , 1.0, 0.0 ) );
935 0 : LatticeExpr<Float> maskinv( iif( (deno) > scalepb , 0.0, 1.0 ) );
936 0 : ratio=LatticeExpr<Float> ( ( (*(residual(tix))) * mask ) / ( deno + maskinv ) );
937 0 : }
938 : else{
939 0 : throw(AipsError("Don't know how to proceed with normtype "+normtype));
940 : }
941 :
942 :
943 :
944 0 : residual(tix)->copyData(ratio);
945 0 : }
946 :
947 0 : if( (residual(tix)->getDefaultMask()=="") && hasPB() && pblimit >= 0.0 )
948 0 : {copyMask(pb(),residual(tix));}
949 :
950 0 : if( pblimit <0.0 && (residual(tix)->getDefaultMask()).matches("mask0") ) removeMask( residual(tix) );
951 :
952 : }
953 0 : }
954 :
955 0 : void SIImageStoreMultiTerm::divideModelByWeight(Float pblimit, const String normtype)
956 : {
957 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","divideModelByWeight",WHERE) );
958 :
959 0 : if( itsUseWeight // only when needed
960 0 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
961 : {
962 :
963 0 : if( normtype=="flatsky") {
964 0 : Array<Float> arrmod;
965 0 : model(0)->get( arrmod, true );
966 :
967 0 : os << "Model is already flat sky with peak flux : " << max(arrmod);
968 0 : os << ". No need to divide before prediction" << LogIO::POST;
969 :
970 0 : return;
971 0 : }
972 :
973 0 : itsPBScaleFactor = getPbMax();
974 :
975 0 : for(uInt tix=0;tix<itsNTerms;tix++)
976 0 : { LatticeExpr<Float> deno;
977 0 : if(normtype=="flatnoise"){
978 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".model")+String::toString(tix);
979 0 : os << " by [ sqrt(weight) / " << itsPBScaleFactor ;
980 0 : os <<" ] to get to flat sky model before prediction" << LogIO::POST;
981 :
982 0 : deno = LatticeExpr<Float> ( sqrt( abs(*(weight(0))) ) / itsPBScaleFactor );
983 : }
984 0 : else if(normtype=="pbsquare"){
985 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".model")+String::toString(tix);
986 0 : os << " by [ (weight) / " << itsPBScaleFactor*itsPBScaleFactor ;
987 0 : os <<" ] to get an optimal sig/noise model before prediction" << LogIO::POST;
988 :
989 0 : deno = LatticeExpr<Float> ( abs(*(weight(0))) / (itsPBScaleFactor*itsPBScaleFactor) );
990 : }
991 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
992 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
993 0 : LatticeExpr<Float> ratio( ( (*(model(tix))) * mask ) / ( deno + maskinv ) );
994 :
995 0 : itsModels[tix]->copyData(ratio);
996 0 : }
997 : }
998 0 : }
999 :
1000 :
1001 0 : void SIImageStoreMultiTerm::multiplyModelByWeight(Float pblimit, const String normtype)
1002 : {
1003 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","multiplyModelByWeight",WHERE) );
1004 :
1005 :
1006 0 : if( itsUseWeight // only when needed
1007 0 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
1008 : {
1009 :
1010 0 : if( normtype=="flatsky") {
1011 0 : os << "Model is already flat sky. No need to multiply back after prediction" << LogIO::POST;
1012 0 : return;
1013 : }
1014 :
1015 0 : itsPBScaleFactor = getPbMax();
1016 :
1017 0 : for(uInt tix=0;tix<itsNTerms;tix++)
1018 : {
1019 0 : LatticeExpr<Float> deno;
1020 0 : if( normtype=="flatnoise") {
1021 0 : os << LogIO::NORMAL1 << "Multiplying " << itsImageName+String(".model")+String::toString(tix);
1022 0 : os << " by [ sqrt(weight) / " << itsPBScaleFactor;
1023 0 : os << " ] to take model back to flat noise with unit pb peak." << LogIO::POST;
1024 :
1025 0 : deno=LatticeExpr<Float> ( sqrt( abs(*(weight(0)) ) ) / itsPBScaleFactor );
1026 : }
1027 0 : else if ( normtype=="pbsquare"){
1028 0 : os << LogIO::NORMAL1 << "Multiplying " << itsImageName+String(".model")+String::toString(tix);
1029 0 : os << " by [ weight / " << itsPBScaleFactor*itsPBScaleFactor;
1030 0 : os << " ] to take model back to optima sig/noise with unit pb peak." << LogIO::POST;
1031 :
1032 0 : deno=LatticeExpr<Float> ( abs(*(weight(0)) ) / (itsPBScaleFactor*itsPBScaleFactor) );
1033 : }
1034 : else{
1035 0 : throw(AipsError("No idea of what to do for "+normtype));
1036 : }
1037 :
1038 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1039 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1040 0 : LatticeExpr<Float> ratio( ( (*(model(tix))) * mask ) * ( deno + maskinv ) );
1041 :
1042 0 : itsModels[tix]->copyData(ratio);
1043 0 : }
1044 : }
1045 0 : }
1046 :
1047 :
1048 0 : void SIImageStoreMultiTerm::restore(GaussianBeam& rbeam, String& usebeam, uInt /*term*/, Float psfcutoff)
1049 : {
1050 :
1051 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","restore",WHERE) );
1052 :
1053 0 : for(uInt tix=0; tix<itsNTerms; tix++)
1054 : {
1055 0 : SIImageStore::restore(rbeam, usebeam, tix, psfcutoff);
1056 : }
1057 :
1058 0 : calculateAlphaBeta("image");
1059 :
1060 0 : }// restore
1061 :
1062 0 : void SIImageStoreMultiTerm::calculateAlphaBeta(String imtype)
1063 : {
1064 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","calculateAlphaBeta",WHERE) );
1065 : try
1066 : {
1067 :
1068 : // Check if this is Stokes I.
1069 0 : Bool isOnlyStokesI=True;
1070 :
1071 0 : Vector<Int> stokes ( (itsParentCoordSys.stokesCoordinate()).stokes() );
1072 0 : AlwaysAssert( stokes.nelements()>0 , AipsError);
1073 0 : if( stokes.nelements()>1 || stokes[0]!=1 ) isOnlyStokesI=False;
1074 :
1075 0 : if( ! isOnlyStokesI )
1076 0 : { os << "Alpha and Beta images will not be computed for images that contain anything other than stokes I in them." << LogIO::POST; }
1077 :
1078 : //Calculate alpha, beta
1079 0 : if( itsNTerms > 1 && isOnlyStokesI)
1080 : {
1081 : // Calculate alpha and beta
1082 0 : LatticeExprNode leMaxRes = max( *( residual(0) ) );
1083 0 : Float maxres = leMaxRes.getFloat();
1084 0 : Float specthreshold = maxres/10.0; //////////// do something better here.....
1085 :
1086 0 : os << "Calculating spectral parameters for Intensity > peakresidual/10 = " << specthreshold << " Jy/beam" << LogIO::POST;
1087 :
1088 : /////////////////////////////////////////////////////////
1089 : /////// Calculate alpha
1090 0 : if(imtype=="image")
1091 : {
1092 0 : LatticeExpr<Float> mask1(iif(((*(image(0))))>(specthreshold),1.0,0.0));
1093 0 : LatticeExpr<Float> mask0(iif(((*(image(0))))>(specthreshold),0.0,1.0));
1094 0 : LatticeExpr<Float> alphacalc( (((*(image(1))))*mask1)/(((*(image(0))))+(mask0)) );
1095 0 : alpha()->copyData(alphacalc);
1096 :
1097 0 : ImageInfo ii = image(0)->imageInfo();
1098 : // Set the restoring beam for alpha
1099 0 : alpha()->setImageInfo(ii);
1100 : //alpha()->table().unmarkForDelete();
1101 :
1102 : // Make a mask for the alpha image
1103 0 : LatticeExpr<Bool> lemask(iif(((*(image(0))) > specthreshold) , True, False));
1104 0 : removeMask( alpha() );
1105 0 : createMask( lemask, (alpha()) );
1106 :
1107 : /////// Calculate alpha error
1108 0 : alphaerror()->set(0.0);
1109 :
1110 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) ) ) );
1111 0 : alphaerror()->copyData(alphacalcerror);
1112 0 : alphaerror()->setImageInfo(ii);
1113 0 : removeMask(alphaerror());
1114 0 : createMask(lemask, alphaerror());
1115 : // alphaerror()->table().unmarkForDelete();
1116 0 : os << "Written Spectral Index Error Image : " << alphaerror()->name() << LogIO::POST;
1117 :
1118 0 : if(itsNTerms>2) // calculate beta too.
1119 : {
1120 0 : beta()->set(0.0);
1121 0 : LatticeExpr<Float> betacalc( (*image(2)*mask1)/((*image(0))+(mask0))-0.5*(*alpha())*((*alpha())-1.0) );
1122 0 : beta()->copyData(betacalc);
1123 0 : beta()->setImageInfo(ii);
1124 : //imbeta.setUnits(Unit("Spectral Curvature"));
1125 0 : removeMask(beta());
1126 0 : createMask(lemask, beta());
1127 0 : os << "Written Spectral Curvature Image : " << beta()->name() << LogIO::POST;
1128 0 : }
1129 :
1130 0 : }
1131 : else
1132 : {
1133 0 : LatticeExpr<Float> mask1(iif(((*(imagepbcor(0))))>(specthreshold),1.0,0.0));
1134 0 : LatticeExpr<Float> mask0(iif(((*(imagepbcor(0))))>(specthreshold),0.0,1.0));
1135 0 : LatticeExpr<Float> alphacalc( (((*(imagepbcor(1))))*mask1)/(((*(imagepbcor(0))))+(mask0)) );
1136 0 : alphapbcor()->copyData(alphacalc);
1137 :
1138 0 : ImageInfo ii = image(0)->imageInfo();
1139 : // Set the restoring beam for alpha
1140 0 : alphapbcor()->setImageInfo(ii);
1141 : //alpha()->table().unmarkForDelete();
1142 :
1143 : // Make a mask for the alpha image
1144 0 : LatticeExpr<Bool> lemask(iif(((*(imagepbcor(0))) > specthreshold) , True, False));
1145 0 : removeMask( alphapbcor() );
1146 0 : createMask( lemask, (alphapbcor()) );
1147 :
1148 : /////////////////////////////////////////////////////////
1149 0 : if(itsNTerms>2) // calculate beta too.
1150 : {
1151 0 : betapbcor()->set(0.0);
1152 0 : LatticeExpr<Float> betacalc( (*imagepbcor(2)*mask1)/((*imagepbcor(0))+(mask0))-0.5*(*alphapbcor())*((*alphapbcor())-1.0) );
1153 0 : betapbcor()->copyData(betacalc);
1154 0 : betapbcor()->setImageInfo(ii);
1155 : //imbeta.setUnits(Unit("Spectral Curvature"));
1156 0 : removeMask(betapbcor());
1157 0 : createMask(lemask, betapbcor());
1158 0 : os << "Written Spectral Curvature Image : " << betapbcor()->name() << LogIO::POST;
1159 0 : }
1160 :
1161 0 : }// pbcor
1162 :
1163 0 : }//if nterms>1
1164 0 : }
1165 0 : catch(AipsError &x)
1166 : {
1167 0 : throw( AipsError("Error in computing Alpha (and Beta) : " + x.getMesg() ) );
1168 0 : }
1169 :
1170 0 : }// calculateAlphaBeta
1171 :
1172 0 : void SIImageStoreMultiTerm::pbcor()
1173 : {
1174 :
1175 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","pbcorPlane",WHERE) );
1176 :
1177 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;
1178 :
1179 0 : for(uInt tix=0; tix<itsNTerms; tix++)
1180 : {
1181 0 : SIImageStore::pbcor(tix);
1182 : }
1183 :
1184 : //calculateAlphaBeta("pbcor");
1185 :
1186 0 : }
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 0 : }
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 0 : }
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 0 : }
1270 0 : }
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 0 : }
1297 : }
1298 0 : return fbw;
1299 0 : }
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 0 : }
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 :
|
