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: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Exceptions/Error.h>
29 : #include <iostream>
30 : #include <sstream>
31 :
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 :
36 : #include <casacore/casa/Logging.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : #include <casacore/casa/Logging/LogMessage.h>
39 : #include <casacore/casa/Logging/LogSink.h>
40 : #include <casacore/casa/Logging/LogMessage.h>
41 : #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 4576 : SIImageStore::SIImageStore()
110 : {
111 :
112 4576 : itsPsf.reset( );
113 4576 : itsModel.reset( );
114 4576 : itsResidual.reset( );
115 4576 : itsWeight.reset( );
116 4576 : itsImage.reset( );
117 4576 : itsMask.reset( );
118 4576 : itsGridWt.reset( );
119 4576 : itsPB.reset( );
120 4576 : itsImagePBcor.reset();
121 4576 : itsTempWorkIm.reset();
122 :
123 4576 : itsSumWt.reset( );
124 4576 : itsOverWrite=False;
125 4576 : itsUseWeight=False;
126 4576 : itsPBScaleFactor=1.0;
127 :
128 4576 : itsNFacets=1;
129 4576 : itsFacetId=0;
130 4576 : itsNChanChunks = 1;
131 4576 : itsChanId = 0;
132 4576 : itsNPolChunks = 1;
133 4576 : itsPolId = 0;
134 :
135 4576 : itsImageShape=IPosition(4,0,0,0,0);
136 4576 : itsImageName=String("");
137 4576 : itsCoordSys=CoordinateSystem();
138 4576 : itsMiscInfo=Record();
139 4576 : init();
140 :
141 :
142 : // validate();
143 :
144 4576 : }
145 :
146 : // Used from SynthesisNormalizer::makeImageStore()
147 741 : 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 741 : const Bool useweightimage)
155 : // TODO : Add parameter to indicate weight image shape.
156 : {
157 2223 : LogIO os( LogOrigin("SIImageStore","Open new Images",WHERE) );
158 :
159 :
160 741 : itsPsf.reset( );
161 741 : itsModel.reset( );
162 741 : itsResidual.reset( );
163 741 : itsWeight.reset( );
164 741 : itsImage.reset( );
165 741 : itsMask.reset( );
166 741 : itsGridWt.reset( );
167 741 : itsPB.reset( );
168 741 : itsImagePBcor.reset( );
169 741 : itsTempWorkIm.reset();
170 :
171 741 : itsSumWt.reset( );
172 741 : itsOverWrite=False; // Hard Coding this. See CAS-6937. overwrite;
173 741 : itsUseWeight=useweightimage;
174 741 : itsPBScaleFactor=1.0;
175 :
176 741 : itsNFacets=1;
177 741 : itsFacetId=0;
178 741 : itsNChanChunks = 1;
179 741 : itsChanId = 0;
180 741 : itsNPolChunks = 1;
181 741 : itsPolId = 0;
182 :
183 741 : itsImageName = imagename;
184 741 : itsCoordSys = imcoordsys;
185 741 : itsImageShape = imshape;
186 741 : itsObjectName = objectname;
187 741 : itsMiscInfo = miscinfo;
188 :
189 741 : init();
190 :
191 741 : validate();
192 741 : }
193 :
194 : // Used from SynthesisNormalizer::makeImageStore()
195 2850 : SIImageStore::SIImageStore(const String &imagename, const Bool ignorefacets, const Bool noRequireSumwt)
196 : {
197 8550 : LogIO os( LogOrigin("SIImageStore","Open existing Images",WHERE) );
198 :
199 :
200 2850 : itsPsf.reset( );
201 2850 : itsModel.reset( );
202 2850 : itsResidual.reset( );
203 2850 : itsWeight.reset( );
204 2850 : itsImage.reset( );
205 2850 : itsMask.reset( );
206 2850 : itsGridWt.reset( );
207 2850 : itsPB.reset( );
208 2850 : itsImagePBcor.reset( );
209 2850 : itsTempWorkIm.reset();
210 2850 : itsMiscInfo=Record();
211 :
212 2850 : itsSumWt.reset( );
213 2850 : itsNFacets=1;
214 2850 : itsFacetId=0;
215 2850 : itsNChanChunks = 1;
216 2850 : itsChanId = 0;
217 2850 : itsNPolChunks = 1;
218 2850 : itsPolId = 0;
219 :
220 2850 : itsOverWrite=False;
221 : //need to to this init now so that imageExts is initialized
222 2850 : init();
223 :
224 2850 : 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 5700 : if( doesImageExist(itsImageName+String(".residual")) ||
230 4132 : doesImageExist(itsImageName+String(".psf")) ||
231 2852 : doesImageExist(itsImageName+String(".model")) ||
232 2850 : doesImageExist(itsImageName+String(".gridwt")) ||
233 9832 : doesImageExist(itsImageName+String(".pb")) ||
234 2850 : doesImageExist(itsImageName+String(".weight"))
235 : )
236 : {
237 0 : std::shared_ptr<ImageInterface<Float> > imptr;
238 2850 : if( doesImageExist(itsImageName+String(".psf")) )
239 : {
240 2719 : buildImage( imptr, (itsImageName+String(".psf")) );
241 : // itsObjectName=imptr->imageInfo().objectName();
242 : // itsMiscInfo=imptr->miscInfo();
243 : }
244 131 : else if ( doesImageExist(itsImageName+String(".residual")) ){
245 129 : 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 2850 : itsObjectName=imptr->imageInfo().objectName();
271 2850 : itsImageShape=imptr->shape();
272 2850 : itsCoordSys = imptr->coordinates();
273 2850 : itsMiscInfo=imptr->miscInfo();
274 :
275 : }
276 : else
277 : {
278 0 : throw( AipsError( "PSF, Residual, Model Image (or sumwt) do not exist. Please create one of them." ) );
279 : }
280 :
281 6982 : if( doesImageExist(itsImageName+String(".residual")) ||
282 4132 : doesImageExist(itsImageName+String(".psf")) )
283 : {
284 2848 : if( doesImageExist(itsImageName+String(".sumwt")) )
285 : {
286 2763 : std::shared_ptr<ImageInterface<Float> > imptr;
287 : //imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt")) );
288 2763 : buildImage( imptr, (itsImageName+String(".sumwt")) );
289 2763 : itsNFacets = imptr->shape()[0];
290 2763 : itsFacetId = 0;
291 2763 : itsUseWeight = getUseWeightImage( *imptr );
292 2763 : 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 2763 : itsCoordSys = imptr->coordinates();
295 2763 : itsMiscInfo=imptr->miscInfo();
296 2763 : 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 : }
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 0 : 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 : }
322 : }
323 : }// if psf or residual exist...
324 :
325 2850 : if( ignorefacets==True ) itsNFacets= 1;
326 :
327 2850 : init();
328 :
329 2850 : validate();
330 2850 : }
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 750 : 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 750 : const Bool useweightimage)
355 : {
356 :
357 :
358 750 : itsPsf=psfim;
359 750 : itsModel=modelim;
360 750 : itsResidual=residim;
361 : /* if(residim){
362 : LatticeLocker lock1(*itsResidual, FileLocker::Read);
363 : cerr << "RESIDMAX " << max(itsResidual->get()) << " " << max(residim->get()) << endl;
364 : }*/
365 750 : itsWeight=weightim;
366 750 : itsImage=restoredim;
367 750 : itsMask=maskim;
368 :
369 750 : itsSumWt=sumwtim;
370 :
371 750 : itsGridWt=gridwtim;
372 750 : itsPB=pbim;
373 750 : itsImagePBcor=restoredpbcorim;
374 750 : itsTempWorkIm=tempworkimage;
375 :
376 750 : itsPBScaleFactor=1.0;// No need to set properly here as it will be computed in makePSF.
377 :
378 750 : itsObjectName = objectname;
379 750 : itsMiscInfo = miscinfo;
380 :
381 750 : itsNFacets = nfacets;
382 750 : itsFacetId = facet;
383 750 : itsNChanChunks = nchanchunks;
384 750 : itsChanId = chan;
385 750 : itsNPolChunks = npolchunks;
386 750 : itsPolId = pol;
387 :
388 750 : itsOverWrite=False;
389 750 : itsUseWeight=useweightimage;
390 :
391 750 : itsParentImageShape = imshape;
392 750 : itsImageShape = imshape;
393 :
394 750 : itsParentCoordSys = csys;
395 750 : itsCoordSys = csys; // Hopefully this doesn't change for a reference image
396 750 : itsImageName = imagename;
397 :
398 : //-----------------------
399 750 : init(); // Connect parent pointers to the images.
400 : //-----------------------
401 :
402 : // Set these to null, to be set later upon first access.
403 750 : itsPsf.reset( );
404 750 : itsModel.reset( );
405 750 : itsResidual.reset( );
406 750 : itsWeight.reset( );
407 750 : itsImage.reset( );
408 750 : itsMask.reset( );
409 750 : itsSumWt.reset( );
410 750 : itsPB.reset( );
411 :
412 750 : validate();
413 750 : }
414 :
415 8600 : void SIImageStore::validate()
416 : {
417 : /// There are two valid states. Check for at least one of them.
418 8600 : Bool state = False;
419 :
420 17200 : stringstream oss;
421 8600 : oss << "shape:" << itsImageShape << " parentimageshape:" << itsParentImageShape
422 8600 : << " nfacets:" << itsNFacets << "x" << itsNFacets << " facetid:" << itsFacetId
423 8600 : << " nchanchunks:" << itsNChanChunks << " chanid:" << itsChanId
424 8600 : << " npolchunks:" << itsNPolChunks << " polid:" << itsPolId
425 8600 : << " coord-dim:" << itsCoordSys.nPixelAxes()
426 8600 : << " psf/res:" << (hasPsf() || hasResidual()) ;
427 8600 : if( hasSumWt() ) oss << " sumwtshape : " << sumwt()->shape() ;
428 8600 : oss << endl;
429 :
430 :
431 : try {
432 :
433 8600 : if( itsCoordSys.nPixelAxes() != 4 ) state=False;
434 :
435 : /// (1) Regular imagestore
436 8600 : if( itsNFacets==1 && itsFacetId==0
437 8408 : && itsNChanChunks==1 && itsChanId==0
438 8408 : && itsNPolChunks==1 && itsPolId==0 ) {
439 8328 : Bool check1 = hasSumWt() && sumwt()->shape()[0]==1;
440 16656 : if( (itsImageShape.isEqual(itsParentImageShape) ) && ( check1 || !hasSumWt() )
441 16656 : && itsParentImageShape.product() > 0 ) state=True;
442 : }
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 : }
463 :
464 : // cout << " SIIM:validate : " << oss.str() << endl;
465 :
466 8600 : if( state == False ) throw(AipsError("Internal Error : Invalid ImageStore state : "+ oss.str()) );
467 :
468 17200 : return;
469 : }
470 :
471 :
472 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
473 :
474 750 : 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 750 : 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 30577 : 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 30577 : std::shared_ptr<ImageInterface<Float> > imPtr;
492 61154 : IPosition useShape( itsParentImageShape );
493 :
494 30577 : if( dosumwt ) // change shape to sumwt image shape.
495 : {
496 4193 : useShape[0] = nfacetsperside;
497 4193 : 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 30577 : Bool dbg=False;
517 30577 : if( doesImageExist( imagenamefull ) )
518 : {
519 : ///// not used since itsOverWrite is hardcoded to FALSE (CAS-6937)
520 24654 : 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 : }
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 : }
541 : }// is table writable
542 : else
543 : {
544 0 : throw( AipsError("Cannot overwrite existing image. : " + x.getMesg() ) );
545 : }
546 : }
547 : }// overwrite existing image
548 : else // open existing image ( Always tries this )
549 : {
550 : if(1) //Table::isWritable( imagenamefull ))
551 : {
552 24654 : if(dbg) cout << "Trying to open existing image : "<< imagenamefull << endl;
553 : try{
554 24654 : buildImage( imPtr, imagenamefull ) ;
555 :
556 24654 : 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 21399 : if( useShape.product()>0 && ! useShape.isEqual( imPtr->shape() ) )
562 : {
563 45 : 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 : }
567 :
568 :
569 :
570 21384 : 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 : }
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 34 : catch (AipsError &x){
600 17 : throw( AipsError("Cannot open existing image : "+imagenamefull+" : " + x.getMesg() ) );
601 : }
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 5923 : if(dbg) cout << "Trying to open new image named : " << imagenamefull << endl;
613 : try{
614 5923 : buildImage(imPtr, useShape, itsParentCoordSys, imagenamefull) ;
615 : // initialize to zeros...
616 : {
617 11846 : LatticeLocker lock1(*imPtr, FileLocker::Write);
618 5923 : imPtr->set(0.0);
619 5923 : imPtr->flush();
620 5923 : imPtr->unlock();
621 : }
622 : }
623 0 : catch (AipsError &x){
624 0 : throw( AipsError("Cannot make new image. : " + x.getMesg() ) );
625 : }
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 61120 : return imPtr;
662 : }
663 :
664 5923 : void SIImageStore::buildImage(std::shared_ptr<ImageInterface<Float> > &imptr, IPosition shape,
665 : CoordinateSystem csys, const String name)
666 : {
667 17769 : LogIO os( LogOrigin("SIImageStore", "Open non-existing image", WHERE) );
668 5923 : os <<"Opening image, name: " << name << LogIO::DEBUG1;
669 :
670 5923 : itsOpened++;
671 : //For some reason cannot open a new paged image with UserNoread directly
672 : {
673 11846 : PagedImage<Float> godot(shape, csys, name);
674 5923 : godot.unlock();
675 : }
676 5923 : 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 5923 : imptr.reset( new PagedImage<Float> (name, locktype) );
682 : {
683 11846 : LatticeLocker lock1(*imptr, FileLocker::Write);
684 5923 : initMetaInfo(imptr, name);
685 5923 : imptr->unlock();
686 : }
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 5923 : }
701 :
702 31599 : void SIImageStore::buildImage(std::shared_ptr<ImageInterface<Float> > &imptr, const String name)
703 : {
704 63198 : LogIO os(LogOrigin("SIImageStore", "Open existing Images", WHERE));
705 31599 : os <<"Opening image, name: " << name << LogIO::DEBUG1;
706 :
707 31599 : itsOpened++;
708 31599 : if(Table::isReadable(name)){
709 31599 : 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 31599 : Table table(name, locktype);
714 31599 : String type = table.tableInfo().type();
715 31599 : if (type != TableInfo::type(TableInfo::PAGEDIMAGE)) {
716 :
717 0 : imptr.reset( new PagedImage<Float>( table ) );
718 0 : imptr->unlock();
719 0 : return;
720 : }
721 : }
722 31599 : LatticeBase* latt =ImageOpener::openImage(name);
723 31599 : if(!latt)
724 : {
725 0 : throw(AipsError("Error in opening Image : "+name));
726 : }
727 31599 : DataType dtype=latt->dataType();
728 31599 : if(dtype==TpFloat)
729 : {
730 31599 : 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 : }
762 :
763 : /**
764 : * Sets ImageInfo and MiscInfo on an image
765 : *
766 : * @param imptr image to initialize
767 : */
768 5925 : 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 11850 : LatticeLocker lock1(*imptr, FileLocker::Write);
774 5925 : if (not itsObjectName.empty()) {
775 11850 : ImageInfo info = imptr->imageInfo();
776 5925 : info.setObjectName(itsObjectName);
777 5925 : imptr->setImageInfo(info);
778 5925 : imptr->setMiscInfo(itsMiscInfo);
779 0 : } 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 : }
786 5925 : imptr->unlock();
787 5925 : }
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 33308 : 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 33308 : Int nx_facets=Int(sqrt(Double(nfacets)));
810 99924 : LogIO os( LogOrigin("SynthesisImager","makeFacet") );
811 : // Make the output image
812 66616 : Slicer imageSlicer;
813 :
814 : // Add checks for all dimensions..
815 33308 : 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 66616 : IPosition imshp=image.shape();
820 66616 : IPosition blc(imshp.nelements(), 0);
821 66616 : IPosition trc=imshp-1;
822 66616 : IPosition inc(imshp.nelements(), 1);
823 :
824 : /// Facets
825 33308 : Int facetx = facet % nx_facets;
826 33308 : Int facety = (facet - facetx) / nx_facets;
827 33308 : Int sizex = imshp(0) / nx_facets;
828 33308 : Int sizey = imshp(1) / nx_facets;
829 33308 : blc(1) = facety * sizey;
830 33308 : trc(1) = blc(1) + sizey - 1;
831 33308 : blc(0) = facetx * sizex;
832 33308 : trc(0) = blc(0) + sizex - 1;
833 :
834 : /// Pol Chunks
835 33308 : Int sizepol = imshp(2) / npolchunks;
836 33308 : blc(2) = pol * sizepol;
837 33308 : trc(2) = blc(2) + sizepol - 1;
838 :
839 : /// Chan Chunks
840 33308 : Int sizechan = imshp(3) / nchanchunks;
841 33308 : blc(3) = chan * sizechan;
842 33308 : trc(3) = blc(3) + sizechan - 1;
843 :
844 33308 : LCBox::verify(blc, trc, inc, imshp);
845 66616 : Slicer imslice(blc, trc, inc, Slicer::endIsLast);
846 :
847 : // Now create the sub image
848 66616 : std::shared_ptr<ImageInterface<Float> > referenceImage( new SubImage<Float>(image, imslice, True) );
849 : {
850 66616 : LatticeLocker lock1(*referenceImage, FileLocker::Write);
851 33308 : referenceImage->setMiscInfo(image.miscInfo());
852 33308 : referenceImage->setUnits(image.units());
853 :
854 : }
855 :
856 : // cout << "Made Ref subimage of shape : " << referenceImage->shape() << endl;
857 :
858 33308 : return referenceImage;
859 :
860 : }
861 :
862 :
863 :
864 : //////////////////////////////////////////////////////////////////////////////////////////////////////
865 :
866 11767 : SIImageStore::~SIImageStore()
867 : {
868 11767 : }
869 :
870 : //////////////////////////////////////////////////////////////////////////////////////////////////////
871 : //////////////////////////////////////////////////////////////////////////////////////////////////////
872 :
873 12050 : Bool SIImageStore::releaseLocks()
874 : {
875 : //LogIO os( LogOrigin("SIImageStore","releaseLocks",WHERE) );
876 :
877 : // String fname( itsImageName+String(".info") );
878 : // makePersistent( fname );
879 :
880 12050 : if( itsPsf ) releaseImage( itsPsf );
881 12050 : if( itsModel ) { releaseImage( itsModel ); }
882 12050 : if( itsResidual ) releaseImage( itsResidual );
883 12050 : if( itsImage ) releaseImage( itsImage );
884 12050 : if( itsWeight ) releaseImage( itsWeight );
885 12050 : if( itsMask ) releaseImage( itsMask );
886 12050 : if( itsSumWt ) releaseImage( itsSumWt );
887 12050 : if( itsGridWt ) releaseImage( itsGridWt );
888 12050 : if( itsPB ) releaseImage( itsPB );
889 12050 : if( itsImagePBcor ) releaseImage( itsImagePBcor );
890 :
891 12050 : 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 25306 : 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 25306 : im->flush();
909 : //os << LogIO::WARN << "clear cache" << LogIO::POST;
910 25306 : im->clearCache();
911 : //os << LogIO::WARN << "unlock" << LogIO::POST;
912 25306 : im->unlock();
913 : //os << LogIO::WARN << "tempClose" << LogIO::POST;
914 : //im->tempClose();
915 : //os << LogIO::WARN << "NULL" << LogIO::POST;
916 25306 : im.reset(); // This was added to allow modification by modules independently
917 25306 : }
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 : }
932 :
933 729 : Long SIImageStore::estimateRAM(){
934 : //right now this is estimated at 2MB for the 2 complex lattices;
935 729 : return Long(2000);
936 : }
937 24 : void SIImageStore::setModelImage( const Vector<String> &modelnames)
938 : {
939 72 : 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 23 : }
948 :
949 :
950 :
951 43 : void SIImageStore::setModelImageOne( const String &modelname , Int nterm)
952 : {
953 88 : 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 43 : 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 9 : TempImage<Float> maskmodel( newmodel->shape(), newmodel->coordinates() );
993 6 : 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 6 : 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 9 : (*newmodel) );
1003 :
1004 : std::shared_ptr<ImageInterface<Float> > submaskmodel=makeSubImage(0,1,
1005 6 : chan, inshape[3],
1006 3 : pol, inshape[2],
1007 6 : maskmodel );
1008 :
1009 6 : ArrayLattice<Bool> pixmask( subim->getMask() );
1010 9 : LatticeExpr<Float> masked( (*subim) * iif(pixmask,1.0,0.0) );
1011 3 : submaskmodel->copyData( masked );
1012 : }
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 5 : 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 0 : } catch (AipsError &x) {
1031 0 : throw(AipsError("Setting masked pixels to zero for input startmodel : "+ x.getMesg()));
1032 : }
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 20 : 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 33 : 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 : }
1051 :
1052 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1053 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1054 :
1055 9338 : IPosition SIImageStore::getShape()
1056 : {
1057 9338 : return itsImageShape;
1058 : }
1059 :
1060 6523 : String SIImageStore::getName()
1061 : {
1062 6523 : return itsImageName;
1063 : }
1064 :
1065 8121 : uInt SIImageStore::getNTaylorTerms(Bool /*dopsf*/)
1066 : {
1067 8121 : 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 125415 : 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 125415 : Bool sw=False;
1087 125415 : if( label.contains(imageExts(SUMWT)) ) sw=True;
1088 :
1089 125415 : if( !ptr )
1090 : {
1091 : //cout << itsNFacets << " " << itsNChanChunks << " " << itsNPolChunks << endl;
1092 31178 : 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 1460 : 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 30362 : if(!parentptr){
1120 : ///coordsys for psf can be different ...shape should be the same.
1121 29715 : ptr = openImage(itsImageName+label , itsOverWrite, sw, 1, !(label.contains(imageExts(PSF)))); }
1122 : else{
1123 647 : 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 125398 : }
1136 :
1137 :
1138 13867 : std::shared_ptr<ImageInterface<Float> > SIImageStore::psf(uInt /*nterm*/)
1139 : {
1140 13867 : accessImage( itsPsf, itsParentPsf, imageExts(PSF) );
1141 :
1142 13866 : return itsPsf;
1143 : }
1144 30411 : std::shared_ptr<ImageInterface<Float> > SIImageStore::residual(uInt /*nterm*/)
1145 : {
1146 30411 : accessImage( itsResidual, itsParentResidual, imageExts(RESIDUAL) );
1147 : // Record mi = itsResidual->miscInfo(); ostringstream oss;mi.print(oss);cout<<"MiscInfo(res) : " << oss.str() << endl;
1148 30411 : return itsResidual;
1149 : }
1150 1635 : std::shared_ptr<ImageInterface<Float> > SIImageStore::weight(uInt /*nterm*/)
1151 : {
1152 1635 : accessImage( itsWeight, itsParentWeight, imageExts(WEIGHT) );
1153 1635 : return itsWeight;
1154 : }
1155 :
1156 9578 : std::shared_ptr<ImageInterface<Float> > SIImageStore::sumwt(uInt /*nterm*/)
1157 : {
1158 :
1159 9578 : accessImage( itsSumWt, itsParentSumWt, imageExts(SUMWT) );
1160 :
1161 9578 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
1162 536 : { itsUseWeight = getUseWeightImage( *itsParentSumWt );}
1163 9578 : setUseWeightImage( *itsSumWt , itsUseWeight); // Sets a flag in the SumWt image.
1164 :
1165 9578 : return itsSumWt;
1166 : }
1167 :
1168 10768 : std::shared_ptr<ImageInterface<Float> > SIImageStore::model(uInt /*nterm*/)
1169 : {
1170 10770 : accessImage( itsModel, itsParentModel, imageExts(MODEL) );
1171 21532 : LatticeLocker lock1(*itsModel, FileLocker::Write);
1172 : // Set up header info the first time.
1173 10766 : itsModel->setUnits("Jy/pixel");
1174 :
1175 21532 : return itsModel;
1176 : }
1177 5912 : std::shared_ptr<ImageInterface<Float> > SIImageStore::image(uInt /*nterm*/)
1178 : {
1179 5912 : accessImage( itsImage, itsParentImage, imageExts(IMAGE) );
1180 11824 : LatticeLocker lock1(*itsImage, FileLocker::Write);
1181 5912 : itsImage->setUnits(Unit("Jy/beam"));
1182 11824 : return itsImage;
1183 : }
1184 :
1185 14697 : std::shared_ptr<ImageInterface<Float> > SIImageStore::mask(uInt /*nterm*/)
1186 : {
1187 14697 : accessImage( itsMask, itsParentMask, imageExts(MASK) );
1188 14693 : 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 : }
1205 2 : itsGridWt.reset(new PagedImage<Float>(newshape, newcoordsys, itsImageName+ imageExts(GRIDWT)));
1206 2 : initMetaInfo(itsGridWt, itsImageName+ imageExts(GRIDWT));
1207 :
1208 : }
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 12080 : std::shared_ptr<ImageInterface<Float> > SIImageStore::pb(uInt /*nterm*/)
1224 : {
1225 12080 : accessImage( itsPB, itsParentPB, imageExts(PB) );
1226 12078 : return itsPB;
1227 : }
1228 261 : std::shared_ptr<ImageInterface<Float> > SIImageStore::imagepbcor(uInt /*nterm*/)
1229 : {
1230 261 : accessImage( itsImagePBcor, itsParentImagePBcor, imageExts(IMAGEPBCOR) );
1231 522 : LatticeLocker lock1(*itsImagePBcor, FileLocker::Write);
1232 261 : itsImagePBcor->setUnits(Unit("Jy/beam"));
1233 522 : return itsImagePBcor;
1234 : }
1235 :
1236 1900 : std::shared_ptr<ImageInterface<Complex> > SIImageStore::forwardGrid(uInt /*nterm*/){
1237 1900 : if( itsForwardGrid ) // && (itsForwardGrid->shape() == itsImageShape))
1238 : {
1239 : // cout << "Forward grid has shape : " << itsForwardGrid->shape() << endl;
1240 1672 : return itsForwardGrid;
1241 : }
1242 456 : Vector<Int> whichStokes(0);
1243 456 : IPosition cimageShape;
1244 228 : cimageShape=itsImageShape;
1245 228 : 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 228 : if(freqframe != MFrequency::LSRK && freqframe!=MFrequency::Undefined && freqframe!=MFrequency::REST)
1248 0 : { itsCoordSys.setSpectralConversion("LSRK"); }
1249 228 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
1250 456 : whichStokes, itsDataPolRep);
1251 228 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
1252 228 : cimageShape(2)=whichStokes.nelements();
1253 : //cout << "Making forward grid of shape : " << cimageShape << " for imshape : " << itsImageShape << endl;
1254 228 : itsForwardGrid.reset( new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()) );
1255 : //if(image())
1256 228 : if(hasRestored()){
1257 15 : LatticeLocker lock1(*itsForwardGrid, FileLocker::Write);
1258 15 : itsForwardGrid->setImageInfo((image())->imageInfo());
1259 :
1260 : }
1261 228 : return itsForwardGrid;
1262 : }
1263 :
1264 2450 : std::shared_ptr<ImageInterface<Complex> > SIImageStore::backwardGrid(uInt /*nterm*/){
1265 2450 : if( itsBackwardGrid ) //&& (itsBackwardGrid->shape() == itsImageShape))
1266 : {
1267 : // cout << "Backward grid has shape : " << itsBackwardGrid->shape() << endl;
1268 1982 : return itsBackwardGrid;
1269 : }
1270 936 : Vector<Int> whichStokes(0);
1271 936 : IPosition cimageShape;
1272 468 : cimageShape=itsImageShape;
1273 468 : 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 468 : if(freqframe != MFrequency::LSRK && freqframe!=MFrequency::Undefined && freqframe!=MFrequency::REST )
1276 0 : { itsCoordSys.setSpectralConversion("LSRK"); }
1277 468 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
1278 936 : whichStokes, itsDataPolRep);
1279 468 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
1280 468 : cimageShape(2)=whichStokes.nelements();
1281 : //cout << "Making backward grid of shape : " << cimageShape << " for imshape : " << itsImageShape << endl;
1282 468 : itsBackwardGrid.reset( new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()) );
1283 : //if(image())
1284 468 : if(hasRestored()){
1285 7 : LatticeLocker lock1(*itsBackwardGrid, FileLocker::Write);
1286 7 : itsBackwardGrid->setImageInfo((image())->imageInfo());
1287 : }
1288 468 : return itsBackwardGrid;
1289 : }
1290 2222 : Double SIImageStore::memoryBeforeLattice(){
1291 : //Calculate how much memory to use per temporary images before disking
1292 2222 : return 1.0; /// in MB
1293 : }
1294 2222 : IPosition SIImageStore::tileShape(){
1295 : //Need to have settable stuff here or algorith to determine this
1296 2222 : 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 45871 : Bool SIImageStore::doesImageExist(String imagename)
1301 : {
1302 137613 : LogIO os( LogOrigin("SIImageStore","doesImageExist",WHERE) );
1303 91742 : Directory image( imagename );
1304 91742 : return image.exists();
1305 : }
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 : }
1315 0 : if( resetresidual ) {
1316 : // removeMask( residual() );
1317 0 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1318 0 : residual()->set(0.0);
1319 : }
1320 0 : if( resetweight && itsWeight ){
1321 0 : LatticeLocker lock1(*(weight()), FileLocker::Write);
1322 0 : weight()->set(0.0);
1323 : }
1324 0 : if( resetweight ){
1325 0 : LatticeLocker lock1(*(sumwt()), FileLocker::Write);
1326 0 : sumwt()->set(0.0);
1327 : }
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 : }
1339 0 : if(addresidual)
1340 : {
1341 0 : LatticeExpr<Float> adderRes( *(residual()) + *(imagestoadd->residual()) );
1342 0 : residual()->copyData(adderRes);
1343 : }
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 : }
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 : }
1364 0 : if(adddensity)
1365 : {
1366 0 : LatticeExpr<Float> adderDensity( *(gridwt()) + *(imagestoadd->gridwt()) );
1367 0 : gridwt()->copyData(adderDensity);
1368 : }
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 : }
1404 :
1405 965 : 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 965 : CountedPtr<ImageInterface<Float> > subim=makeSubImage(0, 1,
1411 965 : chan, itsImageShape[3],
1412 965 : pol, itsImageShape[2],
1413 1930 : *weight(0));
1414 :
1415 1930 : return sqrt(max(subim->get()));
1416 : }
1417 :
1418 104 : void SIImageStore::makePBFromWeight(const Float pblimit)
1419 : {
1420 208 : LogIO os( LogOrigin("SIImageStore","makePBFromWeight",WHERE) );
1421 :
1422 208 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1423 : {
1424 315 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1425 : {
1426 :
1427 211 : itsPBScaleFactor = getPbMax(pol,chan);
1428 :
1429 211 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1430 : else {
1431 :
1432 208 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1433 208 : chan, itsImageShape[3],
1434 208 : pol, itsImageShape[2],
1435 624 : *weight() );
1436 208 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
1437 208 : chan, itsImageShape[3],
1438 208 : pol, itsImageShape[2],
1439 624 : *pb() );
1440 :
1441 :
1442 624 : LatticeExpr<Float> normed( sqrt(abs(*wtsubim)) / itsPBScaleFactor );
1443 624 : LatticeExpr<Float> limited( iif( normed > fabs(pblimit) , normed, 0.0 ) );
1444 208 : pbsubim->copyData( limited );
1445 : }// if not zero
1446 : }
1447 : }
1448 :
1449 104 : 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 160 : LatticeExpr<Bool> pbmask( iif( *pb() > fabs(pblimit) , True , False ) );
1458 : //MSK//
1459 80 : createMask( pbmask, pb() );
1460 80 : pb()->pixelMask().unlock();
1461 : }
1462 :
1463 : }
1464 104 : weight()->unlock();
1465 104 : pb()->unlock();
1466 104 : }
1467 :
1468 572 : void SIImageStore::makePBImage(const Float pblimit)
1469 : {
1470 1144 : LogIO os( LogOrigin("SIImageStore","makePBImage",WHERE) );
1471 :
1472 572 : if( hasPB() )
1473 : {
1474 :
1475 1253 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1476 : {
1477 3308 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1478 : {
1479 :
1480 2627 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
1481 2627 : chan, itsImageShape[3],
1482 2627 : pol, itsImageShape[2],
1483 7881 : *pb() );
1484 :
1485 : // Norm by the max
1486 5254 : LatticeExprNode elmax= max( *pbsubim );
1487 2627 : 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 2627 : if(fmax>1.0)
1491 : {
1492 351 : LatticeExpr<Float> normed( (*pbsubim) / fmax );
1493 351 : LatticeExpr<Float> limited( iif( normed > fabs(pblimit) , normed, 0.0 ) );
1494 117 : pbsubim->copyData( limited );
1495 : }
1496 : else
1497 : {
1498 7530 : LatticeExpr<Float> limited( iif((*pbsubim) > fabs(pblimit) , (*pbsubim), 0.0 ) );
1499 2510 : pbsubim->copyData( limited );
1500 : }
1501 : }
1502 : }
1503 :
1504 572 : if((pb()->getDefaultMask()==""))// && pblimit >= 0.0)
1505 : {
1506 : // removeMask( pb() );
1507 : //MSK//
1508 1112 : LatticeExpr<Bool> pbmask( iif( *pb() > fabs(pblimit) , True , False ) );
1509 : //MSK//
1510 556 : createMask( pbmask, pb() );
1511 556 : pb()->pixelMask().unlock();
1512 : }
1513 :
1514 : }// if hasPB
1515 572 : pb()->unlock();
1516 :
1517 572 : }
1518 :
1519 862 : 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 1724 : LatticeLocker lock1(*outimage, FileLocker::Write);
1526 1724 : ImageRegion outreg = outimage->makeMask("mask0",False,True);
1527 862 : LCRegion& outmask=outreg.asMask();
1528 862 : outmask.copyData(lemask);
1529 862 : outimage->defineRegion("mask0",outreg, RegionHandler::Masks, True);
1530 :
1531 862 : outimage->setDefaultMask("mask0");
1532 :
1533 862 : outimage->unlock();
1534 862 : outimage->tempClose();
1535 :
1536 : // outimage->table().unmarkForDelete();
1537 : }
1538 0 : catch (const AipsError& x) {
1539 0 : throw(AipsError("Error in creating internal T/F mask : " + x.getMesg() ));
1540 : }
1541 :
1542 862 : return True;
1543 : }
1544 :
1545 1695 : 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 1695 : if( (inimage->getDefaultMask()).matches("mask0") ) // input mask exists.
1553 3026 : {LatticeLocker lock1(*outimage, FileLocker::Write);
1554 1513 : 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 1516 : ImageRegion outreg=outimage->makeMask("mask0", False, True);
1563 1510 : LCRegion& outmask=outreg.asMask();
1564 1510 : outmask.copyData(inimage->getRegion("mask0").asLCRegion());
1565 1510 : outimage->defineRegion("mask0",outreg, RegionHandler::Masks,True);
1566 1510 : outimage->setDefaultMask("mask0");
1567 : }
1568 : }
1569 6 : catch (const AipsError& x) {
1570 3 : throw(AipsError("Error in copying internal T/F mask : " + x.getMesg() ));
1571 : }
1572 :
1573 1692 : return True;
1574 : }
1575 :
1576 1825 : 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 3650 : LatticeLocker lock1(*im, FileLocker::Write);
1586 1825 : if (im-> getDefaultMask() != String(""))
1587 : {
1588 748 : String strung=im->getDefaultMask();
1589 374 : im->setDefaultMask("");
1590 374 : im->removeRegion(strung);
1591 : }
1592 1825 : if( im->hasRegion("mask0") )
1593 : {
1594 0 : im->removeRegion("mask0");
1595 : }
1596 : }
1597 0 : catch (const AipsError& x)
1598 : {
1599 0 : throw(AipsError("Error in deleting internal T/F mask : " + x.getMesg() ));
1600 : }
1601 1825 : }
1602 226 : void SIImageStore:: rescaleResolution(Int chan,
1603 : ImageInterface<Float>& image,
1604 : const GaussianBeam& newbeam,
1605 : const GaussianBeam& oldbeam){
1606 :
1607 678 : LogIO os( LogOrigin("SIImageStore","rescaleResolution",WHERE) );
1608 452 : GaussianBeam toBeUsed(Quantity(0.0, "arcsec"),Quantity(0.0, "arcsec"),
1609 904 : 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 2 : 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 : }
1640 :
1641 :
1642 226 : }
1643 :
1644 :
1645 :
1646 :
1647 571 : void SIImageStore::dividePSFByWeight(const Float /* pblimit*/)
1648 : {
1649 1713 : LogIO os( LogOrigin("SIImageStore","dividePSFByWeight",WHERE) );
1650 :
1651 571 : LatticeLocker lock1 (*(psf()), FileLocker::Write);
1652 571 : normPSF();
1653 :
1654 571 : if( itsUseWeight )
1655 : {
1656 :
1657 64 : divideImageByWeightVal( *weight() );
1658 : }
1659 571 : (psf())->unlock();
1660 :
1661 571 : }
1662 :
1663 571 : void SIImageStore::normalizePrimaryBeam(const Float pblimit)
1664 : {
1665 1142 : LogIO os( LogOrigin("SIImageStore","normalizePrimaryBeam",WHERE) );
1666 :
1667 : // cout << "In dividePSFByWeight : itsUseWeight : " << itsUseWeight << endl;
1668 571 : if( itsUseWeight )
1669 : {
1670 :
1671 128 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1672 : {
1673 235 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1674 : {
1675 171 : 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 64 : makePBFromWeight(pblimit);
1680 :
1681 : }//if itsUseWeight
1682 507 : else { makePBImage(pblimit); } // OR... just check that it exists already.
1683 571 : pb()->unlock();
1684 :
1685 571 : }
1686 :
1687 : // Make another for the PSF too.
1688 1296 : void SIImageStore::divideResidualByWeight(Float pblimit, String normtype) {
1689 :
1690 3888 : LogIO os(LogOrigin("SIImageStore", "divideResidualByWeight", WHERE));
1691 1296 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1692 :
1693 328 : auto logTemplate = [&](Int const &chan, Int const &pol, string const &normalizer, string const &result) {
1694 328 : os << LogIO::NORMAL1
1695 656 : << "[C" + String::toString(chan) + ":P" + String::toString(pol) + "] "
1696 656 : << "Dividing " << itsImageName + String(".residual") << " by "
1697 : << "[ " << normalizer << " ] "
1698 984 : << "to get " << result << "." << LogIO::POST;
1699 328 : };
1700 :
1701 : // Normalize by the sumwt, per plane.
1702 1296 : Bool didNorm = divideImageByWeightVal(*residual());
1703 :
1704 1296 : if (itsUseWeight) {
1705 256 : for (Int pol = 0; pol < itsImageShape[2]; pol++) {
1706 465 : for (Int chan = 0; chan < itsImageShape[3]; chan++) {
1707 337 : itsPBScaleFactor = getPbMax(pol, chan);
1708 :
1709 337 : 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 328 : CountedPtr<ImageInterface<Float> > wtsubim = makeSubImage(0, 1,
1716 328 : chan, itsImageShape[3],
1717 328 : pol, itsImageShape[2],
1718 984 : *weight());
1719 328 : CountedPtr<ImageInterface<Float> > ressubim = makeSubImage(0, 1,
1720 328 : chan, itsImageShape[3],
1721 328 : pol, itsImageShape[2],
1722 984 : *residual());
1723 656 : LatticeExpr<Float> ratio;
1724 328 : Float scalepb = 1.0;
1725 :
1726 328 : if (normtype == "flatnoise") {
1727 984 : logTemplate(chan, pol,
1728 656 : "sqrt(weightimage) * " + String::toString(itsPBScaleFactor),
1729 : "flat noise with unit pb peak");
1730 :
1731 656 : LatticeExpr<Float> deno = itsPBScaleFactor * sqrt(abs(LatticeExpr<Float>(*(wtsubim))));
1732 328 : scalepb = fabs(pblimit) * itsPBScaleFactor * itsPBScaleFactor;
1733 328 : ratio = iif(deno > scalepb, (*(ressubim) / deno), 0.0);
1734 :
1735 0 : } 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 : }
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 328 : 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 : } // 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 1296 : 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 1296 : if ((residual()->getDefaultMask() == "") && hasPB() && pblimit >=0.0) {
1779 220 : copyMask(pb(), residual());
1780 : }
1781 :
1782 1296 : if ((pblimit < 0.0) && (residual()->getDefaultMask()).matches("mask0")) {
1783 1 : removeMask(residual());
1784 : }
1785 :
1786 1296 : residual()->unlock();
1787 1296 : }
1788 :
1789 129 : void SIImageStore::divideResidualByWeightSD(Float pblimit) {
1790 :
1791 387 : LogIO os(LogOrigin("SIImageStore", "divideResidualByWeightSD", WHERE));
1792 129 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1793 :
1794 129 : if (itsUseWeight) {
1795 387 : LatticeExpr<Float> deno = LatticeExpr<Float>(*weight());
1796 258 : LatticeExpr<Float> ratio = iif(deno > 0.0, *(residual()) / deno, 0.0);
1797 129 : residual()->copyData(ratio);
1798 : }
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 129 : if ((residual()->getDefaultMask() == "") && hasPB() && pblimit >=0.0) {
1808 0 : copyMask(pb(), residual());
1809 : }
1810 :
1811 129 : if ((pblimit < 0.0) && (residual()->getDefaultMask()).matches("mask0")) {
1812 0 : removeMask(residual());
1813 : }
1814 :
1815 129 : residual()->unlock();
1816 129 : }
1817 :
1818 850 : void SIImageStore::divideModelByWeight(Float pblimit, const String normtype)
1819 : {
1820 1700 : LogIO os( LogOrigin("SIImageStore","divideModelByWeight",WHERE) );
1821 :
1822 :
1823 1700 : if(itsUseWeight // only when needed
1824 850 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
1825 : {
1826 :
1827 84 : 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 : }
1835 84 : else if( normtype=="flatnoise"){
1836 :
1837 168 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1838 : {
1839 313 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1840 : {
1841 :
1842 229 : itsPBScaleFactor = getPbMax(pol,chan);
1843 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1844 229 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1845 : else {
1846 :
1847 203 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1848 203 : chan, itsImageShape[3],
1849 203 : pol, itsImageShape[2],
1850 609 : *weight() );
1851 203 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1852 203 : chan, itsImageShape[3],
1853 203 : pol, itsImageShape[2],
1854 609 : *model() );
1855 203 : os << LogIO::NORMAL1 ;
1856 203 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1857 203 : os << "Dividing " << itsImageName+String(".model") ;
1858 203 : os << " by [ sqrt(weight) / " << itsPBScaleFactor ;
1859 203 : os <<" ] to get to flat sky model before prediction" << LogIO::POST;
1860 :
1861 :
1862 609 : LatticeExpr<Float> deno( sqrt( abs(*(wtsubim)) ) / itsPBScaleFactor );
1863 :
1864 609 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1865 609 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1866 609 : 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 406 : 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 203 : 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 : }// 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 : }// if not zero
1941 : }//chan
1942 : }//pol
1943 :
1944 : }
1945 :
1946 : // storeImg(String("flatmodel.im"), *model());
1947 :
1948 : }
1949 : }
1950 :
1951 784 : void SIImageStore::multiplyModelByWeight(Float pblimit, const String normtype)
1952 : {
1953 1568 : LogIO os( LogOrigin("SIImageStore","multiplyModelByWeight",WHERE) );
1954 :
1955 1568 : if(itsUseWeight // only when needed
1956 784 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
1957 : {
1958 17 : 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 17 : else if( normtype=="flatnoise"){
1963 :
1964 34 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1965 : {
1966 34 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1967 : {
1968 :
1969 17 : itsPBScaleFactor = getPbMax(pol,chan);
1970 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1971 17 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1972 : else {
1973 :
1974 17 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1975 17 : chan, itsImageShape[3],
1976 17 : pol, itsImageShape[2],
1977 51 : *weight() );
1978 17 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1979 17 : chan, itsImageShape[3],
1980 17 : pol, itsImageShape[2],
1981 51 : *model() );
1982 :
1983 :
1984 :
1985 17 : os << LogIO::NORMAL1 ;
1986 17 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1987 17 : os << "Multiplying " << itsImageName+String(".model") ;
1988 17 : os << " by [ sqrt(weight) / " << itsPBScaleFactor;
1989 17 : os << " ] to take model back to flat noise with unit pb peak." << LogIO::POST;
1990 :
1991 51 : LatticeExpr<Float> deno( sqrt( abs(*(wtsubim)) ) / itsPBScaleFactor );
1992 :
1993 51 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1994 51 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1995 51 : 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 17 : modsubim->copyData(ratio);
2000 :
2001 : }// 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 : }// if not zero
2044 : }//chan
2045 : }//pol
2046 : }
2047 : }
2048 : }
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 : }
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 : }
2067 :
2068 0 : return beam;
2069 : }
2070 :
2071 1738 : void SIImageStore::makeImageBeamSet(Float psfcutoff, const Bool forcefit)
2072 : {
2073 1738 : clock_t begin = clock();
2074 :
2075 3476 : LogIO os( LogOrigin("SIImageStore","getPSFGaussian",WHERE) );
2076 : // For all chans/pols, call getPSFGaussian() and put it into ImageBeamSet(chan,pol).
2077 1738 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2078 1738 : uInt nx = itsImageShape[0];
2079 1738 : uInt ny = itsImageShape[1];
2080 1738 : uInt npol = itsImageShape[2];
2081 1738 : uInt nchan = itsImageShape[3];
2082 1738 : ImageInfo ii = psf()->imageInfo();
2083 1738 : ImageBeamSet iibeamset=ii.getBeamSet();
2084 1738 : if(iibeamset.nchan()==nchan && iibeamset.nstokes()==npol && forcefit==False){
2085 1043 : itsPSFBeams=iibeamset;
2086 1043 : itsRestoredBeams=iibeamset;
2087 1043 : return;
2088 : }
2089 695 : itsPSFBeams.resize( nchan, npol );
2090 695 : itsRestoredBeams.resize(nchan, npol);
2091 : // cout << "makeImBeamSet : imshape : " << itsImageShape << endl;
2092 :
2093 1390 : String blankpsfs="";
2094 3383 : for( uInt chanid=0; chanid<nchan;chanid++) {
2095 5591 : for( uInt polid=0; polid<npol; polid++ ) {
2096 5806 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
2097 :
2098 5806 : IPosition substart(4,0,0,polid,chanid);
2099 5806 : IPosition substop(4,nx-1,ny-1,polid,chanid);
2100 5806 : Slicer psfslice(substart, substop,Slicer::endIsLast);
2101 8709 : SubImage<Float> subPsf( *psf() , psfslice, True );
2102 5806 : GaussianBeam beam;
2103 :
2104 2903 : Bool tryfit=True;
2105 :
2106 5806 : LatticeExprNode le( max(subPsf) );
2107 2903 : tryfit = le.getFloat()>0.0;
2108 2903 : if(tryfit)
2109 : {
2110 : try
2111 : {
2112 2799 : StokesImageUtil::FitGaussianPSF( subPsf, beam,psfcutoff );
2113 2799 : itsPSFBeams.setBeam( chanid, polid, beam );
2114 2799 : 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 : }
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 : }// end of pol loop
2130 : }// end of chan loop
2131 :
2132 695 : 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 1390 : 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 695 : if(defaultbeam.getArea("rad2")==0.0){
2146 10 : 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 : }
2150 3383 : for( uInt chanid=0; chanid<nchan;chanid++) {
2151 5591 : for( uInt polid=0; polid<npol; polid++ ) {
2152 2903 : 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 695 : ii.setBeams( itsPSFBeams );
2182 : {
2183 695 : LatticeLocker lock1(*(psf()), FileLocker::Write);
2184 695 : psf()->setImageInfo(ii);
2185 695 : psf()->unlock();
2186 : }
2187 695 : clock_t end = clock();
2188 695 : os << LogIO::NORMAL << "Time to fit Gaussian to PSF " << double(end - begin) / CLOCKS_PER_SEC << LogIO::POST;
2189 : }// 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 : }
2213 264 : void SIImageStore::setBeamSet(const ImageBeamSet& bs){
2214 :
2215 264 : itsPSFBeams=bs;
2216 264 : }
2217 :
2218 1034 : ImageBeamSet SIImageStore::getBeamSet(Float psfcutoff)
2219 : {
2220 2068 : IPosition beamshp = itsPSFBeams.shape();
2221 1034 : AlwaysAssert( beamshp.nelements()==2 , AipsError );
2222 1034 : if( beamshp[0]==0 || beamshp[1]==0 ) {makeImageBeamSet(psfcutoff);}
2223 2068 : return itsPSFBeams;
2224 : }
2225 :
2226 709 : void SIImageStore::printBeamSet(Bool verbose)
2227 : {
2228 2127 : LogIO os( LogOrigin("SIImageStore","printBeamSet",WHERE) );
2229 709 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2230 709 : if( itsImageShape[3] == 1 && itsImageShape[2]==1 )
2231 : {
2232 406 : GaussianBeam beam = itsPSFBeams.getBeam();
2233 406 : os << "Beam : " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2234 : }
2235 : else
2236 : {
2237 : // per CAS-11415, verbose=True when restoringbeam='common'
2238 303 : 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 : }//for chanid
2249 9 : os << oss.str() << LogIO::POST;
2250 : }
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 294 : 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 : }//for chanid
2275 : }// catch
2276 : }
2277 : }
2278 709 : }
2279 :
2280 : /////////////////////////////// Restore all planes.
2281 :
2282 837 : void SIImageStore::restore(GaussianBeam& rbeam, String& usebeam, uInt term, Float psfcutoff)
2283 : {
2284 2511 : LogIO os( LogOrigin("SIImageStore","restore",WHERE) );
2285 : // << ". Optionally, PB-correct too." << LogIO::POST;
2286 :
2287 837 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2288 837 : Int nx = itsImageShape[0];
2289 837 : Int ny = itsImageShape[1];
2290 837 : Int npol = itsImageShape[2];
2291 837 : 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 1674 : IPosition beamset = itsPSFBeams.shape();
2303 837 : AlwaysAssert( beamset.nelements()==2 , AipsError );
2304 837 : if( beamset[0] != nchan || beamset[1] != npol )
2305 : {
2306 :
2307 : // Get PSF Beams....
2308 444 : ImageInfo ii = psf()->imageInfo();
2309 222 : itsPSFBeams = ii.getBeamSet();
2310 222 : itsRestoredBeams=itsPSFBeams;
2311 444 : IPosition beamset2 = itsPSFBeams.shape();
2312 :
2313 222 : AlwaysAssert( beamset2.nelements()==2 , AipsError );
2314 222 : 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 : }
2321 :
2322 : // toggle for printing common beam to the log
2323 837 : 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 837 : 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 837 : 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 34 : 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 : }// 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 837 : Bool emptymodel = isModelEmpty();
2357 837 : 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 1674 : LatticeLocker lock1(*(image(term)), FileLocker::Write);
2361 1674 : LatticeLocker lock2(*(residual(term)), FileLocker::Write);
2362 1674 : LatticeLocker lock3(*(model(term)), FileLocker::Read);
2363 : //// Use beamset to restore
2364 3623 : for( Int chanid=0; chanid<nchan;chanid++) {
2365 5781 : for( Int polid=0; polid<npol; polid++ ) {
2366 :
2367 5990 : IPosition substart(4,0,0,polid,chanid);
2368 5990 : IPosition substop(4,nx-1,ny-1,polid,chanid);
2369 5990 : Slicer imslice(substart, substop,Slicer::endIsLast);
2370 8985 : SubImage<Float> subRestored( *image(term) , imslice, True );
2371 8985 : SubImage<Float> subModel( *model(term) , imslice, False );
2372 8985 : SubImage<Float> subResidual( *residual(term) , imslice, True );
2373 :
2374 5990 : 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 2995 : if(!printcommonbeam) {
2378 2769 : 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 2995 : subRestored.set(0.0);
2385 2995 : if( !emptymodel ) {
2386 : // Copy model into it
2387 2685 : subRestored.copyData( LatticeExpr<Float>( subModel ) );
2388 : // Smooth model by beam
2389 2685 : StokesImageUtil::Convolve( subRestored, beam);
2390 : }
2391 : // Add residual image
2392 2995 : 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 : }
2401 : else// if no need to rescale residuals, just add the residuals.
2402 : {
2403 :
2404 :
2405 2769 : 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 : }
2414 :
2415 : }// end of pol loop
2416 : }// end of chan loop
2417 :
2418 : try
2419 : {
2420 : //MSK//
2421 837 : if( hasPB() )
2422 : {
2423 815 : if( (image(term)->getDefaultMask()).matches("mask0") ) removeMask( image(term) );
2424 815 : copyMask(residual(term),image(term));
2425 : }
2426 :
2427 : // if(hasPB()){copyMask(residual(term),image(term));}
2428 837 : ImageInfo iminf = image(term)->imageInfo();
2429 : //iminf.setBeams( itsRestoredBeams);
2430 :
2431 837 : os << "Beam Set : Single beam : " << itsRestoredBeams.hasSingleBeam() << " Multi-beam : " << itsRestoredBeams.hasMultiBeam() << LogIO::DEBUG2;
2432 :
2433 837 : iminf.removeRestoringBeam();
2434 :
2435 837 : if( itsRestoredBeams.hasSingleBeam() )
2436 581 : { iminf.setRestoringBeam( itsRestoredBeams.getBeam() );}
2437 : else
2438 256 : {iminf.setBeams( itsRestoredBeams);}
2439 :
2440 837 : image(term)->setImageInfo(iminf);
2441 :
2442 : }
2443 0 : catch(AipsError &x)
2444 : {
2445 0 : throw( AipsError("Restoration Error : " + x.getMesg() ) );
2446 : }
2447 :
2448 837 : }// 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 : }//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 : }// 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 41 : void SIImageStore::pbcor(uInt term)
2691 : {
2692 :
2693 82 : LogIO os( LogOrigin("SIImageStore","pbcor",WHERE) );
2694 :
2695 41 : 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 82 : LatticeLocker lock1(*(image(term)), FileLocker::Write);
2702 82 : LatticeLocker lock2(*(pb(term)), FileLocker::Write);
2703 82 : LatticeLocker lock3(*(imagepbcor(term)), FileLocker::Write);
2704 :
2705 82 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2706 : {
2707 156 : for(Int chan=0; chan<itsImageShape[3]; chan++)
2708 : {
2709 :
2710 115 : CountedPtr<ImageInterface<Float> > restoredsubim=makeSubImage(0,1,
2711 115 : chan, itsImageShape[3],
2712 115 : pol, itsImageShape[2],
2713 345 : *image(term) );
2714 115 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
2715 115 : chan, itsImageShape[3],
2716 115 : pol, itsImageShape[2],
2717 345 : *pb() );
2718 :
2719 115 : CountedPtr<ImageInterface<Float> > pbcorsubim=makeSubImage(0,1,
2720 115 : chan, itsImageShape[3],
2721 115 : pol, itsImageShape[2],
2722 345 : *imagepbcor(term) );
2723 :
2724 :
2725 230 : LatticeExprNode pbmax( max( *pbsubim ) );
2726 115 : Float pbmaxval = pbmax.getFloat();
2727 :
2728 115 : 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 345 : LatticeExpr<Float> thepbcor( iif( *(pbsubim) > 0.0 , (*(restoredsubim))/(*(pbsubim)) , 0.0 ) );
2737 115 : pbcorsubim->copyData( thepbcor );
2738 :
2739 : }// if not zero
2740 : }//chan
2741 : }//pol
2742 :
2743 : // Copy over the PB mask.
2744 41 : if((imagepbcor(term)->getDefaultMask()=="") && hasPB())
2745 23 : {copyMask(pb(),imagepbcor(term));}
2746 :
2747 : // Set restoring beam info
2748 41 : ImageInfo iminf = image(term)->imageInfo();
2749 : //iminf.setBeams( itsRestoredBeams );
2750 41 : imagepbcor(term)->setImageInfo(iminf);
2751 :
2752 : }// 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 : }
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 4149 : Bool SIImageStore::getUseWeightImage(ImageInterface<Float>& target)
2778 : {
2779 4149 : Record miscinfo = target.miscInfo();
2780 : Bool useweightimage;
2781 4149 : if( miscinfo.isDefined("useweightimage") && miscinfo.dataType("useweightimage")==TpBool )
2782 4149 : { miscinfo.get( "useweightimage", useweightimage ); }
2783 0 : else { useweightimage = False; }
2784 :
2785 8298 : return useweightimage;
2786 : }
2787 : /*
2788 : Bool SIImageStore::getUseWeightImage()
2789 : {
2790 : if( ! itsParentSumWt )
2791 : return False;
2792 : else
2793 : return getUseWeightImage( *itsParentSumWt );
2794 : }
2795 : */
2796 15298 : void SIImageStore::setUseWeightImage(ImageInterface<Float>& target, Bool useweightimage)
2797 : {
2798 30596 : Record miscinfo = target.miscInfo();
2799 15298 : miscinfo.define("useweightimage", useweightimage);
2800 30596 : LatticeLocker lock1(target, FileLocker::Write);
2801 15298 : target.setMiscInfo( miscinfo );
2802 15298 : }
2803 :
2804 :
2805 :
2806 1833 : Bool SIImageStore::divideImageByWeightVal( ImageInterface<Float>& target )
2807 : {
2808 :
2809 3666 : Array<Float> lsumwt;
2810 1833 : sumwt()->get( lsumwt , False ); // For MT, this will always pick the zeroth sumwt, which it should.
2811 3666 : LatticeLocker lock2(target, FileLocker::Write);
2812 :
2813 1833 : IPosition imshape = target.shape();
2814 :
2815 : //cerr << " SumWt : " << lsumwt << " sumwtshape : " << lsumwt.shape() << " image shape : " << imshape << endl;
2816 :
2817 1833 : AlwaysAssert( lsumwt.shape()[2] == imshape[2] , AipsError ); // polplanes
2818 1833 : AlwaysAssert( lsumwt.shape()[3] == imshape[3] , AipsError ); // chanplanes
2819 :
2820 1833 : Bool div=False; // flag to signal if division actually happened, or weights are all 1.0
2821 :
2822 3890 : for(Int pol=0; pol<lsumwt.shape()[2]; pol++)
2823 : {
2824 8123 : for(Int chan=0; chan<lsumwt.shape()[3]; chan++)
2825 : {
2826 12132 : IPosition pos(4,0,0,pol,chan);
2827 6066 : 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 6066 : chan, lsumwt.shape()[3],
2832 6066 : pol, lsumwt.shape()[2],
2833 12132 : target );
2834 6066 : if ( lsumwt(pos) > 1e-07 ) {
2835 17595 : LatticeExpr<Float> le( (*subim)/lsumwt(pos) );
2836 5865 : subim->copyData( le );
2837 : }
2838 : else {
2839 201 : subim->set(0.0);
2840 : }
2841 6066 : div=True;
2842 : }
2843 : }
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 3666 : return div;
2852 : }
2853 :
2854 880 : void SIImageStore::normPSF(Int term)
2855 : {
2856 :
2857 1881 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2858 : {
2859 4055 : 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 6108 : chan, itsImageShape[3],
2865 3054 : pol, itsImageShape[2],
2866 9162 : (*psf(term)) );
2867 :
2868 : std::shared_ptr<ImageInterface<Float> > subim0=makeSubImage(0,1,
2869 6108 : chan, itsImageShape[3],
2870 3054 : pol, itsImageShape[2],
2871 9162 : (*psf(0)) );
2872 :
2873 6108 : LatticeLocker lock1(*(subim), FileLocker::Write);
2874 :
2875 6108 : LatticeExprNode themax( max(*(subim0)) );
2876 3054 : Float maxim = themax.getFloat();
2877 :
2878 3054 : if ( maxim > 1e-07 )
2879 : {
2880 8859 : LatticeExpr<Float> normed( (*(subim)) / maxim );
2881 2953 : subim->copyData( normed );
2882 : }
2883 : else
2884 : {
2885 101 : subim->set(0.0);
2886 : }
2887 : }//chan
2888 : }//pol
2889 :
2890 880 : }
2891 :
2892 571 : void SIImageStore::calcSensitivity()
2893 : {
2894 1713 : LogIO os( LogOrigin("SIImageStore","calcSensitivity",WHERE) );
2895 :
2896 1142 : Array<Float> lsumwt;
2897 571 : sumwt()->get( lsumwt , False ); // For MT, this will always pick the zeroth sumwt, which it should.
2898 :
2899 1142 : IPosition shp( lsumwt.shape() );
2900 : //cout << "Sumwt shape : " << shp << " : " << lsumwt << endl;
2901 : //AlwaysAssert( shp.nelements()==4 , AipsError );
2902 :
2903 571 : os << "[" << itsImageName << "] Theoretical sensitivity (Jy/bm):" ;
2904 :
2905 1142 : IPosition it(4,0,0,0,0);
2906 1149 : for ( it[0]=0; it[0]<shp[0]; it[0]++)
2907 1182 : for ( it[1]=0; it[1]<shp[1]; it[1]++)
2908 1311 : for ( it[2]=0; it[2]<shp[2]; it[2]++)
2909 3467 : for ( it[3]=0; it[3]<shp[3]; it[3]++)
2910 : {
2911 2760 : if( shp[0]>1 ){os << "f"<< it[0]+(it[1]*shp[0]) << ":" ;}
2912 2760 : if( shp[3]>1 ) { os << "c"<< it[3] << ":"; }
2913 2760 : if( shp[2]>1 ) { os << "p"<< it[2]<< ":" ; }
2914 2760 : if( lsumwt( it ) > 1e-07 )
2915 : {
2916 2639 : os << 1.0/(sqrt(lsumwt(it))) << " " ;
2917 : }
2918 : else
2919 : {
2920 121 : os << "none" << " ";
2921 : }
2922 : }
2923 :
2924 571 : os << LogIO::POST;
2925 :
2926 : // Array<Float> sens = 1.0/sqrtsumwt;
2927 :
2928 :
2929 571 : }
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 8401 : Float SIImageStore::getPeakResidual()
2942 : {
2943 25203 : LogIO os( LogOrigin("SIImageStore","getPeakResidual",WHERE) );
2944 16802 : LatticeLocker lock2 (*(residual()), FileLocker::Read);
2945 25203 : ArrayLattice<Bool> pixelmask(residual()->getMask());
2946 25203 : LatticeExpr<Float> resd(iif(pixelmask,abs(*residual()),0));
2947 : //LatticeExprNode pres( max(abs( *residual() ) ));
2948 8401 : LatticeExprNode pres( max(resd) );
2949 8401 : Float maxresidual = pres.getFloat();
2950 :
2951 16802 : return maxresidual;
2952 : }
2953 :
2954 4883 : Float SIImageStore::getPeakResidualWithinMask()
2955 : {
2956 14649 : LogIO os( LogOrigin("SIImageStore","getPeakResidualWithinMask",WHERE) );
2957 : //Float minresmask, maxresmask, minres, maxres;
2958 : //findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
2959 :
2960 14649 : ArrayLattice<Bool> pixelmask(residual()->getMask());
2961 : //findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
2962 14649 : LatticeExpr<Float> resd(iif(((*mask()) > 0) && pixelmask ,abs((*residual())),0));
2963 4883 : LatticeExprNode pres( max(resd) );
2964 4883 : Float maxresidual = pres.getFloat();
2965 : //Float maxresidual=max( abs(maxresmask), abs(minresmask) );
2966 9766 : return maxresidual;
2967 : }
2968 :
2969 : // Calculate the total model flux
2970 5106 : Float SIImageStore::getModelFlux(uInt term)
2971 : {
2972 : // LogIO os( LogOrigin("SIImageStore","getModelFlux",WHERE) );
2973 10210 : LatticeLocker lock2 (*(model(term)), FileLocker::Read);
2974 10208 : LatticeExprNode mflux( sum( *model(term) ) );
2975 5104 : Float modelflux = mflux.getFloat();
2976 : // Float modelflux = sum( model(term)->get() );
2977 :
2978 10208 : return modelflux;
2979 : }
2980 :
2981 : // Check for non-zero model (this is different from getting model flux, for derived SIIMMT)
2982 1906 : Bool SIImageStore::isModelEmpty()
2983 : {
2984 1906 : if( !itsModel && (! hasModel()) ) return True;
2985 1263 : else return ( fabs( getModelFlux(0) ) < 1e-08 );
2986 : }
2987 :
2988 :
2989 264 : void SIImageStore::setPSFSidelobeLevel(const Float level){
2990 :
2991 264 : itsPSFSideLobeLevel=level;
2992 264 : }
2993 : // Calculate the PSF sidelobe level...
2994 2799 : Float SIImageStore::getPSFSidelobeLevel()
2995 : {
2996 5598 : LogIO os( LogOrigin("SIImageStore","getPSFSidelobeLevel",WHERE) );
2997 : //cerr << "*****PSF sidelobe "<< itsPSFSideLobeLevel << endl;
2998 : /// Calculate only once, store and return for all subsequent calls.
2999 2799 : if( itsPSFSideLobeLevel == 0.0 )
3000 : {
3001 :
3002 1270 : ImageBeamSet thebeams = getBeamSet();
3003 1270 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
3004 :
3005 : //------------------------------------------------------------
3006 1270 : IPosition oneplaneshape( itsImageShape );
3007 635 : AlwaysAssert( oneplaneshape.nelements()==4, AipsError );
3008 635 : oneplaneshape[2]=1; oneplaneshape[3]=1;
3009 635 : TempImage<Float> psfbeam( oneplaneshape, itsCoordSys );
3010 :
3011 : // In a loop through channels, subtract out or mask out the main lobe
3012 635 : Float allmin=0.0, allmax=0.0;
3013 1379 : for(Int pol=0; pol<itsImageShape[2]; pol++)
3014 : {
3015 3365 : for(Int chan=0; chan<itsImageShape[3]; chan++)
3016 : {
3017 : std::shared_ptr<ImageInterface<Float> > onepsf=makeSubImage(0,1,
3018 5242 : chan, itsImageShape[3],
3019 2621 : pol, itsImageShape[2],
3020 7863 : (*psf()) );
3021 :
3022 :
3023 5242 : GaussianBeam beam = thebeams.getBeam( chan, pol );
3024 5242 : Vector<Float> abeam(3); // Holds bmaj, bmin, pa in asec, asec, deg
3025 2621 : abeam[0] = beam.getMajor().get("arcsec").getValue() * C::arcsec;
3026 2621 : abeam[1] = beam.getMinor().get("arcsec").getValue() * C::arcsec;
3027 2621 : abeam[2] = (beam.getPA().get("deg").getValue() + 90.0)* C::degree;
3028 :
3029 : //cout << "Beam : " << abeam << endl;
3030 :
3031 2621 : StokesImageUtil::MakeGaussianPSF( psfbeam, abeam, False);
3032 :
3033 : // storeImg( String("psfbeam.im"), psfbeam );
3034 :
3035 : //Subtract from PSF plane
3036 7863 : 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 5242 : LatticeExprNode minval_le( min( *onepsf ) );
3044 5242 : LatticeExprNode maxval_le( max( delobed ) );
3045 :
3046 2621 : Float minval = minval_le.getFloat();
3047 2621 : Float maxval = maxval_le.getFloat();
3048 :
3049 2621 : if( minval < allmin ) allmin = minval;
3050 2621 : if( maxval > allmax ) allmax = maxval;
3051 :
3052 : //cout << "Chan : " << chan << " minval : " << minval << " maxval : " << maxval << endl;
3053 :
3054 : }//chan
3055 : }//pol
3056 :
3057 : //------------------------------------------------------------
3058 :
3059 635 : itsPSFSideLobeLevel = max( fabs(allmin), fabs(allmax) );
3060 :
3061 : //os << "PSF min : " << allmin << " max : " << allmax << " psfsidelobelevel : " << itsPSFSideLobeLevel << LogIO::POST;
3062 :
3063 : }// if changed.
3064 :
3065 : // LatticeExprNode psfside( min( *psf() ) );
3066 : // itsPSFSideLobeLevel = fabs( psfside.getFloat() );
3067 :
3068 : //cout << "PSF sidelobe level : " << itsPSFSideLobeLevel << endl;
3069 5598 : return itsPSFSideLobeLevel;
3070 : }
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 339 : LogIO os( LogOrigin("SIImageStore","calcRobustRMS",WHERE) );
3089 113 : Record* regionPtr=0;
3090 226 : String LELmask("");
3091 226 : LatticeLocker lockres (*(residual()), FileLocker::Read);
3092 339 : ArrayLattice<Bool> pbmasklat(residual()->shape());
3093 113 : pbmasklat.set(False);
3094 226 : 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 226 : 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 : }
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 : }
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 226 : 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 : }
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 : }
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 : }
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 5680 : Float SIImageStore::getMaskSum()
3193 : {
3194 17040 : LogIO os( LogOrigin("SIImageStore","getMaskSum",WHERE) );
3195 11356 : LatticeLocker lock2 (*(mask()), FileLocker::Read);
3196 11352 : LatticeExprNode msum( sum( *mask() ) );
3197 5676 : Float masksum = msum.getFloat();
3198 :
3199 : // Float masksum = sum( mask()->get() );
3200 :
3201 11352 : return masksum;
3202 : }
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 : }
3252 :
3253 : }
3254 : }
3255 :
3256 0 : return True;
3257 :
3258 :
3259 : }
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 11767 : void SIImageStore::init()
3273 : {
3274 11767 : imageExts.resize(MAX_IMAGE_IDS);
3275 :
3276 11767 : imageExts(MASK)=".mask";
3277 11767 : imageExts(PSF)=".psf";
3278 11767 : imageExts(MODEL)=".model";
3279 11767 : imageExts(RESIDUAL)=".residual";
3280 11767 : imageExts(WEIGHT)=".weight";
3281 11767 : imageExts(IMAGE)=".image";
3282 11767 : imageExts(SUMWT)=".sumwt";
3283 11767 : imageExts(GRIDWT)=".gridwt";
3284 11767 : imageExts(PB)=".pb";
3285 11767 : imageExts(FORWARDGRID)=".forward";
3286 11767 : imageExts(BACKWARDGRID)=".backward";
3287 11767 : imageExts(IMAGEPBCOR)=".image.pbcor";
3288 :
3289 11767 : itsParentPsf = itsPsf;
3290 11767 : itsParentModel=itsModel;
3291 11767 : itsParentResidual=itsResidual;
3292 11767 : itsParentWeight=itsWeight;
3293 11767 : itsParentImage=itsImage;
3294 11767 : itsParentSumWt=itsSumWt;
3295 11767 : itsParentMask=itsMask;
3296 11767 : itsParentImagePBcor=itsImagePBcor;
3297 :
3298 : // cout << "parent shape : " << itsParentImageShape << " shape : " << itsImageShape << endl;
3299 11767 : itsParentImageShape = itsImageShape;
3300 11767 : itsParentCoordSys = itsCoordSys;
3301 :
3302 11767 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 ) { itsImageShape=IPosition(4,0,0,0,0); }
3303 :
3304 11767 : itsOpened=0;
3305 :
3306 11767 : itsPSFSideLobeLevel=0.0;
3307 11767 : itsReadLock=nullptr;
3308 11767 : itsDataPolRep=StokesImageUtil::UNKNOWN; //Should throw an exception if it is not initialized properly
3309 11767 : }
3310 :
3311 :
3312 17 : void SIImageStore::regridToModelImage( ImageInterface<Float> &inputimage, Int term )
3313 : {
3314 : try
3315 : {
3316 :
3317 : //output coords
3318 34 : IPosition outshape = itsImageShape;
3319 34 : CoordinateSystem outcsys = itsCoordSys;
3320 34 : DirectionCoordinate outDirCsys = outcsys.directionCoordinate();
3321 34 : SpectralCoordinate outSpecCsys = outcsys.spectralCoordinate();
3322 :
3323 : // do regrid
3324 34 : IPosition axes(3,0, 1, 2);
3325 34 : IPosition inshape = inputimage.shape();
3326 34 : CoordinateSystem incsys = inputimage.coordinates();
3327 34 : DirectionCoordinate inDirCsys = incsys.directionCoordinate();
3328 34 : SpectralCoordinate inSpecCsys = incsys.spectralCoordinate();
3329 :
3330 34 : 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 19 : ImageRegrid<Float> imregrid;
3337 19 : imregrid.regrid( *(model(term)), Interpolate2D::LINEAR, axes, inputimage );
3338 4 : } catch (AipsError &x) {
3339 2 : throw(AipsError("ImageRegrid error : "+ x.getMesg()));
3340 : }
3341 :
3342 4 : }catch(AipsError &x)
3343 : {
3344 2 : throw(AipsError("Error in regridding input model image to target coordsys : " + x.getMesg()));
3345 : }
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 : }
3467 :
3468 0 : }
3469 : //////////////
3470 15 : Bool SIImageStore::intersectComplexImage(const String& ComplexImageName){
3471 30 : Vector<Int> whichStokes(0);
3472 15 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
3473 30 : whichStokes, itsDataPolRep);
3474 :
3475 :
3476 : //cerr <<"itsDataPolRep " << itsDataPolRep << endl;
3477 :
3478 15 : CountedPtr<PagedImage<Complex> > compliantImage =nullptr;
3479 : {
3480 30 : PagedImage<Complex>inputImage(ComplexImageName);
3481 30 : IPosition theShape=itsImageShape;
3482 15 : theShape(0)=inputImage.shape()(0);
3483 15 : theShape(1)=inputImage.shape()(1);
3484 30 : CoordinateSystem inpcsys=inputImage.coordinates();
3485 30 : 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 45 : String tmpImage=File::newUniqueName(".", "TempImage").absoluteName();
3490 15 : compliantImage=new PagedImage<Complex>(theShape, cimageCoord, tmpImage);
3491 15 : compliantImage->set(0.0);
3492 30 : IPosition iblc(theShape.nelements(),0);
3493 30 : IPosition itrc=theShape-1;
3494 : //cerr << "blc " << iblc << " trc " << itrc << " shape " << theShape << endl;
3495 30 : LCBox lbox(iblc, itrc, theShape);
3496 30 : ImageRegion imagreg(WCBox(lbox, cimageCoord));
3497 :
3498 :
3499 30 : 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 : }
3504 15 : TableUtil::deleteTable(ComplexImageName);
3505 15 : compliantImage->rename(ComplexImageName);
3506 30 : return True;
3507 :
3508 : }
3509 :
3510 :
3511 : //////////////////////////////////////////////////////////////////////////////////////////////////////
3512 : //////////////////////////////////////////////////////////////////////////////////////////////////////
3513 :
3514 : } //# NAMESPACE CASA - END
3515 :
|
