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 774 : 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 774 : const Bool useweightimage)
155 : // TODO : Add parameter to indicate weight image shape.
156 : {
157 1548 : LogIO os( LogOrigin("SIImageStore","Open new Images",WHERE) );
158 :
159 :
160 774 : itsPsf.reset( );
161 774 : itsModel.reset( );
162 774 : itsResidual.reset( );
163 774 : itsWeight.reset( );
164 774 : itsImage.reset( );
165 774 : itsMask.reset( );
166 774 : itsGridWt.reset( );
167 774 : itsPB.reset( );
168 774 : itsImagePBcor.reset( );
169 774 : itsTempWorkIm.reset();
170 :
171 774 : itsSumWt.reset( );
172 774 : itsOverWrite=False; // Hard Coding this. See CAS-6937. overwrite;
173 774 : itsUseWeight=useweightimage;
174 774 : itsPBScaleFactor=1.0;
175 :
176 774 : itsNFacets=1;
177 774 : itsFacetId=0;
178 774 : itsNChanChunks = 1;
179 774 : itsChanId = 0;
180 774 : itsNPolChunks = 1;
181 774 : itsPolId = 0;
182 :
183 774 : itsImageName = imagename;
184 774 : itsCoordSys = imcoordsys;
185 774 : itsImageShape = imshape;
186 774 : itsObjectName = objectname;
187 774 : itsMiscInfo = miscinfo;
188 :
189 774 : init();
190 :
191 774 : validate();
192 774 : }
193 :
194 : // Used from SynthesisNormalizer::makeImageStore()
195 3012 : SIImageStore::SIImageStore(const String &imagename, const Bool ignorefacets, const Bool noRequireSumwt)
196 : {
197 6024 : LogIO os( LogOrigin("SIImageStore","Open existing Images",WHERE) );
198 :
199 :
200 3012 : itsPsf.reset( );
201 3012 : itsModel.reset( );
202 3012 : itsResidual.reset( );
203 3012 : itsWeight.reset( );
204 3012 : itsImage.reset( );
205 3012 : itsMask.reset( );
206 3012 : itsGridWt.reset( );
207 3012 : itsPB.reset( );
208 3012 : itsImagePBcor.reset( );
209 3012 : itsTempWorkIm.reset();
210 3012 : itsMiscInfo=Record();
211 :
212 3012 : itsSumWt.reset( );
213 3012 : itsNFacets=1;
214 3012 : itsFacetId=0;
215 3012 : itsNChanChunks = 1;
216 3012 : itsChanId = 0;
217 3012 : itsNPolChunks = 1;
218 3012 : itsPolId = 0;
219 :
220 3012 : itsOverWrite=False;
221 : //need to to this init now so that imageExts is initialized
222 3012 : init();
223 :
224 3012 : 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 6024 : if( doesImageExist(itsImageName+String(".residual")) ||
230 4347 : doesImageExist(itsImageName+String(".psf")) ||
231 3014 : doesImageExist(itsImageName+String(".model")) ||
232 3012 : doesImageExist(itsImageName+String(".gridwt")) ||
233 10371 : doesImageExist(itsImageName+String(".pb")) ||
234 3012 : doesImageExist(itsImageName+String(".weight"))
235 : )
236 : {
237 3012 : std::shared_ptr<ImageInterface<Float> > imptr;
238 3012 : if( doesImageExist(itsImageName+String(".psf")) )
239 : {
240 2866 : buildImage( imptr, (itsImageName+String(".psf")) );
241 : // itsObjectName=imptr->imageInfo().objectName();
242 : // itsMiscInfo=imptr->miscInfo();
243 : }
244 146 : else if ( doesImageExist(itsImageName+String(".residual")) ){
245 144 : 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 3012 : itsObjectName=imptr->imageInfo().objectName();
271 3012 : itsImageShape=imptr->shape();
272 3012 : itsCoordSys = imptr->coordinates();
273 3012 : itsMiscInfo=imptr->miscInfo();
274 :
275 3012 : }
276 : else
277 : {
278 0 : throw( AipsError( "PSF, Residual, Model Image (or sumwt) do not exist. Please create one of them." ) );
279 : }
280 :
281 7359 : if( doesImageExist(itsImageName+String(".residual")) ||
282 4347 : doesImageExist(itsImageName+String(".psf")) )
283 : {
284 3010 : if( doesImageExist(itsImageName+String(".sumwt")) )
285 : {
286 2925 : std::shared_ptr<ImageInterface<Float> > imptr;
287 : //imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt")) );
288 2925 : buildImage( imptr, (itsImageName+String(".sumwt")) );
289 2925 : itsNFacets = imptr->shape()[0];
290 2925 : itsFacetId = 0;
291 2925 : itsUseWeight = getUseWeightImage( *imptr );
292 2925 : 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 2925 : itsCoordSys = imptr->coordinates();
295 2925 : itsMiscInfo=imptr->miscInfo();
296 2925 : 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 2925 : }
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 3012 : if( ignorefacets==True ) itsNFacets= 1;
326 :
327 3012 : init();
328 :
329 3012 : validate();
330 3012 : }
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 8992 : void SIImageStore::validate()
416 : {
417 : /// There are two valid states. Check for at least one of them.
418 8992 : Bool state = False;
419 :
420 8992 : stringstream oss;
421 8992 : oss << "shape:" << itsImageShape << " parentimageshape:" << itsParentImageShape
422 8992 : << " nfacets:" << itsNFacets << "x" << itsNFacets << " facetid:" << itsFacetId
423 8992 : << " nchanchunks:" << itsNChanChunks << " chanid:" << itsChanId
424 8992 : << " npolchunks:" << itsNPolChunks << " polid:" << itsPolId
425 8992 : << " coord-dim:" << itsCoordSys.nPixelAxes()
426 8992 : << " psf/res:" << (hasPsf() || hasResidual()) ;
427 8992 : if( hasSumWt() ) oss << " sumwtshape : " << sumwt()->shape() ;
428 8992 : oss << endl;
429 :
430 :
431 : try {
432 :
433 8992 : if( itsCoordSys.nPixelAxes() != 4 ) state=False;
434 :
435 : /// (1) Regular imagestore
436 8992 : if( itsNFacets==1 && itsFacetId==0
437 8800 : && itsNChanChunks==1 && itsChanId==0
438 8800 : && itsNPolChunks==1 && itsPolId==0 ) {
439 8720 : Bool check1 = hasSumWt() && sumwt()->shape()[0]==1;
440 17440 : if( (itsImageShape.isEqual(itsParentImageShape) ) && ( check1 || !hasSumWt() )
441 17440 : && itsParentImageShape.product() > 0 ) state=True;
442 8720 : }
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 8992 : if( state == False ) throw(AipsError("Internal Error : Invalid ImageStore state : "+ oss.str()) );
467 :
468 17984 : return;
469 8992 : }
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 31833 : 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 31833 : std::shared_ptr<ImageInterface<Float> > imPtr;
492 31833 : IPosition useShape( itsParentImageShape );
493 :
494 31833 : if( dosumwt ) // change shape to sumwt image shape.
495 : {
496 4351 : useShape[0] = nfacetsperside;
497 4351 : 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 31833 : Bool dbg=False;
517 31833 : if( doesImageExist( imagenamefull ) )
518 : {
519 : ///// not used since itsOverWrite is hardcoded to FALSE (CAS-6937)
520 25691 : 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 25691 : if(dbg) cout << "Trying to open existing image : "<< imagenamefull << endl;
553 : try{
554 25691 : buildImage( imPtr, imagenamefull ) ;
555 :
556 25691 : 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 22311 : 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 22296 : 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 6142 : if(dbg) cout << "Trying to open new image named : " << imagenamefull << endl;
613 : try{
614 6142 : buildImage(imPtr, useShape, itsParentCoordSys, imagenamefull) ;
615 : // initialize to zeros...
616 : {
617 6142 : LatticeLocker lock1(*imPtr, FileLocker::Write);
618 6142 : imPtr->set(0.0);
619 6142 : imPtr->flush();
620 6142 : imPtr->unlock();
621 6142 : }
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 63632 : return imPtr;
662 31850 : }
663 :
664 6142 : void SIImageStore::buildImage(std::shared_ptr<ImageInterface<Float> > &imptr, IPosition shape,
665 : CoordinateSystem csys, const String name)
666 : {
667 12284 : LogIO os( LogOrigin("SIImageStore", "Open non-existing image", WHERE) );
668 6142 : os <<"Opening image, name: " << name << LogIO::DEBUG1;
669 :
670 6142 : itsOpened++;
671 : //For some reason cannot open a new paged image with UserNoread directly
672 : {
673 6142 : PagedImage<Float> godot(shape, csys, name);
674 6142 : godot.unlock();
675 6142 : }
676 6142 : 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 6142 : imptr.reset( new PagedImage<Float> (name, locktype) );
682 : {
683 6142 : LatticeLocker lock1(*imptr, FileLocker::Write);
684 6142 : initMetaInfo(imptr, name);
685 6142 : imptr->unlock();
686 6142 : }
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 6142 : }
701 :
702 33170 : void SIImageStore::buildImage(std::shared_ptr<ImageInterface<Float> > &imptr, const String name)
703 : {
704 66340 : LogIO os(LogOrigin("SIImageStore", "Open existing Images", WHERE));
705 33170 : os <<"Opening image, name: " << name << LogIO::DEBUG1;
706 :
707 33170 : itsOpened++;
708 33170 : if(Table::isReadable(name)){
709 33170 : 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 33170 : Table table(name, locktype);
714 33170 : String type = table.tableInfo().type();
715 33170 : if (type != TableInfo::type(TableInfo::PAGEDIMAGE)) {
716 :
717 0 : imptr.reset( new PagedImage<Float>( table ) );
718 0 : imptr->unlock();
719 0 : return;
720 : }
721 33170 : }
722 33170 : LatticeBase* latt =ImageOpener::openImage(name);
723 33170 : if(!latt)
724 : {
725 0 : throw(AipsError("Error in opening Image : "+name));
726 : }
727 33170 : DataType dtype=latt->dataType();
728 33170 : if(dtype==TpFloat)
729 : {
730 33170 : 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 33170 : }
762 :
763 : /**
764 : * Sets ImageInfo and MiscInfo on an image
765 : *
766 : * @param imptr image to initialize
767 : */
768 6144 : 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 6144 : LatticeLocker lock1(*imptr, FileLocker::Write);
774 6144 : if (not itsObjectName.empty()) {
775 6144 : ImageInfo info = imptr->imageInfo();
776 6144 : info.setObjectName(itsObjectName);
777 6144 : imptr->setImageInfo(info);
778 6144 : imptr->setMiscInfo(itsMiscInfo);
779 6144 : } 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 6144 : imptr->unlock();
787 6144 : }
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 34914 : 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 34914 : Int nx_facets=Int(sqrt(Double(nfacets)));
810 69828 : LogIO os( LogOrigin("SynthesisImager","makeFacet") );
811 : // Make the output image
812 34914 : Slicer imageSlicer;
813 :
814 : // Add checks for all dimensions..
815 34914 : 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 34914 : IPosition imshp=image.shape();
820 34914 : IPosition blc(imshp.nelements(), 0);
821 34914 : IPosition trc=imshp-1;
822 34914 : IPosition inc(imshp.nelements(), 1);
823 :
824 : /// Facets
825 34914 : Int facetx = facet % nx_facets;
826 34914 : Int facety = (facet - facetx) / nx_facets;
827 34914 : Int sizex = imshp(0) / nx_facets;
828 34914 : Int sizey = imshp(1) / nx_facets;
829 34914 : blc(1) = facety * sizey;
830 34914 : trc(1) = blc(1) + sizey - 1;
831 34914 : blc(0) = facetx * sizex;
832 34914 : trc(0) = blc(0) + sizex - 1;
833 :
834 : /// Pol Chunks
835 34914 : Int sizepol = imshp(2) / npolchunks;
836 34914 : blc(2) = pol * sizepol;
837 34914 : trc(2) = blc(2) + sizepol - 1;
838 :
839 : /// Chan Chunks
840 34914 : Int sizechan = imshp(3) / nchanchunks;
841 34914 : blc(3) = chan * sizechan;
842 34914 : trc(3) = blc(3) + sizechan - 1;
843 :
844 34914 : LCBox::verify(blc, trc, inc, imshp);
845 34914 : Slicer imslice(blc, trc, inc, Slicer::endIsLast);
846 :
847 : // Now create the sub image
848 34914 : std::shared_ptr<ImageInterface<Float> > referenceImage( new SubImage<Float>(image, imslice, True) );
849 : {
850 34914 : LatticeLocker lock1(*referenceImage, FileLocker::Write);
851 34914 : referenceImage->setMiscInfo(image.miscInfo());
852 34914 : referenceImage->setUnits(image.units());
853 :
854 34914 : }
855 :
856 : // cout << "Made Ref subimage of shape : " << referenceImage->shape() << endl;
857 :
858 34914 : return referenceImage;
859 :
860 34914 : }
861 :
862 :
863 :
864 : //////////////////////////////////////////////////////////////////////////////////////////////////////
865 :
866 12278 : SIImageStore::~SIImageStore()
867 : {
868 12278 : }
869 :
870 : //////////////////////////////////////////////////////////////////////////////////////////////////////
871 : //////////////////////////////////////////////////////////////////////////////////////////////////////
872 :
873 12335 : Bool SIImageStore::releaseLocks()
874 : {
875 : //LogIO os( LogOrigin("SIImageStore","releaseLocks",WHERE) );
876 :
877 : // String fname( itsImageName+String(".info") );
878 : // makePersistent( fname );
879 :
880 12335 : if( itsPsf ) releaseImage( itsPsf );
881 12335 : if( itsModel ) { releaseImage( itsModel ); }
882 12335 : if( itsResidual ) releaseImage( itsResidual );
883 12335 : if( itsImage ) releaseImage( itsImage );
884 12335 : if( itsWeight ) releaseImage( itsWeight );
885 12335 : if( itsMask ) releaseImage( itsMask );
886 12335 : if( itsSumWt ) releaseImage( itsSumWt );
887 12335 : if( itsGridWt ) releaseImage( itsGridWt );
888 12335 : if( itsPB ) releaseImage( itsPB );
889 12335 : if( itsImagePBcor ) releaseImage( itsImagePBcor );
890 :
891 12335 : 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 25884 : 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 25884 : im->flush();
909 : //os << LogIO::WARN << "clear cache" << LogIO::POST;
910 25884 : im->clearCache();
911 : //os << LogIO::WARN << "unlock" << LogIO::POST;
912 25884 : im->unlock();
913 : //os << LogIO::WARN << "tempClose" << LogIO::POST;
914 : //im->tempClose();
915 : //os << LogIO::WARN << "NULL" << LogIO::POST;
916 25884 : im.reset(); // This was added to allow modification by modules independently
917 25884 : }
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 747 : Long SIImageStore::estimateRAM(){
934 : //right now this is estimated at 2MB for the 2 complex lattices;
935 747 : 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 14498 : IPosition SIImageStore::getShape()
1056 : {
1057 14498 : return itsImageShape;
1058 : }
1059 :
1060 6585 : String SIImageStore::getName()
1061 : {
1062 6585 : return itsImageName;
1063 : }
1064 :
1065 8385 : uInt SIImageStore::getNTaylorTerms(Bool /*dopsf*/)
1066 : {
1067 8385 : 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 130153 : 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 130153 : Bool sw=False;
1087 130153 : if( label.contains(imageExts(SUMWT)) ) sw=True;
1088 :
1089 130153 : if( !ptr )
1090 : {
1091 : //cout << itsNFacets << " " << itsNChanChunks << " " << itsNPolChunks << endl;
1092 32435 : 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 31619 : if(!parentptr){
1120 : ///coordsys for psf can be different ...shape should be the same.
1121 30957 : 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 130136 : }
1136 :
1137 :
1138 14223 : std::shared_ptr<ImageInterface<Float> > SIImageStore::psf(uInt /*nterm*/)
1139 : {
1140 14223 : accessImage( itsPsf, itsParentPsf, imageExts(PSF) );
1141 :
1142 14222 : return itsPsf;
1143 : }
1144 30888 : std::shared_ptr<ImageInterface<Float> > SIImageStore::residual(uInt /*nterm*/)
1145 : {
1146 30888 : accessImage( itsResidual, itsParentResidual, imageExts(RESIDUAL) );
1147 : // Record mi = itsResidual->miscInfo(); ostringstream oss;mi.print(oss);cout<<"MiscInfo(res) : " << oss.str() << endl;
1148 30888 : return itsResidual;
1149 : }
1150 1868 : std::shared_ptr<ImageInterface<Float> > SIImageStore::weight(uInt /*nterm*/)
1151 : {
1152 1868 : accessImage( itsWeight, itsParentWeight, imageExts(WEIGHT) );
1153 1868 : return itsWeight;
1154 : }
1155 :
1156 10097 : std::shared_ptr<ImageInterface<Float> > SIImageStore::sumwt(uInt /*nterm*/)
1157 : {
1158 :
1159 10097 : accessImage( itsSumWt, itsParentSumWt, imageExts(SUMWT) );
1160 :
1161 10097 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
1162 536 : { itsUseWeight = getUseWeightImage( *itsParentSumWt );}
1163 10097 : setUseWeightImage( *itsSumWt , itsUseWeight); // Sets a flag in the SumWt image.
1164 :
1165 10097 : return itsSumWt;
1166 : }
1167 :
1168 11071 : std::shared_ptr<ImageInterface<Float> > SIImageStore::model(uInt /*nterm*/)
1169 : {
1170 11073 : accessImage( itsModel, itsParentModel, imageExts(MODEL) );
1171 11069 : LatticeLocker lock1(*itsModel, FileLocker::Write);
1172 : // Set up header info the first time.
1173 11069 : itsModel->setUnits("Jy/pixel");
1174 :
1175 22138 : return itsModel;
1176 11069 : }
1177 6003 : std::shared_ptr<ImageInterface<Float> > SIImageStore::image(uInt /*nterm*/)
1178 : {
1179 6003 : accessImage( itsImage, itsParentImage, imageExts(IMAGE) );
1180 6003 : LatticeLocker lock1(*itsImage, FileLocker::Write);
1181 6003 : itsImage->setUnits(Unit("Jy/beam"));
1182 12006 : return itsImage;
1183 6003 : }
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 12861 : std::shared_ptr<ImageInterface<Float> > SIImageStore::pb(uInt /*nterm*/)
1224 : {
1225 12861 : accessImage( itsPB, itsParentPB, imageExts(PB) );
1226 12859 : 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 2532 : std::shared_ptr<ImageInterface<Complex> > SIImageStore::backwardGrid(uInt /*nterm*/){
1265 2532 : if( itsBackwardGrid ) //&& (itsBackwardGrid->shape() == itsImageShape))
1266 : {
1267 : // cout << "Backward grid has shape : " << itsBackwardGrid->shape() << endl;
1268 2037 : return itsBackwardGrid;
1269 : }
1270 495 : Vector<Int> whichStokes(0);
1271 495 : IPosition cimageShape;
1272 495 : cimageShape=itsImageShape;
1273 495 : 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 495 : if(freqframe != MFrequency::LSRK && freqframe!=MFrequency::Undefined && freqframe!=MFrequency::REST )
1276 0 : { itsCoordSys.setSpectralConversion("LSRK"); }
1277 495 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
1278 495 : whichStokes, itsDataPolRep);
1279 495 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
1280 495 : cimageShape(2)=whichStokes.nelements();
1281 : //cout << "Making backward grid of shape : " << cimageShape << " for imshape : " << itsImageShape << endl;
1282 495 : itsBackwardGrid.reset( new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()) );
1283 : //if(image())
1284 495 : if(hasRestored()){
1285 8 : LatticeLocker lock1(*itsBackwardGrid, FileLocker::Write);
1286 8 : itsBackwardGrid->setImageInfo((image())->imageInfo());
1287 8 : }
1288 495 : return itsBackwardGrid;
1289 495 : }
1290 2292 : Double SIImageStore::memoryBeforeLattice(){
1291 : //Calculate how much memory to use per temporary images before disking
1292 2292 : return 1.0; /// in MB
1293 : }
1294 2292 : IPosition SIImageStore::tileShape(){
1295 : //Need to have settable stuff here or algorith to determine this
1296 2292 : 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 47610 : Bool SIImageStore::doesImageExist(String imagename)
1301 : {
1302 95220 : LogIO os( LogOrigin("SIImageStore","doesImageExist",WHERE) );
1303 47610 : Directory image( imagename );
1304 95220 : return image.exists();
1305 47610 : }
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 1075 : 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 1075 : CountedPtr<ImageInterface<Float> > subim=makeSubImage(0, 1,
1411 1075 : chan, itsImageShape[3],
1412 1075 : pol, itsImageShape[2],
1413 2150 : *weight(0));
1414 :
1415 2150 : return sqrt(max(subim->get()));
1416 1075 : }
1417 :
1418 113 : void SIImageStore::makePBFromWeight(const Float pblimit)
1419 : {
1420 226 : LogIO os( LogOrigin("SIImageStore","makePBFromWeight",WHERE) );
1421 :
1422 226 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1423 : {
1424 350 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1425 : {
1426 :
1427 237 : itsPBScaleFactor = getPbMax(pol,chan);
1428 :
1429 237 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1430 : else {
1431 :
1432 234 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1433 234 : chan, itsImageShape[3],
1434 234 : pol, itsImageShape[2],
1435 468 : *weight() );
1436 234 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
1437 234 : chan, itsImageShape[3],
1438 234 : pol, itsImageShape[2],
1439 468 : *pb() );
1440 :
1441 :
1442 468 : LatticeExpr<Float> normed( sqrt(abs(*wtsubim)) / itsPBScaleFactor );
1443 468 : LatticeExpr<Float> limited( iif( normed > fabs(pblimit) , normed, 0.0 ) );
1444 234 : pbsubim->copyData( limited );
1445 234 : }// if not zero
1446 : }
1447 : }
1448 :
1449 113 : 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 172 : LatticeExpr<Bool> pbmask( iif( *pb() > fabs(pblimit) , True , False ) );
1458 : //MSK//
1459 86 : createMask( pbmask, pb() );
1460 86 : pb()->pixelMask().unlock();
1461 86 : }
1462 :
1463 : }
1464 113 : weight()->unlock();
1465 113 : pb()->unlock();
1466 113 : }
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 901 : 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 901 : LatticeLocker lock1(*outimage, FileLocker::Write);
1526 901 : ImageRegion outreg = outimage->makeMask("mask0",False,True);
1527 901 : LCRegion& outmask=outreg.asMask();
1528 901 : outmask.copyData(lemask);
1529 901 : outimage->defineRegion("mask0",outreg, RegionHandler::Masks, True);
1530 :
1531 901 : outimage->setDefaultMask("mask0");
1532 :
1533 901 : outimage->unlock();
1534 901 : outimage->tempClose();
1535 :
1536 : // outimage->table().unmarkForDelete();
1537 901 : }
1538 0 : catch (const AipsError& x) {
1539 0 : throw(AipsError("Error in creating internal T/F mask : " + x.getMesg() ));
1540 0 : }
1541 :
1542 901 : return True;
1543 : }
1544 :
1545 1748 : 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 1748 : if( (inimage->getDefaultMask()).matches("mask0") ) // input mask exists.
1553 1558 : {LatticeLocker lock1(*outimage, FileLocker::Write);
1554 1558 : 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 1561 : ImageRegion outreg=outimage->makeMask("mask0", False, True);
1563 1555 : LCRegion& outmask=outreg.asMask();
1564 1555 : outmask.copyData(inimage->getRegion("mask0").asLCRegion());
1565 1555 : outimage->defineRegion("mask0",outreg, RegionHandler::Masks,True);
1566 1555 : outimage->setDefaultMask("mask0");
1567 1558 : }
1568 : }
1569 3 : catch (const AipsError& x) {
1570 3 : throw(AipsError("Error in copying internal T/F mask : " + x.getMesg() ));
1571 3 : }
1572 :
1573 1745 : return True;
1574 : }
1575 :
1576 1893 : 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 1893 : LatticeLocker lock1(*im, FileLocker::Write);
1586 1893 : if (im-> getDefaultMask() != String(""))
1587 : {
1588 388 : String strung=im->getDefaultMask();
1589 388 : im->setDefaultMask("");
1590 388 : im->removeRegion(strung);
1591 388 : }
1592 1893 : if( im->hasRegion("mask0") )
1593 : {
1594 0 : im->removeRegion("mask0");
1595 : }
1596 1893 : }
1597 0 : catch (const AipsError& x)
1598 : {
1599 0 : throw(AipsError("Error in deleting internal T/F mask : " + x.getMesg() ));
1600 0 : }
1601 1893 : }
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 591 : void SIImageStore::dividePSFByWeight(const Float /* pblimit*/)
1648 : {
1649 1182 : LogIO os( LogOrigin("SIImageStore","dividePSFByWeight",WHERE) );
1650 :
1651 591 : LatticeLocker lock1 (*(psf()), FileLocker::Write);
1652 591 : normPSF();
1653 :
1654 591 : if( itsUseWeight )
1655 : {
1656 :
1657 71 : divideImageByWeightVal( *weight() );
1658 : }
1659 591 : (psf())->unlock();
1660 :
1661 591 : }
1662 :
1663 598 : void SIImageStore::normalizePrimaryBeam(const Float pblimit)
1664 : {
1665 1196 : LogIO os( LogOrigin("SIImageStore","normalizePrimaryBeam",WHERE) );
1666 :
1667 : // cout << "In dividePSFByWeight : itsUseWeight : " << itsUseWeight << endl;
1668 598 : if( itsUseWeight )
1669 : {
1670 :
1671 146 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1672 : {
1673 270 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1674 : {
1675 197 : 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 73 : makePBFromWeight(pblimit);
1680 :
1681 : }//if itsUseWeight
1682 525 : else { makePBImage(pblimit); } // OR... just check that it exists already.
1683 598 : pb()->unlock();
1684 :
1685 598 : }
1686 :
1687 : // Make another for the PSF too.
1688 1329 : void SIImageStore::divideResidualByWeight(Float pblimit, String normtype) {
1689 :
1690 2658 : LogIO os(LogOrigin("SIImageStore", "divideResidualByWeight", WHERE));
1691 1329 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1692 :
1693 363 : auto logTemplate = [&](Int const &chan, Int const &pol, string const &normalizer, string const &result) {
1694 363 : os << LogIO::NORMAL1
1695 726 : << "[C" + String::toString(chan) + ":P" + String::toString(pol) + "] "
1696 726 : << "Dividing " << itsImageName + String(".residual") << " by "
1697 : << "[ " << normalizer << " ] "
1698 1089 : << "to get " << result << "." << LogIO::POST;
1699 363 : };
1700 :
1701 : // Normalize by the sumwt, per plane.
1702 1329 : Bool didNorm = divideImageByWeightVal(*residual());
1703 :
1704 1329 : if (itsUseWeight) {
1705 280 : for (Int pol = 0; pol < itsImageShape[2]; pol++) {
1706 512 : for (Int chan = 0; chan < itsImageShape[3]; chan++) {
1707 372 : itsPBScaleFactor = getPbMax(pol, chan);
1708 :
1709 372 : 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 363 : CountedPtr<ImageInterface<Float> > wtsubim = makeSubImage(0, 1,
1716 363 : chan, itsImageShape[3],
1717 363 : pol, itsImageShape[2],
1718 726 : *weight());
1719 363 : CountedPtr<ImageInterface<Float> > ressubim = makeSubImage(0, 1,
1720 363 : chan, itsImageShape[3],
1721 363 : pol, itsImageShape[2],
1722 726 : *residual());
1723 363 : LatticeExpr<Float> ratio;
1724 363 : Float scalepb = 1.0;
1725 :
1726 363 : if (normtype == "flatnoise") {
1727 1089 : logTemplate(chan, pol,
1728 726 : "sqrt(weightimage) * " + String::toString(itsPBScaleFactor),
1729 : "flat noise with unit pb peak");
1730 :
1731 726 : LatticeExpr<Float> deno = itsPBScaleFactor * sqrt(abs(LatticeExpr<Float>(*(wtsubim))));
1732 363 : scalepb = fabs(pblimit) * itsPBScaleFactor * itsPBScaleFactor;
1733 363 : ratio = iif(deno > scalepb, (*(ressubim) / deno), 0.0);
1734 :
1735 363 : } 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 363 : 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 363 : } // 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 1329 : 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 1329 : if ((residual()->getDefaultMask() == "") && hasPB() && pblimit >=0.0) {
1779 230 : copyMask(pb(), residual());
1780 : }
1781 :
1782 1329 : if ((pblimit < 0.0) && (residual()->getDefaultMask()).matches("mask0")) {
1783 1 : removeMask(residual());
1784 : }
1785 :
1786 1329 : residual()->unlock();
1787 1329 : }
1788 :
1789 144 : void SIImageStore::divideResidualByWeightSD(Float pblimit) {
1790 :
1791 288 : LogIO os(LogOrigin("SIImageStore", "divideResidualByWeightSD", WHERE));
1792 144 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1793 :
1794 144 : if (itsUseWeight) {
1795 288 : LatticeExpr<Float> deno = LatticeExpr<Float>(*weight());
1796 288 : LatticeExpr<Float> ratio = iif(deno > 0.0, *(residual()) / deno, 0.0);
1797 144 : residual()->copyData(ratio);
1798 144 : }
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 144 : if ((residual()->getDefaultMask() == "") && hasPB() && pblimit >=0.0) {
1808 0 : copyMask(pb(), residual());
1809 : }
1810 :
1811 144 : if ((pblimit < 0.0) && (residual()->getDefaultMask()).matches("mask0")) {
1812 0 : removeMask(residual());
1813 : }
1814 :
1815 144 : residual()->unlock();
1816 144 : }
1817 :
1818 869 : void SIImageStore::divideModelByWeight(Float pblimit, const String normtype)
1819 : {
1820 1738 : LogIO os( LogOrigin("SIImageStore","divideModelByWeight",WHERE) );
1821 :
1822 :
1823 1738 : if(itsUseWeight // only when needed
1824 869 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
1825 : {
1826 :
1827 93 : 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 93 : else if( normtype=="flatnoise"){
1836 :
1837 186 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1838 : {
1839 341 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1840 : {
1841 :
1842 248 : itsPBScaleFactor = getPbMax(pol,chan);
1843 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1844 248 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1845 : else {
1846 :
1847 222 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1848 222 : chan, itsImageShape[3],
1849 222 : pol, itsImageShape[2],
1850 444 : *weight() );
1851 222 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1852 222 : chan, itsImageShape[3],
1853 222 : pol, itsImageShape[2],
1854 444 : *model() );
1855 222 : os << LogIO::NORMAL1 ;
1856 222 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1857 222 : os << "Dividing " << itsImageName+String(".model") ;
1858 222 : os << " by [ sqrt(weight) / " << itsPBScaleFactor ;
1859 222 : os <<" ] to get to flat sky model before prediction" << LogIO::POST;
1860 :
1861 :
1862 444 : LatticeExpr<Float> deno( sqrt( abs(*(wtsubim)) ) / itsPBScaleFactor );
1863 :
1864 444 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1865 444 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1866 444 : 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 222 : 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 222 : 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 222 : }// 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 869 : }
1950 :
1951 798 : void SIImageStore::multiplyModelByWeight(Float pblimit, const String normtype)
1952 : {
1953 1596 : LogIO os( LogOrigin("SIImageStore","multiplyModelByWeight",WHERE) );
1954 :
1955 1596 : if(itsUseWeight // only when needed
1956 798 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
1957 : {
1958 21 : 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 21 : else if( normtype=="flatnoise"){
1963 :
1964 42 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1965 : {
1966 42 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1967 : {
1968 :
1969 21 : itsPBScaleFactor = getPbMax(pol,chan);
1970 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1971 21 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1972 : else {
1973 :
1974 21 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1975 21 : chan, itsImageShape[3],
1976 21 : pol, itsImageShape[2],
1977 42 : *weight() );
1978 21 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1979 21 : chan, itsImageShape[3],
1980 21 : pol, itsImageShape[2],
1981 42 : *model() );
1982 :
1983 :
1984 :
1985 21 : os << LogIO::NORMAL1 ;
1986 21 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1987 21 : os << "Multiplying " << itsImageName+String(".model") ;
1988 21 : os << " by [ sqrt(weight) / " << itsPBScaleFactor;
1989 21 : os << " ] to take model back to flat noise with unit pb peak." << LogIO::POST;
1990 :
1991 42 : LatticeExpr<Float> deno( sqrt( abs(*(wtsubim)) ) / itsPBScaleFactor );
1992 :
1993 42 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1994 42 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1995 42 : 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 21 : modsubim->copyData(ratio);
2000 :
2001 21 : }// 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 798 : }
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 1783 : void SIImageStore::makeImageBeamSet(Float psfcutoff, const Bool forcefit)
2072 : {
2073 1783 : clock_t begin = clock();
2074 :
2075 3566 : LogIO os( LogOrigin("SIImageStore","getPSFGaussian",WHERE) );
2076 : // For all chans/pols, call getPSFGaussian() and put it into ImageBeamSet(chan,pol).
2077 1783 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2078 1783 : uInt nx = itsImageShape[0];
2079 1783 : uInt ny = itsImageShape[1];
2080 1783 : uInt npol = itsImageShape[2];
2081 1783 : uInt nchan = itsImageShape[3];
2082 1783 : ImageInfo ii = psf()->imageInfo();
2083 1783 : ImageBeamSet iibeamset=ii.getBeamSet();
2084 1783 : if(iibeamset.nchan()==nchan && iibeamset.nstokes()==npol && forcefit==False){
2085 1064 : itsPSFBeams=iibeamset;
2086 1064 : itsRestoredBeams=iibeamset;
2087 1064 : return;
2088 : }
2089 719 : itsPSFBeams.resize( nchan, npol );
2090 719 : itsRestoredBeams.resize(nchan, npol);
2091 : // cout << "makeImBeamSet : imshape : " << itsImageShape << endl;
2092 :
2093 719 : String blankpsfs="";
2094 3471 : for( uInt chanid=0; chanid<nchan;chanid++) {
2095 5719 : for( uInt polid=0; polid<npol; polid++ ) {
2096 2967 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
2097 :
2098 2967 : IPosition substart(4,0,0,polid,chanid);
2099 2967 : IPosition substop(4,nx-1,ny-1,polid,chanid);
2100 2967 : Slicer psfslice(substart, substop,Slicer::endIsLast);
2101 5934 : SubImage<Float> subPsf( *psf() , psfslice, True );
2102 2967 : GaussianBeam beam;
2103 :
2104 2967 : Bool tryfit=True;
2105 :
2106 2967 : LatticeExprNode le( max(subPsf) );
2107 2967 : tryfit = le.getFloat()>0.0;
2108 2967 : if(tryfit)
2109 : {
2110 : try
2111 : {
2112 2863 : StokesImageUtil::FitGaussianPSF( subPsf, beam,psfcutoff );
2113 2863 : itsPSFBeams.setBeam( chanid, polid, beam );
2114 2863 : 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 2967 : }// end of pol loop
2130 : }// end of chan loop
2131 :
2132 719 : 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 719 : 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 719 : 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 3471 : for( uInt chanid=0; chanid<nchan;chanid++) {
2151 5719 : for( uInt polid=0; polid<npol; polid++ ) {
2152 2967 : 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 719 : ii.setBeams( itsPSFBeams );
2182 : {
2183 719 : LatticeLocker lock1(*(psf()), FileLocker::Write);
2184 719 : psf()->setImageInfo(ii);
2185 719 : psf()->unlock();
2186 719 : }
2187 719 : clock_t end = clock();
2188 719 : os << LogIO::NORMAL << "Time to fit Gaussian to PSF " << double(end - begin) / CLOCKS_PER_SEC << LogIO::POST;
2189 3911 : }// 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 1062 : ImageBeamSet SIImageStore::getBeamSet(Float psfcutoff)
2219 : {
2220 1062 : IPosition beamshp = itsPSFBeams.shape();
2221 1062 : AlwaysAssert( beamshp.nelements()==2 , AipsError );
2222 1062 : if( beamshp[0]==0 || beamshp[1]==0 ) {makeImageBeamSet(psfcutoff);}
2223 2124 : return itsPSFBeams;
2224 1062 : }
2225 :
2226 729 : void SIImageStore::printBeamSet(Bool verbose)
2227 : {
2228 1458 : LogIO os( LogOrigin("SIImageStore","printBeamSet",WHERE) );
2229 729 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2230 729 : if( itsImageShape[3] == 1 && itsImageShape[2]==1 )
2231 : {
2232 418 : GaussianBeam beam = itsPSFBeams.getBeam();
2233 418 : os << "Beam : " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2234 418 : }
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 729 : }
2279 :
2280 : /////////////////////////////// Restore all planes.
2281 :
2282 864 : void SIImageStore::restore(GaussianBeam& rbeam, String& usebeam, uInt term, Float psfcutoff)
2283 : {
2284 1728 : LogIO os( LogOrigin("SIImageStore","restore",WHERE) );
2285 : // << ". Optionally, PB-correct too." << LogIO::POST;
2286 :
2287 864 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2288 864 : Int nx = itsImageShape[0];
2289 864 : Int ny = itsImageShape[1];
2290 864 : Int npol = itsImageShape[2];
2291 864 : 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 864 : IPosition beamset = itsPSFBeams.shape();
2303 864 : AlwaysAssert( beamset.nelements()==2 , AipsError );
2304 864 : if( beamset[0] != nchan || beamset[1] != npol )
2305 : {
2306 :
2307 : // Get PSF Beams....
2308 237 : ImageInfo ii = psf()->imageInfo();
2309 237 : itsPSFBeams = ii.getBeamSet();
2310 237 : itsRestoredBeams=itsPSFBeams;
2311 237 : IPosition beamset2 = itsPSFBeams.shape();
2312 :
2313 237 : AlwaysAssert( beamset2.nelements()==2 , AipsError );
2314 237 : 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 237 : }
2321 :
2322 : // toggle for printing common beam to the log
2323 864 : 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 864 : 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 864 : 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 864 : Bool emptymodel = isModelEmpty();
2357 864 : 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 864 : LatticeLocker lock1(*(image(term)), FileLocker::Write);
2361 864 : LatticeLocker lock2(*(residual(term)), FileLocker::Write);
2362 864 : LatticeLocker lock3(*(model(term)), FileLocker::Read);
2363 : //// Use beamset to restore
2364 3686 : for( Int chanid=0; chanid<nchan;chanid++) {
2365 5853 : for( Int polid=0; polid<npol; polid++ ) {
2366 :
2367 3031 : IPosition substart(4,0,0,polid,chanid);
2368 3031 : IPosition substop(4,nx-1,ny-1,polid,chanid);
2369 3031 : Slicer imslice(substart, substop,Slicer::endIsLast);
2370 6062 : SubImage<Float> subRestored( *image(term) , imslice, True );
2371 6062 : SubImage<Float> subModel( *model(term) , imslice, False );
2372 6062 : SubImage<Float> subResidual( *residual(term) , imslice, True );
2373 :
2374 3031 : 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 3031 : if(!printcommonbeam) {
2378 2805 : 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 3031 : subRestored.set(0.0);
2385 3031 : 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 3031 : 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 2805 : 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 3031 : }// end of pol loop
2416 : }// end of chan loop
2417 :
2418 : try
2419 : {
2420 : //MSK//
2421 864 : if( hasPB() )
2422 : {
2423 842 : if( (image(term)->getDefaultMask()).matches("mask0") ) removeMask( image(term) );
2424 842 : copyMask(residual(term),image(term));
2425 : }
2426 :
2427 : // if(hasPB()){copyMask(residual(term),image(term));}
2428 864 : ImageInfo iminf = image(term)->imageInfo();
2429 : //iminf.setBeams( itsRestoredBeams);
2430 :
2431 864 : os << "Beam Set : Single beam : " << itsRestoredBeams.hasSingleBeam() << " Multi-beam : " << itsRestoredBeams.hasMultiBeam() << LogIO::DEBUG2;
2432 :
2433 864 : iminf.removeRestoringBeam();
2434 :
2435 864 : if( itsRestoredBeams.hasSingleBeam() )
2436 607 : { iminf.setRestoringBeam( itsRestoredBeams.getBeam() );}
2437 : else
2438 257 : {iminf.setBeams( itsRestoredBeams);}
2439 :
2440 864 : image(term)->setImageInfo(iminf);
2441 :
2442 864 : }
2443 0 : catch(AipsError &x)
2444 : {
2445 0 : throw( AipsError("Restoration Error : " + x.getMesg() ) );
2446 0 : }
2447 :
2448 864 : }// 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 4416 : Bool SIImageStore::getUseWeightImage(ImageInterface<Float>& target)
2778 : {
2779 4416 : Record miscinfo = target.miscInfo();
2780 : Bool useweightimage;
2781 4416 : if( miscinfo.isDefined("useweightimage") && miscinfo.dataType("useweightimage")==TpBool )
2782 4416 : { miscinfo.get( "useweightimage", useweightimage ); }
2783 0 : else { useweightimage = False; }
2784 :
2785 4416 : return useweightimage;
2786 4416 : }
2787 : /*
2788 : Bool SIImageStore::getUseWeightImage()
2789 : {
2790 : if( ! itsParentSumWt )
2791 : return False;
2792 : else
2793 : return getUseWeightImage( *itsParentSumWt );
2794 : }
2795 : */
2796 15896 : void SIImageStore::setUseWeightImage(ImageInterface<Float>& target, Bool useweightimage)
2797 : {
2798 15896 : Record miscinfo = target.miscInfo();
2799 15896 : miscinfo.define("useweightimage", useweightimage);
2800 15896 : LatticeLocker lock1(target, FileLocker::Write);
2801 15896 : target.setMiscInfo( miscinfo );
2802 15896 : }
2803 :
2804 :
2805 :
2806 1875 : Bool SIImageStore::divideImageByWeightVal( ImageInterface<Float>& target )
2807 : {
2808 :
2809 1875 : Array<Float> lsumwt;
2810 1875 : sumwt()->get( lsumwt , False ); // For MT, this will always pick the zeroth sumwt, which it should.
2811 1875 : LatticeLocker lock2(target, FileLocker::Write);
2812 :
2813 1875 : IPosition imshape = target.shape();
2814 :
2815 : //cerr << " SumWt : " << lsumwt << " sumwtshape : " << lsumwt.shape() << " image shape : " << imshape << endl;
2816 :
2817 1875 : AlwaysAssert( lsumwt.shape()[2] == imshape[2] , AipsError ); // polplanes
2818 1875 : AlwaysAssert( lsumwt.shape()[3] == imshape[3] , AipsError ); // chanplanes
2819 :
2820 1875 : Bool div=False; // flag to signal if division actually happened, or weights are all 1.0
2821 :
2822 3974 : for(Int pol=0; pol<lsumwt.shape()[2]; pol++)
2823 : {
2824 8316 : for(Int chan=0; chan<lsumwt.shape()[3]; chan++)
2825 : {
2826 6217 : IPosition pos(4,0,0,pol,chan);
2827 6217 : 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 6217 : chan, lsumwt.shape()[3],
2832 6217 : pol, lsumwt.shape()[2],
2833 12434 : target );
2834 6217 : if ( lsumwt(pos) > 1e-07 ) {
2835 12032 : LatticeExpr<Float> le( (*subim)/lsumwt(pos) );
2836 6016 : subim->copyData( le );
2837 6016 : }
2838 : else {
2839 201 : subim->set(0.0);
2840 : }
2841 6217 : div=True;
2842 6217 : }
2843 6217 : }
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 1875 : return div;
2852 1875 : }
2853 :
2854 900 : void SIImageStore::normPSF(Int term)
2855 : {
2856 :
2857 1921 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2858 : {
2859 4135 : 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 6228 : chan, itsImageShape[3],
2865 3114 : pol, itsImageShape[2],
2866 6228 : (*psf(term)) );
2867 :
2868 : std::shared_ptr<ImageInterface<Float> > subim0=makeSubImage(0,1,
2869 6228 : chan, itsImageShape[3],
2870 3114 : pol, itsImageShape[2],
2871 6228 : (*psf(0)) );
2872 :
2873 3114 : LatticeLocker lock1(*(subim), FileLocker::Write);
2874 :
2875 3114 : LatticeExprNode themax( max(*(subim0)) );
2876 3114 : Float maxim = themax.getFloat();
2877 :
2878 3114 : if ( maxim > 1e-07 )
2879 : {
2880 6026 : LatticeExpr<Float> normed( (*(subim)) / maxim );
2881 3013 : subim->copyData( normed );
2882 3013 : }
2883 : else
2884 : {
2885 101 : subim->set(0.0);
2886 : }
2887 3114 : }//chan
2888 : }//pol
2889 :
2890 900 : }
2891 :
2892 591 : void SIImageStore::calcSensitivity()
2893 : {
2894 1182 : LogIO os( LogOrigin("SIImageStore","calcSensitivity",WHERE) );
2895 :
2896 591 : Array<Float> lsumwt;
2897 591 : sumwt()->get( lsumwt , False ); // For MT, this will always pick the zeroth sumwt, which it should.
2898 :
2899 591 : IPosition shp( lsumwt.shape() );
2900 : //cout << "Sumwt shape : " << shp << " : " << lsumwt << endl;
2901 : //AlwaysAssert( shp.nelements()==4 , AipsError );
2902 :
2903 591 : os << "[" << itsImageName << "] Theoretical sensitivity (Jy/bm):" ;
2904 :
2905 591 : IPosition it(4,0,0,0,0);
2906 1189 : for ( it[0]=0; it[0]<shp[0]; it[0]++)
2907 1222 : for ( it[1]=0; it[1]<shp[1]; it[1]++)
2908 1351 : for ( it[2]=0; it[2]<shp[2]; it[2]++)
2909 3547 : for ( it[3]=0; it[3]<shp[3]; it[3]++)
2910 : {
2911 2820 : if( shp[0]>1 ){os << "f"<< it[0]+(it[1]*shp[0]) << ":" ;}
2912 2820 : if( shp[3]>1 ) { os << "c"<< it[3] << ":"; }
2913 2820 : if( shp[2]>1 ) { os << "p"<< it[2]<< ":" ; }
2914 2820 : if( lsumwt( it ) > 1e-07 )
2915 : {
2916 2699 : os << 1.0/(sqrt(lsumwt(it))) << " " ;
2917 : }
2918 : else
2919 : {
2920 121 : os << "none" << " ";
2921 : }
2922 : }
2923 :
2924 591 : os << LogIO::POST;
2925 :
2926 : // Array<Float> sens = 1.0/sqrtsumwt;
2927 :
2928 :
2929 591 : }
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 5171 : Float SIImageStore::getModelFlux(uInt term)
2971 : {
2972 : // LogIO os( LogOrigin("SIImageStore","getModelFlux",WHERE) );
2973 5171 : LatticeLocker lock2 (*(model(term)), FileLocker::Read);
2974 10338 : LatticeExprNode mflux( sum( *model(term) ) );
2975 5169 : Float modelflux = mflux.getFloat();
2976 : // Float modelflux = sum( model(term)->get() );
2977 :
2978 5169 : return modelflux;
2979 5169 : }
2980 :
2981 : // Check for non-zero model (this is different from getting model flux, for derived SIIMMT)
2982 1952 : Bool SIImageStore::isModelEmpty()
2983 : {
2984 1952 : if( !itsModel && (! hasModel()) ) return True;
2985 1288 : 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 5291 : Float SIImageStore::getPSFSidelobeLevel()
2995 : {
2996 10582 : LogIO os( LogOrigin("SIImageStore","getPSFSidelobeLevel",WHERE) );
2997 : //cerr << "*****PSF sidelobe "<< itsPSFSideLobeLevel << endl;
2998 : /// Calculate only once, store and return for all subsequent calls.
2999 5291 : if( itsPSFSideLobeLevel == 0.0 )
3000 : {
3001 :
3002 654 : ImageBeamSet thebeams = getBeamSet();
3003 654 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
3004 :
3005 : //------------------------------------------------------------
3006 654 : IPosition oneplaneshape( itsImageShape );
3007 654 : AlwaysAssert( oneplaneshape.nelements()==4, AipsError );
3008 654 : oneplaneshape[2]=1; oneplaneshape[3]=1;
3009 654 : TempImage<Float> psfbeam( oneplaneshape, itsCoordSys );
3010 :
3011 : // In a loop through channels, subtract out or mask out the main lobe
3012 654 : Float allmin=0.0, allmax=0.0;
3013 1417 : for(Int pol=0; pol<itsImageShape[2]; pol++)
3014 : {
3015 3407 : for(Int chan=0; chan<itsImageShape[3]; chan++)
3016 : {
3017 : std::shared_ptr<ImageInterface<Float> > onepsf=makeSubImage(0,1,
3018 5288 : chan, itsImageShape[3],
3019 2644 : pol, itsImageShape[2],
3020 5288 : (*psf()) );
3021 :
3022 :
3023 2644 : GaussianBeam beam = thebeams.getBeam( chan, pol );
3024 2644 : Vector<Float> abeam(3); // Holds bmaj, bmin, pa in asec, asec, deg
3025 2644 : abeam[0] = beam.getMajor().get("arcsec").getValue() * C::arcsec;
3026 2644 : abeam[1] = beam.getMinor().get("arcsec").getValue() * C::arcsec;
3027 2644 : abeam[2] = (beam.getPA().get("deg").getValue() + 90.0)* C::degree;
3028 :
3029 : //cout << "Beam : " << abeam << endl;
3030 :
3031 2644 : StokesImageUtil::MakeGaussianPSF( psfbeam, abeam, False);
3032 :
3033 : // storeImg( String("psfbeam.im"), psfbeam );
3034 :
3035 : //Subtract from PSF plane
3036 5288 : 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 2644 : LatticeExprNode minval_le( min( *onepsf ) );
3044 2644 : LatticeExprNode maxval_le( max( delobed ) );
3045 :
3046 2644 : Float minval = minval_le.getFloat();
3047 2644 : Float maxval = maxval_le.getFloat();
3048 :
3049 2644 : if( minval < allmin ) allmin = minval;
3050 2644 : if( maxval > allmax ) allmax = maxval;
3051 :
3052 : //cout << "Chan : " << chan << " minval : " << minval << " maxval : " << maxval << endl;
3053 :
3054 2644 : }//chan
3055 : }//pol
3056 :
3057 : //------------------------------------------------------------
3058 :
3059 654 : itsPSFSideLobeLevel = max( fabs(allmin), fabs(allmax) );
3060 :
3061 : //os << "PSF min : " << allmin << " max : " << allmax << " psfsidelobelevel : " << itsPSFSideLobeLevel << LogIO::POST;
3062 :
3063 654 : }// if changed.
3064 :
3065 : // LatticeExprNode psfside( min( *psf() ) );
3066 : // itsPSFSideLobeLevel = fabs( psfside.getFloat() );
3067 :
3068 : //cout << "PSF sidelobe level : " << itsPSFSideLobeLevel << endl;
3069 5291 : return itsPSFSideLobeLevel;
3070 5291 : }
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 12278 : void SIImageStore::init()
3273 : {
3274 12278 : imageExts.resize(MAX_IMAGE_IDS);
3275 :
3276 12278 : imageExts(MASK)=".mask";
3277 12278 : imageExts(PSF)=".psf";
3278 12278 : imageExts(MODEL)=".model";
3279 12278 : imageExts(RESIDUAL)=".residual";
3280 12278 : imageExts(WEIGHT)=".weight";
3281 12278 : imageExts(IMAGE)=".image";
3282 12278 : imageExts(SUMWT)=".sumwt";
3283 12278 : imageExts(GRIDWT)=".gridwt";
3284 12278 : imageExts(PB)=".pb";
3285 12278 : imageExts(FORWARDGRID)=".forward";
3286 12278 : imageExts(BACKWARDGRID)=".backward";
3287 12278 : imageExts(IMAGEPBCOR)=".image.pbcor";
3288 :
3289 12278 : itsParentPsf = itsPsf;
3290 12278 : itsParentModel=itsModel;
3291 12278 : itsParentResidual=itsResidual;
3292 12278 : itsParentWeight=itsWeight;
3293 12278 : itsParentImage=itsImage;
3294 12278 : itsParentSumWt=itsSumWt;
3295 12278 : itsParentMask=itsMask;
3296 12278 : itsParentImagePBcor=itsImagePBcor;
3297 :
3298 : // cout << "parent shape : " << itsParentImageShape << " shape : " << itsImageShape << endl;
3299 12278 : itsParentImageShape = itsImageShape;
3300 12278 : itsParentCoordSys = itsCoordSys;
3301 :
3302 12278 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 ) { itsImageShape=IPosition(4,0,0,0,0); }
3303 :
3304 12278 : itsOpened=0;
3305 :
3306 12278 : itsPSFSideLobeLevel=0.0;
3307 12278 : itsReadLock=nullptr;
3308 12278 : itsDataPolRep=StokesImageUtil::UNKNOWN; //Should throw an exception if it is not initialized properly
3309 12278 : }
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 :
|
