Line data Source code
1 : //# SIImageStore.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 : #include <casacore/casa/OS/DirectoryIterator.h>
42 : #include <casacore/casa/OS/File.h>
43 : #include <casacore/casa/OS/Path.h>
44 : #include <casacore/lattices/Lattices/LatticeLocker.h>
45 : #include <casacore/casa/OS/HostInfo.h>
46 : #include <components/ComponentModels/GaussianDeconvolver.h>
47 : #include <casacore/images/Images/TempImage.h>
48 : #include <casacore/images/Images/PagedImage.h>
49 : #include <imageanalysis/ImageAnalysis/CasaImageBeamSet.h>
50 : #include <casacore/lattices/LatticeMath/LatticeMathUtil.h>
51 : #include <casacore/ms/MeasurementSets/MSHistoryHandler.h>
52 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
53 : #include <casacore/tables/Tables/TableUtil.h>
54 :
55 : #include <synthesis/ImagerObjects/SIImageStore.h>
56 : #include <synthesis/ImagerObjects/SDMaskHandler.h>
57 : #include <synthesis/TransformMachines/StokesImageUtil.h>
58 : #include <synthesis/TransformMachines2/Utils.h>
59 : #include <synthesis/ImagerObjects/SynthesisUtilMethods.h>
60 : #include <casacore/images/Images/ImageRegrid.h>
61 : #include <imageanalysis/ImageAnalysis/ImageStatsCalculator.h>
62 :
63 : //#include <imageanalysis/ImageAnalysis/ImageMaskedPixelReplacer.h>
64 :
65 : #include <sys/types.h>
66 : #include <unistd.h>
67 : using namespace std;
68 :
69 : using namespace casacore;
70 :
71 : namespace casa { //# NAMESPACE CASA - BEGIN
72 :
73 : //
74 : //===========================================================================
75 : // Global method that I (SB) could make work in SynthesisUtilsMethods.
76 : //
77 : template <class T>
78 0 : void openImage(const String& imagenamefull,std::shared_ptr<ImageInterface<T> >& imPtr )
79 : {
80 0 : LogIO logIO ( LogOrigin("SynthesisImager","openImage(name)") );
81 : try
82 : {
83 0 : if (Table::isReadable(imagenamefull))
84 0 : imPtr.reset( new PagedImage<T>( imagenamefull ) );
85 : }
86 0 : catch (AipsError &x)
87 : {
88 0 : logIO << "Error in reading image \"" << imagenamefull << "\"" << LogIO::EXCEPTION;
89 : }
90 0 : }
91 : //
92 : //--------------------------------------------------------------
93 : //
94 : template
95 : void openImage(const String& imagenamefull, std::shared_ptr<ImageInterface<Float> >& img );
96 : template
97 : void openImage(const String& imagenamefull, std::shared_ptr<ImageInterface<Complex> >& img );
98 : //
99 : //===========================================================================
100 :
101 :
102 : //////////////////////////////////////////////////////////////////////////////////////////////////////
103 : //////////////////////////////////////////////////////////////////////////////////////////////////////
104 : ///// START of SIIMAGESTORE implementation
105 : //////////////////////////////////////////////////////////////////////////////////////////////////////
106 :
107 : //////////////////////////////////////////////////////////////////////////////////////////////////////
108 :
109 4729 : SIImageStore::SIImageStore()
110 : {
111 :
112 4729 : itsPsf.reset( );
113 4729 : itsModel.reset( );
114 4729 : itsResidual.reset( );
115 4729 : itsWeight.reset( );
116 4729 : itsImage.reset( );
117 4729 : itsMask.reset( );
118 4729 : itsGridWt.reset( );
119 4729 : itsPB.reset( );
120 4729 : itsImagePBcor.reset();
121 4729 : itsTempWorkIm.reset();
122 :
123 4729 : itsSumWt.reset( );
124 4729 : itsOverWrite=False;
125 4729 : itsUseWeight=False;
126 4729 : itsPBScaleFactor=1.0;
127 :
128 4729 : itsNFacets=1;
129 4729 : itsFacetId=0;
130 4729 : itsNChanChunks = 1;
131 4729 : itsChanId = 0;
132 4729 : itsNPolChunks = 1;
133 4729 : itsPolId = 0;
134 :
135 4729 : itsImageShape=IPosition(4,0,0,0,0);
136 4729 : itsImageName=String("");
137 4729 : itsCoordSys=CoordinateSystem();
138 4729 : itsMiscInfo=Record();
139 4729 : init();
140 :
141 :
142 : // validate();
143 :
144 4729 : }
145 :
146 : // Used from SynthesisNormalizer::makeImageStore()
147 766 : SIImageStore::SIImageStore(const String &imagename,
148 : const CoordinateSystem &imcoordsys,
149 : const IPosition &imshape,
150 : const String &objectname,
151 : const Record &miscinfo,
152 : // const Int nfacets,
153 : const Bool /*overwrite*/,
154 766 : const Bool useweightimage)
155 : // TODO : Add parameter to indicate weight image shape.
156 : {
157 1532 : LogIO os( LogOrigin("SIImageStore","Open new Images",WHERE) );
158 :
159 :
160 766 : itsPsf.reset( );
161 766 : itsModel.reset( );
162 766 : itsResidual.reset( );
163 766 : itsWeight.reset( );
164 766 : itsImage.reset( );
165 766 : itsMask.reset( );
166 766 : itsGridWt.reset( );
167 766 : itsPB.reset( );
168 766 : itsImagePBcor.reset( );
169 766 : itsTempWorkIm.reset();
170 :
171 766 : itsSumWt.reset( );
172 766 : itsOverWrite=False; // Hard Coding this. See CAS-6937. overwrite;
173 766 : itsUseWeight=useweightimage;
174 766 : itsPBScaleFactor=1.0;
175 :
176 766 : itsNFacets=1;
177 766 : itsFacetId=0;
178 766 : itsNChanChunks = 1;
179 766 : itsChanId = 0;
180 766 : itsNPolChunks = 1;
181 766 : itsPolId = 0;
182 :
183 766 : itsImageName = imagename;
184 766 : itsCoordSys = imcoordsys;
185 766 : itsImageShape = imshape;
186 766 : itsObjectName = objectname;
187 766 : itsMiscInfo = miscinfo;
188 :
189 766 : init();
190 :
191 766 : validate();
192 766 : }
193 :
194 : // Used from SynthesisNormalizer::makeImageStore()
195 3001 : SIImageStore::SIImageStore(const String &imagename, const Bool ignorefacets, const Bool noRequireSumwt)
196 : {
197 6002 : LogIO os( LogOrigin("SIImageStore","Open existing Images",WHERE) );
198 :
199 :
200 3001 : itsPsf.reset( );
201 3001 : itsModel.reset( );
202 3001 : itsResidual.reset( );
203 3001 : itsWeight.reset( );
204 3001 : itsImage.reset( );
205 3001 : itsMask.reset( );
206 3001 : itsGridWt.reset( );
207 3001 : itsPB.reset( );
208 3001 : itsImagePBcor.reset( );
209 3001 : itsTempWorkIm.reset();
210 3001 : itsMiscInfo=Record();
211 :
212 3001 : itsSumWt.reset( );
213 3001 : itsNFacets=1;
214 3001 : itsFacetId=0;
215 3001 : itsNChanChunks = 1;
216 3001 : itsChanId = 0;
217 3001 : itsNPolChunks = 1;
218 3001 : itsPolId = 0;
219 :
220 3001 : itsOverWrite=False;
221 : //need to to this init now so that imageExts is initialized
222 3001 : init();
223 :
224 3001 : itsImageName = imagename;
225 :
226 : // Since this constructor creates an ImStore from images on disk, it needs at least one of the
227 : // images to actually be present on disk, from which it can retrieve shape and coordsys information.
228 :
229 6002 : if( doesImageExist(itsImageName+String(".residual")) ||
230 4333 : doesImageExist(itsImageName+String(".psf")) ||
231 3003 : doesImageExist(itsImageName+String(".model")) ||
232 3001 : doesImageExist(itsImageName+String(".gridwt")) ||
233 10335 : doesImageExist(itsImageName+String(".pb")) ||
234 3001 : doesImageExist(itsImageName+String(".weight"))
235 : )
236 : {
237 3001 : std::shared_ptr<ImageInterface<Float> > imptr;
238 3001 : if( doesImageExist(itsImageName+String(".psf")) )
239 : {
240 2862 : buildImage( imptr, (itsImageName+String(".psf")) );
241 : // itsObjectName=imptr->imageInfo().objectName();
242 : // itsMiscInfo=imptr->miscInfo();
243 : }
244 139 : else if ( doesImageExist(itsImageName+String(".residual")) ){
245 137 : buildImage( imptr, (itsImageName+String(".residual")) );
246 : // itsObjectName=imptr->imageInfo().objectName();
247 : // itsMiscInfo=imptr->miscInfo();
248 : }
249 2 : else if ( doesImageExist(itsImageName+String(".model")) ){
250 2 : buildImage( imptr, (itsImageName+String(".model")) );
251 : // itsObjectName=imptr->imageInfo().objectName();
252 : // itsMiscInfo=imptr->miscInfo();
253 : }
254 0 : else if ( doesImageExist(itsImageName+String(".pb")) ){
255 0 : buildImage( imptr, (itsImageName+String(".pb")) );
256 : // itsObjectName=imptr->imageInfo().objectName();
257 : // itsMiscInfo=imptr->miscInfo();
258 : }
259 0 : else if ( doesImageExist(itsImageName+String(".weight")) ){
260 0 : buildImage( imptr, (itsImageName+String(".weight")) );
261 : // itsObjectName=imptr->imageInfo().objectName();
262 : // itsMiscInfo=imptr->miscInfo();
263 : }
264 : else
265 : {
266 : // How can this be right ?
267 0 : buildImage( imptr, (itsImageName+String(".gridwt")) );
268 : }
269 :
270 3001 : itsObjectName=imptr->imageInfo().objectName();
271 3001 : itsImageShape=imptr->shape();
272 3001 : itsCoordSys = imptr->coordinates();
273 3001 : itsMiscInfo=imptr->miscInfo();
274 :
275 3001 : }
276 : else
277 : {
278 0 : throw( AipsError( "PSF, Residual, Model Image (or sumwt) do not exist. Please create one of them." ) );
279 : }
280 :
281 7334 : if( doesImageExist(itsImageName+String(".residual")) ||
282 4333 : doesImageExist(itsImageName+String(".psf")) )
283 : {
284 2999 : if( doesImageExist(itsImageName+String(".sumwt")) )
285 : {
286 2914 : std::shared_ptr<ImageInterface<Float> > imptr;
287 : //imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt")) );
288 2914 : buildImage( imptr, (itsImageName+String(".sumwt")) );
289 2914 : itsNFacets = imptr->shape()[0];
290 2914 : itsFacetId = 0;
291 2914 : itsUseWeight = getUseWeightImage( *imptr );
292 2914 : itsPBScaleFactor=1.0; ///// No need to set properly here as it will be calc'd in ()
293 : /////redo this here as psf may have different coordinates
294 2914 : itsCoordSys = imptr->coordinates();
295 2914 : itsMiscInfo=imptr->miscInfo();
296 2914 : if( itsUseWeight && ! doesImageExist(itsImageName+String(".weight")) )
297 : {
298 0 : throw(AipsError("Internal error : Sumwt has a useweightimage=True but the weight image does not exist."));
299 : }
300 2914 : }
301 : else
302 : {
303 85 : if(!noRequireSumwt) // .sumwt image required? -> probably not for just the minor cycle (aka task deconvolve)
304 0 : {throw( AipsError( "SumWt information does not exist. Please create either a PSF or Residual" ) );}
305 : else
306 : {
307 85 : os << "SumWt does not exist. Proceeding only with PSF" << LogIO::POST;
308 85 : std::shared_ptr<ImageInterface<Float> > imptr;
309 : //imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt")) );
310 85 : if( doesImageExist(itsImageName+String(".residual") ) )
311 84 : { buildImage( imptr, (itsImageName+String(".residual")) ); }
312 : else
313 1 : { buildImage( imptr, (itsImageName+String(".psf")) ); }
314 :
315 85 : itsNFacets=1;
316 85 : itsFacetId=0;
317 85 : itsUseWeight=False;
318 85 : itsPBScaleFactor=1.0;
319 85 : itsCoordSys = imptr->coordinates();
320 85 : itsMiscInfo=imptr->miscInfo();
321 85 : }
322 : }
323 : }// if psf or residual exist...
324 :
325 3001 : if( ignorefacets==True ) itsNFacets= 1;
326 :
327 3001 : init();
328 :
329 3001 : validate();
330 3001 : }
331 :
332 : // used from getSubImageStore(), for example when creating the facets list
333 : // this would be safer if it was refactored as a copy constructor of the generic stuff +
334 : // initialization of the facet related parameters
335 751 : SIImageStore::SIImageStore(const std::shared_ptr<ImageInterface<Float> > &modelim,
336 : const std::shared_ptr<ImageInterface<Float> > &residim,
337 : const std::shared_ptr<ImageInterface<Float> > &psfim,
338 : const std::shared_ptr<ImageInterface<Float> > &weightim,
339 : const std::shared_ptr<ImageInterface<Float> > &restoredim,
340 : const std::shared_ptr<ImageInterface<Float> > &maskim,
341 : const std::shared_ptr<ImageInterface<Float> > &sumwtim,
342 : const std::shared_ptr<ImageInterface<Float> > &gridwtim,
343 : const std::shared_ptr<ImageInterface<Float> > &pbim,
344 : const std::shared_ptr<ImageInterface<Float> > &restoredpbcorim,
345 : const std::shared_ptr<ImageInterface<Float> > & tempworkimage,
346 : const CoordinateSystem &csys,
347 : const IPosition &imshape,
348 : const String &imagename,
349 : const String &objectname,
350 : const Record &miscinfo,
351 : const Int facet, const Int nfacets,
352 : const Int chan, const Int nchanchunks,
353 : const Int pol, const Int npolchunks,
354 751 : const Bool useweightimage)
355 : {
356 :
357 :
358 751 : itsPsf=psfim;
359 751 : itsModel=modelim;
360 751 : itsResidual=residim;
361 : /* if(residim){
362 : LatticeLocker lock1(*itsResidual, FileLocker::Read);
363 : cerr << "RESIDMAX " << max(itsResidual->get()) << " " << max(residim->get()) << endl;
364 : }*/
365 751 : itsWeight=weightim;
366 751 : itsImage=restoredim;
367 751 : itsMask=maskim;
368 :
369 751 : itsSumWt=sumwtim;
370 :
371 751 : itsGridWt=gridwtim;
372 751 : itsPB=pbim;
373 751 : itsImagePBcor=restoredpbcorim;
374 751 : itsTempWorkIm=tempworkimage;
375 :
376 751 : itsPBScaleFactor=1.0;// No need to set properly here as it will be computed in makePSF.
377 :
378 751 : itsObjectName = objectname;
379 751 : itsMiscInfo = miscinfo;
380 :
381 751 : itsNFacets = nfacets;
382 751 : itsFacetId = facet;
383 751 : itsNChanChunks = nchanchunks;
384 751 : itsChanId = chan;
385 751 : itsNPolChunks = npolchunks;
386 751 : itsPolId = pol;
387 :
388 751 : itsOverWrite=False;
389 751 : itsUseWeight=useweightimage;
390 :
391 751 : itsParentImageShape = imshape;
392 751 : itsImageShape = imshape;
393 :
394 751 : itsParentCoordSys = csys;
395 751 : itsCoordSys = csys; // Hopefully this doesn't change for a reference image
396 751 : itsImageName = imagename;
397 :
398 : //-----------------------
399 751 : init(); // Connect parent pointers to the images.
400 : //-----------------------
401 :
402 : // Set these to null, to be set later upon first access.
403 751 : itsPsf.reset( );
404 751 : itsModel.reset( );
405 751 : itsResidual.reset( );
406 751 : itsWeight.reset( );
407 751 : itsImage.reset( );
408 751 : itsMask.reset( );
409 751 : itsSumWt.reset( );
410 751 : itsPB.reset( );
411 :
412 751 : validate();
413 751 : }
414 :
415 8963 : void SIImageStore::validate()
416 : {
417 : /// There are two valid states. Check for at least one of them.
418 8963 : Bool state = False;
419 :
420 8963 : stringstream oss;
421 8963 : oss << "shape:" << itsImageShape << " parentimageshape:" << itsParentImageShape
422 8963 : << " nfacets:" << itsNFacets << "x" << itsNFacets << " facetid:" << itsFacetId
423 8963 : << " nchanchunks:" << itsNChanChunks << " chanid:" << itsChanId
424 8963 : << " npolchunks:" << itsNPolChunks << " polid:" << itsPolId
425 8963 : << " coord-dim:" << itsCoordSys.nPixelAxes()
426 8963 : << " psf/res:" << (hasPsf() || hasResidual()) ;
427 8963 : if( hasSumWt() ) oss << " sumwtshape : " << sumwt()->shape() ;
428 8963 : oss << endl;
429 :
430 :
431 : try {
432 :
433 8963 : if( itsCoordSys.nPixelAxes() != 4 ) state=False;
434 :
435 : /// (1) Regular imagestore
436 8963 : if( itsNFacets==1 && itsFacetId==0
437 8771 : && itsNChanChunks==1 && itsChanId==0
438 8771 : && itsNPolChunks==1 && itsPolId==0 ) {
439 8691 : Bool check1 = hasSumWt() && sumwt()->shape()[0]==1;
440 17382 : if( (itsImageShape.isEqual(itsParentImageShape) ) && ( check1 || !hasSumWt() )
441 17382 : && itsParentImageShape.product() > 0 ) state=True;
442 8691 : }
443 : /// (2) Reference Sub Imagestore
444 272 : else if ( ( itsNFacets>1 && itsFacetId >=0 )
445 80 : || ( itsNChanChunks>1 && itsChanId >=0 )
446 80 : || ( itsNPolChunks>1 && itsPolId >=0 ) ) {
447 : // If shape is still unset, even when the first image has been made....
448 272 : Bool check1 = ( itsImageShape.product() > 0 && ( hasPsf() || hasResidual() ) );
449 272 : Bool check2 = ( itsImageShape.isEqual(IPosition(4,0,0,0,0)) && ( !hasPsf() && !hasResidual() ) );
450 272 : Bool check3 = hasSumWt() && sumwt()->shape()[0]==1; // One facet only.
451 :
452 272 : if( ( check1 || check2 ) && ( itsParentImageShape.product()>0 )
453 272 : && ( itsFacetId < itsNFacets*itsNFacets )
454 : // && ( itsChanId <= itsNChanChunks ) // chanchunks can be larger...
455 272 : && ( itsPolId < itsNPolChunks )
456 544 : && ( check3 || !hasSumWt() ) ) state = True;
457 : }
458 :
459 0 : } catch( AipsError &x ) {
460 0 : state = False;
461 0 : oss << " ||||| " << x.getMesg() << endl;
462 0 : }
463 :
464 : // cout << " SIIM:validate : " << oss.str() << endl;
465 :
466 8963 : if( state == False ) throw(AipsError("Internal Error : Invalid ImageStore state : "+ oss.str()) );
467 :
468 17926 : return;
469 8963 : }
470 :
471 :
472 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
473 :
474 751 : std::shared_ptr<SIImageStore> SIImageStore::getSubImageStore(const Int facet, const Int nfacets,
475 : const Int chan, const Int nchanchunks,
476 : const Int pol, const Int npolchunks)
477 : {
478 :
479 751 : return std::make_shared<SIImageStore>(itsModel, itsResidual, itsPsf, itsWeight, itsImage, itsMask, itsSumWt, itsGridWt, itsPB, itsImagePBcor, itsTempWorkIm, itsCoordSys, itsImageShape, itsImageName, itsObjectName, itsMiscInfo, facet, nfacets, chan, nchanchunks, pol, npolchunks, itsUseWeight);
480 : }
481 :
482 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
483 : //// Either read an image from disk, or construct one.
484 :
485 31763 : std::shared_ptr<ImageInterface<Float> > SIImageStore::openImage(const String imagenamefull,
486 : const Bool overwrite,
487 : const Bool dosumwt, const Int nfacetsperside, const Bool checkCoordSys)
488 : {
489 :
490 :
491 31763 : std::shared_ptr<ImageInterface<Float> > imPtr;
492 31763 : IPosition useShape( itsParentImageShape );
493 :
494 31763 : if( dosumwt ) // change shape to sumwt image shape.
495 : {
496 4340 : useShape[0] = nfacetsperside;
497 4340 : useShape[1] = nfacetsperside;
498 : // cout << "openImage : Making sumwt grid : using shape : " << useShape << endl;
499 : }
500 :
501 : // overwrite=False; /// HARD CODING THIS. See CAS-6937.
502 :
503 : // cout << "Open image : " << imagenamefull << " useShape : " << useShape << endl;
504 :
505 : // if image exists
506 : // if overwrite=T
507 : // try to make new image
508 : // if not, try to open existing image
509 : // if cannot, complain.
510 : // if overwrite=F
511 : // try to open existing image
512 : // if cannot, complain
513 : // if image does not exist
514 : // try to make new image
515 : // if cannot, complain
516 31763 : Bool dbg=False;
517 31763 : if( doesImageExist( imagenamefull ) )
518 : {
519 : ///// not used since itsOverWrite is hardcoded to FALSE (CAS-6937)
520 25656 : if (overwrite) //overwrite and make new image
521 : {
522 0 : if(dbg) cout << "Trying to overwrite and open new image named : " << imagenamefull << " ow:"<< overwrite << endl;
523 : try{
524 0 : buildImage(imPtr, useShape, itsParentCoordSys, imagenamefull) ;
525 0 : LatticeLocker lock1(*imPtr, FileLocker::Write);
526 : // initialize to zeros...
527 0 : imPtr->set(0.0);
528 0 : }
529 0 : catch (AipsError &x){
530 0 : if(dbg)cout << "Cannot overwrite : " << x.getMesg() << endl;
531 0 : if(dbg)cout << "Open already ? : " << Table::isOpened( imagenamefull ) << " Writable ? : " << Table::isWritable( imagenamefull ) << endl;
532 0 : if(Table::isWritable( imagenamefull ))
533 : {
534 0 : if(dbg) cout << "--- Trying to open existing image : "<< imagenamefull << endl;
535 : try{
536 0 : buildImage( imPtr, imagenamefull );
537 : }
538 0 : catch (AipsError &x){
539 0 : throw( AipsError("Writable table exists, but cannot open : " + x.getMesg() ) );
540 0 : }
541 : }// is table writable
542 : else
543 : {
544 0 : throw( AipsError("Cannot overwrite existing image. : " + x.getMesg() ) );
545 : }
546 0 : }
547 : }// overwrite existing image
548 : else // open existing image ( Always tries this )
549 : {
550 : if(1) //Table::isWritable( imagenamefull ))
551 : {
552 25656 : if(dbg) cout << "Trying to open existing image : "<< imagenamefull << endl;
553 : try{
554 25656 : buildImage( imPtr, imagenamefull ) ;
555 :
556 25656 : if( !dosumwt)
557 : {
558 : //cout << useShape << " and " << imPtr->shape() << " ---- " << useShape.product() << " : " << itsCoordSys.nCoordinates() << endl;
559 : // Check if coordsys and shape of this image are consistent with current ones (if filled)
560 :
561 22279 : if( useShape.product()>0 && ! useShape.isEqual( imPtr->shape() ) )
562 : {
563 15 : ostringstream oo1,oo2;
564 15 : oo1 << useShape; oo2 << imPtr->shape();
565 15 : throw( AipsError( "There is a shape mismatch between existing images ("+oo2.str()+") and current parameters ("+oo1.str()+"). If you are attempting to restart a run with a new image shape, please change imagename and supply the old model or mask as inputs (via the startmodel or mask parameters) so that they can be regridded to the new shape before continuing." ) );
566 30 : }
567 :
568 :
569 :
570 22264 : if( itsParentCoordSys.nCoordinates()>0 && checkCoordSys && ! itsParentCoordSys.near( imPtr->coordinates()) )
571 : {
572 : /// Implement an exception to get CAS-9977 to work. Modify the current coordsys to match the parent.
573 : /// "The DirectionCoordinates have differing latpoles"
574 2 : if( itsParentCoordSys.errorMessage().contains("differing latpoles") )
575 : {
576 0 : DirectionCoordinate dcord = itsParentCoordSys.directionCoordinate();
577 : //cout << "Latpole from parent : " << dcord.longLatPoles() << endl;
578 :
579 0 : DirectionCoordinate dcord2 = imPtr->coordinates().directionCoordinate();
580 : //cout << "Latpole from imPtr : " << dcord2.longLatPoles() << endl;
581 :
582 0 : LogIO os( LogOrigin("SIImageStore","Open existing Images",WHERE) );
583 0 : os << " Mismatch in Csys latpoles between existing image on disk (" << dcord2.longLatPoles() << ") and current imaging run (" << dcord.longLatPoles() << ") : " << itsParentCoordSys.errorMessage() << " -- Resetting to match image on disk" << LogIO::WARN << LogIO::POST;
584 : // cout << "Set the coordinate from the Parent value" << endl;
585 : //imPtr->setCoordinateInfo( itsParentCoordSys );
586 0 : itsParentCoordSys = imPtr->coordinates();
587 0 : }
588 :
589 : // Check again.
590 2 : if ( ! itsParentCoordSys.near( imPtr->coordinates() ) )
591 : {
592 : //cout << " Still different..." << endl;
593 2 : throw( AipsError( "There is a coordinate system mismatch between existing images on disk and current parameters ("+itsParentCoordSys.errorMessage()+"). If you are attempting to restart a run, please change imagename and supply the old model or mask as inputs (via the startmodel or mask parameters) so that they can be regridded to the new coordinate system before continuing. " ) );
594 : }
595 : //}
596 : }
597 : }// not dosumwt
598 : }
599 17 : catch (AipsError &x){
600 17 : throw( AipsError("Cannot open existing image : "+imagenamefull+" : " + x.getMesg() ) );
601 17 : }
602 : }// is table writable
603 : else // table exists but not writeable
604 : {
605 : if(dbg)cout << "Table exists but not writeable : " << imagenamefull << " --- Open : " << Table::isOpened( imagenamefull ) << endl;
606 : throw( AipsError("Image exists but not able to open for writes :"+imagenamefull+ ". Opened elsewhere : " + String::toString(Table::isOpened(imagenamefull))) );
607 : }
608 : }// open existing image
609 : }// if image exists
610 : else // image doesn't exist. make new one
611 : {
612 6107 : if(dbg) cout << "Trying to open new image named : " << imagenamefull << endl;
613 : try{
614 6107 : buildImage(imPtr, useShape, itsParentCoordSys, imagenamefull) ;
615 : // initialize to zeros...
616 : {
617 6107 : LatticeLocker lock1(*imPtr, FileLocker::Write);
618 6107 : imPtr->set(0.0);
619 6107 : imPtr->flush();
620 6107 : imPtr->unlock();
621 6107 : }
622 : }
623 0 : catch (AipsError &x){
624 0 : throw( AipsError("Cannot make new image. : " + x.getMesg() ) );
625 0 : }
626 : }
627 :
628 :
629 : //////////////////////////////////////
630 : /*
631 : if( overwrite || !Table::isWritable( imagenamefull ) )
632 : {
633 : cout << "Trying to open new image named : " << imagenamefull << " ow:"<< overwrite << endl;
634 : imPtr.reset( new PagedImage<Float> (useShape, itsCoordSys, imagenamefull) );
635 : // initialize to zeros...
636 : imPtr->set(0.0);
637 : }
638 : else
639 : {
640 : if(Table::isWritable( imagenamefull ))
641 : {
642 : cout << "Trying to open existing image : "<< imagenamefull << endl;
643 : try{
644 : imPtr.reset( new PagedImage<Float>( imagenamefull ) );
645 : }
646 : catch (AipsError &x){
647 : cerr << "Writable table exists, but cannot open. Creating temp image. : " << x.getMesg() << endl;
648 : imPtr.reset( new TempImage<Float> (useShape, itsCoordSys) );
649 : // imPtr=new PagedImage<Float> (useShape, itsCoordSys, imagenamefull);
650 : imPtr->set(0.0);
651 : }
652 : }
653 : else
654 : {
655 : cerr << "Table " << imagenamefull << " is not writeable. Creating temp image." << endl;
656 : imPtr.reset( new TempImage<Float> (useShape, itsCoordSys) );
657 : imPtr->set(0.0);
658 : }
659 : }
660 : */
661 63492 : return imPtr;
662 31780 : }
663 :
664 6107 : void SIImageStore::buildImage(std::shared_ptr<ImageInterface<Float> > &imptr, IPosition shape,
665 : CoordinateSystem csys, const String name)
666 : {
667 12214 : LogIO os( LogOrigin("SIImageStore", "Open non-existing image", WHERE) );
668 6107 : os <<"Opening image, name: " << name << LogIO::DEBUG1;
669 :
670 6107 : itsOpened++;
671 : //For some reason cannot open a new paged image with UserNoread directly
672 : {
673 6107 : PagedImage<Float> godot(shape, csys, name);
674 6107 : godot.unlock();
675 6107 : }
676 6107 : TableLock::LockOption locktype=TableLock::AutoNoReadLocking;
677 : //TableLock::LockOption locktype=TableLock::UserLocking;
678 : /*if((name.contains(imageExts(PSF)) && !name.contains(imageExts(PSF)+".tt"))|| (name.contains(imageExts(RESIDUAL))&& !name.contains(imageExts(RESIDUAL)+".tt")) || (name.contains(imageExts(SUMWT)) && !name.contains(imageExts(SUMWT)+".tt"))){
679 : locktype=TableLock::UserNoReadLocking;
680 : }*/
681 6107 : imptr.reset( new PagedImage<Float> (name, locktype) );
682 : {
683 6107 : LatticeLocker lock1(*imptr, FileLocker::Write);
684 6107 : initMetaInfo(imptr, name);
685 6107 : imptr->unlock();
686 6107 : }
687 : /*
688 : Int MEMFACTOR = 18;
689 : Double memoryMB=HostInfo::memoryTotal(True)/1024/(MEMFACTOR*itsOpened);
690 :
691 :
692 : TempImage<Float> *tptr = new TempImage( TiledShape(shape, tileShape()), csys, memoryBeforeLattice() ) ;
693 :
694 : tptr->setMaximumCacheSize(shape.product());
695 : tptr->cleanCache();
696 :
697 : imptr.reset( tptr );
698 :
699 : */
700 6107 : }
701 :
702 33113 : void SIImageStore::buildImage(std::shared_ptr<ImageInterface<Float> > &imptr, const String name)
703 : {
704 66226 : LogIO os(LogOrigin("SIImageStore", "Open existing Images", WHERE));
705 33113 : os <<"Opening image, name: " << name << LogIO::DEBUG1;
706 :
707 33113 : itsOpened++;
708 33113 : if(Table::isReadable(name)){
709 33113 : TableLock::LockOption locktype=TableLock::AutoNoReadLocking;
710 : /*if((name.contains(imageExts(PSF)) && !name.contains(imageExts(PSF)+".tt"))|| (name.contains(imageExts(RESIDUAL))&& !name.contains(imageExts(RESIDUAL)+".tt")) || (name.contains(imageExts(SUMWT)) && !name.contains(imageExts(SUMWT)+".tt"))){
711 : locktype=TableLock::UserNoReadLocking;
712 : }*/
713 33113 : Table table(name, locktype);
714 33113 : String type = table.tableInfo().type();
715 33113 : if (type != TableInfo::type(TableInfo::PAGEDIMAGE)) {
716 :
717 0 : imptr.reset( new PagedImage<Float>( table ) );
718 0 : imptr->unlock();
719 0 : return;
720 : }
721 33113 : }
722 33113 : LatticeBase* latt =ImageOpener::openImage(name);
723 33113 : if(!latt)
724 : {
725 0 : throw(AipsError("Error in opening Image : "+name));
726 : }
727 33113 : DataType dtype=latt->dataType();
728 33113 : if(dtype==TpFloat)
729 : {
730 33113 : imptr.reset(dynamic_cast<ImageInterface<Float>* >(latt));
731 : }
732 : else
733 : {
734 0 : throw AipsError( "Need image to have float values : "+name);
735 : }
736 :
737 : /*
738 : std::shared_ptr<casacore::ImageInterface<Float> > fim;
739 : std::shared_ptr<casacore::ImageInterface<Complex> > cim;
740 :
741 : std::tie(fim , cim)=ImageFactory::fromFile(name);
742 : if(fim)
743 : {
744 : imptr.reset( dynamic_cast<std::shared_ptr<casacore::ImageInterface<Float> > >(*fim) );
745 : }
746 : else
747 : {
748 : throw( AipsError("Cannot open with ImageFactory : "+name));
749 : }
750 : */
751 :
752 :
753 :
754 :
755 : /*
756 : IPosition cimageShape;
757 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
758 : whichStokes, itsDataPolRep);
759 : */
760 :
761 33113 : }
762 :
763 : /**
764 : * Sets ImageInfo and MiscInfo on an image
765 : *
766 : * @param imptr image to initialize
767 : */
768 6109 : void SIImageStore::initMetaInfo(std::shared_ptr<ImageInterface<Float> > &imptr,
769 : const String name)
770 : {
771 : // Check objectname, as one of the mandatory fields. What this is meant to check is -
772 : // has the metainfo been initialized? If not, grab info from associated PSF
773 6109 : LatticeLocker lock1(*imptr, FileLocker::Write);
774 6109 : if (not itsObjectName.empty()) {
775 6109 : ImageInfo info = imptr->imageInfo();
776 6109 : info.setObjectName(itsObjectName);
777 6109 : imptr->setImageInfo(info);
778 6109 : imptr->setMiscInfo(itsMiscInfo);
779 6109 : } else if (std::string::npos == name.find(imageExts(PSF))) {
780 0 : auto srcImg = psf(0);
781 0 : if (nullptr != srcImg) {
782 0 : imptr->setImageInfo(srcImg->imageInfo());
783 0 : imptr->setMiscInfo(srcImg->miscInfo());
784 : }
785 0 : }
786 6109 : imptr->unlock();
787 6109 : }
788 :
789 :
790 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
791 0 : void SIImageStore::setMiscInfo(const Record miscinfo)
792 : {
793 0 : itsMiscInfo = miscinfo;
794 0 : }
795 :
796 0 : void SIImageStore::setObjectName(const String name)
797 : {
798 0 : itsObjectName = name;
799 0 : }
800 :
801 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
802 :
803 34883 : std::shared_ptr<ImageInterface<Float> > SIImageStore::makeSubImage(const Int facet, const Int nfacets,
804 : const Int chan, const Int nchanchunks,
805 : const Int pol, const Int npolchunks,
806 : ImageInterface<Float>& image)
807 : {
808 : //assuming n x n facets
809 34883 : Int nx_facets=Int(sqrt(Double(nfacets)));
810 69766 : LogIO os( LogOrigin("SynthesisImager","makeFacet") );
811 : // Make the output image
812 34883 : Slicer imageSlicer;
813 :
814 : // Add checks for all dimensions..
815 34883 : if((facet>(nfacets-1))||(facet<0)) {
816 0 : os << LogIO::SEVERE << "Illegal facet " << facet << LogIO::POST;
817 0 : return std::shared_ptr<ImageInterface<Float> >();
818 : }
819 34883 : IPosition imshp=image.shape();
820 34883 : IPosition blc(imshp.nelements(), 0);
821 34883 : IPosition trc=imshp-1;
822 34883 : IPosition inc(imshp.nelements(), 1);
823 :
824 : /// Facets
825 34883 : Int facetx = facet % nx_facets;
826 34883 : Int facety = (facet - facetx) / nx_facets;
827 34883 : Int sizex = imshp(0) / nx_facets;
828 34883 : Int sizey = imshp(1) / nx_facets;
829 34883 : blc(1) = facety * sizey;
830 34883 : trc(1) = blc(1) + sizey - 1;
831 34883 : blc(0) = facetx * sizex;
832 34883 : trc(0) = blc(0) + sizex - 1;
833 :
834 : /// Pol Chunks
835 34883 : Int sizepol = imshp(2) / npolchunks;
836 34883 : blc(2) = pol * sizepol;
837 34883 : trc(2) = blc(2) + sizepol - 1;
838 :
839 : /// Chan Chunks
840 34883 : Int sizechan = imshp(3) / nchanchunks;
841 34883 : blc(3) = chan * sizechan;
842 34883 : trc(3) = blc(3) + sizechan - 1;
843 :
844 34883 : LCBox::verify(blc, trc, inc, imshp);
845 34883 : Slicer imslice(blc, trc, inc, Slicer::endIsLast);
846 :
847 : // Now create the sub image
848 34883 : std::shared_ptr<ImageInterface<Float> > referenceImage( new SubImage<Float>(image, imslice, True) );
849 : {
850 34883 : LatticeLocker lock1(*referenceImage, FileLocker::Write);
851 34883 : referenceImage->setMiscInfo(image.miscInfo());
852 34883 : referenceImage->setUnits(image.units());
853 :
854 34883 : }
855 :
856 : // cout << "Made Ref subimage of shape : " << referenceImage->shape() << endl;
857 :
858 34883 : return referenceImage;
859 :
860 34883 : }
861 :
862 :
863 :
864 : //////////////////////////////////////////////////////////////////////////////////////////////////////
865 :
866 12248 : SIImageStore::~SIImageStore()
867 : {
868 12248 : }
869 :
870 : //////////////////////////////////////////////////////////////////////////////////////////////////////
871 : //////////////////////////////////////////////////////////////////////////////////////////////////////
872 :
873 12295 : Bool SIImageStore::releaseLocks()
874 : {
875 : //LogIO os( LogOrigin("SIImageStore","releaseLocks",WHERE) );
876 :
877 : // String fname( itsImageName+String(".info") );
878 : // makePersistent( fname );
879 :
880 12295 : if( itsPsf ) releaseImage( itsPsf );
881 12295 : if( itsModel ) { releaseImage( itsModel ); }
882 12295 : if( itsResidual ) releaseImage( itsResidual );
883 12295 : if( itsImage ) releaseImage( itsImage );
884 12295 : if( itsWeight ) releaseImage( itsWeight );
885 12295 : if( itsMask ) releaseImage( itsMask );
886 12295 : if( itsSumWt ) releaseImage( itsSumWt );
887 12295 : if( itsGridWt ) releaseImage( itsGridWt );
888 12295 : if( itsPB ) releaseImage( itsPB );
889 12295 : if( itsImagePBcor ) releaseImage( itsImagePBcor );
890 :
891 12295 : return True; // do something more intelligent here.
892 : }
893 :
894 0 : Bool SIImageStore::releaseComplexGrids()
895 : {
896 : //LogIO os( LogOrigin("SIImageStore","releaseComplexGrids",WHERE) );
897 :
898 0 : if( itsForwardGrid ) releaseImage( itsForwardGrid );
899 0 : if( itsBackwardGrid ) releaseImage( itsBackwardGrid );
900 :
901 0 : return True; // do something more intelligent here.
902 : }
903 :
904 25840 : void SIImageStore::releaseImage( std::shared_ptr<ImageInterface<Float> > &im )
905 : {
906 : //cerr << this << " releaseimage " << im.get() << endl;
907 : //LogIO os( LogOrigin("SIImageStore","releaseLocks",WHERE) );
908 25840 : im->flush();
909 : //os << LogIO::WARN << "clear cache" << LogIO::POST;
910 25840 : im->clearCache();
911 : //os << LogIO::WARN << "unlock" << LogIO::POST;
912 25840 : im->unlock();
913 : //os << LogIO::WARN << "tempClose" << LogIO::POST;
914 : //im->tempClose();
915 : //os << LogIO::WARN << "NULL" << LogIO::POST;
916 25840 : im.reset(); // This was added to allow modification by modules independently
917 25840 : }
918 :
919 0 : void SIImageStore::releaseImage( std::shared_ptr<ImageInterface<Complex> > &im )
920 : {
921 0 : im->tempClose();
922 0 : }
923 :
924 : //////////////////////////////////////////////////////////////////////////////////////////////////////
925 : //////////////////////////////////////////////////////////////////////////////////////////////////////
926 0 : Vector<String> SIImageStore::getModelImageName()
927 : {
928 0 : Vector<String> mods(1);
929 0 : mods[0]=itsImageName + imageExts(MODEL);
930 0 : return mods;
931 0 : }
932 :
933 746 : Long SIImageStore::estimateRAM(){
934 : //right now this is estimated at 2MB for the 2 complex lattices;
935 746 : return Long(2000);
936 : }
937 24 : void SIImageStore::setModelImage( const Vector<String> &modelnames)
938 : {
939 48 : LogIO os( LogOrigin("SIImageStore","setModelImage",WHERE) );
940 :
941 24 : if( modelnames.nelements() > 1 )
942 0 : {throw( AipsError("Multiple starting model images are currently not supported. Please merge them before supplying as input to startmodel"));
943 : /// If needed, THIS is the place to add code to support lists of model images... perhaps regrid separately and add them up or some such thing.
944 : }
945 :
946 24 : if( modelnames.nelements()==1 ) { setModelImageOne( modelnames[0] ); }
947 24 : }
948 :
949 :
950 :
951 43 : void SIImageStore::setModelImageOne( const String &modelname , Int nterm)
952 : {
953 86 : LogIO os( LogOrigin("SIImageStore","setModelImageOne",WHERE) );
954 :
955 43 : if(modelname==String("")) return;
956 :
957 41 : Bool multiterm=False;
958 41 : if(nterm>-1)multiterm=True;
959 41 : if(nterm==-1)nterm=0;
960 :
961 41 : Directory immodel( modelname ); // +String(".model") );
962 41 : if( !immodel.exists() )
963 : {
964 0 : os << "Starting model image " << modelname << " does not exist. No initial prediction will be done" << ((multiterm)?String(" for term")+String::toString(nterm) :"") << LogIO::POST;
965 0 : return;
966 : }
967 :
968 : // master merge 2019.01.08 - leaving in the commnets for now but clean up after it is verified
969 : //SHARED_PTR<PagedImage<Float> > newmodel( new PagedImage<Float>( modelname ) ); //+String(".model") ) );
970 : //SHARED_PTR<ImageInterface<Float> > newmodel;
971 41 : std::shared_ptr<ImageInterface<Float> > newmodel;
972 41 : buildImage(newmodel, modelname);
973 : // in master
974 : //std::shared_ptr<PagedImage<Float> > newmodel( new PagedImage<Float>( modelname ) ); //+String(".model") ) );
975 :
976 41 : Bool hasMask = newmodel->isMasked(); /// || newmodel->hasPixelMask() ;
977 :
978 41 : if( hasMask )
979 : {
980 :
981 3 : os << "Input startmodel has an internal mask. Setting masked pixels to zero" << LogIO::POST;
982 :
983 : try {
984 :
985 : ////// Neat function to replace masked pixels, but it will do it in-place.
986 : ////// We need to not modify the input model on disk, but apply the mask only OTF before
987 : ////// regridding to the target image,.
988 : // ImageMaskedPixelReplacer<Float> impr( newmodel, 0, "" );
989 : // impr.replace("0", False, False );
990 :
991 :
992 6 : TempImage<Float> maskmodel( newmodel->shape(), newmodel->coordinates() );
993 3 : IPosition inshape = newmodel->shape();
994 6 : for(Int pol=0; pol<inshape[2]; pol++)
995 : {
996 6 : for(Int chan=0; chan<inshape[3]; chan++)
997 : {
998 3 : IPosition pos(4,0,0,pol,chan);
999 : std::shared_ptr<ImageInterface<Float> > subim=makeSubImage(0,1,
1000 6 : chan, inshape[3],
1001 3 : pol, inshape[2],
1002 6 : (*newmodel) );
1003 :
1004 : std::shared_ptr<ImageInterface<Float> > submaskmodel=makeSubImage(0,1,
1005 6 : chan, inshape[3],
1006 3 : pol, inshape[2],
1007 3 : maskmodel );
1008 :
1009 3 : ArrayLattice<Bool> pixmask( subim->getMask() );
1010 6 : LatticeExpr<Float> masked( (*subim) * iif(pixmask,1.0,0.0) );
1011 3 : submaskmodel->copyData( masked );
1012 3 : }
1013 : }
1014 :
1015 :
1016 :
1017 : // Check shapes, coordsys with those of other images. If different, try to re-grid here.
1018 3 : if( (newmodel->shape() != model(nterm)->shape()) || (! itsCoordSys.near(newmodel->coordinates() )) )
1019 : {
1020 0 : os << "Regridding input model " << modelname << " to target coordinate system for " << itsImageName << ".model" << ((multiterm)?".tt"+String::toString(nterm) :"") << LogIO::POST;
1021 0 : regridToModelImage( maskmodel, nterm );
1022 : }
1023 : else
1024 : {
1025 3 : os << "Copying input model " << modelname << " to " << itsImageName << ".model" << ((multiterm)?".tt"+String::toString(nterm) :"") << LogIO::POST;
1026 3 : model(nterm)->copyData( LatticeExpr<Float> (maskmodel) );
1027 : }
1028 :
1029 :
1030 3 : } catch (AipsError &x) {
1031 0 : throw(AipsError("Setting masked pixels to zero for input startmodel : "+ x.getMesg()));
1032 0 : }
1033 :
1034 : }// hasMask
1035 : else // nomask
1036 : {
1037 :
1038 : // Check shapes, coordsys with those of other images. If different, try to re-grid here.
1039 38 : if( (newmodel->shape() != model(nterm)->shape()) || (! itsCoordSys.near(newmodel->coordinates() )) )
1040 : {
1041 17 : os << "Regridding input model " << modelname << " to target coordinate system for " << itsImageName << ".model" << ((multiterm)?".tt"+String::toString(nterm) :"") << LogIO::POST;
1042 17 : regridToModelImage( *newmodel, nterm );
1043 : }
1044 : else
1045 : {
1046 21 : os << "Copying input model " << modelname << " to " << itsImageName << ".model" << ((multiterm)?".tt"+String::toString(nterm) :"") << LogIO::POST;
1047 21 : model(nterm)->copyData( LatticeExpr<Float> (*newmodel) );
1048 : }
1049 : }//nomask
1050 47 : }
1051 :
1052 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1053 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1054 :
1055 9556 : IPosition SIImageStore::getShape()
1056 : {
1057 9556 : return itsImageShape;
1058 : }
1059 :
1060 6584 : String SIImageStore::getName()
1061 : {
1062 6584 : return itsImageName;
1063 : }
1064 :
1065 8367 : uInt SIImageStore::getNTaylorTerms(Bool /*dopsf*/)
1066 : {
1067 8367 : return 1;
1068 : }
1069 :
1070 : /*
1071 : void SIImageStore::checkRef( std::shared_ptr<ImageInterface<Float> > ptr, const String label )
1072 : {
1073 : if( ! ptr && itsImageName==String("reference") )
1074 : {throw(AipsError("Internal Error : Attempt to access null subImageStore "+label + " by reference."));}
1075 : }
1076 : */
1077 :
1078 129887 : void SIImageStore::accessImage( std::shared_ptr<ImageInterface<Float> > &ptr,
1079 : std::shared_ptr<ImageInterface<Float> > &parentptr,
1080 : const String label )
1081 : {
1082 : // if ptr is not null, assume it's OK. Perhaps add more checks.
1083 :
1084 : //cerr << "ACCIM itsptr" << ptr.get() << " parent " << parentptr.get() << " label " << label << endl;
1085 :
1086 129887 : Bool sw=False;
1087 129887 : if( label.contains(imageExts(SUMWT)) ) sw=True;
1088 :
1089 129887 : if( !ptr )
1090 : {
1091 : //cout << itsNFacets << " " << itsNChanChunks << " " << itsNPolChunks << endl;
1092 32365 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
1093 : {
1094 816 : if( ! parentptr )
1095 : {
1096 : //cout << "Making parent : " << itsImageName+label << " sw : " << sw << endl;
1097 636 : parentptr = openImage(itsImageName+label , itsOverWrite, sw, itsNFacets );
1098 636 : if( sw) {setUseWeightImage( *parentptr, itsUseWeight ); }
1099 : }
1100 : // cout << "Making facet " << itsFacetId << " out of " << itsNFacets << endl;
1101 : //cout << "Making chunk " << itsChanId << " out of " << itsNChanChunks << endl;
1102 : //ptr = makeFacet( itsFacetId, itsNFacets*itsNFacets, *parentptr );
1103 1632 : ptr = makeSubImage( itsFacetId, itsNFacets*itsNFacets,
1104 : itsChanId, itsNChanChunks,
1105 : itsPolId, itsNPolChunks,
1106 1632 : *parentptr );
1107 816 : if( ! sw )
1108 : {
1109 644 : itsImageShape = ptr->shape(); // over and over again.... FIX.
1110 644 : itsCoordSys = ptr->coordinates();
1111 644 : itsMiscInfo = ptr->miscInfo();
1112 : }
1113 :
1114 : //cout << "accessImage : " << label << " : sumwt : " << sw << " : shape : " << itsImageShape << endl;
1115 :
1116 : }
1117 : else
1118 : { //no facets of chanchunks
1119 31549 : if(!parentptr){
1120 : ///coordsys for psf can be different ...shape should be the same.
1121 30887 : ptr = openImage(itsImageName+label , itsOverWrite, sw, 1, !(label.contains(imageExts(PSF)))); }
1122 : else{
1123 662 : ptr=parentptr;
1124 : }
1125 : //cout << "Opening image : " << itsImageName+label << " of shape " << ptr->shape() << endl;
1126 : }
1127 : }
1128 :
1129 : //cerr << "ACCIM2 ptr " << ptr.get() << " lock " << itsReadLock.get() << endl;
1130 : //itsReadLock.reset(new LatticeLocker(*ptr, FileLocker::Read));
1131 :
1132 : // DirectionCoordinate dcord = ptr->coordinates().directionCoordinate();
1133 : // cout << "Latpole from " << label << " : " << dcord.longLatPoles() << endl;
1134 :
1135 129870 : }
1136 :
1137 :
1138 14174 : std::shared_ptr<ImageInterface<Float> > SIImageStore::psf(uInt /*nterm*/)
1139 : {
1140 14174 : accessImage( itsPsf, itsParentPsf, imageExts(PSF) );
1141 :
1142 14173 : return itsPsf;
1143 : }
1144 30828 : std::shared_ptr<ImageInterface<Float> > SIImageStore::residual(uInt /*nterm*/)
1145 : {
1146 30828 : accessImage( itsResidual, itsParentResidual, imageExts(RESIDUAL) );
1147 : // Record mi = itsResidual->miscInfo(); ostringstream oss;mi.print(oss);cout<<"MiscInfo(res) : " << oss.str() << endl;
1148 30828 : return itsResidual;
1149 : }
1150 1821 : std::shared_ptr<ImageInterface<Float> > SIImageStore::weight(uInt /*nterm*/)
1151 : {
1152 1821 : accessImage( itsWeight, itsParentWeight, imageExts(WEIGHT) );
1153 1821 : return itsWeight;
1154 : }
1155 :
1156 10010 : std::shared_ptr<ImageInterface<Float> > SIImageStore::sumwt(uInt /*nterm*/)
1157 : {
1158 :
1159 10010 : accessImage( itsSumWt, itsParentSumWt, imageExts(SUMWT) );
1160 :
1161 10010 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
1162 536 : { itsUseWeight = getUseWeightImage( *itsParentSumWt );}
1163 10010 : setUseWeightImage( *itsSumWt , itsUseWeight); // Sets a flag in the SumWt image.
1164 :
1165 10010 : return itsSumWt;
1166 : }
1167 :
1168 11063 : std::shared_ptr<ImageInterface<Float> > SIImageStore::model(uInt /*nterm*/)
1169 : {
1170 11065 : accessImage( itsModel, itsParentModel, imageExts(MODEL) );
1171 11061 : LatticeLocker lock1(*itsModel, FileLocker::Write);
1172 : // Set up header info the first time.
1173 11061 : itsModel->setUnits("Jy/pixel");
1174 :
1175 22122 : return itsModel;
1176 11061 : }
1177 5997 : std::shared_ptr<ImageInterface<Float> > SIImageStore::image(uInt /*nterm*/)
1178 : {
1179 5997 : accessImage( itsImage, itsParentImage, imageExts(IMAGE) );
1180 5997 : LatticeLocker lock1(*itsImage, FileLocker::Write);
1181 5997 : itsImage->setUnits(Unit("Jy/beam"));
1182 11994 : return itsImage;
1183 5997 : }
1184 :
1185 14919 : std::shared_ptr<ImageInterface<Float> > SIImageStore::mask(uInt /*nterm*/)
1186 : {
1187 14919 : accessImage( itsMask, itsParentMask, imageExts(MASK) );
1188 14915 : return itsMask;
1189 : }
1190 6 : std::shared_ptr<ImageInterface<Float> > SIImageStore::gridwt(IPosition newshape)
1191 :
1192 : {
1193 6 : if(newshape.empty()){
1194 2 : accessImage( itsGridWt, itsParentGridWt, imageExts(GRIDWT) );
1195 : }
1196 : else{
1197 4 : if(!itsGridWt || (itsGridWt && (itsGridWt->shape() != newshape))){
1198 2 : itsGridWt.reset(); // deassign previous one hopefully it'll close it
1199 2 : CoordinateSystem newcoordsys=itsCoordSys;
1200 2 : if(newshape.nelements() > 4){
1201 0 : Matrix<Double> pc(1,1); pc = 1.0;
1202 0 : LinearCoordinate zc(Vector<String>(1, "FIELD_ORDER"), Vector<String>(1, ""), Vector<Double>(1, 0.0), Vector<Double>(1,1.0), pc, Vector<Double>(1,0.0));
1203 0 : newcoordsys.addCoordinate(zc);
1204 0 : }
1205 2 : itsGridWt.reset(new PagedImage<Float>(newshape, newcoordsys, itsImageName+ imageExts(GRIDWT)));
1206 2 : initMetaInfo(itsGridWt, itsImageName+ imageExts(GRIDWT));
1207 :
1208 2 : }
1209 : }
1210 : /// change the coordinate system here, to uv.
1211 6 : return itsGridWt;
1212 : }
1213 :
1214 31 : std::shared_ptr<ImageInterface<Float> > SIImageStore::tempworkimage(uInt /*term*/){
1215 31 : if(itsTempWorkIm) return itsTempWorkIm;
1216 31 : itsTempWorkIm.reset(new PagedImage<Float>(itsImageShape, itsCoordSys, itsImageName+ ".work.temp"));
1217 31 : static_cast<PagedImage<Float>* > (itsTempWorkIm.get())->set(0.0);
1218 31 : itsTempWorkIm->flush();
1219 31 : static_cast<PagedImage<Float>* > (itsTempWorkIm.get())->table().markForDelete();
1220 31 : return itsTempWorkIm;
1221 : }
1222 :
1223 12852 : std::shared_ptr<ImageInterface<Float> > SIImageStore::pb(uInt /*nterm*/)
1224 : {
1225 12852 : accessImage( itsPB, itsParentPB, imageExts(PB) );
1226 12850 : return itsPB;
1227 : }
1228 261 : std::shared_ptr<ImageInterface<Float> > SIImageStore::imagepbcor(uInt /*nterm*/)
1229 : {
1230 261 : accessImage( itsImagePBcor, itsParentImagePBcor, imageExts(IMAGEPBCOR) );
1231 261 : LatticeLocker lock1(*itsImagePBcor, FileLocker::Write);
1232 261 : itsImagePBcor->setUnits(Unit("Jy/beam"));
1233 522 : return itsImagePBcor;
1234 261 : }
1235 :
1236 1916 : std::shared_ptr<ImageInterface<Complex> > SIImageStore::forwardGrid(uInt /*nterm*/){
1237 1916 : if( itsForwardGrid ) // && (itsForwardGrid->shape() == itsImageShape))
1238 : {
1239 : // cout << "Forward grid has shape : " << itsForwardGrid->shape() << endl;
1240 1684 : return itsForwardGrid;
1241 : }
1242 232 : Vector<Int> whichStokes(0);
1243 232 : IPosition cimageShape;
1244 232 : cimageShape=itsImageShape;
1245 232 : MFrequency::Types freqframe = itsCoordSys.spectralCoordinate(itsCoordSys.findCoordinate(Coordinate::SPECTRAL)).frequencySystem(True);
1246 : // No need to set a conversion layer if image is in LSRK already or it is 'Undefined'
1247 232 : if(freqframe != MFrequency::LSRK && freqframe!=MFrequency::Undefined && freqframe!=MFrequency::REST)
1248 0 : { itsCoordSys.setSpectralConversion("LSRK"); }
1249 232 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
1250 232 : whichStokes, itsDataPolRep);
1251 232 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
1252 232 : cimageShape(2)=whichStokes.nelements();
1253 : //cout << "Making forward grid of shape : " << cimageShape << " for imshape : " << itsImageShape << endl;
1254 232 : itsForwardGrid.reset( new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()) );
1255 : //if(image())
1256 232 : if(hasRestored()){
1257 16 : LatticeLocker lock1(*itsForwardGrid, FileLocker::Write);
1258 16 : itsForwardGrid->setImageInfo((image())->imageInfo());
1259 :
1260 16 : }
1261 232 : return itsForwardGrid;
1262 232 : }
1263 :
1264 2514 : std::shared_ptr<ImageInterface<Complex> > SIImageStore::backwardGrid(uInt /*nterm*/){
1265 2514 : if( itsBackwardGrid ) //&& (itsBackwardGrid->shape() == itsImageShape))
1266 : {
1267 : // cout << "Backward grid has shape : " << itsBackwardGrid->shape() << endl;
1268 2027 : return itsBackwardGrid;
1269 : }
1270 487 : Vector<Int> whichStokes(0);
1271 487 : IPosition cimageShape;
1272 487 : cimageShape=itsImageShape;
1273 487 : MFrequency::Types freqframe = itsCoordSys.spectralCoordinate(itsCoordSys.findCoordinate(Coordinate::SPECTRAL)).frequencySystem(True);
1274 : // No need to set a conversion layer if image is in LSRK already or it is 'Undefined'
1275 487 : if(freqframe != MFrequency::LSRK && freqframe!=MFrequency::Undefined && freqframe!=MFrequency::REST )
1276 0 : { itsCoordSys.setSpectralConversion("LSRK"); }
1277 487 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
1278 487 : whichStokes, itsDataPolRep);
1279 487 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
1280 487 : cimageShape(2)=whichStokes.nelements();
1281 : //cout << "Making backward grid of shape : " << cimageShape << " for imshape : " << itsImageShape << endl;
1282 487 : itsBackwardGrid.reset( new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()) );
1283 : //if(image())
1284 487 : if(hasRestored()){
1285 8 : LatticeLocker lock1(*itsBackwardGrid, FileLocker::Write);
1286 8 : itsBackwardGrid->setImageInfo((image())->imageInfo());
1287 8 : }
1288 487 : return itsBackwardGrid;
1289 487 : }
1290 2284 : Double SIImageStore::memoryBeforeLattice(){
1291 : //Calculate how much memory to use per temporary images before disking
1292 2284 : return 1.0; /// in MB
1293 : }
1294 2284 : IPosition SIImageStore::tileShape(){
1295 : //Need to have settable stuff here or algorith to determine this
1296 2284 : return IPosition(4, min(itsImageShape[0],1000), min(itsImageShape[1],1000), 1, 1);
1297 : }
1298 :
1299 : // TODO : Move to an image-wrapper class ? Same function exists in SynthesisDeconvolver.
1300 47452 : Bool SIImageStore::doesImageExist(String imagename)
1301 : {
1302 94904 : LogIO os( LogOrigin("SIImageStore","doesImageExist",WHERE) );
1303 47452 : Directory image( imagename );
1304 94904 : return image.exists();
1305 47452 : }
1306 :
1307 :
1308 0 : void SIImageStore::resetImages( Bool resetpsf, Bool resetresidual, Bool resetweight )
1309 : {
1310 0 : if( resetpsf ){
1311 0 : LatticeLocker lock1(*(psf()), FileLocker::Write);
1312 0 : psf()->set(0.0);
1313 :
1314 0 : }
1315 0 : if( resetresidual ) {
1316 : // removeMask( residual() );
1317 0 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1318 0 : residual()->set(0.0);
1319 0 : }
1320 0 : if( resetweight && itsWeight ){
1321 0 : LatticeLocker lock1(*(weight()), FileLocker::Write);
1322 0 : weight()->set(0.0);
1323 0 : }
1324 0 : if( resetweight ){
1325 0 : LatticeLocker lock1(*(sumwt()), FileLocker::Write);
1326 0 : sumwt()->set(0.0);
1327 0 : }
1328 0 : }
1329 :
1330 0 : void SIImageStore::addImages( std::shared_ptr<SIImageStore> imagestoadd,
1331 : Bool addpsf, Bool addresidual, Bool addweight, Bool adddensity)
1332 : {
1333 :
1334 0 : if(addpsf)
1335 : {
1336 0 : LatticeExpr<Float> adderPsf( *(psf()) + *(imagestoadd->psf()) );
1337 0 : psf()->copyData(adderPsf);
1338 0 : }
1339 0 : if(addresidual)
1340 : {
1341 0 : LatticeExpr<Float> adderRes( *(residual()) + *(imagestoadd->residual()) );
1342 0 : residual()->copyData(adderRes);
1343 0 : }
1344 0 : if(addweight)
1345 : {
1346 0 : if( getUseWeightImage( *(imagestoadd->sumwt()) ) ) // Access and add weight only if it is needed.
1347 : {
1348 0 : LatticeExpr<Float> adderWeight( *(weight()) + *(imagestoadd->weight()) );
1349 0 : weight()->copyData(adderWeight);
1350 0 : }
1351 :
1352 : /*
1353 : Array<Float> qqq, www;
1354 : imagestoadd->sumwt()->get(qqq,True);
1355 : sumwt()->get(www,True);
1356 : cout << "SUMWT : Adding : " << qqq << " to " << www << endl;
1357 : */
1358 :
1359 0 : LatticeExpr<Float> adderSumWt( *(sumwt()) + *(imagestoadd->sumwt()) );
1360 0 : sumwt()->copyData(adderSumWt);
1361 0 : setUseWeightImage( *sumwt(), getUseWeightImage(*(imagestoadd->sumwt()) ) );
1362 :
1363 0 : }
1364 0 : if(adddensity)
1365 : {
1366 0 : LatticeExpr<Float> adderDensity( *(gridwt()) + *(imagestoadd->gridwt()) );
1367 0 : gridwt()->copyData(adderDensity);
1368 0 : }
1369 :
1370 0 : }
1371 :
1372 0 : void SIImageStore::setWeightDensity( std::shared_ptr<SIImageStore> imagetoset )
1373 : {
1374 0 : LogIO os( LogOrigin("SIImageStore","setWeightDensity",WHERE) );
1375 : //gridwt()->copyData( LatticeExpr<Float> ( *(imagetoset->gridwt()) ) );
1376 : //looks like you have to call them before hand or it crashes in parallel sometimes
1377 0 : IPosition shape=(imagetoset->gridwt())->shape();
1378 0 : os << LogIO::DEBUG2 << "SHAPES " << imagetoset->gridwt()->shape() << " " << gridwt()->shape() << endl;
1379 0 : (gridwt(shape))->copyData( LatticeExpr<Float> ( *(imagetoset->gridwt()) ) );
1380 :
1381 0 : }
1382 :
1383 : // TODO
1384 : // cout << "WARN : get getPbMax right for cube.... weight is indexed on chan and pol." << endl;
1385 192 : Double SIImageStore::getPbMax()
1386 : {
1387 :
1388 : //// Don't do any extra norm. Minor cycle will operate with native PB.
1389 : //return 1.0;
1390 :
1391 : //// Normalize PB to 1 at the center of the image OF CHANNEL ZERO
1392 :
1393 : // IPosition shp = weight(0)->shape();
1394 : // IPosition center(4, shp[0]/2, shp[1]/2,0,0);
1395 : // return sqrt( weight(0)->getAt( center ) );
1396 :
1397 :
1398 : //// Normalize PB to 1 at the location of the maximum (across all chans..)
1399 :
1400 384 : LatticeExprNode le( sqrt(max( *(weight(0)) )) );
1401 384 : return le.getFloat();
1402 :
1403 192 : }
1404 :
1405 1070 : Double SIImageStore::getPbMax(Int pol,Int chan)
1406 : {
1407 :
1408 : //// Normalize PB to 1 at the location of the maximum (per pol,chan)
1409 :
1410 1070 : CountedPtr<ImageInterface<Float> > subim=makeSubImage(0, 1,
1411 1070 : chan, itsImageShape[3],
1412 1070 : pol, itsImageShape[2],
1413 2140 : *weight(0));
1414 :
1415 2140 : return sqrt(max(subim->get()));
1416 1070 : }
1417 :
1418 112 : void SIImageStore::makePBFromWeight(const Float pblimit)
1419 : {
1420 224 : LogIO os( LogOrigin("SIImageStore","makePBFromWeight",WHERE) );
1421 :
1422 224 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1423 : {
1424 348 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1425 : {
1426 :
1427 236 : itsPBScaleFactor = getPbMax(pol,chan);
1428 :
1429 236 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1430 : else {
1431 :
1432 233 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1433 233 : chan, itsImageShape[3],
1434 233 : pol, itsImageShape[2],
1435 466 : *weight() );
1436 233 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
1437 233 : chan, itsImageShape[3],
1438 233 : pol, itsImageShape[2],
1439 466 : *pb() );
1440 :
1441 :
1442 466 : LatticeExpr<Float> normed( sqrt(abs(*wtsubim)) / itsPBScaleFactor );
1443 466 : LatticeExpr<Float> limited( iif( normed > fabs(pblimit) , normed, 0.0 ) );
1444 233 : pbsubim->copyData( limited );
1445 233 : }// if not zero
1446 : }
1447 : }
1448 :
1449 112 : if((pb()->getDefaultMask()==""))
1450 : {
1451 : //Remove the old mask as it is no longer valid
1452 : //removeMask( pb() );
1453 :
1454 : // if( pblimit >= 0.0 )
1455 : {
1456 : //MSK//
1457 170 : LatticeExpr<Bool> pbmask( iif( *pb() > fabs(pblimit) , True , False ) );
1458 : //MSK//
1459 85 : createMask( pbmask, pb() );
1460 85 : pb()->pixelMask().unlock();
1461 85 : }
1462 :
1463 : }
1464 112 : weight()->unlock();
1465 112 : pb()->unlock();
1466 112 : }
1467 :
1468 590 : void SIImageStore::makePBImage(const Float pblimit)
1469 : {
1470 1180 : LogIO os( LogOrigin("SIImageStore","makePBImage",WHERE) );
1471 :
1472 590 : if( hasPB() )
1473 : {
1474 :
1475 1289 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1476 : {
1477 3398 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1478 : {
1479 :
1480 2699 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
1481 2699 : chan, itsImageShape[3],
1482 2699 : pol, itsImageShape[2],
1483 5398 : *pb() );
1484 :
1485 : // Norm by the max
1486 2699 : LatticeExprNode elmax= max( *pbsubim );
1487 2699 : Float fmax = abs(elmax.getFloat());
1488 : //If there are multiple overlap of beam such that the peak is larger than 1 then normalize
1489 : //otherwise leave as is
1490 2699 : if(fmax>1.0)
1491 : {
1492 238 : LatticeExpr<Float> normed( (*pbsubim) / fmax );
1493 238 : LatticeExpr<Float> limited( iif( normed > fabs(pblimit) , normed, 0.0 ) );
1494 119 : pbsubim->copyData( limited );
1495 119 : }
1496 : else
1497 : {
1498 5160 : LatticeExpr<Float> limited( iif((*pbsubim) > fabs(pblimit) , (*pbsubim), 0.0 ) );
1499 2580 : pbsubim->copyData( limited );
1500 2580 : }
1501 2699 : }
1502 : }
1503 :
1504 590 : if((pb()->getDefaultMask()==""))// && pblimit >= 0.0)
1505 : {
1506 : // removeMask( pb() );
1507 : //MSK//
1508 1136 : LatticeExpr<Bool> pbmask( iif( *pb() > fabs(pblimit) , True , False ) );
1509 : //MSK//
1510 568 : createMask( pbmask, pb() );
1511 568 : pb()->pixelMask().unlock();
1512 568 : }
1513 :
1514 : }// if hasPB
1515 590 : pb()->unlock();
1516 :
1517 590 : }
1518 :
1519 900 : Bool SIImageStore::createMask(LatticeExpr<Bool> &lemask,
1520 : CountedPtr<ImageInterface<Float> > outimage)
1521 : {
1522 : //cout << "Calling makeMask for mask0 for " << outimage->name() << endl;
1523 : try
1524 : {
1525 900 : LatticeLocker lock1(*outimage, FileLocker::Write);
1526 900 : ImageRegion outreg = outimage->makeMask("mask0",False,True);
1527 900 : LCRegion& outmask=outreg.asMask();
1528 900 : outmask.copyData(lemask);
1529 900 : outimage->defineRegion("mask0",outreg, RegionHandler::Masks, True);
1530 :
1531 900 : outimage->setDefaultMask("mask0");
1532 :
1533 900 : outimage->unlock();
1534 900 : outimage->tempClose();
1535 :
1536 : // outimage->table().unmarkForDelete();
1537 900 : }
1538 0 : catch (const AipsError& x) {
1539 0 : throw(AipsError("Error in creating internal T/F mask : " + x.getMesg() ));
1540 0 : }
1541 :
1542 900 : return True;
1543 : }
1544 :
1545 1746 : Bool SIImageStore::copyMask(CountedPtr<ImageInterface<Float> > inimage,
1546 : CountedPtr<ImageInterface<Float> > outimage)
1547 : {
1548 : // cout << "In copyMask for " << outimage->name() << endl;
1549 :
1550 : try
1551 : {
1552 1746 : if( (inimage->getDefaultMask()).matches("mask0") ) // input mask exists.
1553 1556 : {LatticeLocker lock1(*outimage, FileLocker::Write);
1554 1556 : removeMask(outimage);
1555 :
1556 : // // clear output image mask
1557 : // if( (outimage->getDefaultMask()).matches("mask0") )
1558 : // {outimage->setDefaultMask("");
1559 : // outimage->removeRegion("mask0");}
1560 : // get mask from input image
1561 :
1562 1559 : ImageRegion outreg=outimage->makeMask("mask0", False, True);
1563 1553 : LCRegion& outmask=outreg.asMask();
1564 1553 : outmask.copyData(inimage->getRegion("mask0").asLCRegion());
1565 1553 : outimage->defineRegion("mask0",outreg, RegionHandler::Masks,True);
1566 1553 : outimage->setDefaultMask("mask0");
1567 1556 : }
1568 : }
1569 3 : catch (const AipsError& x) {
1570 3 : throw(AipsError("Error in copying internal T/F mask : " + x.getMesg() ));
1571 3 : }
1572 :
1573 1743 : return True;
1574 : }
1575 :
1576 1891 : void SIImageStore::removeMask(CountedPtr<ImageInterface<Float> > im)
1577 : {
1578 : try
1579 : {
1580 : // // clear output image mask
1581 : // if( (im->getDefaultMask()).matches("mask0") )
1582 : // {im->setDefaultMask("");
1583 : // im->removeRegion("mask0");}
1584 : ///Remove the old mask as it is no longer valid
1585 1891 : LatticeLocker lock1(*im, FileLocker::Write);
1586 1891 : if (im-> getDefaultMask() != String(""))
1587 : {
1588 388 : String strung=im->getDefaultMask();
1589 388 : im->setDefaultMask("");
1590 388 : im->removeRegion(strung);
1591 388 : }
1592 1891 : if( im->hasRegion("mask0") )
1593 : {
1594 0 : im->removeRegion("mask0");
1595 : }
1596 1891 : }
1597 0 : catch (const AipsError& x)
1598 : {
1599 0 : throw(AipsError("Error in deleting internal T/F mask : " + x.getMesg() ));
1600 0 : }
1601 1891 : }
1602 226 : void SIImageStore:: rescaleResolution(Int chan,
1603 : ImageInterface<Float>& image,
1604 : const GaussianBeam& newbeam,
1605 : const GaussianBeam& oldbeam){
1606 :
1607 452 : LogIO os( LogOrigin("SIImageStore","rescaleResolution",WHERE) );
1608 452 : GaussianBeam toBeUsed(Quantity(0.0, "arcsec"),Quantity(0.0, "arcsec"),
1609 678 : Quantity(0.0, "deg")) ;
1610 : try {
1611 226 : Bool samesize = GaussianDeconvolver::deconvolve(toBeUsed, newbeam, oldbeam);
1612 :
1613 : /*
1614 : os << LogIO::NORMAL2 << "Chan : " << chan << " : Input beam : : " << oldbeam.getMajor(Unit("arcsec")) << " arcsec, " << oldbeam.getMinor(Unit("arcsec"))<< " arcsec, " << oldbeam.getPA(Unit("deg")) << " deg" << LogIO::POST;
1615 : os << LogIO::NORMAL2 << "Target beam : " << newbeam.getMajor(Unit("arcsec")) << " arcsec, " << newbeam.getMinor(Unit("arcsec"))<< " arcsec, " << newbeam.getPA(Unit("deg")) << " deg" << LogIO::POST;
1616 : os << LogIO::NORMAL2 << "Beam to be used : " << toBeUsed.getMajor(Unit("arcsec")) << " arcsec, " << toBeUsed.getMinor(Unit("arcsec"))<< " arcsec, " << toBeUsed.getPA(Unit("deg")) << " deg" << LogIO::POST;
1617 : os << LogIO::NORMAL2 << "SameSize ? " << samesize << endl;
1618 : */
1619 :
1620 225 : if( samesize )
1621 : {
1622 8 : os << LogIO::NORMAL2 << "Input and output beam sizes are the same for Channel : " << chan << ". Not convolving residuals." << LogIO::POST;
1623 : }
1624 : else
1625 : {
1626 217 : Double pixwidth=sqrt(image.coordinates().increment()(0)*image.coordinates().increment()(0)+image.coordinates().increment()(1)*image.coordinates().increment()(1));
1627 :
1628 217 : if(toBeUsed.getMinor(image.coordinates().worldAxisUnits()[0]) > pixwidth)
1629 : {
1630 : //cerr << "old beam area " << oldbeam.getArea("rad2") << " new beam " << newbeam.getArea("rad2") << endl;
1631 211 : StokesImageUtil::Convolve(image, toBeUsed, True);
1632 211 : image.copyData(LatticeExpr<Float>(image*newbeam.getArea("rad2")/ oldbeam.getArea("rad2")));
1633 : }
1634 : }
1635 : }
1636 1 : catch (const AipsError& x) {
1637 : //throw(AipsError("Cannot convolve to new beam: may be smaller than old beam : " + x.getMesg() ));
1638 1 : os << LogIO::WARN << "Cannot convolve to new beam for Channel : " << chan << " : " << x.getMesg() << LogIO::POST;
1639 1 : }
1640 :
1641 :
1642 226 : }
1643 :
1644 :
1645 :
1646 :
1647 590 : void SIImageStore::dividePSFByWeight(const Float /* pblimit*/)
1648 : {
1649 1180 : LogIO os( LogOrigin("SIImageStore","dividePSFByWeight",WHERE) );
1650 :
1651 590 : LatticeLocker lock1 (*(psf()), FileLocker::Write);
1652 590 : normPSF();
1653 :
1654 590 : if( itsUseWeight )
1655 : {
1656 :
1657 70 : divideImageByWeightVal( *weight() );
1658 : }
1659 590 : (psf())->unlock();
1660 :
1661 590 : }
1662 :
1663 597 : void SIImageStore::normalizePrimaryBeam(const Float pblimit)
1664 : {
1665 1194 : LogIO os( LogOrigin("SIImageStore","normalizePrimaryBeam",WHERE) );
1666 :
1667 : // cout << "In dividePSFByWeight : itsUseWeight : " << itsUseWeight << endl;
1668 597 : if( itsUseWeight )
1669 : {
1670 :
1671 144 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1672 : {
1673 268 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1674 : {
1675 196 : os << LogIO::NORMAL1 << "Scale factor for [C" +String::toString(chan) + ":P" + String::toString(pol) + "] to keep the model image w.r.to a PB of max=1 is " << getPbMax(pol,chan) << LogIO::POST;
1676 : }//chan
1677 : }//pol
1678 :
1679 72 : makePBFromWeight(pblimit);
1680 :
1681 : }//if itsUseWeight
1682 525 : else { makePBImage(pblimit); } // OR... just check that it exists already.
1683 597 : pb()->unlock();
1684 :
1685 597 : }
1686 :
1687 : // Make another for the PSF too.
1688 1328 : void SIImageStore::divideResidualByWeight(Float pblimit, String normtype) {
1689 :
1690 2656 : LogIO os(LogOrigin("SIImageStore", "divideResidualByWeight", WHERE));
1691 1328 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1692 :
1693 362 : auto logTemplate = [&](Int const &chan, Int const &pol, string const &normalizer, string const &result) {
1694 362 : os << LogIO::NORMAL1
1695 724 : << "[C" + String::toString(chan) + ":P" + String::toString(pol) + "] "
1696 724 : << "Dividing " << itsImageName + String(".residual") << " by "
1697 : << "[ " << normalizer << " ] "
1698 1086 : << "to get " << result << "." << LogIO::POST;
1699 362 : };
1700 :
1701 : // Normalize by the sumwt, per plane.
1702 1328 : Bool didNorm = divideImageByWeightVal(*residual());
1703 :
1704 1328 : if (itsUseWeight) {
1705 278 : for (Int pol = 0; pol < itsImageShape[2]; pol++) {
1706 510 : for (Int chan = 0; chan < itsImageShape[3]; chan++) {
1707 371 : itsPBScaleFactor = getPbMax(pol, chan);
1708 :
1709 371 : if (itsPBScaleFactor <= 0) {
1710 : os << LogIO::NORMAL1
1711 : << "Skipping normalization for C:" << chan << " P:" << pol
1712 9 : << " because pb max is zero " << LogIO::POST;
1713 :
1714 : } else {
1715 362 : CountedPtr<ImageInterface<Float> > wtsubim = makeSubImage(0, 1,
1716 362 : chan, itsImageShape[3],
1717 362 : pol, itsImageShape[2],
1718 724 : *weight());
1719 362 : CountedPtr<ImageInterface<Float> > ressubim = makeSubImage(0, 1,
1720 362 : chan, itsImageShape[3],
1721 362 : pol, itsImageShape[2],
1722 724 : *residual());
1723 362 : LatticeExpr<Float> ratio;
1724 362 : Float scalepb = 1.0;
1725 :
1726 362 : if (normtype == "flatnoise") {
1727 1086 : logTemplate(chan, pol,
1728 724 : "sqrt(weightimage) * " + String::toString(itsPBScaleFactor),
1729 : "flat noise with unit pb peak");
1730 :
1731 724 : LatticeExpr<Float> deno = itsPBScaleFactor * sqrt(abs(LatticeExpr<Float>(*(wtsubim))));
1732 362 : scalepb = fabs(pblimit) * itsPBScaleFactor * itsPBScaleFactor;
1733 362 : ratio = iif(deno > scalepb, (*(ressubim) / deno), 0.0);
1734 :
1735 362 : } else if (normtype == "pbsquare") {
1736 0 : logTemplate(chan, pol,
1737 0 : String::toString(itsPBScaleFactor),
1738 : "optimal noise with unit pb peak");
1739 :
1740 0 : Float deno = itsPBScaleFactor * itsPBScaleFactor;
1741 0 : ratio = (*(ressubim) / deno);
1742 :
1743 0 : } else if (normtype == "flatsky") {
1744 0 : logTemplate(chan, pol, "weight", "flat sky");
1745 :
1746 0 : LatticeExpr<Float> deno = LatticeExpr<Float>(*(wtsubim));
1747 0 : scalepb = fabs(pblimit * pblimit) * itsPBScaleFactor * itsPBScaleFactor;
1748 0 : ratio = iif(deno > scalepb, (*(ressubim) / deno), 0.0);
1749 :
1750 0 : }
1751 :
1752 : //IPosition ip(4, itsImageShape[0] / 2, itsImageShape[1]/2, 0, 0);
1753 : //Float resval = ressubim->getAt(ip);
1754 : //LatticeExpr<Float> mask(iif((deno) > scalepb, 1.0, 0.0));
1755 : //LatticeExpr<Float> maskinv(iif((deno) > scalepb, 0.0, 1.0));
1756 : //LatticeExpr<Float> ratio(((*(ressubim)) * mask) / (deno + maskinv));
1757 :
1758 : //above blocks all sources outside minpb but visible with weight coverage
1759 : //which could be cleaned out...one could use below for that
1760 : //LatticeExpr<Float> ratio(iif(deno > scalepb, (*(ressubim)) / deno, *ressubim));
1761 :
1762 362 : ressubim->copyData(ratio);
1763 :
1764 : //cout << "Val of residual before|after normalizing at center for pol " << pol << " chan " << chan << " : " << resval << "|" << ressubim->getAt(ip) << " weight : " << wtsubim->getAt(ip) << endl;
1765 362 : } // if not zero
1766 : } //chan
1767 : } //pol
1768 : }
1769 :
1770 : // If no normalization happened, print a warning. The user must check if it's right or not.
1771 : // Or... later if we get a gridder that does pre-norms, this warning can go.
1772 1328 : if ((didNorm | itsUseWeight) != True) {
1773 0 : os << LogIO::WARN << "No normalization done to residual" << LogIO::POST;
1774 : }
1775 :
1776 : ///// A T/F mask in the residual will confuse users looking at the interactive clean
1777 : ///// window
1778 1328 : if ((residual()->getDefaultMask() == "") && hasPB() && pblimit >=0.0) {
1779 229 : copyMask(pb(), residual());
1780 : }
1781 :
1782 1328 : if ((pblimit < 0.0) && (residual()->getDefaultMask()).matches("mask0")) {
1783 1 : removeMask(residual());
1784 : }
1785 :
1786 1328 : residual()->unlock();
1787 1328 : }
1788 :
1789 137 : void SIImageStore::divideResidualByWeightSD(Float pblimit) {
1790 :
1791 274 : LogIO os(LogOrigin("SIImageStore", "divideResidualByWeightSD", WHERE));
1792 137 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1793 :
1794 137 : if (itsUseWeight) {
1795 274 : LatticeExpr<Float> deno = LatticeExpr<Float>(*weight());
1796 274 : LatticeExpr<Float> ratio = iif(deno > 0.0, *(residual()) / deno, 0.0);
1797 137 : residual()->copyData(ratio);
1798 137 : }
1799 : else {
1800 : // If no normalization happened, print a warning. The user must check if it's right or not.
1801 : // Or... later if we get a gridder that does pre-norms, this warning can go.
1802 0 : os << LogIO::WARN << "No normalization done to residual" << LogIO::POST;
1803 : }
1804 :
1805 : ///// A T/F mask in the residual will confuse users looking at the interactive clean
1806 : ///// window
1807 137 : if ((residual()->getDefaultMask() == "") && hasPB() && pblimit >=0.0) {
1808 0 : copyMask(pb(), residual());
1809 : }
1810 :
1811 137 : if ((pblimit < 0.0) && (residual()->getDefaultMask()).matches("mask0")) {
1812 0 : removeMask(residual());
1813 : }
1814 :
1815 137 : residual()->unlock();
1816 137 : }
1817 :
1818 868 : void SIImageStore::divideModelByWeight(Float pblimit, const String normtype)
1819 : {
1820 1736 : LogIO os( LogOrigin("SIImageStore","divideModelByWeight",WHERE) );
1821 :
1822 :
1823 1736 : if(itsUseWeight // only when needed
1824 868 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
1825 : {
1826 :
1827 92 : if( normtype=="flatsky") {
1828 :
1829 0 : LatticeExprNode LEN = max( *model());
1830 0 : os << LogIO::NORMAL1 << "Model is already flat sky with peak flux : " << LEN.getFloat();
1831 0 : os << ". No need to divide before prediction" << LogIO::POST;
1832 :
1833 0 : return;
1834 0 : }
1835 92 : else if( normtype=="flatnoise"){
1836 :
1837 184 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1838 : {
1839 339 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1840 : {
1841 :
1842 247 : itsPBScaleFactor = getPbMax(pol,chan);
1843 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1844 247 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1845 : else {
1846 :
1847 221 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1848 221 : chan, itsImageShape[3],
1849 221 : pol, itsImageShape[2],
1850 442 : *weight() );
1851 221 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1852 221 : chan, itsImageShape[3],
1853 221 : pol, itsImageShape[2],
1854 442 : *model() );
1855 221 : os << LogIO::NORMAL1 ;
1856 221 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1857 221 : os << "Dividing " << itsImageName+String(".model") ;
1858 221 : os << " by [ sqrt(weight) / " << itsPBScaleFactor ;
1859 221 : os <<" ] to get to flat sky model before prediction" << LogIO::POST;
1860 :
1861 :
1862 442 : LatticeExpr<Float> deno( sqrt( abs(*(wtsubim)) ) / itsPBScaleFactor );
1863 :
1864 442 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1865 442 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1866 442 : LatticeExpr<Float> ratio( ( (*(modsubim)) * mask ) / ( deno + maskinv ) );
1867 :
1868 : //
1869 : //The above has a problem...mask is cutting out clean components found
1870 : // outside pblimit ...use below if this is what is wanted
1871 : //LatticeExpr<Float> ratio(iif(sqrt(abs(*(wtsubim))) < (fabs(pblimit)*itsPBScaleFactor), 0.0, (*(modsubim))/(sqrt( abs(*(wtsubim))) / itsPBScaleFactor)));
1872 : // LatticeExpr<Float> ratio(iif((sqrt(abs(*(wtsubim))) ==0.0), 0.0, ((*(modsubim))*itsPBScaleFactor)/sqrt( abs(*(wtsubim))) ));
1873 :
1874 221 : IPosition ip(4,itsImageShape[0]/2,itsImageShape[1]/2,0,0);
1875 : /// Float modval = modsubim->getAt(ip);
1876 : //LatticeExprNode aminval( min(*modsubim) );
1877 : //LatticeExprNode amaxval( max(*modsubim) );
1878 : //cout << "Before ---- min : " << aminval.getFloat() << " max : " << amaxval.getFloat() << endl;
1879 :
1880 221 : modsubim->copyData(ratio);
1881 : // cout << "Val of model before|after flattening at center for pol " << pol << " chan " << chan << " : " << modval << "|" << modsubim->getAt(ip) << " weight : " << wtsubim->getAt(ip) << endl;
1882 : //LatticeExprNode minval( min(*modsubim) );
1883 : //LatticeExprNode maxval( max(*modsubim) );
1884 : //cout << "After ---- min : " << minval.getFloat() << " max : " << maxval.getFloat() << endl;
1885 221 : }// if not zero
1886 : }//chan
1887 : }//pol
1888 :
1889 : }
1890 0 : else if( normtype=="pbsquare"){
1891 :
1892 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1893 : {
1894 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1895 : {
1896 :
1897 0 : itsPBScaleFactor = getPbMax(pol,chan);
1898 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1899 0 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1900 : else {
1901 :
1902 0 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1903 0 : chan, itsImageShape[3],
1904 0 : pol, itsImageShape[2],
1905 0 : *weight() );
1906 0 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1907 0 : chan, itsImageShape[3],
1908 0 : pol, itsImageShape[2],
1909 0 : *model() );
1910 0 : os << LogIO::NORMAL1 ;
1911 0 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1912 0 : os << "Dividing " << itsImageName+String(".model") ;
1913 0 : os << " by [ (weight) / " << itsPBScaleFactor ;
1914 0 : os <<" ] to restore to optimal sky model before prediction" << LogIO::POST;
1915 :
1916 :
1917 0 : LatticeExpr<Float> deno( abs(*(wtsubim)) / (itsPBScaleFactor*itsPBScaleFactor) );
1918 :
1919 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1920 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1921 0 : LatticeExpr<Float> ratio( ( (*(modsubim)) * mask ) / ( deno + maskinv ) );
1922 :
1923 : //
1924 : //The above has a problem...mask is cutting out clean components found
1925 : // outside pblimit ...use below if this is what is wanted
1926 : // LatticeExpr<Float> ratio(iif(abs(*(wtsubim)) == 0.0, *modsubim, (*(modsubim))/( abs(*(wtsubim)) / (itsPBScaleFactor*itsPBScaleFactor))));
1927 :
1928 0 : IPosition ip(4,itsImageShape[0]/2,itsImageShape[1]/2,0,0);
1929 : /// Float modval = modsubim->getAt(ip);
1930 : //LatticeExprNode aminval( min(*modsubim) );
1931 : //LatticeExprNode amaxval( max(*modsubim) );
1932 : //cout << "Before ---- min : " << aminval.getFloat() << " max : " << amaxval.getFloat() << endl;
1933 :
1934 0 : modsubim->copyData(ratio);
1935 :
1936 : // cout << "Val of model before|after flattening at center for pol " << pol << " chan " << chan << " : " << modval << "|" << modsubim->getAt(ip) << " weight : " << wtsubim->getAt(ip) << endl;
1937 : //LatticeExprNode minval( min(*modsubim) );
1938 : //LatticeExprNode maxval( max(*modsubim) );
1939 : //cout << "After ---- min : " << minval.getFloat() << " max : " << maxval.getFloat() << endl;
1940 0 : }// if not zero
1941 : }//chan
1942 : }//pol
1943 :
1944 : }
1945 :
1946 : // storeImg(String("flatmodel.im"), *model());
1947 :
1948 : }
1949 868 : }
1950 :
1951 797 : void SIImageStore::multiplyModelByWeight(Float pblimit, const String normtype)
1952 : {
1953 1594 : LogIO os( LogOrigin("SIImageStore","multiplyModelByWeight",WHERE) );
1954 :
1955 1594 : if(itsUseWeight // only when needed
1956 797 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
1957 : {
1958 20 : if( normtype=="flatsky") {
1959 0 : os << "Model is already flat sky. No need to multiply back after prediction" << LogIO::POST;
1960 0 : return;
1961 : }
1962 20 : else if( normtype=="flatnoise"){
1963 :
1964 40 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1965 : {
1966 40 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1967 : {
1968 :
1969 20 : itsPBScaleFactor = getPbMax(pol,chan);
1970 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1971 20 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1972 : else {
1973 :
1974 20 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1975 20 : chan, itsImageShape[3],
1976 20 : pol, itsImageShape[2],
1977 40 : *weight() );
1978 20 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1979 20 : chan, itsImageShape[3],
1980 20 : pol, itsImageShape[2],
1981 40 : *model() );
1982 :
1983 :
1984 :
1985 20 : os << LogIO::NORMAL1 ;
1986 20 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1987 20 : os << "Multiplying " << itsImageName+String(".model") ;
1988 20 : os << " by [ sqrt(weight) / " << itsPBScaleFactor;
1989 20 : os << " ] to take model back to flat noise with unit pb peak." << LogIO::POST;
1990 :
1991 40 : LatticeExpr<Float> deno( sqrt( abs(*(wtsubim)) ) / itsPBScaleFactor );
1992 :
1993 40 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1994 40 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1995 40 : LatticeExpr<Float> ratio( ( (*(modsubim)) * mask ) * ( deno + maskinv ) );
1996 :
1997 : /////See comment in divmodel and divresidual for below usage
1998 : //LatticeExpr<Float> ratio ( (*(modsubim)) * sqrt( abs(*(wtsubim))) / itsPBScaleFactor );
1999 20 : modsubim->copyData(ratio);
2000 :
2001 20 : }// if not zero
2002 : }//chan
2003 : }//pol
2004 : }
2005 0 : else if( normtype=="pbsquare"){
2006 :
2007 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2008 : {
2009 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
2010 : {
2011 :
2012 0 : itsPBScaleFactor = getPbMax(pol,chan);
2013 : // cout << " pbscale : " << itsPBScaleFactor << endl;
2014 0 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
2015 : else {
2016 :
2017 0 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
2018 0 : chan, itsImageShape[3],
2019 0 : pol, itsImageShape[2],
2020 0 : *weight() );
2021 0 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
2022 0 : chan, itsImageShape[3],
2023 0 : pol, itsImageShape[2],
2024 0 : *model() );
2025 :
2026 :
2027 :
2028 0 : os << LogIO::NORMAL1 ;
2029 0 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
2030 0 : os << "Multiplying " << itsImageName+String(".model") ;
2031 0 : os << " by [ weight / " << itsPBScaleFactor*itsPBScaleFactor;
2032 0 : os << " ] to take model back to flat noise with unit pb peak." << LogIO::POST;
2033 :
2034 0 : LatticeExpr<Float> deno( abs(*(wtsubim)) / (itsPBScaleFactor*itsPBScaleFactor) );
2035 :
2036 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
2037 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
2038 0 : LatticeExpr<Float> ratio( ( (*(modsubim)) * mask ) * ( deno + maskinv ) );
2039 :
2040 : /////See comment in divmodel and divresidual for below usage
2041 : //LatticeExpr<Float> ratio ( (*(modsubim)) * sqrt( abs(*(wtsubim)) / itsPBScaleFactor) );
2042 0 : modsubim->copyData(ratio);
2043 0 : }// if not zero
2044 : }//chan
2045 : }//pol
2046 : }
2047 : }
2048 797 : }
2049 :
2050 0 : GaussianBeam SIImageStore::getPSFGaussian(Float psfcutoff)
2051 : {
2052 0 : GaussianBeam beam;
2053 : try
2054 : {
2055 0 : if( psf()->ndim() > 0 )
2056 : {
2057 0 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
2058 0 : StokesImageUtil::FitGaussianPSF( *(psf()), beam, psfcutoff);
2059 0 : }
2060 : }
2061 0 : catch(AipsError &x)
2062 : {
2063 : // LogIO os( LogOrigin("SIImageStore","getPSFGaussian",WHERE) );
2064 : // os << "Error in fitting a Gaussian to the PSF : " << x.getMesg() << LogIO::POST;
2065 0 : throw( AipsError("Error in fitting a Gaussian to the PSF : " + x.getMesg()) );
2066 0 : }
2067 :
2068 0 : return beam;
2069 0 : }
2070 :
2071 1782 : void SIImageStore::makeImageBeamSet(Float psfcutoff, const Bool forcefit)
2072 : {
2073 1782 : clock_t begin = clock();
2074 :
2075 3564 : LogIO os( LogOrigin("SIImageStore","getPSFGaussian",WHERE) );
2076 : // For all chans/pols, call getPSFGaussian() and put it into ImageBeamSet(chan,pol).
2077 1782 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2078 1782 : uInt nx = itsImageShape[0];
2079 1782 : uInt ny = itsImageShape[1];
2080 1782 : uInt npol = itsImageShape[2];
2081 1782 : uInt nchan = itsImageShape[3];
2082 1782 : ImageInfo ii = psf()->imageInfo();
2083 1782 : ImageBeamSet iibeamset=ii.getBeamSet();
2084 1782 : if(iibeamset.nchan()==nchan && iibeamset.nstokes()==npol && forcefit==False){
2085 1064 : itsPSFBeams=iibeamset;
2086 1064 : itsRestoredBeams=iibeamset;
2087 1064 : return;
2088 : }
2089 718 : itsPSFBeams.resize( nchan, npol );
2090 718 : itsRestoredBeams.resize(nchan, npol);
2091 : // cout << "makeImBeamSet : imshape : " << itsImageShape << endl;
2092 :
2093 718 : String blankpsfs="";
2094 3469 : for( uInt chanid=0; chanid<nchan;chanid++) {
2095 5717 : for( uInt polid=0; polid<npol; polid++ ) {
2096 2966 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
2097 :
2098 2966 : IPosition substart(4,0,0,polid,chanid);
2099 2966 : IPosition substop(4,nx-1,ny-1,polid,chanid);
2100 2966 : Slicer psfslice(substart, substop,Slicer::endIsLast);
2101 5932 : SubImage<Float> subPsf( *psf() , psfslice, True );
2102 2966 : GaussianBeam beam;
2103 :
2104 2966 : Bool tryfit=True;
2105 :
2106 2966 : LatticeExprNode le( max(subPsf) );
2107 2966 : tryfit = le.getFloat()>0.0;
2108 2966 : if(tryfit)
2109 : {
2110 : try
2111 : {
2112 2862 : StokesImageUtil::FitGaussianPSF( subPsf, beam,psfcutoff );
2113 2862 : itsPSFBeams.setBeam( chanid, polid, beam );
2114 2862 : itsRestoredBeams.setBeam(chanid, polid, beam);
2115 : }
2116 0 : catch(AipsError &x)
2117 : {
2118 0 : os << LogIO::WARN << "[Chan" << chanid << ":Pol" << polid << "] Error Gaussian fit to PSF : " << x.getMesg() ;
2119 : // os << LogIO::POST;
2120 0 : os << " : Setting restoring beam to largest valid beam." << LogIO::POST;
2121 0 : }
2122 : }
2123 : else
2124 : {
2125 : // os << LogIO::WARN << "[Chan" << chanid << ":Pol" << polid << "] PSF is blank. Setting null restoring beam." << LogIO::POST ;
2126 104 : blankpsfs += "[C" +String::toString(chanid) + ":P" + String::toString(polid) + "] ";
2127 : }
2128 :
2129 2966 : }// end of pol loop
2130 : }// end of chan loop
2131 :
2132 718 : if( blankpsfs.length() >0 )
2133 39 : os << LogIO::WARN << "PSF is blank for" << blankpsfs << LogIO::POST;
2134 :
2135 : //// Replace null (and bad) beams with the good one.
2136 : ////GaussianBeam maxbeam = findGoodBeam(True);
2137 :
2138 : //// Replace null beams by a tiny tiny beam, just to get past the ImageInfo restriction that
2139 : //// all planes must have non-null beams.
2140 :
2141 718 : GaussianBeam defaultbeam=itsPSFBeams.getMaxAreaBeam();
2142 : ///many of the unittests in tsdimaging seem to depend on having 0 area beams
2143 : ///which throws and exception when it is stored in the image
2144 : ///(so setting them to some small number)!
2145 718 : if(defaultbeam.getArea("rad2")==0.0){
2146 5 : Quantity majax(1e-06,"arcsec"),minax(1e-06,"arcsec"),pa(0.0,"deg");
2147 5 : defaultbeam.setMajorMinor(majax,minax);
2148 5 : defaultbeam.setPA(pa);
2149 5 : }
2150 3469 : for( uInt chanid=0; chanid<nchan;chanid++) {
2151 5717 : for( uInt polid=0; polid<npol; polid++ ) {
2152 2966 : if( (itsPSFBeams.getBeam(chanid, polid)).isNull() )
2153 104 : { itsPSFBeams.setBeam( chanid, polid, defaultbeam );
2154 104 : itsRestoredBeams.setBeam( chanid, polid, defaultbeam );
2155 : }
2156 : }// end of pol loop
2157 : }// end of chan loop
2158 :
2159 : /*
2160 : //// Fill in gaps if there are any --- with the MAX Area beam.
2161 : ///// GaussianBeam maxbeam = itsPSFBeams.getMaxAreaBeam();
2162 : if( maxbeam.isNull() ) {
2163 : os << LogIO::WARN << "No planes have non-zero restoring beams. Forcing artificial 1.0arcsec beam." << LogIO::POST;
2164 : Quantity majax(1.0,"arcsec"),minax(1.0,"arcsec"),pa(0.0,"deg");
2165 : maxbeam.setMajorMinor(majax,minax);
2166 : maxbeam.setPA(pa);
2167 : }
2168 : else {
2169 : for( Int chanid=0; chanid<nchan;chanid++) {
2170 : for( Int polid=0; polid<npol; polid++ ) {
2171 : if( (itsPSFBeams.getBeam(chanid, polid)).isNull() )
2172 : { itsPSFBeams.setBeam( chanid, polid, maxbeam ); }
2173 : }// end of pol loop
2174 : }// end of chan loop
2175 : }
2176 : */
2177 :
2178 :
2179 : /// For lack of a better place, store this inside the PSF image. To be read later and used to restore
2180 :
2181 718 : ii.setBeams( itsPSFBeams );
2182 : {
2183 718 : LatticeLocker lock1(*(psf()), FileLocker::Write);
2184 718 : psf()->setImageInfo(ii);
2185 718 : psf()->unlock();
2186 718 : }
2187 718 : clock_t end = clock();
2188 718 : os << LogIO::NORMAL << "Time to fit Gaussian to PSF " << double(end - begin) / CLOCKS_PER_SEC << LogIO::POST;
2189 3910 : }// end of make beam set
2190 :
2191 :
2192 0 : ImageBeamSet SIImageStore::getChannelBeamSet(const Int channel){
2193 :
2194 0 : return getChannelSliceBeamSet(channel, channel);
2195 : }
2196 :
2197 :
2198 0 : ImageBeamSet SIImageStore::getChannelSliceBeamSet(const Int begChan, const Int endChan){
2199 0 : ImageBeamSet bs=getBeamSet();
2200 0 : if(bs.shape()[0]==1)
2201 0 : return bs;
2202 0 : if(begChan > endChan || begChan <0)
2203 0 : throw(AipsError("Inconsistent slice of beam in channel requested"));
2204 0 : if(bs.shape()[0] < (endChan-1))
2205 0 : throw(AipsError("beam of channel "+String::toString(endChan)+" does not exist"));
2206 0 : IPosition blc(2, begChan, 0);
2207 0 : IPosition trc(2, endChan, bs.shape()[1]-1);
2208 0 : Matrix<GaussianBeam> sliceBeam=bs.getBeams()(blc, trc);
2209 0 : ImageBeamSet subBeamSet(sliceBeam);
2210 0 : return subBeamSet;
2211 :
2212 0 : }
2213 263 : void SIImageStore::setBeamSet(const ImageBeamSet& bs){
2214 :
2215 263 : itsPSFBeams=bs;
2216 263 : }
2217 :
2218 1052 : ImageBeamSet SIImageStore::getBeamSet(Float psfcutoff)
2219 : {
2220 1052 : IPosition beamshp = itsPSFBeams.shape();
2221 1052 : AlwaysAssert( beamshp.nelements()==2 , AipsError );
2222 1052 : if( beamshp[0]==0 || beamshp[1]==0 ) {makeImageBeamSet(psfcutoff);}
2223 2104 : return itsPSFBeams;
2224 1052 : }
2225 :
2226 728 : void SIImageStore::printBeamSet(Bool verbose)
2227 : {
2228 1456 : LogIO os( LogOrigin("SIImageStore","printBeamSet",WHERE) );
2229 728 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2230 728 : if( itsImageShape[3] == 1 && itsImageShape[2]==1 )
2231 : {
2232 417 : GaussianBeam beam = itsPSFBeams.getBeam();
2233 417 : os << "Beam : " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2234 417 : }
2235 : else
2236 : {
2237 : // per CAS-11415, verbose=True when restoringbeam='common'
2238 311 : if( verbose )
2239 : {
2240 9 : ostringstream oss;
2241 9 : oss<<"Beam";
2242 9 : Int nchan = itsImageShape[3];
2243 189 : for( Int chanid=0; chanid<nchan;chanid++) {
2244 180 : Int polid=0;
2245 : // for( Int polid=0; polid<npol; polid++ ) {
2246 180 : GaussianBeam beam = itsPSFBeams.getBeam( chanid, polid );
2247 180 : oss << " [C" << chanid << "]: " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg";
2248 180 : }//for chanid
2249 9 : os << oss.str() << LogIO::POST;
2250 9 : }
2251 : else
2252 : {
2253 : // TODO : Enable this, when this function doesn't complain about 0 rest freq.
2254 : // or when rest freq is never zero !
2255 : try{
2256 302 : itsPSFBeams.summarize( os, False, itsCoordSys );
2257 : // per CAS-11415 request, not turn on this one (it prints out per-channel beam in each line in the logger)
2258 : //itsPSFBeams.summarize( os, verbose, itsCoordSys );
2259 : }
2260 0 : catch(AipsError &x)
2261 : {
2262 0 : os << LogIO::WARN << "Error while printing (compact) restoring beam summary : " << x.getMesg() << LogIO::POST;
2263 0 : os << "Printing long summary" << LogIO::POST;
2264 :
2265 0 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2266 : //Int npol = itsImageShape[2];
2267 0 : Int nchan = itsImageShape[3];
2268 0 : for( Int chanid=0; chanid<nchan;chanid++) {
2269 0 : Int polid=0;
2270 : // for( Int polid=0; polid<npol; polid++ ) {
2271 0 : GaussianBeam beam = itsPSFBeams.getBeam( chanid, polid );
2272 0 : os << "Beam [C" << chanid << "]: " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2273 : //}//for polid
2274 0 : }//for chanid
2275 0 : }// catch
2276 : }
2277 : }
2278 728 : }
2279 :
2280 : /////////////////////////////// Restore all planes.
2281 :
2282 863 : void SIImageStore::restore(GaussianBeam& rbeam, String& usebeam, uInt term, Float psfcutoff)
2283 : {
2284 1726 : LogIO os( LogOrigin("SIImageStore","restore",WHERE) );
2285 : // << ". Optionally, PB-correct too." << LogIO::POST;
2286 :
2287 863 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2288 863 : Int nx = itsImageShape[0];
2289 863 : Int ny = itsImageShape[1];
2290 863 : Int npol = itsImageShape[2];
2291 863 : Int nchan = itsImageShape[3];
2292 :
2293 : /* if( !hasResidualImage() )
2294 : {
2295 : // Cannot restore without residual/dirty image.
2296 : os << "Cannot restore without residual image" << LogIO::POST;
2297 : return;
2298 : }
2299 : */
2300 :
2301 : //// Get/fill the beamset
2302 863 : IPosition beamset = itsPSFBeams.shape();
2303 863 : AlwaysAssert( beamset.nelements()==2 , AipsError );
2304 863 : if( beamset[0] != nchan || beamset[1] != npol )
2305 : {
2306 :
2307 : // Get PSF Beams....
2308 236 : ImageInfo ii = psf()->imageInfo();
2309 236 : itsPSFBeams = ii.getBeamSet();
2310 236 : itsRestoredBeams=itsPSFBeams;
2311 236 : IPosition beamset2 = itsPSFBeams.shape();
2312 :
2313 236 : AlwaysAssert( beamset2.nelements()==2 , AipsError );
2314 236 : if( beamset2[0] != nchan || beamset2[1] != npol )
2315 : {
2316 : // Make new beams.
2317 0 : os << LogIO::WARN << "Couldn't find pre-computed restoring beams. Re-fitting." << LogIO::POST;
2318 0 : makeImageBeamSet(psfcutoff);
2319 : }
2320 236 : }
2321 :
2322 : // toggle for printing common beam to the log
2323 863 : Bool printcommonbeam(False);
2324 : //// Modify the beamset if needed
2325 : //// if ( rbeam is Null and usebeam=="" ) Don't do anything.
2326 : //// If rbeam is Null but usebeam=='common', calculate a common beam and set 'rbeam'
2327 : //// If rbeam is given (or exists due to 'common'), just use it.
2328 863 : if( rbeam.isNull() && usebeam=="common") {
2329 12 : os << "Getting common beam" << LogIO::POST;
2330 12 : rbeam = CasaImageBeamSet(itsPSFBeams).getCommonBeam();
2331 12 : os << "Common Beam : " << rbeam.getMajor(Unit("arcsec")) << " arcsec, " << rbeam.getMinor(Unit("arcsec"))<< " arcsec, " << rbeam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2332 12 : printcommonbeam=True;
2333 : }
2334 863 : if( !rbeam.isNull() ) {
2335 : /*for( Int chanid=0; chanid<nchan;chanid++) {
2336 : for( Int polid=0; polid<npol; polid++ ) {
2337 : itsPSFBeams.setBeam( chanid, polid, rbeam );
2338 : /// Still need to add the 'common beam' option.
2339 : }//for chanid
2340 : }//for polid
2341 : */
2342 17 : itsRestoredBeams=ImageBeamSet(rbeam);
2343 17 : GaussianBeam beam = itsRestoredBeams.getBeam();
2344 :
2345 : //if commonbeam has not printed in the log
2346 17 : if (!printcommonbeam) {
2347 5 : os << "Common Beam : " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2348 5 : printcommonbeam=True; // to avoid duplicate the info is printed to th elog
2349 : }
2350 17 : }// if rbeam not NULL
2351 : //// Done modifying beamset if needed
2352 :
2353 :
2354 : /// Before restoring, check for an empty model image and don't convolve (but still smooth residuals)
2355 : /// (for CAS-8271 : make output restored image even if .model is zero... )
2356 863 : Bool emptymodel = isModelEmpty();
2357 863 : if(emptymodel) os << LogIO::WARN << "Restoring with an empty model image. Only residuals will be processed to form the output restored image." << LogIO::POST;
2358 :
2359 :
2360 863 : LatticeLocker lock1(*(image(term)), FileLocker::Write);
2361 863 : LatticeLocker lock2(*(residual(term)), FileLocker::Write);
2362 863 : LatticeLocker lock3(*(model(term)), FileLocker::Read);
2363 : //// Use beamset to restore
2364 3684 : for( Int chanid=0; chanid<nchan;chanid++) {
2365 5851 : for( Int polid=0; polid<npol; polid++ ) {
2366 :
2367 3030 : IPosition substart(4,0,0,polid,chanid);
2368 3030 : IPosition substop(4,nx-1,ny-1,polid,chanid);
2369 3030 : Slicer imslice(substart, substop,Slicer::endIsLast);
2370 6060 : SubImage<Float> subRestored( *image(term) , imslice, True );
2371 6060 : SubImage<Float> subModel( *model(term) , imslice, False );
2372 6060 : SubImage<Float> subResidual( *residual(term) , imslice, True );
2373 :
2374 3030 : GaussianBeam beam = itsRestoredBeams.getBeam( chanid, polid );;
2375 : //os << "Common Beam for chan : " << chanid << " : " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2376 : // only print per-chan beam if the common beam is not used for restoring beam
2377 3030 : if(!printcommonbeam) {
2378 2804 : os << "Beam for chan : " << chanid << " : " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2379 : }
2380 :
2381 : try
2382 : {
2383 : // Initialize restored image
2384 3030 : subRestored.set(0.0);
2385 3030 : if( !emptymodel ) {
2386 : // Copy model into it
2387 2702 : subRestored.copyData( LatticeExpr<Float>( subModel ) );
2388 : // Smooth model by beam
2389 2702 : StokesImageUtil::Convolve( subRestored, beam);
2390 : }
2391 : // Add residual image
2392 3030 : if( !rbeam.isNull() || usebeam == "common") // If need to rescale residuals, make a copy and do it.
2393 : {
2394 : // rescaleResolution(chanid, subResidual, beam, itsPSFBeams.getBeam(chanid, polid));
2395 452 : TempImage<Float> tmpSubResidualCopy( IPosition(4,nx,ny,1,1), subResidual.coordinates());
2396 226 : tmpSubResidualCopy.copyData( subResidual );
2397 :
2398 226 : rescaleResolution(chanid, tmpSubResidualCopy, beam, itsPSFBeams.getBeam(chanid, polid));
2399 226 : subRestored.copyData( LatticeExpr<Float>( subRestored + tmpSubResidualCopy ) );
2400 226 : }
2401 : else// if no need to rescale residuals, just add the residuals.
2402 : {
2403 :
2404 :
2405 2804 : subRestored.copyData( LatticeExpr<Float>( subRestored + subResidual ) );
2406 :
2407 : }
2408 :
2409 : }
2410 0 : catch(AipsError &x)
2411 : {
2412 0 : throw( AipsError("Restoration Error in chan" + String::toString(chanid) + ":pol" + String::toString(polid) + " : " + x.getMesg() ) );
2413 0 : }
2414 :
2415 3030 : }// end of pol loop
2416 : }// end of chan loop
2417 :
2418 : try
2419 : {
2420 : //MSK//
2421 863 : if( hasPB() )
2422 : {
2423 841 : if( (image(term)->getDefaultMask()).matches("mask0") ) removeMask( image(term) );
2424 841 : copyMask(residual(term),image(term));
2425 : }
2426 :
2427 : // if(hasPB()){copyMask(residual(term),image(term));}
2428 863 : ImageInfo iminf = image(term)->imageInfo();
2429 : //iminf.setBeams( itsRestoredBeams);
2430 :
2431 863 : os << "Beam Set : Single beam : " << itsRestoredBeams.hasSingleBeam() << " Multi-beam : " << itsRestoredBeams.hasMultiBeam() << LogIO::DEBUG2;
2432 :
2433 863 : iminf.removeRestoringBeam();
2434 :
2435 863 : if( itsRestoredBeams.hasSingleBeam() )
2436 606 : { iminf.setRestoringBeam( itsRestoredBeams.getBeam() );}
2437 : else
2438 257 : {iminf.setBeams( itsRestoredBeams);}
2439 :
2440 863 : image(term)->setImageInfo(iminf);
2441 :
2442 863 : }
2443 0 : catch(AipsError &x)
2444 : {
2445 0 : throw( AipsError("Restoration Error : " + x.getMesg() ) );
2446 0 : }
2447 :
2448 863 : }// end of restore()
2449 :
2450 0 : GaussianBeam SIImageStore::findGoodBeam()
2451 : {
2452 0 : LogIO os( LogOrigin("SIImageStore","findGoodBeam",WHERE) );
2453 0 : IPosition beamshp = itsPSFBeams.shape();
2454 0 : AlwaysAssert( beamshp.nelements()==2 , AipsError );
2455 :
2456 : /*
2457 : if( beamshp[0] != nchan || beamshp[1] != npol )
2458 : {
2459 : // Make new beams.
2460 : os << LogIO::WARN << "Couldn't find pre-computed restoring beams. Re-fitting." << LogIO::POST;
2461 : makeImageBeamSet();
2462 : }
2463 : */
2464 :
2465 0 : Vector<Float> areas(beamshp[0]*beamshp[1]);
2466 0 : Vector<Float> axrat(beamshp[0]*beamshp[1]);
2467 0 : areas=0.0; axrat=1.0;
2468 0 : Vector<Bool> flags( areas.nelements() );
2469 0 : flags=False;
2470 :
2471 0 : Int cnt=0;
2472 0 : for( Int chanid=0; chanid<beamshp[0];chanid++) {
2473 0 : for( Int polid=0; polid<beamshp[1]; polid++ ) {
2474 0 : GaussianBeam beam = itsPSFBeams(chanid, polid);
2475 0 : if( !beam.isNull() && beam.getMajor(Unit("arcsec"))>1.1e-06 ) // larger than default filler beam.
2476 : {
2477 0 : areas[cnt] = beam.getArea( Unit("arcsec2") );
2478 0 : axrat[cnt] = beam.getMajor( Unit("arcsec") ) / beam.getMinor( Unit("arcsec") );
2479 : }
2480 : else {
2481 0 : flags[cnt] = True;
2482 : }
2483 0 : cnt++;
2484 0 : }//for chanid
2485 : }//for polid
2486 :
2487 0 : Vector<Float> fit( areas.nelements() );
2488 0 : Vector<Float> fitaxr( areas.nelements() );
2489 0 : for (Int loop=0;loop<5;loop++) {
2490 : /// Filter on outliers in PSF beam area
2491 0 : lineFit( areas, flags, fit, 0, areas.nelements()-1 );
2492 0 : Float sd = calcStd( areas , flags, fit );
2493 0 : for (uInt i=0;i<areas.nelements();i++) {
2494 0 : if( fabs( areas[i] - fit[i] ) > 3*sd ) flags[i]=True;
2495 : }
2496 : /// Filter on outliers in PSF axial ratio
2497 0 : lineFit( axrat, flags, fitaxr, 0, areas.nelements()-1 );
2498 0 : Float sdaxr = calcStd( axrat , flags, fitaxr );
2499 0 : for (uInt i=0;i<areas.nelements();i++) {
2500 0 : if( fabs( axrat[i] - fitaxr[i] ) > 3*sdaxr ) flags[i]=True;
2501 : }
2502 : }
2503 : // cout << "Original areas : " << areas << endl;
2504 : // cout << "Original axrats : " << axrat << endl;
2505 : // cout << "Flags : " << flags << endl;
2506 :
2507 : // Find max area good beam.
2508 0 : GaussianBeam goodbeam;
2509 0 : Int cid=0,pid=0;
2510 0 : Float maxval=0.0;
2511 0 : cnt=0;
2512 0 : for( Int chanid=0; chanid<beamshp[0];chanid++) {
2513 0 : for( Int polid=0; polid<beamshp[1]; polid++ ) {
2514 0 : if( flags[cnt] == False ){
2515 0 : if( areas[cnt] > maxval ) {maxval = areas[cnt]; goodbeam = itsPSFBeams.getBeam(chanid,polid);cid=chanid;pid=polid;}
2516 : }
2517 0 : cnt++;
2518 : }//polid
2519 : }//chanid
2520 :
2521 0 : os << "Picking common beam from C"<<cid<<":P"<<pid<<" : " << goodbeam.getMajor(Unit("arcsec")) << " arcsec, " << goodbeam.getMinor(Unit("arcsec"))<< " arcsec, " << goodbeam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2522 :
2523 0 : Bool badbeam=False;
2524 0 : for(uInt i=0;i<flags.nelements();i++){if(flags[i]==True) badbeam=True;}
2525 :
2526 0 : if( badbeam == True )
2527 : {
2528 0 : os << "(Ignored beams from :";
2529 0 : cnt=0;
2530 0 : for( Int chanid=0; chanid<beamshp[0];chanid++) {
2531 0 : for( Int polid=0; polid<beamshp[1]; polid++ ) {
2532 0 : if( flags[cnt] == True ){
2533 0 : os << " C"<<chanid<<":P"<<polid;
2534 : }
2535 0 : cnt++;
2536 : }//polid
2537 : }//chanid
2538 0 : os << " as outliers either by area or by axial ratio)" << LogIO::POST;
2539 : }
2540 :
2541 :
2542 : /*
2543 : // Replace 'bad' psfs with the chosen one.
2544 : if( goodbeam.isNull() ) {
2545 : os << LogIO::WARN << "No planes have non-zero restoring beams. Forcing artificial 1.0arcsec beam." << LogIO::POST;
2546 : Quantity majax(1.0,"arcsec"),minax(1.0,"arcsec"),pa(0.0,"deg");
2547 : goodbeam.setMajorMinor(majax,minax);
2548 : goodbeam.setPA(pa);
2549 : }
2550 : else {
2551 : cnt=0;
2552 : for( Int chanid=0; chanid<nchan;chanid++) {
2553 : for( Int polid=0; polid<npol; polid++ ) {
2554 : if( flags[cnt]==True )
2555 : { itsPSFBeams.setBeam( chanid, polid, goodbeam ); }
2556 : cnt++;
2557 : }// end of pol loop
2558 : }// end of chan loop
2559 : }
2560 : */
2561 :
2562 0 : return goodbeam;
2563 0 : }// end of findGoodBeam
2564 :
2565 : ///////////////////////// Funcs to calculate robust mean and fit, taking into account 'flagged' points.
2566 0 : void SIImageStore :: lineFit(Vector<Float> &data, Vector<Bool> &flag, Vector<Float> &fit, uInt lim1, uInt lim2)
2567 : {
2568 0 : float Sx = 0, Sy = 0, Sxx = 0, Sxy = 0, S = 0, a, b, sd, mn;
2569 :
2570 0 : mn = calcMean(data, flag);
2571 0 : sd = calcStd (data, flag, mn);
2572 :
2573 0 : for (uInt i = lim1; i <= lim2; i++)
2574 : {
2575 0 : if (flag[i] == False) // if unflagged
2576 : {
2577 0 : S += 1 / (sd * sd);
2578 0 : Sx += i / (sd * sd);
2579 0 : Sy += data[i] / (sd * sd);
2580 0 : Sxx += (i * i) / (sd * sd);
2581 0 : Sxy += (i * data[i]) / (sd * sd);
2582 : }
2583 : }
2584 0 : a = (Sxx * Sy - Sx * Sxy) / (S * Sxx - Sx * Sx);
2585 0 : b = (S * Sxy - Sx * Sy) / (S * Sxx - Sx * Sx);
2586 :
2587 0 : for (uInt i = lim1; i <= lim2; i++)
2588 0 : fit[i] = a + b * i;
2589 :
2590 0 : }
2591 : /* Calculate the MEAN of 'vect' ignoring values flagged in 'flag' */
2592 0 : Float SIImageStore :: calcMean(Vector<Float> &vect, Vector<Bool> &flag)
2593 : {
2594 0 : Float mean=0;
2595 0 : Int cnt=0;
2596 0 : for(uInt i=0;i<vect.nelements();i++)
2597 0 : if(flag[i]==False)
2598 : {
2599 0 : mean += vect[i];
2600 0 : cnt++;
2601 : }
2602 0 : if(cnt==0) cnt=1;
2603 0 : return mean/cnt;
2604 : }
2605 0 : Float SIImageStore :: calcStd(Vector<Float> &vect, Vector<Bool> &flag, Vector<Float> &fit)
2606 : {
2607 0 : Float std=0;
2608 0 : uInt n=0,cnt=0;
2609 0 : n = vect.nelements() < fit.nelements() ? vect.nelements() : fit.nelements();
2610 0 : for(uInt i=0;i<n;i++)
2611 0 : if(flag[i]==False)
2612 : {
2613 0 : cnt++;
2614 0 : std += (vect[i]-fit[i])*(vect[i]-fit[i]);
2615 : }
2616 0 : if(cnt==0) cnt=1;
2617 0 : return sqrt(std/cnt);
2618 :
2619 : }
2620 0 : Float SIImageStore :: calcStd(Vector<Float> &vect, Vector<Bool> &flag, Float mean)
2621 : {
2622 0 : Float std=0;
2623 0 : uInt cnt=0;
2624 0 : for(uInt i=0;i<vect.nelements();i++)
2625 0 : if(flag[i]==False)
2626 : {
2627 0 : cnt++;
2628 0 : std += (vect[i]-mean)*(vect[i]-mean);
2629 : }
2630 0 : return sqrt(std/cnt);
2631 : }
2632 :
2633 : ///////////////////////// End of Funcs to calculate robust mean and fit.
2634 :
2635 :
2636 :
2637 : /*
2638 : GaussianBeam SIImageStore::restorePlane()
2639 : {
2640 :
2641 : LogIO os( LogOrigin("SIImageStore","restorePlane",WHERE) );
2642 : // << ". Optionally, PB-correct too." << LogIO::POST;
2643 :
2644 : Bool validbeam=False;
2645 : GaussianBeam beam;
2646 : try
2647 : {
2648 : // Fit a Gaussian to the PSF.
2649 : beam = getPSFGaussian();
2650 : validbeam = True;
2651 : }
2652 : catch(AipsError &x)
2653 : {
2654 : os << LogIO::WARN << "Beam fit error : " + x.getMesg() << LogIO::POST;
2655 : }
2656 :
2657 : try
2658 : {
2659 : if( validbeam==True )
2660 : {
2661 : //os << "[" << itsImageName << "] " ; // Add when parent image name is available.
2662 : //os << "Restore with beam : " << beam.getMajor(Unit("arcmin")) << " arcmin, " << beam.getMinor(Unit("arcmin"))<< " arcmin, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2663 :
2664 : // Initialize restored image
2665 : image()->set(0.0);
2666 : // Copy model into it
2667 : image()->copyData( LatticeExpr<Float>( *(model()) ) );
2668 : // Smooth model by beam
2669 : StokesImageUtil::Convolve( *(image()), beam);
2670 : // Add residual image
2671 : image()->copyData( LatticeExpr<Float>( *(image()) + *(residual()) ) );
2672 :
2673 : // Set restoring beam into the image
2674 : ImageInfo ii = image()->imageInfo();
2675 : //ii.setRestoringBeam(beam);
2676 : ii.setBeams(beam);
2677 : image()->setImageInfo(ii);
2678 : }
2679 : }
2680 : catch(AipsError &x)
2681 : {
2682 : throw( AipsError("Restoration Error : " + x.getMesg() ) );
2683 : }
2684 :
2685 : return beam;
2686 :
2687 : }
2688 : */
2689 :
2690 45 : void SIImageStore::pbcor(uInt term)
2691 : {
2692 :
2693 90 : LogIO os( LogOrigin("SIImageStore","pbcor",WHERE) );
2694 :
2695 45 : if( !hasRestored() || !hasPB() )
2696 : {
2697 : // Cannot pbcor without restored image and pb
2698 0 : os << LogIO::WARN << "Cannot pbcor without restored image and pb" << LogIO::POST;
2699 0 : return;
2700 : }
2701 45 : LatticeLocker lock1(*(image(term)), FileLocker::Write);
2702 45 : LatticeLocker lock2(*(pb(term)), FileLocker::Write);
2703 45 : LatticeLocker lock3(*(imagepbcor(term)), FileLocker::Write);
2704 :
2705 90 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2706 : {
2707 164 : for(Int chan=0; chan<itsImageShape[3]; chan++)
2708 : {
2709 :
2710 119 : CountedPtr<ImageInterface<Float> > restoredsubim=makeSubImage(0,1,
2711 119 : chan, itsImageShape[3],
2712 119 : pol, itsImageShape[2],
2713 238 : *image(term) );
2714 119 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
2715 119 : chan, itsImageShape[3],
2716 119 : pol, itsImageShape[2],
2717 238 : *pb() );
2718 :
2719 119 : CountedPtr<ImageInterface<Float> > pbcorsubim=makeSubImage(0,1,
2720 119 : chan, itsImageShape[3],
2721 119 : pol, itsImageShape[2],
2722 238 : *imagepbcor(term) );
2723 :
2724 :
2725 119 : LatticeExprNode pbmax( max( *pbsubim ) );
2726 119 : Float pbmaxval = pbmax.getFloat();
2727 :
2728 119 : if( pbmaxval<=0.0 )
2729 : {
2730 0 : os << LogIO::WARN << "Skipping PBCOR for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;
2731 0 : pbcorsubim->set(0.0);
2732 : }
2733 : else
2734 : {
2735 :
2736 238 : LatticeExpr<Float> thepbcor( iif( *(pbsubim) > 0.0 , (*(restoredsubim))/(*(pbsubim)) , 0.0 ) );
2737 119 : pbcorsubim->copyData( thepbcor );
2738 :
2739 119 : }// if not zero
2740 119 : }//chan
2741 : }//pol
2742 :
2743 : // Copy over the PB mask.
2744 45 : if((imagepbcor(term)->getDefaultMask()=="") && hasPB())
2745 27 : {copyMask(pb(),imagepbcor(term));}
2746 :
2747 : // Set restoring beam info
2748 45 : ImageInfo iminf = image(term)->imageInfo();
2749 : //iminf.setBeams( itsRestoredBeams );
2750 45 : imagepbcor(term)->setImageInfo(iminf);
2751 :
2752 45 : }// end pbcor
2753 :
2754 0 : Matrix<Float> SIImageStore::getSumWt(ImageInterface<Float>& target)
2755 : {
2756 0 : Record miscinfo = target.miscInfo();
2757 :
2758 0 : Matrix<Float> sumwt;
2759 0 : sumwt.resize();
2760 0 : if( miscinfo.isDefined("sumwt")
2761 0 : && (miscinfo.dataType("sumwt")==TpArrayFloat || miscinfo.dataType("sumwt")==TpArrayDouble ) )
2762 0 : { miscinfo.get( "sumwt" , sumwt ); }
2763 0 : else { sumwt.resize( IPosition(2, target.shape()[2], target.shape()[3] ) ); sumwt = 1.0; }
2764 :
2765 0 : return sumwt;
2766 0 : }
2767 :
2768 0 : void SIImageStore::setSumWt(ImageInterface<Float>& target, Matrix<Float>& sumwt)
2769 : {
2770 0 : Record miscinfo = target.miscInfo();
2771 0 : miscinfo.define("sumwt", sumwt);
2772 0 : LatticeLocker lock1(target, FileLocker::Write);
2773 0 : target.setMiscInfo( miscinfo );
2774 0 : }
2775 :
2776 :
2777 4405 : Bool SIImageStore::getUseWeightImage(ImageInterface<Float>& target)
2778 : {
2779 4405 : Record miscinfo = target.miscInfo();
2780 : Bool useweightimage;
2781 4405 : if( miscinfo.isDefined("useweightimage") && miscinfo.dataType("useweightimage")==TpBool )
2782 4405 : { miscinfo.get( "useweightimage", useweightimage ); }
2783 0 : else { useweightimage = False; }
2784 :
2785 4405 : return useweightimage;
2786 4405 : }
2787 : /*
2788 : Bool SIImageStore::getUseWeightImage()
2789 : {
2790 : if( ! itsParentSumWt )
2791 : return False;
2792 : else
2793 : return getUseWeightImage( *itsParentSumWt );
2794 : }
2795 : */
2796 15809 : void SIImageStore::setUseWeightImage(ImageInterface<Float>& target, Bool useweightimage)
2797 : {
2798 15809 : Record miscinfo = target.miscInfo();
2799 15809 : miscinfo.define("useweightimage", useweightimage);
2800 15809 : LatticeLocker lock1(target, FileLocker::Write);
2801 15809 : target.setMiscInfo( miscinfo );
2802 15809 : }
2803 :
2804 :
2805 :
2806 1873 : Bool SIImageStore::divideImageByWeightVal( ImageInterface<Float>& target )
2807 : {
2808 :
2809 1873 : Array<Float> lsumwt;
2810 1873 : sumwt()->get( lsumwt , False ); // For MT, this will always pick the zeroth sumwt, which it should.
2811 1873 : LatticeLocker lock2(target, FileLocker::Write);
2812 :
2813 1873 : IPosition imshape = target.shape();
2814 :
2815 : //cerr << " SumWt : " << lsumwt << " sumwtshape : " << lsumwt.shape() << " image shape : " << imshape << endl;
2816 :
2817 1873 : AlwaysAssert( lsumwt.shape()[2] == imshape[2] , AipsError ); // polplanes
2818 1873 : AlwaysAssert( lsumwt.shape()[3] == imshape[3] , AipsError ); // chanplanes
2819 :
2820 1873 : Bool div=False; // flag to signal if division actually happened, or weights are all 1.0
2821 :
2822 3970 : for(Int pol=0; pol<lsumwt.shape()[2]; pol++)
2823 : {
2824 8312 : for(Int chan=0; chan<lsumwt.shape()[3]; chan++)
2825 : {
2826 6215 : IPosition pos(4,0,0,pol,chan);
2827 6215 : if( lsumwt(pos) != 1.0 )
2828 : {
2829 : // SubImage<Float>* subim=makePlane( chan, True ,pol, True, target );
2830 : std::shared_ptr<ImageInterface<Float> > subim=makeSubImage(0,1,
2831 6215 : chan, lsumwt.shape()[3],
2832 6215 : pol, lsumwt.shape()[2],
2833 12430 : target );
2834 6215 : if ( lsumwt(pos) > 1e-07 ) {
2835 12028 : LatticeExpr<Float> le( (*subim)/lsumwt(pos) );
2836 6014 : subim->copyData( le );
2837 6014 : }
2838 : else {
2839 201 : subim->set(0.0);
2840 : }
2841 6215 : div=True;
2842 6215 : }
2843 6215 : }
2844 : }
2845 :
2846 : // if( div==True ) cout << "Div image by sumwt : " << lsumwt << endl;
2847 : // else cout << "Already normalized" << endl;
2848 :
2849 : // lsumwt = 1.0; setSumWt( target , lsumwt );
2850 :
2851 1873 : return div;
2852 1873 : }
2853 :
2854 899 : void SIImageStore::normPSF(Int term)
2855 : {
2856 :
2857 1919 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2858 : {
2859 4133 : for(Int chan=0; chan<itsImageShape[3]; chan++)
2860 : {
2861 : /// IPosition center(4,itsImageShape[0]/2,itsImageShape[1]/2,pol,chan);
2862 :
2863 : std::shared_ptr<ImageInterface<Float> > subim=makeSubImage(0,1,
2864 6226 : chan, itsImageShape[3],
2865 3113 : pol, itsImageShape[2],
2866 6226 : (*psf(term)) );
2867 :
2868 : std::shared_ptr<ImageInterface<Float> > subim0=makeSubImage(0,1,
2869 6226 : chan, itsImageShape[3],
2870 3113 : pol, itsImageShape[2],
2871 6226 : (*psf(0)) );
2872 :
2873 3113 : LatticeLocker lock1(*(subim), FileLocker::Write);
2874 :
2875 3113 : LatticeExprNode themax( max(*(subim0)) );
2876 3113 : Float maxim = themax.getFloat();
2877 :
2878 3113 : if ( maxim > 1e-07 )
2879 : {
2880 6024 : LatticeExpr<Float> normed( (*(subim)) / maxim );
2881 3012 : subim->copyData( normed );
2882 3012 : }
2883 : else
2884 : {
2885 101 : subim->set(0.0);
2886 : }
2887 3113 : }//chan
2888 : }//pol
2889 :
2890 899 : }
2891 :
2892 590 : void SIImageStore::calcSensitivity()
2893 : {
2894 1180 : LogIO os( LogOrigin("SIImageStore","calcSensitivity",WHERE) );
2895 :
2896 590 : Array<Float> lsumwt;
2897 590 : sumwt()->get( lsumwt , False ); // For MT, this will always pick the zeroth sumwt, which it should.
2898 :
2899 590 : IPosition shp( lsumwt.shape() );
2900 : //cout << "Sumwt shape : " << shp << " : " << lsumwt << endl;
2901 : //AlwaysAssert( shp.nelements()==4 , AipsError );
2902 :
2903 590 : os << "[" << itsImageName << "] Theoretical sensitivity (Jy/bm):" ;
2904 :
2905 590 : IPosition it(4,0,0,0,0);
2906 1187 : for ( it[0]=0; it[0]<shp[0]; it[0]++)
2907 1220 : for ( it[1]=0; it[1]<shp[1]; it[1]++)
2908 1349 : for ( it[2]=0; it[2]<shp[2]; it[2]++)
2909 3545 : for ( it[3]=0; it[3]<shp[3]; it[3]++)
2910 : {
2911 2819 : if( shp[0]>1 ){os << "f"<< it[0]+(it[1]*shp[0]) << ":" ;}
2912 2819 : if( shp[3]>1 ) { os << "c"<< it[3] << ":"; }
2913 2819 : if( shp[2]>1 ) { os << "p"<< it[2]<< ":" ; }
2914 2819 : if( lsumwt( it ) > 1e-07 )
2915 : {
2916 2698 : os << 1.0/(sqrt(lsumwt(it))) << " " ;
2917 : }
2918 : else
2919 : {
2920 121 : os << "none" << " ";
2921 : }
2922 : }
2923 :
2924 590 : os << LogIO::POST;
2925 :
2926 : // Array<Float> sens = 1.0/sqrtsumwt;
2927 :
2928 :
2929 590 : }
2930 :
2931 0 : Double SIImageStore::calcFractionalBandwidth()
2932 : {
2933 0 : throw(AipsError("calcFractionalBandwidth is not implemented for SIImageStore. Only SIImageStoreMultiTerm"));
2934 : }
2935 :
2936 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2937 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2938 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2939 : /////////////////// Utility Functions to gather statistics on the imagestore.
2940 :
2941 8452 : Float SIImageStore::getPeakResidual()
2942 : {
2943 16904 : LogIO os( LogOrigin("SIImageStore","getPeakResidual",WHERE) );
2944 8452 : LatticeLocker lock2 (*(residual()), FileLocker::Read);
2945 16904 : ArrayLattice<Bool> pixelmask(residual()->getMask());
2946 16904 : LatticeExpr<Float> resd(iif(pixelmask,abs(*residual()),0));
2947 : //LatticeExprNode pres( max(abs( *residual() ) ));
2948 8452 : LatticeExprNode pres( max(resd) );
2949 8452 : Float maxresidual = pres.getFloat();
2950 :
2951 8452 : return maxresidual;
2952 8452 : }
2953 :
2954 4910 : Float SIImageStore::getPeakResidualWithinMask()
2955 : {
2956 9820 : LogIO os( LogOrigin("SIImageStore","getPeakResidualWithinMask",WHERE) );
2957 : //Float minresmask, maxresmask, minres, maxres;
2958 : //findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
2959 :
2960 9820 : ArrayLattice<Bool> pixelmask(residual()->getMask());
2961 : //findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
2962 9820 : LatticeExpr<Float> resd(iif(((*mask()) > 0) && pixelmask ,abs((*residual())),0));
2963 4910 : LatticeExprNode pres( max(resd) );
2964 4910 : Float maxresidual = pres.getFloat();
2965 : //Float maxresidual=max( abs(maxresmask), abs(minresmask) );
2966 4910 : return maxresidual;
2967 4910 : }
2968 :
2969 : // Calculate the total model flux
2970 5169 : Float SIImageStore::getModelFlux(uInt term)
2971 : {
2972 : // LogIO os( LogOrigin("SIImageStore","getModelFlux",WHERE) );
2973 5169 : LatticeLocker lock2 (*(model(term)), FileLocker::Read);
2974 10334 : LatticeExprNode mflux( sum( *model(term) ) );
2975 5167 : Float modelflux = mflux.getFloat();
2976 : // Float modelflux = sum( model(term)->get() );
2977 :
2978 5167 : return modelflux;
2979 5167 : }
2980 :
2981 : // Check for non-zero model (this is different from getting model flux, for derived SIIMMT)
2982 1950 : Bool SIImageStore::isModelEmpty()
2983 : {
2984 1950 : if( !itsModel && (! hasModel()) ) return True;
2985 1286 : else return ( fabs( getModelFlux(0) ) < 1e-08 );
2986 : }
2987 :
2988 :
2989 263 : void SIImageStore::setPSFSidelobeLevel(const Float level){
2990 :
2991 263 : itsPSFSideLobeLevel=level;
2992 263 : }
2993 : // Calculate the PSF sidelobe level...
2994 2829 : Float SIImageStore::getPSFSidelobeLevel()
2995 : {
2996 5658 : LogIO os( LogOrigin("SIImageStore","getPSFSidelobeLevel",WHERE) );
2997 : //cerr << "*****PSF sidelobe "<< itsPSFSideLobeLevel << endl;
2998 : /// Calculate only once, store and return for all subsequent calls.
2999 2829 : if( itsPSFSideLobeLevel == 0.0 )
3000 : {
3001 :
3002 644 : ImageBeamSet thebeams = getBeamSet();
3003 644 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
3004 :
3005 : //------------------------------------------------------------
3006 644 : IPosition oneplaneshape( itsImageShape );
3007 644 : AlwaysAssert( oneplaneshape.nelements()==4, AipsError );
3008 644 : oneplaneshape[2]=1; oneplaneshape[3]=1;
3009 644 : TempImage<Float> psfbeam( oneplaneshape, itsCoordSys );
3010 :
3011 : // In a loop through channels, subtract out or mask out the main lobe
3012 644 : Float allmin=0.0, allmax=0.0;
3013 1397 : for(Int pol=0; pol<itsImageShape[2]; pol++)
3014 : {
3015 3383 : for(Int chan=0; chan<itsImageShape[3]; chan++)
3016 : {
3017 : std::shared_ptr<ImageInterface<Float> > onepsf=makeSubImage(0,1,
3018 5260 : chan, itsImageShape[3],
3019 2630 : pol, itsImageShape[2],
3020 5260 : (*psf()) );
3021 :
3022 :
3023 2630 : GaussianBeam beam = thebeams.getBeam( chan, pol );
3024 2630 : Vector<Float> abeam(3); // Holds bmaj, bmin, pa in asec, asec, deg
3025 2630 : abeam[0] = beam.getMajor().get("arcsec").getValue() * C::arcsec;
3026 2630 : abeam[1] = beam.getMinor().get("arcsec").getValue() * C::arcsec;
3027 2630 : abeam[2] = (beam.getPA().get("deg").getValue() + 90.0)* C::degree;
3028 :
3029 : //cout << "Beam : " << abeam << endl;
3030 :
3031 2630 : StokesImageUtil::MakeGaussianPSF( psfbeam, abeam, False);
3032 :
3033 : // storeImg( String("psfbeam.im"), psfbeam );
3034 :
3035 : //Subtract from PSF plane
3036 5260 : LatticeExpr<Float> delobed( (*onepsf) - psfbeam );
3037 :
3038 : // For debugging
3039 : //onepsf->copyData( delobed );
3040 :
3041 : //Calc max and min and accumulate across channels.
3042 :
3043 2630 : LatticeExprNode minval_le( min( *onepsf ) );
3044 2630 : LatticeExprNode maxval_le( max( delobed ) );
3045 :
3046 2630 : Float minval = minval_le.getFloat();
3047 2630 : Float maxval = maxval_le.getFloat();
3048 :
3049 2630 : if( minval < allmin ) allmin = minval;
3050 2630 : if( maxval > allmax ) allmax = maxval;
3051 :
3052 : //cout << "Chan : " << chan << " minval : " << minval << " maxval : " << maxval << endl;
3053 :
3054 2630 : }//chan
3055 : }//pol
3056 :
3057 : //------------------------------------------------------------
3058 :
3059 644 : itsPSFSideLobeLevel = max( fabs(allmin), fabs(allmax) );
3060 :
3061 : //os << "PSF min : " << allmin << " max : " << allmax << " psfsidelobelevel : " << itsPSFSideLobeLevel << LogIO::POST;
3062 :
3063 644 : }// if changed.
3064 :
3065 : // LatticeExprNode psfside( min( *psf() ) );
3066 : // itsPSFSideLobeLevel = fabs( psfside.getFloat() );
3067 :
3068 : //cout << "PSF sidelobe level : " << itsPSFSideLobeLevel << endl;
3069 2829 : return itsPSFSideLobeLevel;
3070 2829 : }
3071 :
3072 0 : void SIImageStore::findMinMax(const Array<Float>& lattice,
3073 : const Array<Float>& mask,
3074 : Float& minVal, Float& maxVal,
3075 : Float& minValMask, Float& maxValMask)
3076 : {
3077 0 : IPosition posmin(lattice.shape().nelements(), 0);
3078 0 : IPosition posmax(lattice.shape().nelements(), 0);
3079 :
3080 0 : if( sum(mask) <1e-06 ) {minValMask=0.0; maxValMask=0.0;}
3081 0 : else { minMaxMasked(minValMask, maxValMask, posmin, posmax, lattice,mask); }
3082 :
3083 0 : minMax( minVal, maxVal, posmin, posmax, lattice );
3084 0 : }
3085 :
3086 113 : Array<Double> SIImageStore::calcRobustRMS(Array<Double>& mdns, const Float pbmasklevel, const Bool fastcalc)
3087 : {
3088 226 : LogIO os( LogOrigin("SIImageStore","calcRobustRMS",WHERE) );
3089 113 : Record* regionPtr=0;
3090 113 : String LELmask("");
3091 113 : LatticeLocker lockres (*(residual()), FileLocker::Read);
3092 226 : ArrayLattice<Bool> pbmasklat(residual()->shape());
3093 113 : pbmasklat.set(False);
3094 113 : LatticeExpr<Bool> pbmask(pbmasklat);
3095 113 : if (hasPB()) {
3096 : // set bool mask: False = masked
3097 113 : pbmask = LatticeExpr<Bool> (iif(*pb() > pbmasklevel, True, False));
3098 113 : os << LogIO::DEBUG1 << "pbmask to be attached===> nfalse(pbmask.getMask())="<<nfalse(pbmask.getMask())<<endl;
3099 : }
3100 :
3101 :
3102 113 : Record thestats;
3103 113 : if (fastcalc) { // older calculation
3104 : // need to apply pbmask if present to be consistent between fastcalc = true and false.
3105 : //TempImage<Float>* tempRes = new TempImage<Float>(residual()->shape(), residual()->coordinates(), SDMaskHandler::memoryToUse());
3106 222 : auto tempRes = std::shared_ptr<TempImage<Float>>(new TempImage<Float>(residual()->shape(), residual()->coordinates(), SDMaskHandler::memoryToUse()));
3107 111 : tempRes->copyData(*residual());
3108 111 : tempRes->attachMask(pbmask);
3109 : //thestats = SDMaskHandler::calcImageStatistics(*residual(), LELmask, regionPtr, True);
3110 111 : thestats = SDMaskHandler::calcImageStatistics(*tempRes, LELmask, regionPtr, True);
3111 111 : }
3112 : else { // older way to calculate
3113 : // use the new statistic calculation algorithm
3114 2 : Vector<Bool> dummyvec;
3115 : // TT: 2018.08.01 using revised version (the older version of this is renameed to calcRobustImageStatisticsOld)
3116 2 : thestats = SDMaskHandler::calcRobustImageStatistics(*residual(), *mask(), pbmask, LELmask, regionPtr, True, dummyvec);
3117 2 : }
3118 :
3119 :
3120 : /***
3121 : ImageStatsCalculator imcalc( residual(), regionPtr, LELmask, False);
3122 :
3123 : Vector<Int> axes(2);
3124 : axes[0] = 0;
3125 : axes[1] = 1;
3126 : imcalc.setAxes(axes);
3127 : imcalc.setRobust(True);
3128 : Record thestats = imcalc.statistics();
3129 : //cout<<"thestats="<<thestats<<endl;
3130 : ***/
3131 :
3132 : //Array<Double>rmss, mads, mdns;
3133 113 : Array<Double>rmss, mads;
3134 : //thestats.get(RecordFieldId("max"), maxs);
3135 113 : thestats.get(RecordFieldId("rms"), rmss);
3136 113 : thestats.get(RecordFieldId("medabsdevmed"), mads);
3137 113 : thestats.get(RecordFieldId("median"), mdns);
3138 :
3139 : //os << LogIO::DEBUG1 << "Max : " << maxs << LogIO::POST;
3140 113 : os << LogIO::DEBUG1 << "RMS : " << rmss << LogIO::POST;
3141 113 : os << LogIO::DEBUG1 << "MAD : " << mads << LogIO::POST;
3142 :
3143 : // this for the new noise calc
3144 : //return mdns+mads*1.4826;
3145 : // this is the old noise calc
3146 226 : return mads*1.4826;
3147 113 : }
3148 :
3149 0 : void SIImageStore::printImageStats()
3150 : {
3151 0 : LogIO os( LogOrigin("SIImageStore","printImageStats",WHERE) );
3152 0 : Float minresmask=0, maxresmask=0, minres=0, maxres=0;
3153 : // findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
3154 0 : LatticeLocker lock1 (*(residual()), FileLocker::Read);
3155 0 : ArrayLattice<Bool> pixelmask(residual()->getMask());
3156 0 : if(hasMask())
3157 : {
3158 0 : LatticeLocker lock2 (*(mask()), FileLocker::Read);
3159 0 : findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
3160 0 : }
3161 : else
3162 : {
3163 : //LatticeExprNode pres( max( *residual() ) );
3164 0 : LatticeExprNode pres( max( iif(pixelmask,*residual(),0) ) );
3165 0 : maxres = pres.getFloat();
3166 : //LatticeExprNode pres2( min( *residual() ) );
3167 0 : LatticeExprNode pres2( min( iif(pixelmask,*residual(),0) ) );
3168 0 : minres = pres2.getFloat();
3169 0 : }
3170 :
3171 0 : os << "[" << itsImageName << "]" ;
3172 0 : os << " Peak residual (max,min) " ;
3173 0 : if( minresmask!=0.0 || maxresmask!=0.0 )
3174 0 : { os << "within mask : (" << maxresmask << "," << minresmask << ") "; }
3175 0 : os << "over full image : (" << maxres << "," << minres << ")" << LogIO::POST;
3176 :
3177 0 : os << "[" << itsImageName << "] Total Model Flux : " << getModelFlux() << LogIO::POST;
3178 :
3179 :
3180 0 : Record* regionPtr=0;
3181 0 : String LELmask("");
3182 0 : Record thestats = SDMaskHandler::calcImageStatistics(*residual(), LELmask, regionPtr, True);
3183 0 : Array<Double> maxs, mins;
3184 0 : thestats.get(RecordFieldId("max"), maxs);
3185 0 : thestats.get(RecordFieldId("min"), mins);
3186 : //os << LogIO::DEBUG1 << "Max : " << maxs << LogIO::POST;
3187 : //os << LogIO::DEBUG1 << "Min : " << mins << LogIO::POST;
3188 :
3189 0 : }
3190 :
3191 : // Calculate the total model flux
3192 5717 : Float SIImageStore::getMaskSum()
3193 : {
3194 11434 : LogIO os( LogOrigin("SIImageStore","getMaskSum",WHERE) );
3195 5717 : LatticeLocker lock2 (*(mask()), FileLocker::Read);
3196 11426 : LatticeExprNode msum( sum( *mask() ) );
3197 5713 : Float masksum = msum.getFloat();
3198 :
3199 : // Float masksum = sum( mask()->get() );
3200 :
3201 5713 : return masksum;
3202 5717 : }
3203 :
3204 0 : Bool SIImageStore::findMinMaxLattice(const Lattice<Float>& lattice,
3205 : const Lattice<Float>& mask,
3206 : const Lattice<Bool>& pixelmask,
3207 : Float& maxAbs, Float& maxAbsMask,
3208 : Float& minAbs, Float& minAbsMask )
3209 : {
3210 :
3211 : //FOR DEGUG
3212 : //LogIO os( LogOrigin("SIImageStore","findMinMaxLattice",WHERE) );
3213 :
3214 0 : maxAbs=0.0;maxAbsMask=0.0;
3215 0 : minAbs=1e+10;minAbsMask=1e+10;
3216 0 : LatticeLocker lock1 (const_cast<Lattice<Float>& > (lattice), FileLocker::Read);
3217 0 : LatticeLocker lock2 (const_cast<Lattice<Float>& >(mask), FileLocker::Read);
3218 0 : LatticeLocker lock3 (const_cast<Lattice<Bool>& >(pixelmask), FileLocker::Read);
3219 0 : const IPosition tileShape = lattice.niceCursorShape();
3220 0 : TiledLineStepper ls(lattice.shape(), tileShape, 0);
3221 : {
3222 0 : RO_LatticeIterator<Float> li(lattice, ls);
3223 0 : RO_LatticeIterator<Float> mi(mask, ls);
3224 0 : RO_LatticeIterator<Bool> pmi(pixelmask, ls);
3225 0 : for(li.reset(),mi.reset(),pmi.reset();!li.atEnd();li++, mi++, pmi++) {
3226 0 : IPosition posMax=li.position();
3227 0 : IPosition posMin=li.position();
3228 0 : IPosition posMaxMask=li.position();
3229 0 : IPosition posMinMask=li.position();
3230 0 : Float maxVal=0.0;
3231 0 : Float minVal=0.0;
3232 0 : Float maxValMask=0.0;
3233 0 : Float minValMask=0.0;
3234 :
3235 :
3236 : // skip if lattice chunk is masked entirely.
3237 0 : if(ntrue(pmi.cursor()) > 0 ) {
3238 0 : const MaskedArray<Float> marr(li.cursor(), pmi.cursor());
3239 0 : const MaskedArray<Float> marrinmask(li.cursor() * mi.cursor(), pmi.cursor());
3240 : //minMax( minVal, maxVal, posMin, posMax, li.cursor() );
3241 0 : minMax( minVal, maxVal, posMin, posMax, marr );
3242 : //minMaxMasked(minValMask, maxValMask, posMin, posMax, li.cursor(), mi.cursor());
3243 0 : minMax(minValMask, maxValMask, posMin, posMax, marrinmask);
3244 :
3245 : //os<<"DONE minMax"<<LogIO::POST;
3246 0 : if( (maxVal) > (maxAbs) ) maxAbs = maxVal;
3247 0 : if( (maxValMask) > (maxAbsMask) ) maxAbsMask = maxValMask;
3248 :
3249 0 : if( (minVal) < (minAbs) ) minAbs = minVal;
3250 0 : if( (minValMask) < (minAbsMask) ) minAbsMask = minValMask;
3251 0 : }
3252 :
3253 0 : }
3254 0 : }
3255 :
3256 0 : return True;
3257 :
3258 :
3259 0 : }
3260 :
3261 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3262 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3263 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3264 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3265 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3266 :
3267 : //
3268 : //-------------------------------------------------------------
3269 : // Initialize the internals of the object. Perhaps other such
3270 : // initializations in the constructors can be moved here too.
3271 : //
3272 12248 : void SIImageStore::init()
3273 : {
3274 12248 : imageExts.resize(MAX_IMAGE_IDS);
3275 :
3276 12248 : imageExts(MASK)=".mask";
3277 12248 : imageExts(PSF)=".psf";
3278 12248 : imageExts(MODEL)=".model";
3279 12248 : imageExts(RESIDUAL)=".residual";
3280 12248 : imageExts(WEIGHT)=".weight";
3281 12248 : imageExts(IMAGE)=".image";
3282 12248 : imageExts(SUMWT)=".sumwt";
3283 12248 : imageExts(GRIDWT)=".gridwt";
3284 12248 : imageExts(PB)=".pb";
3285 12248 : imageExts(FORWARDGRID)=".forward";
3286 12248 : imageExts(BACKWARDGRID)=".backward";
3287 12248 : imageExts(IMAGEPBCOR)=".image.pbcor";
3288 :
3289 12248 : itsParentPsf = itsPsf;
3290 12248 : itsParentModel=itsModel;
3291 12248 : itsParentResidual=itsResidual;
3292 12248 : itsParentWeight=itsWeight;
3293 12248 : itsParentImage=itsImage;
3294 12248 : itsParentSumWt=itsSumWt;
3295 12248 : itsParentMask=itsMask;
3296 12248 : itsParentImagePBcor=itsImagePBcor;
3297 :
3298 : // cout << "parent shape : " << itsParentImageShape << " shape : " << itsImageShape << endl;
3299 12248 : itsParentImageShape = itsImageShape;
3300 12248 : itsParentCoordSys = itsCoordSys;
3301 :
3302 12248 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 ) { itsImageShape=IPosition(4,0,0,0,0); }
3303 :
3304 12248 : itsOpened=0;
3305 :
3306 12248 : itsPSFSideLobeLevel=0.0;
3307 12248 : itsReadLock=nullptr;
3308 12248 : itsDataPolRep=StokesImageUtil::UNKNOWN; //Should throw an exception if it is not initialized properly
3309 12248 : }
3310 :
3311 :
3312 17 : void SIImageStore::regridToModelImage( ImageInterface<Float> &inputimage, Int term )
3313 : {
3314 : try
3315 : {
3316 :
3317 : //output coords
3318 17 : IPosition outshape = itsImageShape;
3319 17 : CoordinateSystem outcsys = itsCoordSys;
3320 17 : DirectionCoordinate outDirCsys = outcsys.directionCoordinate();
3321 17 : SpectralCoordinate outSpecCsys = outcsys.spectralCoordinate();
3322 :
3323 : // do regrid
3324 17 : IPosition axes(3,0, 1, 2);
3325 17 : IPosition inshape = inputimage.shape();
3326 17 : CoordinateSystem incsys = inputimage.coordinates();
3327 17 : DirectionCoordinate inDirCsys = incsys.directionCoordinate();
3328 17 : SpectralCoordinate inSpecCsys = incsys.spectralCoordinate();
3329 :
3330 17 : Vector<Int> dirAxes = CoordinateUtil::findDirectionAxes(incsys);
3331 17 : axes(0) = dirAxes(0);
3332 17 : axes(1) = dirAxes(1);
3333 17 : axes(2) = CoordinateUtil::findSpectralAxis(incsys);
3334 :
3335 : try {
3336 17 : ImageRegrid<Float> imregrid;
3337 19 : imregrid.regrid( *(model(term)), Interpolate2D::LINEAR, axes, inputimage );
3338 19 : } catch (AipsError &x) {
3339 2 : throw(AipsError("ImageRegrid error : "+ x.getMesg()));
3340 2 : }
3341 :
3342 37 : }catch(AipsError &x)
3343 : {
3344 2 : throw(AipsError("Error in regridding input model image to target coordsys : " + x.getMesg()));
3345 2 : }
3346 15 : }
3347 :
3348 : //
3349 : //---------------------------------------------------------------
3350 : //
3351 0 : void SIImageStore::makePersistent(String& fileName)
3352 : {
3353 0 : LogIO logIO(LogOrigin("SIImageStore", "makePersistent"));
3354 0 : ofstream outFile; outFile.open(fileName.c_str(),std::ofstream::out);
3355 0 : if (!outFile) logIO << "Failed to open filed \"" << fileName << "\"" << LogIO::EXCEPTION;
3356 : // String itsImageName;
3357 0 : outFile << "itsImageNameBase: " << itsImageName << endl;
3358 :
3359 : //IPosition itsImageShape;
3360 0 : outFile << "itsImageShape: " << itsImageShape.nelements() << " ";
3361 0 : for (uInt i=0;i<itsImageShape.nelements(); i++) outFile << itsImageShape(i) << " "; outFile << endl;
3362 :
3363 : // Don't know what to do with this. Looks like this gets
3364 : // filled-in from one of the images. So load this from one of the
3365 : // images if they exist?
3366 : //CoordinateSystem itsCoordSys;
3367 :
3368 : // Int itsNFacets;
3369 0 : outFile << "itsNFacets: " << itsNFacets << endl;
3370 0 : outFile << "itsUseWeight: " << itsUseWeight << endl;
3371 :
3372 :
3373 : // Misc Information to go into the header.
3374 : // Record itsMiscInfo;
3375 0 : itsMiscInfo.print(outFile);
3376 :
3377 : // std::shared_ptr<ImageInterface<Float> > itsMask, itsPsf, itsModel, itsResidual, itsWeight, itsImage, itsSumWt;
3378 : // std::shared_ptr<ImageInterface<Complex> > itsForwardGrid, itsBackwardGrid;
3379 :
3380 0 : Vector<Bool> ImageExists(MAX_IMAGE_IDS);
3381 0 : if ( ! itsMask ) ImageExists(MASK)=False;
3382 0 : if ( ! itsPsf ) ImageExists(PSF)=False;
3383 0 : if ( ! itsModel ) ImageExists(MODEL)=False;
3384 0 : if ( ! itsResidual ) ImageExists(RESIDUAL)=False;
3385 0 : if ( ! itsWeight ) ImageExists(WEIGHT)=False;
3386 0 : if ( ! itsImage ) ImageExists(IMAGE)=False;
3387 0 : if ( ! itsSumWt ) ImageExists(SUMWT)=False;
3388 0 : if ( ! itsGridWt ) ImageExists(GRIDWT)=False;
3389 0 : if ( ! itsPB ) ImageExists(PB)=False;
3390 :
3391 0 : if ( ! itsForwardGrid ) ImageExists(FORWARDGRID)=False;
3392 0 : if ( ! itsBackwardGrid ) ImageExists(BACKWARDGRID)=False;
3393 :
3394 0 : outFile << "ImagesExist: " << ImageExists << endl;
3395 0 : }
3396 : //
3397 : //---------------------------------------------------------------
3398 : //
3399 0 : void SIImageStore::recreate(String& fileName)
3400 : {
3401 0 : LogIO logIO(LogOrigin("SIImageStore", "recreate"));
3402 0 : ifstream inFile; inFile.open(fileName.c_str(),std::ofstream::out);
3403 0 : if (!inFile) logIO << "Failed to open filed \"" << fileName << "\"" << LogIO::EXCEPTION;
3404 :
3405 0 : String token;
3406 0 : inFile >> token; if (token == "itsImageNameBase:") inFile >> itsImageName;
3407 :
3408 0 : inFile >> token;
3409 0 : if (token=="itsImageShape:")
3410 : {
3411 : Int n;
3412 0 : inFile >> n;
3413 0 : itsImageShape.resize(n);
3414 0 : for (Int i=0;i<n; i++) inFile >> itsImageShape(i);
3415 : }
3416 :
3417 : // Int itsNFacets;
3418 0 : inFile >> token; if (token=="itsNFacets:") inFile >> itsNFacets;
3419 0 : inFile >> token; if (token=="itsUseWeight:") inFile >> itsUseWeight;
3420 :
3421 0 : Bool coordSysLoaded=False;
3422 0 : String itsName;
3423 : try
3424 : {
3425 0 : itsName=itsImageName+imageExts(MASK);casa::openImage(itsName, itsMask);
3426 0 : if (coordSysLoaded==False) {itsCoordSys=itsMask->coordinates(); itsMiscInfo=itsImage->miscInfo();coordSysLoaded=True;}
3427 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3428 : try
3429 : {
3430 0 : itsName=itsImageName+imageExts(MODEL);casa::openImage(itsName, itsModel);
3431 0 : if (coordSysLoaded==False) {itsCoordSys=itsModel->coordinates(); itsMiscInfo=itsModel->miscInfo();coordSysLoaded=True;}
3432 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3433 : try
3434 : {
3435 0 : itsName=itsImageName+imageExts(RESIDUAL);casa::openImage(itsName, itsResidual);
3436 0 : if (coordSysLoaded==False) {itsCoordSys=itsResidual->coordinates(); itsMiscInfo=itsResidual->miscInfo();coordSysLoaded=True;}
3437 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3438 : try
3439 : {
3440 0 : itsName=itsImageName+imageExts(WEIGHT);casa::openImage(itsName, itsWeight);
3441 0 : if (coordSysLoaded==False) {itsCoordSys=itsWeight->coordinates(); itsMiscInfo=itsWeight->miscInfo();coordSysLoaded=True;}
3442 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3443 : try
3444 : {
3445 0 : itsName=itsImageName+imageExts(IMAGE);casa::openImage(itsName, itsImage);
3446 0 : if (coordSysLoaded==False) {itsCoordSys=itsImage->coordinates(); itsMiscInfo=itsImage->miscInfo();coordSysLoaded=True;}
3447 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3448 : try
3449 : {
3450 0 : itsName=itsImageName+imageExts(SUMWT);casa::openImage(itsName, itsSumWt);
3451 0 : if (coordSysLoaded==False) {itsCoordSys=itsSumWt->coordinates(); itsMiscInfo=itsSumWt->miscInfo();coordSysLoaded=True;}
3452 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3453 : try
3454 : {
3455 0 : itsName=itsImageName+imageExts(PSF);casa::openImage(itsName, itsPsf);
3456 0 : if (coordSysLoaded==False) {itsCoordSys=itsPsf->coordinates(); itsMiscInfo=itsPsf->miscInfo();coordSysLoaded=True;}
3457 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3458 : try
3459 : {
3460 0 : casa::openImage(itsImageName+imageExts(FORWARDGRID), itsForwardGrid);
3461 0 : casa::openImage(itsImageName+imageExts(BACKWARDGRID), itsBackwardGrid);
3462 : }
3463 0 : catch (AipsError& x)
3464 : {
3465 0 : logIO << "Did not find forward and/or backward grid. Just say'n..." << LogIO::POST;
3466 0 : }
3467 :
3468 0 : }
3469 : //////////////
3470 15 : Bool SIImageStore::intersectComplexImage(const String& ComplexImageName){
3471 15 : Vector<Int> whichStokes(0);
3472 15 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
3473 15 : whichStokes, itsDataPolRep);
3474 :
3475 :
3476 : //cerr <<"itsDataPolRep " << itsDataPolRep << endl;
3477 :
3478 15 : CountedPtr<PagedImage<Complex> > compliantImage =nullptr;
3479 : {
3480 15 : PagedImage<Complex>inputImage(ComplexImageName);
3481 15 : IPosition theShape=itsImageShape;
3482 15 : theShape(0)=inputImage.shape()(0);
3483 15 : theShape(1)=inputImage.shape()(1);
3484 15 : CoordinateSystem inpcsys=inputImage.coordinates();
3485 15 : Vector<Double> refpix=cimageCoord.referencePixel();
3486 15 : refpix(0)+=(theShape(0)-itsImageShape(0))/2.0;
3487 15 : refpix(1)+=(theShape(1)-itsImageShape(1))/2.0;
3488 15 : cimageCoord.setReferencePixel(refpix);
3489 30 : String tmpImage=File::newUniqueName(".", "TempImage").absoluteName();
3490 15 : compliantImage=new PagedImage<Complex>(theShape, cimageCoord, tmpImage);
3491 15 : compliantImage->set(0.0);
3492 15 : IPosition iblc(theShape.nelements(),0);
3493 15 : IPosition itrc=theShape-1;
3494 : //cerr << "blc " << iblc << " trc " << itrc << " shape " << theShape << endl;
3495 15 : LCBox lbox(iblc, itrc, theShape);
3496 15 : ImageRegion imagreg(WCBox(lbox, cimageCoord));
3497 :
3498 :
3499 15 : SubImage<Complex> subim(inputImage, imagreg, false);
3500 : //cerr << "shapes " << inputImage.shape() << " sub " << subim.shape() << " compl " << compliantImage->shape() << endl;
3501 15 : compliantImage->copyData(subim);
3502 :
3503 15 : }
3504 15 : TableUtil::deleteTable(ComplexImageName);
3505 15 : compliantImage->rename(ComplexImageName);
3506 15 : return True;
3507 :
3508 15 : }
3509 :
3510 :
3511 : //////////////////////////////////////////////////////////////////////////////////////////////////////
3512 : //////////////////////////////////////////////////////////////////////////////////////////////////////
3513 :
3514 : } //# NAMESPACE CASA - END
3515 :
|
