Line data Source code
1 : //# SIImageStore.cc: Implementation of Imager.h
2 : //# Copyright (C) 1997-2008
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Exceptions/Error.h>
29 : #include <iostream>
30 : #include <sstream>
31 :
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 :
36 : #include <casacore/casa/Logging.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : #include <casacore/casa/Logging/LogMessage.h>
39 : #include <casacore/casa/Logging/LogSink.h>
40 : #include <casacore/casa/Logging/LogMessage.h>
41 : #include <casacore/casa/OS/DirectoryIterator.h>
42 : #include <casacore/casa/OS/File.h>
43 : #include <casacore/casa/OS/Path.h>
44 : #include <casacore/lattices/Lattices/LatticeLocker.h>
45 : #include <casacore/casa/OS/HostInfo.h>
46 : #include <components/ComponentModels/GaussianDeconvolver.h>
47 : #include <casacore/images/Images/TempImage.h>
48 : #include <casacore/images/Images/PagedImage.h>
49 : #include <imageanalysis/ImageAnalysis/CasaImageBeamSet.h>
50 : #include <casacore/lattices/LatticeMath/LatticeMathUtil.h>
51 : #include <casacore/ms/MeasurementSets/MSHistoryHandler.h>
52 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
53 : #include <casacore/tables/Tables/TableUtil.h>
54 :
55 : #include <synthesis/ImagerObjects/SIImageStore.h>
56 : #include <synthesis/ImagerObjects/SDMaskHandler.h>
57 : #include <synthesis/TransformMachines/StokesImageUtil.h>
58 : #include <synthesis/TransformMachines2/Utils.h>
59 : #include <synthesis/ImagerObjects/SynthesisUtilMethods.h>
60 : #include <casacore/images/Images/ImageRegrid.h>
61 : #include <imageanalysis/ImageAnalysis/ImageStatsCalculator.h>
62 :
63 : //#include <imageanalysis/ImageAnalysis/ImageMaskedPixelReplacer.h>
64 :
65 : #include <sys/types.h>
66 : #include <unistd.h>
67 : using namespace std;
68 :
69 : using namespace casacore;
70 :
71 : namespace casa { //# NAMESPACE CASA - BEGIN
72 :
73 : //
74 : //===========================================================================
75 : // Global method that I (SB) could make work in SynthesisUtilsMethods.
76 : //
77 : template <class T>
78 0 : void openImage(const String& imagenamefull,std::shared_ptr<ImageInterface<T> >& imPtr )
79 : {
80 0 : LogIO logIO ( LogOrigin("SynthesisImager","openImage(name)") );
81 : try
82 : {
83 0 : if (Table::isReadable(imagenamefull))
84 0 : imPtr.reset( new PagedImage<T>( imagenamefull ) );
85 : }
86 0 : catch (AipsError &x)
87 : {
88 0 : logIO << "Error in reading image \"" << imagenamefull << "\"" << LogIO::EXCEPTION;
89 : }
90 0 : }
91 : //
92 : //--------------------------------------------------------------
93 : //
94 : template
95 : void openImage(const String& imagenamefull, std::shared_ptr<ImageInterface<Float> >& img );
96 : template
97 : void openImage(const String& imagenamefull, std::shared_ptr<ImageInterface<Complex> >& img );
98 : //
99 : //===========================================================================
100 :
101 :
102 : //////////////////////////////////////////////////////////////////////////////////////////////////////
103 : //////////////////////////////////////////////////////////////////////////////////////////////////////
104 : ///// START of SIIMAGESTORE implementation
105 : //////////////////////////////////////////////////////////////////////////////////////////////////////
106 :
107 : //////////////////////////////////////////////////////////////////////////////////////////////////////
108 :
109 0 : SIImageStore::SIImageStore()
110 : {
111 :
112 0 : itsPsf.reset( );
113 0 : itsModel.reset( );
114 0 : itsResidual.reset( );
115 0 : itsWeight.reset( );
116 0 : itsImage.reset( );
117 0 : itsMask.reset( );
118 0 : itsGridWt.reset( );
119 0 : itsPB.reset( );
120 0 : itsImagePBcor.reset();
121 0 : itsTempWorkIm.reset();
122 :
123 0 : itsSumWt.reset( );
124 0 : itsOverWrite=False;
125 0 : itsUseWeight=False;
126 0 : itsPBScaleFactor=1.0;
127 :
128 0 : itsNFacets=1;
129 0 : itsFacetId=0;
130 0 : itsNChanChunks = 1;
131 0 : itsChanId = 0;
132 0 : itsNPolChunks = 1;
133 0 : itsPolId = 0;
134 :
135 0 : itsImageShape=IPosition(4,0,0,0,0);
136 0 : itsImageName=String("");
137 0 : itsCoordSys=CoordinateSystem();
138 0 : itsMiscInfo=Record();
139 0 : init();
140 :
141 :
142 : // validate();
143 :
144 0 : }
145 :
146 : // Used from SynthesisNormalizer::makeImageStore()
147 0 : 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 0 : const Bool useweightimage)
155 : // TODO : Add parameter to indicate weight image shape.
156 : {
157 0 : LogIO os( LogOrigin("SIImageStore","Open new Images",WHERE) );
158 :
159 :
160 0 : itsPsf.reset( );
161 0 : itsModel.reset( );
162 0 : itsResidual.reset( );
163 0 : itsWeight.reset( );
164 0 : itsImage.reset( );
165 0 : itsMask.reset( );
166 0 : itsGridWt.reset( );
167 0 : itsPB.reset( );
168 0 : itsImagePBcor.reset( );
169 0 : itsTempWorkIm.reset();
170 :
171 0 : itsSumWt.reset( );
172 0 : itsOverWrite=False; // Hard Coding this. See CAS-6937. overwrite;
173 0 : itsUseWeight=useweightimage;
174 0 : itsPBScaleFactor=1.0;
175 :
176 0 : itsNFacets=1;
177 0 : itsFacetId=0;
178 0 : itsNChanChunks = 1;
179 0 : itsChanId = 0;
180 0 : itsNPolChunks = 1;
181 0 : itsPolId = 0;
182 :
183 0 : itsImageName = imagename;
184 0 : itsCoordSys = imcoordsys;
185 0 : itsImageShape = imshape;
186 0 : itsObjectName = objectname;
187 0 : itsMiscInfo = miscinfo;
188 :
189 0 : init();
190 :
191 0 : validate();
192 0 : }
193 :
194 : // Used from SynthesisNormalizer::makeImageStore()
195 0 : SIImageStore::SIImageStore(const String &imagename, const Bool ignorefacets, const Bool noRequireSumwt)
196 : {
197 0 : LogIO os( LogOrigin("SIImageStore","Open existing Images",WHERE) );
198 :
199 :
200 0 : itsPsf.reset( );
201 0 : itsModel.reset( );
202 0 : itsResidual.reset( );
203 0 : itsWeight.reset( );
204 0 : itsImage.reset( );
205 0 : itsMask.reset( );
206 0 : itsGridWt.reset( );
207 0 : itsPB.reset( );
208 0 : itsImagePBcor.reset( );
209 0 : itsTempWorkIm.reset();
210 0 : itsMiscInfo=Record();
211 :
212 0 : itsSumWt.reset( );
213 0 : itsNFacets=1;
214 0 : itsFacetId=0;
215 0 : itsNChanChunks = 1;
216 0 : itsChanId = 0;
217 0 : itsNPolChunks = 1;
218 0 : itsPolId = 0;
219 :
220 0 : itsOverWrite=False;
221 : //need to to this init now so that imageExts is initialized
222 0 : init();
223 :
224 0 : 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 0 : if( doesImageExist(itsImageName+String(".residual")) ||
230 0 : doesImageExist(itsImageName+String(".psf")) ||
231 0 : doesImageExist(itsImageName+String(".model")) ||
232 0 : doesImageExist(itsImageName+String(".gridwt")) ||
233 0 : doesImageExist(itsImageName+String(".pb")) ||
234 0 : doesImageExist(itsImageName+String(".weight"))
235 : )
236 : {
237 0 : std::shared_ptr<ImageInterface<Float> > imptr;
238 0 : if( doesImageExist(itsImageName+String(".psf")) )
239 : {
240 0 : buildImage( imptr, (itsImageName+String(".psf")) );
241 : // itsObjectName=imptr->imageInfo().objectName();
242 : // itsMiscInfo=imptr->miscInfo();
243 : }
244 0 : else if ( doesImageExist(itsImageName+String(".residual")) ){
245 0 : buildImage( imptr, (itsImageName+String(".residual")) );
246 : // itsObjectName=imptr->imageInfo().objectName();
247 : // itsMiscInfo=imptr->miscInfo();
248 : }
249 0 : else if ( doesImageExist(itsImageName+String(".model")) ){
250 0 : 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 0 : itsObjectName=imptr->imageInfo().objectName();
271 0 : itsImageShape=imptr->shape();
272 0 : itsCoordSys = imptr->coordinates();
273 0 : itsMiscInfo=imptr->miscInfo();
274 :
275 0 : }
276 : else
277 : {
278 0 : throw( AipsError( "PSF, Residual, Model Image (or sumwt) do not exist. Please create one of them." ) );
279 : }
280 :
281 0 : if( doesImageExist(itsImageName+String(".residual")) ||
282 0 : doesImageExist(itsImageName+String(".psf")) )
283 : {
284 0 : if( doesImageExist(itsImageName+String(".sumwt")) )
285 : {
286 0 : std::shared_ptr<ImageInterface<Float> > imptr;
287 : //imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt")) );
288 0 : buildImage( imptr, (itsImageName+String(".sumwt")) );
289 0 : itsNFacets = imptr->shape()[0];
290 0 : itsFacetId = 0;
291 0 : itsUseWeight = getUseWeightImage( *imptr );
292 0 : 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 0 : itsCoordSys = imptr->coordinates();
295 0 : itsMiscInfo=imptr->miscInfo();
296 0 : 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 0 : }
301 : else
302 : {
303 0 : 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 0 : 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 0 : if( doesImageExist(itsImageName+String(".residual") ) )
311 0 : { buildImage( imptr, (itsImageName+String(".residual")) ); }
312 : else
313 0 : { buildImage( imptr, (itsImageName+String(".psf")) ); }
314 :
315 0 : itsNFacets=1;
316 0 : itsFacetId=0;
317 0 : itsUseWeight=False;
318 0 : itsPBScaleFactor=1.0;
319 0 : itsCoordSys = imptr->coordinates();
320 0 : itsMiscInfo=imptr->miscInfo();
321 0 : }
322 : }
323 : }// if psf or residual exist...
324 :
325 0 : if( ignorefacets==True ) itsNFacets= 1;
326 :
327 0 : init();
328 :
329 0 : validate();
330 0 : }
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 0 : 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 0 : const Bool useweightimage)
355 : {
356 :
357 :
358 0 : itsPsf=psfim;
359 0 : itsModel=modelim;
360 0 : itsResidual=residim;
361 : /* if(residim){
362 : LatticeLocker lock1(*itsResidual, FileLocker::Read);
363 : cerr << "RESIDMAX " << max(itsResidual->get()) << " " << max(residim->get()) << endl;
364 : }*/
365 0 : itsWeight=weightim;
366 0 : itsImage=restoredim;
367 0 : itsMask=maskim;
368 :
369 0 : itsSumWt=sumwtim;
370 :
371 0 : itsGridWt=gridwtim;
372 0 : itsPB=pbim;
373 0 : itsImagePBcor=restoredpbcorim;
374 0 : itsTempWorkIm=tempworkimage;
375 :
376 0 : itsPBScaleFactor=1.0;// No need to set properly here as it will be computed in makePSF.
377 :
378 0 : itsObjectName = objectname;
379 0 : itsMiscInfo = miscinfo;
380 :
381 0 : itsNFacets = nfacets;
382 0 : itsFacetId = facet;
383 0 : itsNChanChunks = nchanchunks;
384 0 : itsChanId = chan;
385 0 : itsNPolChunks = npolchunks;
386 0 : itsPolId = pol;
387 :
388 0 : itsOverWrite=False;
389 0 : itsUseWeight=useweightimage;
390 :
391 0 : itsParentImageShape = imshape;
392 0 : itsImageShape = imshape;
393 :
394 0 : itsParentCoordSys = csys;
395 0 : itsCoordSys = csys; // Hopefully this doesn't change for a reference image
396 0 : itsImageName = imagename;
397 :
398 : //-----------------------
399 0 : init(); // Connect parent pointers to the images.
400 : //-----------------------
401 :
402 : // Set these to null, to be set later upon first access.
403 0 : itsPsf.reset( );
404 0 : itsModel.reset( );
405 0 : itsResidual.reset( );
406 0 : itsWeight.reset( );
407 0 : itsImage.reset( );
408 0 : itsMask.reset( );
409 0 : itsSumWt.reset( );
410 0 : itsPB.reset( );
411 :
412 0 : validate();
413 0 : }
414 :
415 0 : void SIImageStore::validate()
416 : {
417 : /// There are two valid states. Check for at least one of them.
418 0 : Bool state = False;
419 :
420 0 : stringstream oss;
421 0 : oss << "shape:" << itsImageShape << " parentimageshape:" << itsParentImageShape
422 0 : << " nfacets:" << itsNFacets << "x" << itsNFacets << " facetid:" << itsFacetId
423 0 : << " nchanchunks:" << itsNChanChunks << " chanid:" << itsChanId
424 0 : << " npolchunks:" << itsNPolChunks << " polid:" << itsPolId
425 0 : << " coord-dim:" << itsCoordSys.nPixelAxes()
426 0 : << " psf/res:" << (hasPsf() || hasResidual()) ;
427 0 : if( hasSumWt() ) oss << " sumwtshape : " << sumwt()->shape() ;
428 0 : oss << endl;
429 :
430 :
431 : try {
432 :
433 0 : if( itsCoordSys.nPixelAxes() != 4 ) state=False;
434 :
435 : /// (1) Regular imagestore
436 0 : if( itsNFacets==1 && itsFacetId==0
437 0 : && itsNChanChunks==1 && itsChanId==0
438 0 : && itsNPolChunks==1 && itsPolId==0 ) {
439 0 : Bool check1 = hasSumWt() && sumwt()->shape()[0]==1;
440 0 : if( (itsImageShape.isEqual(itsParentImageShape) ) && ( check1 || !hasSumWt() )
441 0 : && itsParentImageShape.product() > 0 ) state=True;
442 0 : }
443 : /// (2) Reference Sub Imagestore
444 0 : else if ( ( itsNFacets>1 && itsFacetId >=0 )
445 0 : || ( itsNChanChunks>1 && itsChanId >=0 )
446 0 : || ( itsNPolChunks>1 && itsPolId >=0 ) ) {
447 : // If shape is still unset, even when the first image has been made....
448 0 : Bool check1 = ( itsImageShape.product() > 0 && ( hasPsf() || hasResidual() ) );
449 0 : Bool check2 = ( itsImageShape.isEqual(IPosition(4,0,0,0,0)) && ( !hasPsf() && !hasResidual() ) );
450 0 : Bool check3 = hasSumWt() && sumwt()->shape()[0]==1; // One facet only.
451 :
452 0 : if( ( check1 || check2 ) && ( itsParentImageShape.product()>0 )
453 0 : && ( itsFacetId < itsNFacets*itsNFacets )
454 : // && ( itsChanId <= itsNChanChunks ) // chanchunks can be larger...
455 0 : && ( itsPolId < itsNPolChunks )
456 0 : && ( check3 || !hasSumWt() ) ) state = True;
457 : }
458 :
459 0 : } catch( AipsError &x ) {
460 0 : state = False;
461 0 : oss << " ||||| " << x.getMesg() << endl;
462 0 : }
463 :
464 : // cout << " SIIM:validate : " << oss.str() << endl;
465 :
466 0 : if( state == False ) throw(AipsError("Internal Error : Invalid ImageStore state : "+ oss.str()) );
467 :
468 0 : return;
469 0 : }
470 :
471 :
472 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
473 :
474 0 : 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 0 : 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 0 : 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 0 : std::shared_ptr<ImageInterface<Float> > imPtr;
492 0 : IPosition useShape( itsParentImageShape );
493 :
494 0 : if( dosumwt ) // change shape to sumwt image shape.
495 : {
496 0 : useShape[0] = nfacetsperside;
497 0 : 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 0 : Bool dbg=False;
517 0 : if( doesImageExist( imagenamefull ) )
518 : {
519 : ///// not used since itsOverWrite is hardcoded to FALSE (CAS-6937)
520 0 : if (overwrite) //overwrite and make new image
521 : {
522 0 : if(dbg) cout << "Trying to overwrite and open new image named : " << imagenamefull << " ow:"<< overwrite << endl;
523 : try{
524 0 : buildImage(imPtr, useShape, itsParentCoordSys, imagenamefull) ;
525 0 : LatticeLocker lock1(*imPtr, FileLocker::Write);
526 : // initialize to zeros...
527 0 : imPtr->set(0.0);
528 0 : }
529 0 : catch (AipsError &x){
530 0 : if(dbg)cout << "Cannot overwrite : " << x.getMesg() << endl;
531 0 : if(dbg)cout << "Open already ? : " << Table::isOpened( imagenamefull ) << " Writable ? : " << Table::isWritable( imagenamefull ) << endl;
532 0 : if(Table::isWritable( imagenamefull ))
533 : {
534 0 : if(dbg) cout << "--- Trying to open existing image : "<< imagenamefull << endl;
535 : try{
536 0 : buildImage( imPtr, imagenamefull );
537 : }
538 0 : catch (AipsError &x){
539 0 : throw( AipsError("Writable table exists, but cannot open : " + x.getMesg() ) );
540 0 : }
541 : }// is table writable
542 : else
543 : {
544 0 : throw( AipsError("Cannot overwrite existing image. : " + x.getMesg() ) );
545 : }
546 0 : }
547 : }// overwrite existing image
548 : else // open existing image ( Always tries this )
549 : {
550 : if(1) //Table::isWritable( imagenamefull ))
551 : {
552 0 : if(dbg) cout << "Trying to open existing image : "<< imagenamefull << endl;
553 : try{
554 0 : buildImage( imPtr, imagenamefull ) ;
555 :
556 0 : 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 0 : if( useShape.product()>0 && ! useShape.isEqual( imPtr->shape() ) )
562 : {
563 0 : ostringstream oo1,oo2;
564 0 : oo1 << useShape; oo2 << imPtr->shape();
565 0 : 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 0 : }
567 :
568 :
569 :
570 0 : 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 0 : if( itsParentCoordSys.errorMessage().contains("differing latpoles") )
575 : {
576 0 : DirectionCoordinate dcord = itsParentCoordSys.directionCoordinate();
577 : //cout << "Latpole from parent : " << dcord.longLatPoles() << endl;
578 :
579 0 : DirectionCoordinate dcord2 = imPtr->coordinates().directionCoordinate();
580 : //cout << "Latpole from imPtr : " << dcord2.longLatPoles() << endl;
581 :
582 0 : LogIO os( LogOrigin("SIImageStore","Open existing Images",WHERE) );
583 0 : os << " Mismatch in Csys latpoles between existing image on disk (" << dcord2.longLatPoles() << ") and current imaging run (" << dcord.longLatPoles() << ") : " << itsParentCoordSys.errorMessage() << " -- Resetting to match image on disk" << LogIO::WARN << LogIO::POST;
584 : // cout << "Set the coordinate from the Parent value" << endl;
585 : //imPtr->setCoordinateInfo( itsParentCoordSys );
586 0 : itsParentCoordSys = imPtr->coordinates();
587 0 : }
588 :
589 : // Check again.
590 0 : if ( ! itsParentCoordSys.near( imPtr->coordinates() ) )
591 : {
592 : //cout << " Still different..." << endl;
593 0 : 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 0 : catch (AipsError &x){
600 0 : throw( AipsError("Cannot open existing image : "+imagenamefull+" : " + x.getMesg() ) );
601 0 : }
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 0 : if(dbg) cout << "Trying to open new image named : " << imagenamefull << endl;
613 : try{
614 0 : buildImage(imPtr, useShape, itsParentCoordSys, imagenamefull) ;
615 : // initialize to zeros...
616 : {
617 0 : LatticeLocker lock1(*imPtr, FileLocker::Write);
618 0 : imPtr->set(0.0);
619 0 : imPtr->flush();
620 0 : imPtr->unlock();
621 0 : }
622 : }
623 0 : catch (AipsError &x){
624 0 : throw( AipsError("Cannot make new image. : " + x.getMesg() ) );
625 0 : }
626 : }
627 :
628 :
629 : //////////////////////////////////////
630 : /*
631 : if( overwrite || !Table::isWritable( imagenamefull ) )
632 : {
633 : cout << "Trying to open new image named : " << imagenamefull << " ow:"<< overwrite << endl;
634 : imPtr.reset( new PagedImage<Float> (useShape, itsCoordSys, imagenamefull) );
635 : // initialize to zeros...
636 : imPtr->set(0.0);
637 : }
638 : else
639 : {
640 : if(Table::isWritable( imagenamefull ))
641 : {
642 : cout << "Trying to open existing image : "<< imagenamefull << endl;
643 : try{
644 : imPtr.reset( new PagedImage<Float>( imagenamefull ) );
645 : }
646 : catch (AipsError &x){
647 : cerr << "Writable table exists, but cannot open. Creating temp image. : " << x.getMesg() << endl;
648 : imPtr.reset( new TempImage<Float> (useShape, itsCoordSys) );
649 : // imPtr=new PagedImage<Float> (useShape, itsCoordSys, imagenamefull);
650 : imPtr->set(0.0);
651 : }
652 : }
653 : else
654 : {
655 : cerr << "Table " << imagenamefull << " is not writeable. Creating temp image." << endl;
656 : imPtr.reset( new TempImage<Float> (useShape, itsCoordSys) );
657 : imPtr->set(0.0);
658 : }
659 : }
660 : */
661 0 : return imPtr;
662 0 : }
663 :
664 0 : void SIImageStore::buildImage(std::shared_ptr<ImageInterface<Float> > &imptr, IPosition shape,
665 : CoordinateSystem csys, const String name)
666 : {
667 0 : LogIO os( LogOrigin("SIImageStore", "Open non-existing image", WHERE) );
668 0 : os <<"Opening image, name: " << name << LogIO::DEBUG1;
669 :
670 0 : itsOpened++;
671 : //For some reason cannot open a new paged image with UserNoread directly
672 : {
673 0 : PagedImage<Float> godot(shape, csys, name);
674 0 : godot.unlock();
675 0 : }
676 0 : 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 0 : imptr.reset( new PagedImage<Float> (name, locktype) );
682 : {
683 0 : LatticeLocker lock1(*imptr, FileLocker::Write);
684 0 : initMetaInfo(imptr, name);
685 0 : imptr->unlock();
686 0 : }
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 0 : }
701 :
702 0 : void SIImageStore::buildImage(std::shared_ptr<ImageInterface<Float> > &imptr, const String name)
703 : {
704 0 : LogIO os(LogOrigin("SIImageStore", "Open existing Images", WHERE));
705 0 : os <<"Opening image, name: " << name << LogIO::DEBUG1;
706 :
707 0 : itsOpened++;
708 0 : if(Table::isReadable(name)){
709 0 : 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 0 : Table table(name, locktype);
714 0 : String type = table.tableInfo().type();
715 0 : if (type != TableInfo::type(TableInfo::PAGEDIMAGE)) {
716 :
717 0 : imptr.reset( new PagedImage<Float>( table ) );
718 0 : imptr->unlock();
719 0 : return;
720 : }
721 0 : }
722 0 : LatticeBase* latt =ImageOpener::openImage(name);
723 0 : if(!latt)
724 : {
725 0 : throw(AipsError("Error in opening Image : "+name));
726 : }
727 0 : DataType dtype=latt->dataType();
728 0 : if(dtype==TpFloat)
729 : {
730 0 : 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 0 : }
762 :
763 : /**
764 : * Sets ImageInfo and MiscInfo on an image
765 : *
766 : * @param imptr image to initialize
767 : */
768 0 : 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 0 : LatticeLocker lock1(*imptr, FileLocker::Write);
774 0 : if (not itsObjectName.empty()) {
775 0 : ImageInfo info = imptr->imageInfo();
776 0 : info.setObjectName(itsObjectName);
777 0 : imptr->setImageInfo(info);
778 0 : 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 0 : }
786 0 : imptr->unlock();
787 0 : }
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 0 : 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 0 : Int nx_facets=Int(sqrt(Double(nfacets)));
810 0 : LogIO os( LogOrigin("SynthesisImager","makeFacet") );
811 : // Make the output image
812 0 : Slicer imageSlicer;
813 :
814 : // Add checks for all dimensions..
815 0 : 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 0 : IPosition imshp=image.shape();
820 0 : IPosition blc(imshp.nelements(), 0);
821 0 : IPosition trc=imshp-1;
822 0 : IPosition inc(imshp.nelements(), 1);
823 :
824 : /// Facets
825 0 : Int facetx = facet % nx_facets;
826 0 : Int facety = (facet - facetx) / nx_facets;
827 0 : Int sizex = imshp(0) / nx_facets;
828 0 : Int sizey = imshp(1) / nx_facets;
829 0 : blc(1) = facety * sizey;
830 0 : trc(1) = blc(1) + sizey - 1;
831 0 : blc(0) = facetx * sizex;
832 0 : trc(0) = blc(0) + sizex - 1;
833 :
834 : /// Pol Chunks
835 0 : Int sizepol = imshp(2) / npolchunks;
836 0 : blc(2) = pol * sizepol;
837 0 : trc(2) = blc(2) + sizepol - 1;
838 :
839 : /// Chan Chunks
840 0 : Int sizechan = imshp(3) / nchanchunks;
841 0 : blc(3) = chan * sizechan;
842 0 : trc(3) = blc(3) + sizechan - 1;
843 :
844 0 : LCBox::verify(blc, trc, inc, imshp);
845 0 : Slicer imslice(blc, trc, inc, Slicer::endIsLast);
846 :
847 : // Now create the sub image
848 0 : std::shared_ptr<ImageInterface<Float> > referenceImage( new SubImage<Float>(image, imslice, True) );
849 : {
850 0 : LatticeLocker lock1(*referenceImage, FileLocker::Write);
851 0 : referenceImage->setMiscInfo(image.miscInfo());
852 0 : referenceImage->setUnits(image.units());
853 :
854 0 : }
855 :
856 : // cout << "Made Ref subimage of shape : " << referenceImage->shape() << endl;
857 :
858 0 : return referenceImage;
859 :
860 0 : }
861 :
862 :
863 :
864 : //////////////////////////////////////////////////////////////////////////////////////////////////////
865 :
866 0 : SIImageStore::~SIImageStore()
867 : {
868 0 : }
869 :
870 : //////////////////////////////////////////////////////////////////////////////////////////////////////
871 : //////////////////////////////////////////////////////////////////////////////////////////////////////
872 :
873 0 : Bool SIImageStore::releaseLocks()
874 : {
875 : //LogIO os( LogOrigin("SIImageStore","releaseLocks",WHERE) );
876 :
877 : // String fname( itsImageName+String(".info") );
878 : // makePersistent( fname );
879 :
880 0 : if( itsPsf ) releaseImage( itsPsf );
881 0 : if( itsModel ) { releaseImage( itsModel ); }
882 0 : if( itsResidual ) releaseImage( itsResidual );
883 0 : if( itsImage ) releaseImage( itsImage );
884 0 : if( itsWeight ) releaseImage( itsWeight );
885 0 : if( itsMask ) releaseImage( itsMask );
886 0 : if( itsSumWt ) releaseImage( itsSumWt );
887 0 : if( itsGridWt ) releaseImage( itsGridWt );
888 0 : if( itsPB ) releaseImage( itsPB );
889 0 : if( itsImagePBcor ) releaseImage( itsImagePBcor );
890 :
891 0 : 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 0 : 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 0 : im->flush();
909 : //os << LogIO::WARN << "clear cache" << LogIO::POST;
910 0 : im->clearCache();
911 : //os << LogIO::WARN << "unlock" << LogIO::POST;
912 0 : im->unlock();
913 : //os << LogIO::WARN << "tempClose" << LogIO::POST;
914 : //im->tempClose();
915 : //os << LogIO::WARN << "NULL" << LogIO::POST;
916 0 : im.reset(); // This was added to allow modification by modules independently
917 0 : }
918 :
919 0 : void SIImageStore::releaseImage( std::shared_ptr<ImageInterface<Complex> > &im )
920 : {
921 0 : im->tempClose();
922 0 : }
923 :
924 : //////////////////////////////////////////////////////////////////////////////////////////////////////
925 : //////////////////////////////////////////////////////////////////////////////////////////////////////
926 0 : Vector<String> SIImageStore::getModelImageName()
927 : {
928 0 : Vector<String> mods(1);
929 0 : mods[0]=itsImageName + imageExts(MODEL);
930 0 : return mods;
931 0 : }
932 :
933 0 : Long SIImageStore::estimateRAM(){
934 : //right now this is estimated at 2MB for the 2 complex lattices;
935 0 : return Long(2000);
936 : }
937 0 : void SIImageStore::setModelImage( const Vector<String> &modelnames)
938 : {
939 0 : LogIO os( LogOrigin("SIImageStore","setModelImage",WHERE) );
940 :
941 0 : 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 0 : if( modelnames.nelements()==1 ) { setModelImageOne( modelnames[0] ); }
947 0 : }
948 :
949 :
950 :
951 0 : void SIImageStore::setModelImageOne( const String &modelname , Int nterm)
952 : {
953 0 : LogIO os( LogOrigin("SIImageStore","setModelImageOne",WHERE) );
954 :
955 0 : if(modelname==String("")) return;
956 :
957 0 : Bool multiterm=False;
958 0 : if(nterm>-1)multiterm=True;
959 0 : if(nterm==-1)nterm=0;
960 :
961 0 : Directory immodel( modelname ); // +String(".model") );
962 0 : 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 0 : std::shared_ptr<ImageInterface<Float> > newmodel;
972 0 : buildImage(newmodel, modelname);
973 : // in master
974 : //std::shared_ptr<PagedImage<Float> > newmodel( new PagedImage<Float>( modelname ) ); //+String(".model") ) );
975 :
976 0 : Bool hasMask = newmodel->isMasked(); /// || newmodel->hasPixelMask() ;
977 :
978 0 : if( hasMask )
979 : {
980 :
981 0 : 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 0 : TempImage<Float> maskmodel( newmodel->shape(), newmodel->coordinates() );
993 0 : IPosition inshape = newmodel->shape();
994 0 : for(Int pol=0; pol<inshape[2]; pol++)
995 : {
996 0 : for(Int chan=0; chan<inshape[3]; chan++)
997 : {
998 0 : IPosition pos(4,0,0,pol,chan);
999 : std::shared_ptr<ImageInterface<Float> > subim=makeSubImage(0,1,
1000 0 : chan, inshape[3],
1001 0 : pol, inshape[2],
1002 0 : (*newmodel) );
1003 :
1004 : std::shared_ptr<ImageInterface<Float> > submaskmodel=makeSubImage(0,1,
1005 0 : chan, inshape[3],
1006 0 : pol, inshape[2],
1007 0 : maskmodel );
1008 :
1009 0 : ArrayLattice<Bool> pixmask( subim->getMask() );
1010 0 : LatticeExpr<Float> masked( (*subim) * iif(pixmask,1.0,0.0) );
1011 0 : submaskmodel->copyData( masked );
1012 0 : }
1013 : }
1014 :
1015 :
1016 :
1017 : // Check shapes, coordsys with those of other images. If different, try to re-grid here.
1018 0 : 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 0 : os << "Copying input model " << modelname << " to " << itsImageName << ".model" << ((multiterm)?".tt"+String::toString(nterm) :"") << LogIO::POST;
1026 0 : 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 0 : }
1033 :
1034 : }// hasMask
1035 : else // nomask
1036 : {
1037 :
1038 : // Check shapes, coordsys with those of other images. If different, try to re-grid here.
1039 0 : if( (newmodel->shape() != model(nterm)->shape()) || (! itsCoordSys.near(newmodel->coordinates() )) )
1040 : {
1041 0 : os << "Regridding input model " << modelname << " to target coordinate system for " << itsImageName << ".model" << ((multiterm)?".tt"+String::toString(nterm) :"") << LogIO::POST;
1042 0 : regridToModelImage( *newmodel, nterm );
1043 : }
1044 : else
1045 : {
1046 0 : os << "Copying input model " << modelname << " to " << itsImageName << ".model" << ((multiterm)?".tt"+String::toString(nterm) :"") << LogIO::POST;
1047 0 : model(nterm)->copyData( LatticeExpr<Float> (*newmodel) );
1048 : }
1049 : }//nomask
1050 0 : }
1051 :
1052 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1053 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1054 :
1055 0 : IPosition SIImageStore::getShape()
1056 : {
1057 0 : return itsImageShape;
1058 : }
1059 :
1060 0 : String SIImageStore::getName()
1061 : {
1062 0 : return itsImageName;
1063 : }
1064 :
1065 0 : uInt SIImageStore::getNTaylorTerms(Bool /*dopsf*/)
1066 : {
1067 0 : 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 0 : 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 0 : Bool sw=False;
1087 0 : if( label.contains(imageExts(SUMWT)) ) sw=True;
1088 :
1089 0 : if( !ptr )
1090 : {
1091 : //cout << itsNFacets << " " << itsNChanChunks << " " << itsNPolChunks << endl;
1092 0 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
1093 : {
1094 0 : if( ! parentptr )
1095 : {
1096 : //cout << "Making parent : " << itsImageName+label << " sw : " << sw << endl;
1097 0 : parentptr = openImage(itsImageName+label , itsOverWrite, sw, itsNFacets );
1098 0 : 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 0 : ptr = makeSubImage( itsFacetId, itsNFacets*itsNFacets,
1104 : itsChanId, itsNChanChunks,
1105 : itsPolId, itsNPolChunks,
1106 0 : *parentptr );
1107 0 : if( ! sw )
1108 : {
1109 0 : itsImageShape = ptr->shape(); // over and over again.... FIX.
1110 0 : itsCoordSys = ptr->coordinates();
1111 0 : itsMiscInfo = ptr->miscInfo();
1112 : }
1113 :
1114 : //cout << "accessImage : " << label << " : sumwt : " << sw << " : shape : " << itsImageShape << endl;
1115 :
1116 : }
1117 : else
1118 : { //no facets of chanchunks
1119 0 : if(!parentptr){
1120 : ///coordsys for psf can be different ...shape should be the same.
1121 0 : ptr = openImage(itsImageName+label , itsOverWrite, sw, 1, !(label.contains(imageExts(PSF)))); }
1122 : else{
1123 0 : 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 0 : }
1136 :
1137 :
1138 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::psf(uInt /*nterm*/)
1139 : {
1140 0 : accessImage( itsPsf, itsParentPsf, imageExts(PSF) );
1141 :
1142 0 : return itsPsf;
1143 : }
1144 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::residual(uInt /*nterm*/)
1145 : {
1146 0 : accessImage( itsResidual, itsParentResidual, imageExts(RESIDUAL) );
1147 : // Record mi = itsResidual->miscInfo(); ostringstream oss;mi.print(oss);cout<<"MiscInfo(res) : " << oss.str() << endl;
1148 0 : return itsResidual;
1149 : }
1150 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::weight(uInt /*nterm*/)
1151 : {
1152 0 : accessImage( itsWeight, itsParentWeight, imageExts(WEIGHT) );
1153 0 : return itsWeight;
1154 : }
1155 :
1156 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::sumwt(uInt /*nterm*/)
1157 : {
1158 :
1159 0 : accessImage( itsSumWt, itsParentSumWt, imageExts(SUMWT) );
1160 :
1161 0 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
1162 0 : { itsUseWeight = getUseWeightImage( *itsParentSumWt );}
1163 0 : setUseWeightImage( *itsSumWt , itsUseWeight); // Sets a flag in the SumWt image.
1164 :
1165 0 : return itsSumWt;
1166 : }
1167 :
1168 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::model(uInt /*nterm*/)
1169 : {
1170 0 : accessImage( itsModel, itsParentModel, imageExts(MODEL) );
1171 0 : LatticeLocker lock1(*itsModel, FileLocker::Write);
1172 : // Set up header info the first time.
1173 0 : itsModel->setUnits("Jy/pixel");
1174 :
1175 0 : return itsModel;
1176 0 : }
1177 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::image(uInt /*nterm*/)
1178 : {
1179 0 : accessImage( itsImage, itsParentImage, imageExts(IMAGE) );
1180 0 : LatticeLocker lock1(*itsImage, FileLocker::Write);
1181 0 : itsImage->setUnits(Unit("Jy/beam"));
1182 0 : return itsImage;
1183 0 : }
1184 :
1185 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::mask(uInt /*nterm*/)
1186 : {
1187 0 : accessImage( itsMask, itsParentMask, imageExts(MASK) );
1188 0 : return itsMask;
1189 : }
1190 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::gridwt(IPosition newshape)
1191 :
1192 : {
1193 0 : if(newshape.empty()){
1194 0 : accessImage( itsGridWt, itsParentGridWt, imageExts(GRIDWT) );
1195 : }
1196 : else{
1197 0 : if(!itsGridWt || (itsGridWt && (itsGridWt->shape() != newshape))){
1198 0 : itsGridWt.reset(); // deassign previous one hopefully it'll close it
1199 0 : CoordinateSystem newcoordsys=itsCoordSys;
1200 0 : if(newshape.nelements() > 4){
1201 0 : Matrix<Double> pc(1,1); pc = 1.0;
1202 0 : LinearCoordinate zc(Vector<String>(1, "FIELD_ORDER"), Vector<String>(1, ""), Vector<Double>(1, 0.0), Vector<Double>(1,1.0), pc, Vector<Double>(1,0.0));
1203 0 : newcoordsys.addCoordinate(zc);
1204 0 : }
1205 0 : itsGridWt.reset(new PagedImage<Float>(newshape, newcoordsys, itsImageName+ imageExts(GRIDWT)));
1206 0 : initMetaInfo(itsGridWt, itsImageName+ imageExts(GRIDWT));
1207 :
1208 0 : }
1209 : }
1210 : /// change the coordinate system here, to uv.
1211 0 : return itsGridWt;
1212 : }
1213 :
1214 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::tempworkimage(uInt /*term*/){
1215 0 : if(itsTempWorkIm) return itsTempWorkIm;
1216 0 : itsTempWorkIm.reset(new PagedImage<Float>(itsImageShape, itsCoordSys, itsImageName+ ".work.temp"));
1217 0 : static_cast<PagedImage<Float>* > (itsTempWorkIm.get())->set(0.0);
1218 0 : itsTempWorkIm->flush();
1219 0 : static_cast<PagedImage<Float>* > (itsTempWorkIm.get())->table().markForDelete();
1220 0 : return itsTempWorkIm;
1221 : }
1222 :
1223 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::pb(uInt /*nterm*/)
1224 : {
1225 0 : accessImage( itsPB, itsParentPB, imageExts(PB) );
1226 0 : return itsPB;
1227 : }
1228 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::imagepbcor(uInt /*nterm*/)
1229 : {
1230 0 : accessImage( itsImagePBcor, itsParentImagePBcor, imageExts(IMAGEPBCOR) );
1231 0 : LatticeLocker lock1(*itsImagePBcor, FileLocker::Write);
1232 0 : itsImagePBcor->setUnits(Unit("Jy/beam"));
1233 0 : return itsImagePBcor;
1234 0 : }
1235 :
1236 0 : std::shared_ptr<ImageInterface<Complex> > SIImageStore::forwardGrid(uInt /*nterm*/){
1237 0 : if( itsForwardGrid ) // && (itsForwardGrid->shape() == itsImageShape))
1238 : {
1239 : // cout << "Forward grid has shape : " << itsForwardGrid->shape() << endl;
1240 0 : return itsForwardGrid;
1241 : }
1242 0 : Vector<Int> whichStokes(0);
1243 0 : IPosition cimageShape;
1244 0 : cimageShape=itsImageShape;
1245 0 : 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 0 : if(freqframe != MFrequency::LSRK && freqframe!=MFrequency::Undefined && freqframe!=MFrequency::REST)
1248 0 : { itsCoordSys.setSpectralConversion("LSRK"); }
1249 0 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
1250 0 : whichStokes, itsDataPolRep);
1251 0 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
1252 0 : cimageShape(2)=whichStokes.nelements();
1253 : //cout << "Making forward grid of shape : " << cimageShape << " for imshape : " << itsImageShape << endl;
1254 0 : itsForwardGrid.reset( new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()) );
1255 : //if(image())
1256 0 : if(hasRestored()){
1257 0 : LatticeLocker lock1(*itsForwardGrid, FileLocker::Write);
1258 0 : itsForwardGrid->setImageInfo((image())->imageInfo());
1259 :
1260 0 : }
1261 0 : return itsForwardGrid;
1262 0 : }
1263 :
1264 0 : std::shared_ptr<ImageInterface<Complex> > SIImageStore::backwardGrid(uInt /*nterm*/){
1265 0 : if( itsBackwardGrid ) //&& (itsBackwardGrid->shape() == itsImageShape))
1266 : {
1267 : // cout << "Backward grid has shape : " << itsBackwardGrid->shape() << endl;
1268 0 : return itsBackwardGrid;
1269 : }
1270 0 : Vector<Int> whichStokes(0);
1271 0 : IPosition cimageShape;
1272 0 : cimageShape=itsImageShape;
1273 0 : 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 0 : if(freqframe != MFrequency::LSRK && freqframe!=MFrequency::Undefined && freqframe!=MFrequency::REST )
1276 0 : { itsCoordSys.setSpectralConversion("LSRK"); }
1277 0 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
1278 0 : whichStokes, itsDataPolRep);
1279 0 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
1280 0 : cimageShape(2)=whichStokes.nelements();
1281 : //cout << "Making backward grid of shape : " << cimageShape << " for imshape : " << itsImageShape << endl;
1282 0 : itsBackwardGrid.reset( new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()) );
1283 : //if(image())
1284 0 : if(hasRestored()){
1285 0 : LatticeLocker lock1(*itsBackwardGrid, FileLocker::Write);
1286 0 : itsBackwardGrid->setImageInfo((image())->imageInfo());
1287 0 : }
1288 0 : return itsBackwardGrid;
1289 0 : }
1290 0 : Double SIImageStore::memoryBeforeLattice(){
1291 : //Calculate how much memory to use per temporary images before disking
1292 0 : return 1.0; /// in MB
1293 : }
1294 0 : IPosition SIImageStore::tileShape(){
1295 : //Need to have settable stuff here or algorith to determine this
1296 0 : 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 0 : Bool SIImageStore::doesImageExist(String imagename)
1301 : {
1302 0 : LogIO os( LogOrigin("SIImageStore","doesImageExist",WHERE) );
1303 0 : Directory image( imagename );
1304 0 : return image.exists();
1305 0 : }
1306 :
1307 :
1308 0 : void SIImageStore::resetImages( Bool resetpsf, Bool resetresidual, Bool resetweight )
1309 : {
1310 0 : if( resetpsf ){
1311 0 : LatticeLocker lock1(*(psf()), FileLocker::Write);
1312 0 : psf()->set(0.0);
1313 :
1314 0 : }
1315 0 : if( resetresidual ) {
1316 : // removeMask( residual() );
1317 0 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1318 0 : residual()->set(0.0);
1319 0 : }
1320 0 : if( resetweight && itsWeight ){
1321 0 : LatticeLocker lock1(*(weight()), FileLocker::Write);
1322 0 : weight()->set(0.0);
1323 0 : }
1324 0 : if( resetweight ){
1325 0 : LatticeLocker lock1(*(sumwt()), FileLocker::Write);
1326 0 : sumwt()->set(0.0);
1327 0 : }
1328 0 : }
1329 :
1330 0 : void SIImageStore::addImages( std::shared_ptr<SIImageStore> imagestoadd,
1331 : Bool addpsf, Bool addresidual, Bool addweight, Bool adddensity)
1332 : {
1333 :
1334 0 : if(addpsf)
1335 : {
1336 0 : LatticeExpr<Float> adderPsf( *(psf()) + *(imagestoadd->psf()) );
1337 0 : psf()->copyData(adderPsf);
1338 0 : }
1339 0 : if(addresidual)
1340 : {
1341 0 : LatticeExpr<Float> adderRes( *(residual()) + *(imagestoadd->residual()) );
1342 0 : residual()->copyData(adderRes);
1343 0 : }
1344 0 : if(addweight)
1345 : {
1346 0 : if( getUseWeightImage( *(imagestoadd->sumwt()) ) ) // Access and add weight only if it is needed.
1347 : {
1348 0 : LatticeExpr<Float> adderWeight( *(weight()) + *(imagestoadd->weight()) );
1349 0 : weight()->copyData(adderWeight);
1350 0 : }
1351 :
1352 : /*
1353 : Array<Float> qqq, www;
1354 : imagestoadd->sumwt()->get(qqq,True);
1355 : sumwt()->get(www,True);
1356 : cout << "SUMWT : Adding : " << qqq << " to " << www << endl;
1357 : */
1358 :
1359 0 : LatticeExpr<Float> adderSumWt( *(sumwt()) + *(imagestoadd->sumwt()) );
1360 0 : sumwt()->copyData(adderSumWt);
1361 0 : setUseWeightImage( *sumwt(), getUseWeightImage(*(imagestoadd->sumwt()) ) );
1362 :
1363 0 : }
1364 0 : if(adddensity)
1365 : {
1366 0 : LatticeExpr<Float> adderDensity( *(gridwt()) + *(imagestoadd->gridwt()) );
1367 0 : gridwt()->copyData(adderDensity);
1368 0 : }
1369 :
1370 0 : }
1371 :
1372 0 : void SIImageStore::setWeightDensity( std::shared_ptr<SIImageStore> imagetoset )
1373 : {
1374 0 : LogIO os( LogOrigin("SIImageStore","setWeightDensity",WHERE) );
1375 : //gridwt()->copyData( LatticeExpr<Float> ( *(imagetoset->gridwt()) ) );
1376 : //looks like you have to call them before hand or it crashes in parallel sometimes
1377 0 : IPosition shape=(imagetoset->gridwt())->shape();
1378 0 : os << LogIO::DEBUG2 << "SHAPES " << imagetoset->gridwt()->shape() << " " << gridwt()->shape() << endl;
1379 0 : (gridwt(shape))->copyData( LatticeExpr<Float> ( *(imagetoset->gridwt()) ) );
1380 :
1381 0 : }
1382 :
1383 : // TODO
1384 : // cout << "WARN : get getPbMax right for cube.... weight is indexed on chan and pol." << endl;
1385 0 : 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 0 : LatticeExprNode le( sqrt(max( *(weight(0)) )) );
1401 0 : return le.getFloat();
1402 :
1403 0 : }
1404 :
1405 0 : 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 0 : CountedPtr<ImageInterface<Float> > subim=makeSubImage(0, 1,
1411 0 : chan, itsImageShape[3],
1412 0 : pol, itsImageShape[2],
1413 0 : *weight(0));
1414 :
1415 0 : return sqrt(max(subim->get()));
1416 0 : }
1417 :
1418 0 : void SIImageStore::makePBFromWeight(const Float pblimit)
1419 : {
1420 0 : LogIO os( LogOrigin("SIImageStore","makePBFromWeight",WHERE) );
1421 :
1422 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1423 : {
1424 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1425 : {
1426 :
1427 0 : itsPBScaleFactor = getPbMax(pol,chan);
1428 :
1429 0 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1430 : else {
1431 :
1432 0 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1433 0 : chan, itsImageShape[3],
1434 0 : pol, itsImageShape[2],
1435 0 : *weight() );
1436 0 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
1437 0 : chan, itsImageShape[3],
1438 0 : pol, itsImageShape[2],
1439 0 : *pb() );
1440 :
1441 :
1442 0 : LatticeExpr<Float> normed( sqrt(abs(*wtsubim)) / itsPBScaleFactor );
1443 0 : LatticeExpr<Float> limited( iif( normed > fabs(pblimit) , normed, 0.0 ) );
1444 0 : pbsubim->copyData( limited );
1445 0 : }// if not zero
1446 : }
1447 : }
1448 :
1449 0 : 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 0 : LatticeExpr<Bool> pbmask( iif( *pb() > fabs(pblimit) , True , False ) );
1458 : //MSK//
1459 0 : createMask( pbmask, pb() );
1460 0 : pb()->pixelMask().unlock();
1461 0 : }
1462 :
1463 : }
1464 0 : weight()->unlock();
1465 0 : pb()->unlock();
1466 0 : }
1467 :
1468 0 : void SIImageStore::makePBImage(const Float pblimit)
1469 : {
1470 0 : LogIO os( LogOrigin("SIImageStore","makePBImage",WHERE) );
1471 :
1472 0 : if( hasPB() )
1473 : {
1474 :
1475 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1476 : {
1477 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1478 : {
1479 :
1480 0 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
1481 0 : chan, itsImageShape[3],
1482 0 : pol, itsImageShape[2],
1483 0 : *pb() );
1484 :
1485 : // Norm by the max
1486 0 : LatticeExprNode elmax= max( *pbsubim );
1487 0 : 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 0 : if(fmax>1.0)
1491 : {
1492 0 : LatticeExpr<Float> normed( (*pbsubim) / fmax );
1493 0 : LatticeExpr<Float> limited( iif( normed > fabs(pblimit) , normed, 0.0 ) );
1494 0 : pbsubim->copyData( limited );
1495 0 : }
1496 : else
1497 : {
1498 0 : LatticeExpr<Float> limited( iif((*pbsubim) > fabs(pblimit) , (*pbsubim), 0.0 ) );
1499 0 : pbsubim->copyData( limited );
1500 0 : }
1501 0 : }
1502 : }
1503 :
1504 0 : if((pb()->getDefaultMask()==""))// && pblimit >= 0.0)
1505 : {
1506 : // removeMask( pb() );
1507 : //MSK//
1508 0 : LatticeExpr<Bool> pbmask( iif( *pb() > fabs(pblimit) , True , False ) );
1509 : //MSK//
1510 0 : createMask( pbmask, pb() );
1511 0 : pb()->pixelMask().unlock();
1512 0 : }
1513 :
1514 : }// if hasPB
1515 0 : pb()->unlock();
1516 :
1517 0 : }
1518 :
1519 0 : 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 0 : LatticeLocker lock1(*outimage, FileLocker::Write);
1526 0 : ImageRegion outreg = outimage->makeMask("mask0",False,True);
1527 0 : LCRegion& outmask=outreg.asMask();
1528 0 : outmask.copyData(lemask);
1529 0 : outimage->defineRegion("mask0",outreg, RegionHandler::Masks, True);
1530 :
1531 0 : outimage->setDefaultMask("mask0");
1532 :
1533 0 : outimage->unlock();
1534 0 : outimage->tempClose();
1535 :
1536 : // outimage->table().unmarkForDelete();
1537 0 : }
1538 0 : catch (const AipsError& x) {
1539 0 : throw(AipsError("Error in creating internal T/F mask : " + x.getMesg() ));
1540 0 : }
1541 :
1542 0 : return True;
1543 : }
1544 :
1545 0 : 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 0 : if( (inimage->getDefaultMask()).matches("mask0") ) // input mask exists.
1553 0 : {LatticeLocker lock1(*outimage, FileLocker::Write);
1554 0 : 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 0 : ImageRegion outreg=outimage->makeMask("mask0", False, True);
1563 0 : LCRegion& outmask=outreg.asMask();
1564 0 : outmask.copyData(inimage->getRegion("mask0").asLCRegion());
1565 0 : outimage->defineRegion("mask0",outreg, RegionHandler::Masks,True);
1566 0 : outimage->setDefaultMask("mask0");
1567 0 : }
1568 : }
1569 0 : catch (const AipsError& x) {
1570 0 : throw(AipsError("Error in copying internal T/F mask : " + x.getMesg() ));
1571 0 : }
1572 :
1573 0 : return True;
1574 : }
1575 :
1576 0 : 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 0 : LatticeLocker lock1(*im, FileLocker::Write);
1586 0 : if (im-> getDefaultMask() != String(""))
1587 : {
1588 0 : String strung=im->getDefaultMask();
1589 0 : im->setDefaultMask("");
1590 0 : im->removeRegion(strung);
1591 0 : }
1592 0 : if( im->hasRegion("mask0") )
1593 : {
1594 0 : im->removeRegion("mask0");
1595 : }
1596 0 : }
1597 0 : catch (const AipsError& x)
1598 : {
1599 0 : throw(AipsError("Error in deleting internal T/F mask : " + x.getMesg() ));
1600 0 : }
1601 0 : }
1602 0 : void SIImageStore:: rescaleResolution(Int chan,
1603 : ImageInterface<Float>& image,
1604 : const GaussianBeam& newbeam,
1605 : const GaussianBeam& oldbeam){
1606 :
1607 0 : LogIO os( LogOrigin("SIImageStore","rescaleResolution",WHERE) );
1608 0 : GaussianBeam toBeUsed(Quantity(0.0, "arcsec"),Quantity(0.0, "arcsec"),
1609 0 : Quantity(0.0, "deg")) ;
1610 : try {
1611 0 : 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 0 : if( samesize )
1621 : {
1622 0 : os << LogIO::NORMAL2 << "Input and output beam sizes are the same for Channel : " << chan << ". Not convolving residuals." << LogIO::POST;
1623 : }
1624 : else
1625 : {
1626 0 : Double pixwidth=sqrt(image.coordinates().increment()(0)*image.coordinates().increment()(0)+image.coordinates().increment()(1)*image.coordinates().increment()(1));
1627 :
1628 0 : if(toBeUsed.getMinor(image.coordinates().worldAxisUnits()[0]) > pixwidth)
1629 : {
1630 : //cerr << "old beam area " << oldbeam.getArea("rad2") << " new beam " << newbeam.getArea("rad2") << endl;
1631 0 : StokesImageUtil::Convolve(image, toBeUsed, True);
1632 0 : image.copyData(LatticeExpr<Float>(image*newbeam.getArea("rad2")/ oldbeam.getArea("rad2")));
1633 : }
1634 : }
1635 : }
1636 0 : catch (const AipsError& x) {
1637 : //throw(AipsError("Cannot convolve to new beam: may be smaller than old beam : " + x.getMesg() ));
1638 0 : os << LogIO::WARN << "Cannot convolve to new beam for Channel : " << chan << " : " << x.getMesg() << LogIO::POST;
1639 0 : }
1640 :
1641 :
1642 0 : }
1643 :
1644 :
1645 :
1646 :
1647 0 : void SIImageStore::dividePSFByWeight(const Float /* pblimit*/)
1648 : {
1649 0 : LogIO os( LogOrigin("SIImageStore","dividePSFByWeight",WHERE) );
1650 :
1651 0 : LatticeLocker lock1 (*(psf()), FileLocker::Write);
1652 0 : normPSF();
1653 :
1654 0 : if( itsUseWeight )
1655 : {
1656 :
1657 0 : divideImageByWeightVal( *weight() );
1658 : }
1659 0 : (psf())->unlock();
1660 :
1661 0 : }
1662 :
1663 0 : void SIImageStore::normalizePrimaryBeam(const Float pblimit)
1664 : {
1665 0 : LogIO os( LogOrigin("SIImageStore","normalizePrimaryBeam",WHERE) );
1666 :
1667 : // cout << "In dividePSFByWeight : itsUseWeight : " << itsUseWeight << endl;
1668 0 : if( itsUseWeight )
1669 : {
1670 :
1671 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1672 : {
1673 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1674 : {
1675 0 : 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 0 : makePBFromWeight(pblimit);
1680 :
1681 : }//if itsUseWeight
1682 0 : else { makePBImage(pblimit); } // OR... just check that it exists already.
1683 0 : pb()->unlock();
1684 :
1685 0 : }
1686 :
1687 : // Make another for the PSF too.
1688 0 : void SIImageStore::divideResidualByWeight(Float pblimit, String normtype) {
1689 :
1690 0 : LogIO os(LogOrigin("SIImageStore", "divideResidualByWeight", WHERE));
1691 0 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1692 :
1693 0 : auto logTemplate = [&](Int const &chan, Int const &pol, string const &normalizer, string const &result) {
1694 0 : os << LogIO::NORMAL1
1695 0 : << "[C" + String::toString(chan) + ":P" + String::toString(pol) + "] "
1696 0 : << "Dividing " << itsImageName + String(".residual") << " by "
1697 : << "[ " << normalizer << " ] "
1698 0 : << "to get " << result << "." << LogIO::POST;
1699 0 : };
1700 :
1701 : // Normalize by the sumwt, per plane.
1702 0 : Bool didNorm = divideImageByWeightVal(*residual());
1703 :
1704 0 : if (itsUseWeight) {
1705 0 : for (Int pol = 0; pol < itsImageShape[2]; pol++) {
1706 0 : for (Int chan = 0; chan < itsImageShape[3]; chan++) {
1707 0 : itsPBScaleFactor = getPbMax(pol, chan);
1708 :
1709 0 : if (itsPBScaleFactor <= 0) {
1710 : os << LogIO::NORMAL1
1711 : << "Skipping normalization for C:" << chan << " P:" << pol
1712 0 : << " because pb max is zero " << LogIO::POST;
1713 :
1714 : } else {
1715 0 : CountedPtr<ImageInterface<Float> > wtsubim = makeSubImage(0, 1,
1716 0 : chan, itsImageShape[3],
1717 0 : pol, itsImageShape[2],
1718 0 : *weight());
1719 0 : CountedPtr<ImageInterface<Float> > ressubim = makeSubImage(0, 1,
1720 0 : chan, itsImageShape[3],
1721 0 : pol, itsImageShape[2],
1722 0 : *residual());
1723 0 : LatticeExpr<Float> ratio;
1724 0 : Float scalepb = 1.0;
1725 :
1726 0 : if (normtype == "flatnoise") {
1727 0 : logTemplate(chan, pol,
1728 0 : "sqrt(weightimage) * " + String::toString(itsPBScaleFactor),
1729 : "flat noise with unit pb peak");
1730 :
1731 0 : LatticeExpr<Float> deno = itsPBScaleFactor * sqrt(abs(LatticeExpr<Float>(*(wtsubim))));
1732 0 : scalepb = fabs(pblimit) * itsPBScaleFactor * itsPBScaleFactor;
1733 0 : 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 0 : }
1751 :
1752 : //IPosition ip(4, itsImageShape[0] / 2, itsImageShape[1]/2, 0, 0);
1753 : //Float resval = ressubim->getAt(ip);
1754 : //LatticeExpr<Float> mask(iif((deno) > scalepb, 1.0, 0.0));
1755 : //LatticeExpr<Float> maskinv(iif((deno) > scalepb, 0.0, 1.0));
1756 : //LatticeExpr<Float> ratio(((*(ressubim)) * mask) / (deno + maskinv));
1757 :
1758 : //above blocks all sources outside minpb but visible with weight coverage
1759 : //which could be cleaned out...one could use below for that
1760 : //LatticeExpr<Float> ratio(iif(deno > scalepb, (*(ressubim)) / deno, *ressubim));
1761 :
1762 0 : 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 0 : } // 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 0 : 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 0 : if ((residual()->getDefaultMask() == "") && hasPB() && pblimit >=0.0) {
1779 0 : copyMask(pb(), residual());
1780 : }
1781 :
1782 0 : if ((pblimit < 0.0) && (residual()->getDefaultMask()).matches("mask0")) {
1783 0 : removeMask(residual());
1784 : }
1785 :
1786 0 : residual()->unlock();
1787 0 : }
1788 :
1789 0 : void SIImageStore::divideResidualByWeightSD(Float pblimit) {
1790 :
1791 0 : LogIO os(LogOrigin("SIImageStore", "divideResidualByWeightSD", WHERE));
1792 0 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1793 :
1794 0 : if (itsUseWeight) {
1795 0 : LatticeExpr<Float> deno = LatticeExpr<Float>(*weight());
1796 0 : LatticeExpr<Float> ratio = iif(deno > 0.0, *(residual()) / deno, 0.0);
1797 0 : residual()->copyData(ratio);
1798 0 : }
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 0 : if ((residual()->getDefaultMask() == "") && hasPB() && pblimit >=0.0) {
1808 0 : copyMask(pb(), residual());
1809 : }
1810 :
1811 0 : if ((pblimit < 0.0) && (residual()->getDefaultMask()).matches("mask0")) {
1812 0 : removeMask(residual());
1813 : }
1814 :
1815 0 : residual()->unlock();
1816 0 : }
1817 :
1818 0 : void SIImageStore::divideModelByWeight(Float pblimit, const String normtype)
1819 : {
1820 0 : LogIO os( LogOrigin("SIImageStore","divideModelByWeight",WHERE) );
1821 :
1822 :
1823 0 : if(itsUseWeight // only when needed
1824 0 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
1825 : {
1826 :
1827 0 : if( normtype=="flatsky") {
1828 :
1829 0 : LatticeExprNode LEN = max( *model());
1830 0 : os << LogIO::NORMAL1 << "Model is already flat sky with peak flux : " << LEN.getFloat();
1831 0 : os << ". No need to divide before prediction" << LogIO::POST;
1832 :
1833 0 : return;
1834 0 : }
1835 0 : else if( normtype=="flatnoise"){
1836 :
1837 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1838 : {
1839 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1840 : {
1841 :
1842 0 : itsPBScaleFactor = getPbMax(pol,chan);
1843 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1844 0 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1845 : else {
1846 :
1847 0 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1848 0 : chan, itsImageShape[3],
1849 0 : pol, itsImageShape[2],
1850 0 : *weight() );
1851 0 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1852 0 : chan, itsImageShape[3],
1853 0 : pol, itsImageShape[2],
1854 0 : *model() );
1855 0 : os << LogIO::NORMAL1 ;
1856 0 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1857 0 : os << "Dividing " << itsImageName+String(".model") ;
1858 0 : os << " by [ sqrt(weight) / " << itsPBScaleFactor ;
1859 0 : os <<" ] to get to flat sky model before prediction" << LogIO::POST;
1860 :
1861 :
1862 0 : LatticeExpr<Float> deno( sqrt( abs(*(wtsubim)) ) / itsPBScaleFactor );
1863 :
1864 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1865 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1866 0 : 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 0 : 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 0 : 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 0 : }// if not zero
1886 : }//chan
1887 : }//pol
1888 :
1889 : }
1890 0 : else if( normtype=="pbsquare"){
1891 :
1892 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1893 : {
1894 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1895 : {
1896 :
1897 0 : itsPBScaleFactor = getPbMax(pol,chan);
1898 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1899 0 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1900 : else {
1901 :
1902 0 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1903 0 : chan, itsImageShape[3],
1904 0 : pol, itsImageShape[2],
1905 0 : *weight() );
1906 0 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1907 0 : chan, itsImageShape[3],
1908 0 : pol, itsImageShape[2],
1909 0 : *model() );
1910 0 : os << LogIO::NORMAL1 ;
1911 0 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1912 0 : os << "Dividing " << itsImageName+String(".model") ;
1913 0 : os << " by [ (weight) / " << itsPBScaleFactor ;
1914 0 : os <<" ] to restore to optimal sky model before prediction" << LogIO::POST;
1915 :
1916 :
1917 0 : LatticeExpr<Float> deno( abs(*(wtsubim)) / (itsPBScaleFactor*itsPBScaleFactor) );
1918 :
1919 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1920 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1921 0 : LatticeExpr<Float> ratio( ( (*(modsubim)) * mask ) / ( deno + maskinv ) );
1922 :
1923 : //
1924 : //The above has a problem...mask is cutting out clean components found
1925 : // outside pblimit ...use below if this is what is wanted
1926 : // LatticeExpr<Float> ratio(iif(abs(*(wtsubim)) == 0.0, *modsubim, (*(modsubim))/( abs(*(wtsubim)) / (itsPBScaleFactor*itsPBScaleFactor))));
1927 :
1928 0 : IPosition ip(4,itsImageShape[0]/2,itsImageShape[1]/2,0,0);
1929 : /// Float modval = modsubim->getAt(ip);
1930 : //LatticeExprNode aminval( min(*modsubim) );
1931 : //LatticeExprNode amaxval( max(*modsubim) );
1932 : //cout << "Before ---- min : " << aminval.getFloat() << " max : " << amaxval.getFloat() << endl;
1933 :
1934 0 : modsubim->copyData(ratio);
1935 :
1936 : // cout << "Val of model before|after flattening at center for pol " << pol << " chan " << chan << " : " << modval << "|" << modsubim->getAt(ip) << " weight : " << wtsubim->getAt(ip) << endl;
1937 : //LatticeExprNode minval( min(*modsubim) );
1938 : //LatticeExprNode maxval( max(*modsubim) );
1939 : //cout << "After ---- min : " << minval.getFloat() << " max : " << maxval.getFloat() << endl;
1940 0 : }// if not zero
1941 : }//chan
1942 : }//pol
1943 :
1944 : }
1945 :
1946 : // storeImg(String("flatmodel.im"), *model());
1947 :
1948 : }
1949 0 : }
1950 :
1951 0 : void SIImageStore::multiplyModelByWeight(Float pblimit, const String normtype)
1952 : {
1953 0 : LogIO os( LogOrigin("SIImageStore","multiplyModelByWeight",WHERE) );
1954 :
1955 0 : if(itsUseWeight // only when needed
1956 0 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
1957 : {
1958 0 : 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 0 : else if( normtype=="flatnoise"){
1963 :
1964 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1965 : {
1966 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1967 : {
1968 :
1969 0 : itsPBScaleFactor = getPbMax(pol,chan);
1970 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1971 0 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1972 : else {
1973 :
1974 0 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1975 0 : chan, itsImageShape[3],
1976 0 : pol, itsImageShape[2],
1977 0 : *weight() );
1978 0 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1979 0 : chan, itsImageShape[3],
1980 0 : pol, itsImageShape[2],
1981 0 : *model() );
1982 :
1983 :
1984 :
1985 0 : os << LogIO::NORMAL1 ;
1986 0 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1987 0 : os << "Multiplying " << itsImageName+String(".model") ;
1988 0 : os << " by [ sqrt(weight) / " << itsPBScaleFactor;
1989 0 : os << " ] to take model back to flat noise with unit pb peak." << LogIO::POST;
1990 :
1991 0 : LatticeExpr<Float> deno( sqrt( abs(*(wtsubim)) ) / itsPBScaleFactor );
1992 :
1993 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1994 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1995 0 : 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 0 : modsubim->copyData(ratio);
2000 :
2001 0 : }// if not zero
2002 : }//chan
2003 : }//pol
2004 : }
2005 0 : else if( normtype=="pbsquare"){
2006 :
2007 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2008 : {
2009 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
2010 : {
2011 :
2012 0 : itsPBScaleFactor = getPbMax(pol,chan);
2013 : // cout << " pbscale : " << itsPBScaleFactor << endl;
2014 0 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
2015 : else {
2016 :
2017 0 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
2018 0 : chan, itsImageShape[3],
2019 0 : pol, itsImageShape[2],
2020 0 : *weight() );
2021 0 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
2022 0 : chan, itsImageShape[3],
2023 0 : pol, itsImageShape[2],
2024 0 : *model() );
2025 :
2026 :
2027 :
2028 0 : os << LogIO::NORMAL1 ;
2029 0 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
2030 0 : os << "Multiplying " << itsImageName+String(".model") ;
2031 0 : os << " by [ weight / " << itsPBScaleFactor*itsPBScaleFactor;
2032 0 : os << " ] to take model back to flat noise with unit pb peak." << LogIO::POST;
2033 :
2034 0 : LatticeExpr<Float> deno( abs(*(wtsubim)) / (itsPBScaleFactor*itsPBScaleFactor) );
2035 :
2036 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
2037 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
2038 0 : LatticeExpr<Float> ratio( ( (*(modsubim)) * mask ) * ( deno + maskinv ) );
2039 :
2040 : /////See comment in divmodel and divresidual for below usage
2041 : //LatticeExpr<Float> ratio ( (*(modsubim)) * sqrt( abs(*(wtsubim)) / itsPBScaleFactor) );
2042 0 : modsubim->copyData(ratio);
2043 0 : }// if not zero
2044 : }//chan
2045 : }//pol
2046 : }
2047 : }
2048 0 : }
2049 :
2050 0 : GaussianBeam SIImageStore::getPSFGaussian(Float psfcutoff)
2051 : {
2052 0 : GaussianBeam beam;
2053 : try
2054 : {
2055 0 : if( psf()->ndim() > 0 )
2056 : {
2057 0 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
2058 0 : StokesImageUtil::FitGaussianPSF( *(psf()), beam, psfcutoff);
2059 0 : }
2060 : }
2061 0 : catch(AipsError &x)
2062 : {
2063 : // LogIO os( LogOrigin("SIImageStore","getPSFGaussian",WHERE) );
2064 : // os << "Error in fitting a Gaussian to the PSF : " << x.getMesg() << LogIO::POST;
2065 0 : throw( AipsError("Error in fitting a Gaussian to the PSF : " + x.getMesg()) );
2066 0 : }
2067 :
2068 0 : return beam;
2069 0 : }
2070 :
2071 0 : void SIImageStore::makeImageBeamSet(Float psfcutoff, const Bool forcefit)
2072 : {
2073 0 : clock_t begin = clock();
2074 :
2075 0 : LogIO os( LogOrigin("SIImageStore","getPSFGaussian",WHERE) );
2076 : // For all chans/pols, call getPSFGaussian() and put it into ImageBeamSet(chan,pol).
2077 0 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2078 0 : uInt nx = itsImageShape[0];
2079 0 : uInt ny = itsImageShape[1];
2080 0 : uInt npol = itsImageShape[2];
2081 0 : uInt nchan = itsImageShape[3];
2082 0 : ImageInfo ii = psf()->imageInfo();
2083 0 : ImageBeamSet iibeamset=ii.getBeamSet();
2084 0 : if(iibeamset.nchan()==nchan && iibeamset.nstokes()==npol && forcefit==False){
2085 0 : itsPSFBeams=iibeamset;
2086 0 : itsRestoredBeams=iibeamset;
2087 0 : return;
2088 : }
2089 0 : itsPSFBeams.resize( nchan, npol );
2090 0 : itsRestoredBeams.resize(nchan, npol);
2091 : // cout << "makeImBeamSet : imshape : " << itsImageShape << endl;
2092 :
2093 0 : String blankpsfs="";
2094 0 : for( uInt chanid=0; chanid<nchan;chanid++) {
2095 0 : for( uInt polid=0; polid<npol; polid++ ) {
2096 0 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
2097 :
2098 0 : IPosition substart(4,0,0,polid,chanid);
2099 0 : IPosition substop(4,nx-1,ny-1,polid,chanid);
2100 0 : Slicer psfslice(substart, substop,Slicer::endIsLast);
2101 0 : SubImage<Float> subPsf( *psf() , psfslice, True );
2102 0 : GaussianBeam beam;
2103 :
2104 0 : Bool tryfit=True;
2105 :
2106 0 : LatticeExprNode le( max(subPsf) );
2107 0 : tryfit = le.getFloat()>0.0;
2108 0 : if(tryfit)
2109 : {
2110 : try
2111 : {
2112 0 : StokesImageUtil::FitGaussianPSF( subPsf, beam,psfcutoff );
2113 0 : itsPSFBeams.setBeam( chanid, polid, beam );
2114 0 : itsRestoredBeams.setBeam(chanid, polid, beam);
2115 : }
2116 0 : catch(AipsError &x)
2117 : {
2118 0 : os << LogIO::WARN << "[Chan" << chanid << ":Pol" << polid << "] Error Gaussian fit to PSF : " << x.getMesg() ;
2119 : // os << LogIO::POST;
2120 0 : os << " : Setting restoring beam to largest valid beam." << LogIO::POST;
2121 0 : }
2122 : }
2123 : else
2124 : {
2125 : // os << LogIO::WARN << "[Chan" << chanid << ":Pol" << polid << "] PSF is blank. Setting null restoring beam." << LogIO::POST ;
2126 0 : blankpsfs += "[C" +String::toString(chanid) + ":P" + String::toString(polid) + "] ";
2127 : }
2128 :
2129 0 : }// end of pol loop
2130 : }// end of chan loop
2131 :
2132 0 : if( blankpsfs.length() >0 )
2133 0 : 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 0 : 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 0 : if(defaultbeam.getArea("rad2")==0.0){
2146 0 : Quantity majax(1e-06,"arcsec"),minax(1e-06,"arcsec"),pa(0.0,"deg");
2147 0 : defaultbeam.setMajorMinor(majax,minax);
2148 0 : defaultbeam.setPA(pa);
2149 0 : }
2150 0 : for( uInt chanid=0; chanid<nchan;chanid++) {
2151 0 : for( uInt polid=0; polid<npol; polid++ ) {
2152 0 : if( (itsPSFBeams.getBeam(chanid, polid)).isNull() )
2153 0 : { itsPSFBeams.setBeam( chanid, polid, defaultbeam );
2154 0 : 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 0 : ii.setBeams( itsPSFBeams );
2182 : {
2183 0 : LatticeLocker lock1(*(psf()), FileLocker::Write);
2184 0 : psf()->setImageInfo(ii);
2185 0 : psf()->unlock();
2186 0 : }
2187 0 : clock_t end = clock();
2188 0 : os << LogIO::NORMAL << "Time to fit Gaussian to PSF " << double(end - begin) / CLOCKS_PER_SEC << LogIO::POST;
2189 0 : }// end of make beam set
2190 :
2191 :
2192 0 : ImageBeamSet SIImageStore::getChannelBeamSet(const Int channel){
2193 :
2194 0 : return getChannelSliceBeamSet(channel, channel);
2195 : }
2196 :
2197 :
2198 0 : ImageBeamSet SIImageStore::getChannelSliceBeamSet(const Int begChan, const Int endChan){
2199 0 : ImageBeamSet bs=getBeamSet();
2200 0 : if(bs.shape()[0]==1)
2201 0 : return bs;
2202 0 : if(begChan > endChan || begChan <0)
2203 0 : throw(AipsError("Inconsistent slice of beam in channel requested"));
2204 0 : if(bs.shape()[0] < (endChan-1))
2205 0 : throw(AipsError("beam of channel "+String::toString(endChan)+" does not exist"));
2206 0 : IPosition blc(2, begChan, 0);
2207 0 : IPosition trc(2, endChan, bs.shape()[1]-1);
2208 0 : Matrix<GaussianBeam> sliceBeam=bs.getBeams()(blc, trc);
2209 0 : ImageBeamSet subBeamSet(sliceBeam);
2210 0 : return subBeamSet;
2211 :
2212 0 : }
2213 0 : void SIImageStore::setBeamSet(const ImageBeamSet& bs){
2214 :
2215 0 : itsPSFBeams=bs;
2216 0 : }
2217 :
2218 0 : ImageBeamSet SIImageStore::getBeamSet(Float psfcutoff)
2219 : {
2220 0 : IPosition beamshp = itsPSFBeams.shape();
2221 0 : AlwaysAssert( beamshp.nelements()==2 , AipsError );
2222 0 : if( beamshp[0]==0 || beamshp[1]==0 ) {makeImageBeamSet(psfcutoff);}
2223 0 : return itsPSFBeams;
2224 0 : }
2225 :
2226 0 : void SIImageStore::printBeamSet(Bool verbose)
2227 : {
2228 0 : LogIO os( LogOrigin("SIImageStore","printBeamSet",WHERE) );
2229 0 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2230 0 : if( itsImageShape[3] == 1 && itsImageShape[2]==1 )
2231 : {
2232 0 : GaussianBeam beam = itsPSFBeams.getBeam();
2233 0 : os << "Beam : " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2234 0 : }
2235 : else
2236 : {
2237 : // per CAS-11415, verbose=True when restoringbeam='common'
2238 0 : if( verbose )
2239 : {
2240 0 : ostringstream oss;
2241 0 : oss<<"Beam";
2242 0 : Int nchan = itsImageShape[3];
2243 0 : for( Int chanid=0; chanid<nchan;chanid++) {
2244 0 : Int polid=0;
2245 : // for( Int polid=0; polid<npol; polid++ ) {
2246 0 : GaussianBeam beam = itsPSFBeams.getBeam( chanid, polid );
2247 0 : oss << " [C" << chanid << "]: " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg";
2248 0 : }//for chanid
2249 0 : os << oss.str() << LogIO::POST;
2250 0 : }
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 0 : itsPSFBeams.summarize( os, False, itsCoordSys );
2257 : // per CAS-11415 request, not turn on this one (it prints out per-channel beam in each line in the logger)
2258 : //itsPSFBeams.summarize( os, verbose, itsCoordSys );
2259 : }
2260 0 : catch(AipsError &x)
2261 : {
2262 0 : os << LogIO::WARN << "Error while printing (compact) restoring beam summary : " << x.getMesg() << LogIO::POST;
2263 0 : os << "Printing long summary" << LogIO::POST;
2264 :
2265 0 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2266 : //Int npol = itsImageShape[2];
2267 0 : Int nchan = itsImageShape[3];
2268 0 : for( Int chanid=0; chanid<nchan;chanid++) {
2269 0 : Int polid=0;
2270 : // for( Int polid=0; polid<npol; polid++ ) {
2271 0 : GaussianBeam beam = itsPSFBeams.getBeam( chanid, polid );
2272 0 : os << "Beam [C" << chanid << "]: " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2273 : //}//for polid
2274 0 : }//for chanid
2275 0 : }// catch
2276 : }
2277 : }
2278 0 : }
2279 :
2280 : /////////////////////////////// Restore all planes.
2281 :
2282 0 : void SIImageStore::restore(GaussianBeam& rbeam, String& usebeam, uInt term, Float psfcutoff)
2283 : {
2284 0 : LogIO os( LogOrigin("SIImageStore","restore",WHERE) );
2285 : // << ". Optionally, PB-correct too." << LogIO::POST;
2286 :
2287 0 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2288 0 : Int nx = itsImageShape[0];
2289 0 : Int ny = itsImageShape[1];
2290 0 : Int npol = itsImageShape[2];
2291 0 : 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 0 : IPosition beamset = itsPSFBeams.shape();
2303 0 : AlwaysAssert( beamset.nelements()==2 , AipsError );
2304 0 : if( beamset[0] != nchan || beamset[1] != npol )
2305 : {
2306 :
2307 : // Get PSF Beams....
2308 0 : ImageInfo ii = psf()->imageInfo();
2309 0 : itsPSFBeams = ii.getBeamSet();
2310 0 : itsRestoredBeams=itsPSFBeams;
2311 0 : IPosition beamset2 = itsPSFBeams.shape();
2312 :
2313 0 : AlwaysAssert( beamset2.nelements()==2 , AipsError );
2314 0 : 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 0 : }
2321 :
2322 : // toggle for printing common beam to the log
2323 0 : 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 0 : if( rbeam.isNull() && usebeam=="common") {
2329 0 : os << "Getting common beam" << LogIO::POST;
2330 0 : rbeam = CasaImageBeamSet(itsPSFBeams).getCommonBeam();
2331 0 : os << "Common Beam : " << rbeam.getMajor(Unit("arcsec")) << " arcsec, " << rbeam.getMinor(Unit("arcsec"))<< " arcsec, " << rbeam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2332 0 : printcommonbeam=True;
2333 : }
2334 0 : 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 0 : itsRestoredBeams=ImageBeamSet(rbeam);
2343 0 : GaussianBeam beam = itsRestoredBeams.getBeam();
2344 :
2345 : //if commonbeam has not printed in the log
2346 0 : if (!printcommonbeam) {
2347 0 : os << "Common Beam : " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2348 0 : printcommonbeam=True; // to avoid duplicate the info is printed to th elog
2349 : }
2350 0 : }// 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 0 : Bool emptymodel = isModelEmpty();
2357 0 : 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 0 : LatticeLocker lock1(*(image(term)), FileLocker::Write);
2361 0 : LatticeLocker lock2(*(residual(term)), FileLocker::Write);
2362 0 : LatticeLocker lock3(*(model(term)), FileLocker::Read);
2363 : //// Use beamset to restore
2364 0 : for( Int chanid=0; chanid<nchan;chanid++) {
2365 0 : for( Int polid=0; polid<npol; polid++ ) {
2366 :
2367 0 : IPosition substart(4,0,0,polid,chanid);
2368 0 : IPosition substop(4,nx-1,ny-1,polid,chanid);
2369 0 : Slicer imslice(substart, substop,Slicer::endIsLast);
2370 0 : SubImage<Float> subRestored( *image(term) , imslice, True );
2371 0 : SubImage<Float> subModel( *model(term) , imslice, False );
2372 0 : SubImage<Float> subResidual( *residual(term) , imslice, True );
2373 :
2374 0 : 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 0 : if(!printcommonbeam) {
2378 0 : 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 0 : subRestored.set(0.0);
2385 0 : if( !emptymodel ) {
2386 : // Copy model into it
2387 0 : subRestored.copyData( LatticeExpr<Float>( subModel ) );
2388 : // Smooth model by beam
2389 0 : StokesImageUtil::Convolve( subRestored, beam);
2390 : }
2391 : // Add residual image
2392 0 : 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 0 : TempImage<Float> tmpSubResidualCopy( IPosition(4,nx,ny,1,1), subResidual.coordinates());
2396 0 : tmpSubResidualCopy.copyData( subResidual );
2397 :
2398 0 : rescaleResolution(chanid, tmpSubResidualCopy, beam, itsPSFBeams.getBeam(chanid, polid));
2399 0 : subRestored.copyData( LatticeExpr<Float>( subRestored + tmpSubResidualCopy ) );
2400 0 : }
2401 : else// if no need to rescale residuals, just add the residuals.
2402 : {
2403 :
2404 :
2405 0 : subRestored.copyData( LatticeExpr<Float>( subRestored + subResidual ) );
2406 :
2407 : }
2408 :
2409 : }
2410 0 : catch(AipsError &x)
2411 : {
2412 0 : throw( AipsError("Restoration Error in chan" + String::toString(chanid) + ":pol" + String::toString(polid) + " : " + x.getMesg() ) );
2413 0 : }
2414 :
2415 0 : }// end of pol loop
2416 : }// end of chan loop
2417 :
2418 : try
2419 : {
2420 : //MSK//
2421 0 : if( hasPB() )
2422 : {
2423 0 : if( (image(term)->getDefaultMask()).matches("mask0") ) removeMask( image(term) );
2424 0 : copyMask(residual(term),image(term));
2425 : }
2426 :
2427 : // if(hasPB()){copyMask(residual(term),image(term));}
2428 0 : ImageInfo iminf = image(term)->imageInfo();
2429 : //iminf.setBeams( itsRestoredBeams);
2430 :
2431 0 : os << "Beam Set : Single beam : " << itsRestoredBeams.hasSingleBeam() << " Multi-beam : " << itsRestoredBeams.hasMultiBeam() << LogIO::DEBUG2;
2432 :
2433 0 : iminf.removeRestoringBeam();
2434 :
2435 0 : if( itsRestoredBeams.hasSingleBeam() )
2436 0 : { iminf.setRestoringBeam( itsRestoredBeams.getBeam() );}
2437 : else
2438 0 : {iminf.setBeams( itsRestoredBeams);}
2439 :
2440 0 : image(term)->setImageInfo(iminf);
2441 :
2442 0 : }
2443 0 : catch(AipsError &x)
2444 : {
2445 0 : throw( AipsError("Restoration Error : " + x.getMesg() ) );
2446 0 : }
2447 :
2448 0 : }// end of restore()
2449 :
2450 0 : GaussianBeam SIImageStore::findGoodBeam()
2451 : {
2452 0 : LogIO os( LogOrigin("SIImageStore","findGoodBeam",WHERE) );
2453 0 : IPosition beamshp = itsPSFBeams.shape();
2454 0 : AlwaysAssert( beamshp.nelements()==2 , AipsError );
2455 :
2456 : /*
2457 : if( beamshp[0] != nchan || beamshp[1] != npol )
2458 : {
2459 : // Make new beams.
2460 : os << LogIO::WARN << "Couldn't find pre-computed restoring beams. Re-fitting." << LogIO::POST;
2461 : makeImageBeamSet();
2462 : }
2463 : */
2464 :
2465 0 : Vector<Float> areas(beamshp[0]*beamshp[1]);
2466 0 : Vector<Float> axrat(beamshp[0]*beamshp[1]);
2467 0 : areas=0.0; axrat=1.0;
2468 0 : Vector<Bool> flags( areas.nelements() );
2469 0 : flags=False;
2470 :
2471 0 : Int cnt=0;
2472 0 : for( Int chanid=0; chanid<beamshp[0];chanid++) {
2473 0 : for( Int polid=0; polid<beamshp[1]; polid++ ) {
2474 0 : GaussianBeam beam = itsPSFBeams(chanid, polid);
2475 0 : if( !beam.isNull() && beam.getMajor(Unit("arcsec"))>1.1e-06 ) // larger than default filler beam.
2476 : {
2477 0 : areas[cnt] = beam.getArea( Unit("arcsec2") );
2478 0 : axrat[cnt] = beam.getMajor( Unit("arcsec") ) / beam.getMinor( Unit("arcsec") );
2479 : }
2480 : else {
2481 0 : flags[cnt] = True;
2482 : }
2483 0 : cnt++;
2484 0 : }//for chanid
2485 : }//for polid
2486 :
2487 0 : Vector<Float> fit( areas.nelements() );
2488 0 : Vector<Float> fitaxr( areas.nelements() );
2489 0 : for (Int loop=0;loop<5;loop++) {
2490 : /// Filter on outliers in PSF beam area
2491 0 : lineFit( areas, flags, fit, 0, areas.nelements()-1 );
2492 0 : Float sd = calcStd( areas , flags, fit );
2493 0 : for (uInt i=0;i<areas.nelements();i++) {
2494 0 : if( fabs( areas[i] - fit[i] ) > 3*sd ) flags[i]=True;
2495 : }
2496 : /// Filter on outliers in PSF axial ratio
2497 0 : lineFit( axrat, flags, fitaxr, 0, areas.nelements()-1 );
2498 0 : Float sdaxr = calcStd( axrat , flags, fitaxr );
2499 0 : for (uInt i=0;i<areas.nelements();i++) {
2500 0 : if( fabs( axrat[i] - fitaxr[i] ) > 3*sdaxr ) flags[i]=True;
2501 : }
2502 : }
2503 : // cout << "Original areas : " << areas << endl;
2504 : // cout << "Original axrats : " << axrat << endl;
2505 : // cout << "Flags : " << flags << endl;
2506 :
2507 : // Find max area good beam.
2508 0 : GaussianBeam goodbeam;
2509 0 : Int cid=0,pid=0;
2510 0 : Float maxval=0.0;
2511 0 : cnt=0;
2512 0 : for( Int chanid=0; chanid<beamshp[0];chanid++) {
2513 0 : for( Int polid=0; polid<beamshp[1]; polid++ ) {
2514 0 : if( flags[cnt] == False ){
2515 0 : if( areas[cnt] > maxval ) {maxval = areas[cnt]; goodbeam = itsPSFBeams.getBeam(chanid,polid);cid=chanid;pid=polid;}
2516 : }
2517 0 : cnt++;
2518 : }//polid
2519 : }//chanid
2520 :
2521 0 : os << "Picking common beam from C"<<cid<<":P"<<pid<<" : " << goodbeam.getMajor(Unit("arcsec")) << " arcsec, " << goodbeam.getMinor(Unit("arcsec"))<< " arcsec, " << goodbeam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2522 :
2523 0 : Bool badbeam=False;
2524 0 : for(uInt i=0;i<flags.nelements();i++){if(flags[i]==True) badbeam=True;}
2525 :
2526 0 : if( badbeam == True )
2527 : {
2528 0 : os << "(Ignored beams from :";
2529 0 : cnt=0;
2530 0 : for( Int chanid=0; chanid<beamshp[0];chanid++) {
2531 0 : for( Int polid=0; polid<beamshp[1]; polid++ ) {
2532 0 : if( flags[cnt] == True ){
2533 0 : os << " C"<<chanid<<":P"<<polid;
2534 : }
2535 0 : cnt++;
2536 : }//polid
2537 : }//chanid
2538 0 : os << " as outliers either by area or by axial ratio)" << LogIO::POST;
2539 : }
2540 :
2541 :
2542 : /*
2543 : // Replace 'bad' psfs with the chosen one.
2544 : if( goodbeam.isNull() ) {
2545 : os << LogIO::WARN << "No planes have non-zero restoring beams. Forcing artificial 1.0arcsec beam." << LogIO::POST;
2546 : Quantity majax(1.0,"arcsec"),minax(1.0,"arcsec"),pa(0.0,"deg");
2547 : goodbeam.setMajorMinor(majax,minax);
2548 : goodbeam.setPA(pa);
2549 : }
2550 : else {
2551 : cnt=0;
2552 : for( Int chanid=0; chanid<nchan;chanid++) {
2553 : for( Int polid=0; polid<npol; polid++ ) {
2554 : if( flags[cnt]==True )
2555 : { itsPSFBeams.setBeam( chanid, polid, goodbeam ); }
2556 : cnt++;
2557 : }// end of pol loop
2558 : }// end of chan loop
2559 : }
2560 : */
2561 :
2562 0 : return goodbeam;
2563 0 : }// end of findGoodBeam
2564 :
2565 : ///////////////////////// Funcs to calculate robust mean and fit, taking into account 'flagged' points.
2566 0 : void SIImageStore :: lineFit(Vector<Float> &data, Vector<Bool> &flag, Vector<Float> &fit, uInt lim1, uInt lim2)
2567 : {
2568 0 : float Sx = 0, Sy = 0, Sxx = 0, Sxy = 0, S = 0, a, b, sd, mn;
2569 :
2570 0 : mn = calcMean(data, flag);
2571 0 : sd = calcStd (data, flag, mn);
2572 :
2573 0 : for (uInt i = lim1; i <= lim2; i++)
2574 : {
2575 0 : if (flag[i] == False) // if unflagged
2576 : {
2577 0 : S += 1 / (sd * sd);
2578 0 : Sx += i / (sd * sd);
2579 0 : Sy += data[i] / (sd * sd);
2580 0 : Sxx += (i * i) / (sd * sd);
2581 0 : Sxy += (i * data[i]) / (sd * sd);
2582 : }
2583 : }
2584 0 : a = (Sxx * Sy - Sx * Sxy) / (S * Sxx - Sx * Sx);
2585 0 : b = (S * Sxy - Sx * Sy) / (S * Sxx - Sx * Sx);
2586 :
2587 0 : for (uInt i = lim1; i <= lim2; i++)
2588 0 : fit[i] = a + b * i;
2589 :
2590 0 : }
2591 : /* Calculate the MEAN of 'vect' ignoring values flagged in 'flag' */
2592 0 : Float SIImageStore :: calcMean(Vector<Float> &vect, Vector<Bool> &flag)
2593 : {
2594 0 : Float mean=0;
2595 0 : Int cnt=0;
2596 0 : for(uInt i=0;i<vect.nelements();i++)
2597 0 : if(flag[i]==False)
2598 : {
2599 0 : mean += vect[i];
2600 0 : cnt++;
2601 : }
2602 0 : if(cnt==0) cnt=1;
2603 0 : return mean/cnt;
2604 : }
2605 0 : Float SIImageStore :: calcStd(Vector<Float> &vect, Vector<Bool> &flag, Vector<Float> &fit)
2606 : {
2607 0 : Float std=0;
2608 0 : uInt n=0,cnt=0;
2609 0 : n = vect.nelements() < fit.nelements() ? vect.nelements() : fit.nelements();
2610 0 : for(uInt i=0;i<n;i++)
2611 0 : if(flag[i]==False)
2612 : {
2613 0 : cnt++;
2614 0 : std += (vect[i]-fit[i])*(vect[i]-fit[i]);
2615 : }
2616 0 : if(cnt==0) cnt=1;
2617 0 : return sqrt(std/cnt);
2618 :
2619 : }
2620 0 : Float SIImageStore :: calcStd(Vector<Float> &vect, Vector<Bool> &flag, Float mean)
2621 : {
2622 0 : Float std=0;
2623 0 : uInt cnt=0;
2624 0 : for(uInt i=0;i<vect.nelements();i++)
2625 0 : if(flag[i]==False)
2626 : {
2627 0 : cnt++;
2628 0 : std += (vect[i]-mean)*(vect[i]-mean);
2629 : }
2630 0 : return sqrt(std/cnt);
2631 : }
2632 :
2633 : ///////////////////////// End of Funcs to calculate robust mean and fit.
2634 :
2635 :
2636 :
2637 : /*
2638 : GaussianBeam SIImageStore::restorePlane()
2639 : {
2640 :
2641 : LogIO os( LogOrigin("SIImageStore","restorePlane",WHERE) );
2642 : // << ". Optionally, PB-correct too." << LogIO::POST;
2643 :
2644 : Bool validbeam=False;
2645 : GaussianBeam beam;
2646 : try
2647 : {
2648 : // Fit a Gaussian to the PSF.
2649 : beam = getPSFGaussian();
2650 : validbeam = True;
2651 : }
2652 : catch(AipsError &x)
2653 : {
2654 : os << LogIO::WARN << "Beam fit error : " + x.getMesg() << LogIO::POST;
2655 : }
2656 :
2657 : try
2658 : {
2659 : if( validbeam==True )
2660 : {
2661 : //os << "[" << itsImageName << "] " ; // Add when parent image name is available.
2662 : //os << "Restore with beam : " << beam.getMajor(Unit("arcmin")) << " arcmin, " << beam.getMinor(Unit("arcmin"))<< " arcmin, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2663 :
2664 : // Initialize restored image
2665 : image()->set(0.0);
2666 : // Copy model into it
2667 : image()->copyData( LatticeExpr<Float>( *(model()) ) );
2668 : // Smooth model by beam
2669 : StokesImageUtil::Convolve( *(image()), beam);
2670 : // Add residual image
2671 : image()->copyData( LatticeExpr<Float>( *(image()) + *(residual()) ) );
2672 :
2673 : // Set restoring beam into the image
2674 : ImageInfo ii = image()->imageInfo();
2675 : //ii.setRestoringBeam(beam);
2676 : ii.setBeams(beam);
2677 : image()->setImageInfo(ii);
2678 : }
2679 : }
2680 : catch(AipsError &x)
2681 : {
2682 : throw( AipsError("Restoration Error : " + x.getMesg() ) );
2683 : }
2684 :
2685 : return beam;
2686 :
2687 : }
2688 : */
2689 :
2690 0 : void SIImageStore::pbcor(uInt term)
2691 : {
2692 :
2693 0 : LogIO os( LogOrigin("SIImageStore","pbcor",WHERE) );
2694 :
2695 0 : 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 0 : LatticeLocker lock1(*(image(term)), FileLocker::Write);
2702 0 : LatticeLocker lock2(*(pb(term)), FileLocker::Write);
2703 0 : LatticeLocker lock3(*(imagepbcor(term)), FileLocker::Write);
2704 :
2705 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2706 : {
2707 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
2708 : {
2709 :
2710 0 : CountedPtr<ImageInterface<Float> > restoredsubim=makeSubImage(0,1,
2711 0 : chan, itsImageShape[3],
2712 0 : pol, itsImageShape[2],
2713 0 : *image(term) );
2714 0 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
2715 0 : chan, itsImageShape[3],
2716 0 : pol, itsImageShape[2],
2717 0 : *pb() );
2718 :
2719 0 : CountedPtr<ImageInterface<Float> > pbcorsubim=makeSubImage(0,1,
2720 0 : chan, itsImageShape[3],
2721 0 : pol, itsImageShape[2],
2722 0 : *imagepbcor(term) );
2723 :
2724 :
2725 0 : LatticeExprNode pbmax( max( *pbsubim ) );
2726 0 : Float pbmaxval = pbmax.getFloat();
2727 :
2728 0 : 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 0 : LatticeExpr<Float> thepbcor( iif( *(pbsubim) > 0.0 , (*(restoredsubim))/(*(pbsubim)) , 0.0 ) );
2737 0 : pbcorsubim->copyData( thepbcor );
2738 :
2739 0 : }// if not zero
2740 0 : }//chan
2741 : }//pol
2742 :
2743 : // Copy over the PB mask.
2744 0 : if((imagepbcor(term)->getDefaultMask()=="") && hasPB())
2745 0 : {copyMask(pb(),imagepbcor(term));}
2746 :
2747 : // Set restoring beam info
2748 0 : ImageInfo iminf = image(term)->imageInfo();
2749 : //iminf.setBeams( itsRestoredBeams );
2750 0 : imagepbcor(term)->setImageInfo(iminf);
2751 :
2752 0 : }// end pbcor
2753 :
2754 0 : Matrix<Float> SIImageStore::getSumWt(ImageInterface<Float>& target)
2755 : {
2756 0 : Record miscinfo = target.miscInfo();
2757 :
2758 0 : Matrix<Float> sumwt;
2759 0 : sumwt.resize();
2760 0 : if( miscinfo.isDefined("sumwt")
2761 0 : && (miscinfo.dataType("sumwt")==TpArrayFloat || miscinfo.dataType("sumwt")==TpArrayDouble ) )
2762 0 : { miscinfo.get( "sumwt" , sumwt ); }
2763 0 : else { sumwt.resize( IPosition(2, target.shape()[2], target.shape()[3] ) ); sumwt = 1.0; }
2764 :
2765 0 : return sumwt;
2766 0 : }
2767 :
2768 0 : void SIImageStore::setSumWt(ImageInterface<Float>& target, Matrix<Float>& sumwt)
2769 : {
2770 0 : Record miscinfo = target.miscInfo();
2771 0 : miscinfo.define("sumwt", sumwt);
2772 0 : LatticeLocker lock1(target, FileLocker::Write);
2773 0 : target.setMiscInfo( miscinfo );
2774 0 : }
2775 :
2776 :
2777 0 : Bool SIImageStore::getUseWeightImage(ImageInterface<Float>& target)
2778 : {
2779 0 : Record miscinfo = target.miscInfo();
2780 : Bool useweightimage;
2781 0 : if( miscinfo.isDefined("useweightimage") && miscinfo.dataType("useweightimage")==TpBool )
2782 0 : { miscinfo.get( "useweightimage", useweightimage ); }
2783 0 : else { useweightimage = False; }
2784 :
2785 0 : return useweightimage;
2786 0 : }
2787 : /*
2788 : Bool SIImageStore::getUseWeightImage()
2789 : {
2790 : if( ! itsParentSumWt )
2791 : return False;
2792 : else
2793 : return getUseWeightImage( *itsParentSumWt );
2794 : }
2795 : */
2796 0 : void SIImageStore::setUseWeightImage(ImageInterface<Float>& target, Bool useweightimage)
2797 : {
2798 0 : Record miscinfo = target.miscInfo();
2799 0 : miscinfo.define("useweightimage", useweightimage);
2800 0 : LatticeLocker lock1(target, FileLocker::Write);
2801 0 : target.setMiscInfo( miscinfo );
2802 0 : }
2803 :
2804 :
2805 :
2806 0 : Bool SIImageStore::divideImageByWeightVal( ImageInterface<Float>& target )
2807 : {
2808 :
2809 0 : Array<Float> lsumwt;
2810 0 : sumwt()->get( lsumwt , False ); // For MT, this will always pick the zeroth sumwt, which it should.
2811 0 : LatticeLocker lock2(target, FileLocker::Write);
2812 :
2813 0 : IPosition imshape = target.shape();
2814 :
2815 : //cerr << " SumWt : " << lsumwt << " sumwtshape : " << lsumwt.shape() << " image shape : " << imshape << endl;
2816 :
2817 0 : AlwaysAssert( lsumwt.shape()[2] == imshape[2] , AipsError ); // polplanes
2818 0 : AlwaysAssert( lsumwt.shape()[3] == imshape[3] , AipsError ); // chanplanes
2819 :
2820 0 : Bool div=False; // flag to signal if division actually happened, or weights are all 1.0
2821 :
2822 0 : for(Int pol=0; pol<lsumwt.shape()[2]; pol++)
2823 : {
2824 0 : for(Int chan=0; chan<lsumwt.shape()[3]; chan++)
2825 : {
2826 0 : IPosition pos(4,0,0,pol,chan);
2827 0 : 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 0 : chan, lsumwt.shape()[3],
2832 0 : pol, lsumwt.shape()[2],
2833 0 : target );
2834 0 : if ( lsumwt(pos) > 1e-07 ) {
2835 0 : LatticeExpr<Float> le( (*subim)/lsumwt(pos) );
2836 0 : subim->copyData( le );
2837 0 : }
2838 : else {
2839 0 : subim->set(0.0);
2840 : }
2841 0 : div=True;
2842 0 : }
2843 0 : }
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 0 : return div;
2852 0 : }
2853 :
2854 0 : void SIImageStore::normPSF(Int term)
2855 : {
2856 :
2857 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2858 : {
2859 0 : 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 0 : chan, itsImageShape[3],
2865 0 : pol, itsImageShape[2],
2866 0 : (*psf(term)) );
2867 :
2868 : std::shared_ptr<ImageInterface<Float> > subim0=makeSubImage(0,1,
2869 0 : chan, itsImageShape[3],
2870 0 : pol, itsImageShape[2],
2871 0 : (*psf(0)) );
2872 :
2873 0 : LatticeLocker lock1(*(subim), FileLocker::Write);
2874 :
2875 0 : LatticeExprNode themax( max(*(subim0)) );
2876 0 : Float maxim = themax.getFloat();
2877 :
2878 0 : if ( maxim > 1e-07 )
2879 : {
2880 0 : LatticeExpr<Float> normed( (*(subim)) / maxim );
2881 0 : subim->copyData( normed );
2882 0 : }
2883 : else
2884 : {
2885 0 : subim->set(0.0);
2886 : }
2887 0 : }//chan
2888 : }//pol
2889 :
2890 0 : }
2891 :
2892 0 : void SIImageStore::calcSensitivity()
2893 : {
2894 0 : LogIO os( LogOrigin("SIImageStore","calcSensitivity",WHERE) );
2895 :
2896 0 : Array<Float> lsumwt;
2897 0 : sumwt()->get( lsumwt , False ); // For MT, this will always pick the zeroth sumwt, which it should.
2898 :
2899 0 : IPosition shp( lsumwt.shape() );
2900 : //cout << "Sumwt shape : " << shp << " : " << lsumwt << endl;
2901 : //AlwaysAssert( shp.nelements()==4 , AipsError );
2902 :
2903 0 : os << "[" << itsImageName << "] Theoretical sensitivity (Jy/bm):" ;
2904 :
2905 0 : IPosition it(4,0,0,0,0);
2906 0 : for ( it[0]=0; it[0]<shp[0]; it[0]++)
2907 0 : for ( it[1]=0; it[1]<shp[1]; it[1]++)
2908 0 : for ( it[2]=0; it[2]<shp[2]; it[2]++)
2909 0 : for ( it[3]=0; it[3]<shp[3]; it[3]++)
2910 : {
2911 0 : if( shp[0]>1 ){os << "f"<< it[0]+(it[1]*shp[0]) << ":" ;}
2912 0 : if( shp[3]>1 ) { os << "c"<< it[3] << ":"; }
2913 0 : if( shp[2]>1 ) { os << "p"<< it[2]<< ":" ; }
2914 0 : if( lsumwt( it ) > 1e-07 )
2915 : {
2916 0 : os << 1.0/(sqrt(lsumwt(it))) << " " ;
2917 : }
2918 : else
2919 : {
2920 0 : os << "none" << " ";
2921 : }
2922 : }
2923 :
2924 0 : os << LogIO::POST;
2925 :
2926 : // Array<Float> sens = 1.0/sqrtsumwt;
2927 :
2928 :
2929 0 : }
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 0 : Float SIImageStore::getPeakResidual()
2942 : {
2943 0 : LogIO os( LogOrigin("SIImageStore","getPeakResidual",WHERE) );
2944 0 : LatticeLocker lock2 (*(residual()), FileLocker::Read);
2945 0 : ArrayLattice<Bool> pixelmask(residual()->getMask());
2946 0 : LatticeExpr<Float> resd(iif(pixelmask,abs(*residual()),0));
2947 : //LatticeExprNode pres( max(abs( *residual() ) ));
2948 0 : LatticeExprNode pres( max(resd) );
2949 0 : Float maxresidual = pres.getFloat();
2950 :
2951 0 : return maxresidual;
2952 0 : }
2953 :
2954 0 : Float SIImageStore::getPeakResidualWithinMask()
2955 : {
2956 0 : LogIO os( LogOrigin("SIImageStore","getPeakResidualWithinMask",WHERE) );
2957 : //Float minresmask, maxresmask, minres, maxres;
2958 : //findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
2959 :
2960 0 : ArrayLattice<Bool> pixelmask(residual()->getMask());
2961 : //findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
2962 0 : LatticeExpr<Float> resd(iif(((*mask()) > 0) && pixelmask ,abs((*residual())),0));
2963 0 : LatticeExprNode pres( max(resd) );
2964 0 : Float maxresidual = pres.getFloat();
2965 : //Float maxresidual=max( abs(maxresmask), abs(minresmask) );
2966 0 : return maxresidual;
2967 0 : }
2968 :
2969 : // Calculate the total model flux
2970 0 : Float SIImageStore::getModelFlux(uInt term)
2971 : {
2972 : // LogIO os( LogOrigin("SIImageStore","getModelFlux",WHERE) );
2973 0 : LatticeLocker lock2 (*(model(term)), FileLocker::Read);
2974 0 : LatticeExprNode mflux( sum( *model(term) ) );
2975 0 : Float modelflux = mflux.getFloat();
2976 : // Float modelflux = sum( model(term)->get() );
2977 :
2978 0 : return modelflux;
2979 0 : }
2980 :
2981 : // Check for non-zero model (this is different from getting model flux, for derived SIIMMT)
2982 0 : Bool SIImageStore::isModelEmpty()
2983 : {
2984 0 : if( !itsModel && (! hasModel()) ) return True;
2985 0 : else return ( fabs( getModelFlux(0) ) < 1e-08 );
2986 : }
2987 :
2988 :
2989 0 : void SIImageStore::setPSFSidelobeLevel(const Float level){
2990 :
2991 0 : itsPSFSideLobeLevel=level;
2992 0 : }
2993 : // Calculate the PSF sidelobe level...
2994 0 : Float SIImageStore::getPSFSidelobeLevel()
2995 : {
2996 0 : LogIO os( LogOrigin("SIImageStore","getPSFSidelobeLevel",WHERE) );
2997 : //cerr << "*****PSF sidelobe "<< itsPSFSideLobeLevel << endl;
2998 : /// Calculate only once, store and return for all subsequent calls.
2999 0 : if( itsPSFSideLobeLevel == 0.0 )
3000 : {
3001 :
3002 0 : ImageBeamSet thebeams = getBeamSet();
3003 0 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
3004 :
3005 : //------------------------------------------------------------
3006 0 : IPosition oneplaneshape( itsImageShape );
3007 0 : AlwaysAssert( oneplaneshape.nelements()==4, AipsError );
3008 0 : oneplaneshape[2]=1; oneplaneshape[3]=1;
3009 0 : TempImage<Float> psfbeam( oneplaneshape, itsCoordSys );
3010 :
3011 : // In a loop through channels, subtract out or mask out the main lobe
3012 0 : Float allmin=0.0, allmax=0.0;
3013 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
3014 : {
3015 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
3016 : {
3017 : std::shared_ptr<ImageInterface<Float> > onepsf=makeSubImage(0,1,
3018 0 : chan, itsImageShape[3],
3019 0 : pol, itsImageShape[2],
3020 0 : (*psf()) );
3021 :
3022 :
3023 0 : GaussianBeam beam = thebeams.getBeam( chan, pol );
3024 0 : Vector<Float> abeam(3); // Holds bmaj, bmin, pa in asec, asec, deg
3025 0 : abeam[0] = beam.getMajor().get("arcsec").getValue() * C::arcsec;
3026 0 : abeam[1] = beam.getMinor().get("arcsec").getValue() * C::arcsec;
3027 0 : abeam[2] = (beam.getPA().get("deg").getValue() + 90.0)* C::degree;
3028 :
3029 : //cout << "Beam : " << abeam << endl;
3030 :
3031 0 : StokesImageUtil::MakeGaussianPSF( psfbeam, abeam, False);
3032 :
3033 : // storeImg( String("psfbeam.im"), psfbeam );
3034 :
3035 : //Subtract from PSF plane
3036 0 : 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 0 : LatticeExprNode minval_le( min( *onepsf ) );
3044 0 : LatticeExprNode maxval_le( max( delobed ) );
3045 :
3046 0 : Float minval = minval_le.getFloat();
3047 0 : Float maxval = maxval_le.getFloat();
3048 :
3049 0 : if( minval < allmin ) allmin = minval;
3050 0 : if( maxval > allmax ) allmax = maxval;
3051 :
3052 : //cout << "Chan : " << chan << " minval : " << minval << " maxval : " << maxval << endl;
3053 :
3054 0 : }//chan
3055 : }//pol
3056 :
3057 : //------------------------------------------------------------
3058 :
3059 0 : itsPSFSideLobeLevel = max( fabs(allmin), fabs(allmax) );
3060 :
3061 : //os << "PSF min : " << allmin << " max : " << allmax << " psfsidelobelevel : " << itsPSFSideLobeLevel << LogIO::POST;
3062 :
3063 0 : }// if changed.
3064 :
3065 : // LatticeExprNode psfside( min( *psf() ) );
3066 : // itsPSFSideLobeLevel = fabs( psfside.getFloat() );
3067 :
3068 : //cout << "PSF sidelobe level : " << itsPSFSideLobeLevel << endl;
3069 0 : return itsPSFSideLobeLevel;
3070 0 : }
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 0 : Array<Double> SIImageStore::calcRobustRMS(Array<Double>& mdns, const Float pbmasklevel, const Bool fastcalc)
3087 : {
3088 0 : LogIO os( LogOrigin("SIImageStore","calcRobustRMS",WHERE) );
3089 0 : Record* regionPtr=0;
3090 0 : String LELmask("");
3091 0 : LatticeLocker lockres (*(residual()), FileLocker::Read);
3092 0 : ArrayLattice<Bool> pbmasklat(residual()->shape());
3093 0 : pbmasklat.set(False);
3094 0 : LatticeExpr<Bool> pbmask(pbmasklat);
3095 0 : if (hasPB()) {
3096 : // set bool mask: False = masked
3097 0 : pbmask = LatticeExpr<Bool> (iif(*pb() > pbmasklevel, True, False));
3098 0 : os << LogIO::DEBUG1 << "pbmask to be attached===> nfalse(pbmask.getMask())="<<nfalse(pbmask.getMask())<<endl;
3099 : }
3100 :
3101 :
3102 0 : Record thestats;
3103 0 : 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 0 : auto tempRes = std::shared_ptr<TempImage<Float>>(new TempImage<Float>(residual()->shape(), residual()->coordinates(), SDMaskHandler::memoryToUse()));
3107 0 : tempRes->copyData(*residual());
3108 0 : tempRes->attachMask(pbmask);
3109 : //thestats = SDMaskHandler::calcImageStatistics(*residual(), LELmask, regionPtr, True);
3110 0 : thestats = SDMaskHandler::calcImageStatistics(*tempRes, LELmask, regionPtr, True);
3111 0 : }
3112 : else { // older way to calculate
3113 : // use the new statistic calculation algorithm
3114 0 : Vector<Bool> dummyvec;
3115 : // TT: 2018.08.01 using revised version (the older version of this is renameed to calcRobustImageStatisticsOld)
3116 0 : thestats = SDMaskHandler::calcRobustImageStatistics(*residual(), *mask(), pbmask, LELmask, regionPtr, True, dummyvec);
3117 0 : }
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 0 : Array<Double>rmss, mads;
3134 : //thestats.get(RecordFieldId("max"), maxs);
3135 0 : thestats.get(RecordFieldId("rms"), rmss);
3136 0 : thestats.get(RecordFieldId("medabsdevmed"), mads);
3137 0 : thestats.get(RecordFieldId("median"), mdns);
3138 :
3139 : //os << LogIO::DEBUG1 << "Max : " << maxs << LogIO::POST;
3140 0 : os << LogIO::DEBUG1 << "RMS : " << rmss << LogIO::POST;
3141 0 : 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 0 : return mads*1.4826;
3147 0 : }
3148 :
3149 0 : void SIImageStore::printImageStats()
3150 : {
3151 0 : LogIO os( LogOrigin("SIImageStore","printImageStats",WHERE) );
3152 0 : Float minresmask=0, maxresmask=0, minres=0, maxres=0;
3153 : // findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
3154 0 : LatticeLocker lock1 (*(residual()), FileLocker::Read);
3155 0 : ArrayLattice<Bool> pixelmask(residual()->getMask());
3156 0 : if(hasMask())
3157 : {
3158 0 : LatticeLocker lock2 (*(mask()), FileLocker::Read);
3159 0 : findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
3160 0 : }
3161 : else
3162 : {
3163 : //LatticeExprNode pres( max( *residual() ) );
3164 0 : LatticeExprNode pres( max( iif(pixelmask,*residual(),0) ) );
3165 0 : maxres = pres.getFloat();
3166 : //LatticeExprNode pres2( min( *residual() ) );
3167 0 : LatticeExprNode pres2( min( iif(pixelmask,*residual(),0) ) );
3168 0 : minres = pres2.getFloat();
3169 0 : }
3170 :
3171 0 : os << "[" << itsImageName << "]" ;
3172 0 : os << " Peak residual (max,min) " ;
3173 0 : if( minresmask!=0.0 || maxresmask!=0.0 )
3174 0 : { os << "within mask : (" << maxresmask << "," << minresmask << ") "; }
3175 0 : os << "over full image : (" << maxres << "," << minres << ")" << LogIO::POST;
3176 :
3177 0 : os << "[" << itsImageName << "] Total Model Flux : " << getModelFlux() << LogIO::POST;
3178 :
3179 :
3180 0 : Record* regionPtr=0;
3181 0 : String LELmask("");
3182 0 : Record thestats = SDMaskHandler::calcImageStatistics(*residual(), LELmask, regionPtr, True);
3183 0 : Array<Double> maxs, mins;
3184 0 : thestats.get(RecordFieldId("max"), maxs);
3185 0 : thestats.get(RecordFieldId("min"), mins);
3186 : //os << LogIO::DEBUG1 << "Max : " << maxs << LogIO::POST;
3187 : //os << LogIO::DEBUG1 << "Min : " << mins << LogIO::POST;
3188 :
3189 0 : }
3190 :
3191 : // Calculate the total model flux
3192 0 : Float SIImageStore::getMaskSum()
3193 : {
3194 0 : LogIO os( LogOrigin("SIImageStore","getMaskSum",WHERE) );
3195 0 : LatticeLocker lock2 (*(mask()), FileLocker::Read);
3196 0 : LatticeExprNode msum( sum( *mask() ) );
3197 0 : Float masksum = msum.getFloat();
3198 :
3199 : // Float masksum = sum( mask()->get() );
3200 :
3201 0 : return masksum;
3202 0 : }
3203 :
3204 0 : Bool SIImageStore::findMinMaxLattice(const Lattice<Float>& lattice,
3205 : const Lattice<Float>& mask,
3206 : const Lattice<Bool>& pixelmask,
3207 : Float& maxAbs, Float& maxAbsMask,
3208 : Float& minAbs, Float& minAbsMask )
3209 : {
3210 :
3211 : //FOR DEGUG
3212 : //LogIO os( LogOrigin("SIImageStore","findMinMaxLattice",WHERE) );
3213 :
3214 0 : maxAbs=0.0;maxAbsMask=0.0;
3215 0 : minAbs=1e+10;minAbsMask=1e+10;
3216 0 : LatticeLocker lock1 (const_cast<Lattice<Float>& > (lattice), FileLocker::Read);
3217 0 : LatticeLocker lock2 (const_cast<Lattice<Float>& >(mask), FileLocker::Read);
3218 0 : LatticeLocker lock3 (const_cast<Lattice<Bool>& >(pixelmask), FileLocker::Read);
3219 0 : const IPosition tileShape = lattice.niceCursorShape();
3220 0 : TiledLineStepper ls(lattice.shape(), tileShape, 0);
3221 : {
3222 0 : RO_LatticeIterator<Float> li(lattice, ls);
3223 0 : RO_LatticeIterator<Float> mi(mask, ls);
3224 0 : RO_LatticeIterator<Bool> pmi(pixelmask, ls);
3225 0 : for(li.reset(),mi.reset(),pmi.reset();!li.atEnd();li++, mi++, pmi++) {
3226 0 : IPosition posMax=li.position();
3227 0 : IPosition posMin=li.position();
3228 0 : IPosition posMaxMask=li.position();
3229 0 : IPosition posMinMask=li.position();
3230 0 : Float maxVal=0.0;
3231 0 : Float minVal=0.0;
3232 0 : Float maxValMask=0.0;
3233 0 : Float minValMask=0.0;
3234 :
3235 :
3236 : // skip if lattice chunk is masked entirely.
3237 0 : if(ntrue(pmi.cursor()) > 0 ) {
3238 0 : const MaskedArray<Float> marr(li.cursor(), pmi.cursor());
3239 0 : const MaskedArray<Float> marrinmask(li.cursor() * mi.cursor(), pmi.cursor());
3240 : //minMax( minVal, maxVal, posMin, posMax, li.cursor() );
3241 0 : minMax( minVal, maxVal, posMin, posMax, marr );
3242 : //minMaxMasked(minValMask, maxValMask, posMin, posMax, li.cursor(), mi.cursor());
3243 0 : minMax(minValMask, maxValMask, posMin, posMax, marrinmask);
3244 :
3245 : //os<<"DONE minMax"<<LogIO::POST;
3246 0 : if( (maxVal) > (maxAbs) ) maxAbs = maxVal;
3247 0 : if( (maxValMask) > (maxAbsMask) ) maxAbsMask = maxValMask;
3248 :
3249 0 : if( (minVal) < (minAbs) ) minAbs = minVal;
3250 0 : if( (minValMask) < (minAbsMask) ) minAbsMask = minValMask;
3251 0 : }
3252 :
3253 0 : }
3254 0 : }
3255 :
3256 0 : return True;
3257 :
3258 :
3259 0 : }
3260 :
3261 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3262 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3263 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3264 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3265 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3266 :
3267 : //
3268 : //-------------------------------------------------------------
3269 : // Initialize the internals of the object. Perhaps other such
3270 : // initializations in the constructors can be moved here too.
3271 : //
3272 0 : void SIImageStore::init()
3273 : {
3274 0 : imageExts.resize(MAX_IMAGE_IDS);
3275 :
3276 0 : imageExts(MASK)=".mask";
3277 0 : imageExts(PSF)=".psf";
3278 0 : imageExts(MODEL)=".model";
3279 0 : imageExts(RESIDUAL)=".residual";
3280 0 : imageExts(WEIGHT)=".weight";
3281 0 : imageExts(IMAGE)=".image";
3282 0 : imageExts(SUMWT)=".sumwt";
3283 0 : imageExts(GRIDWT)=".gridwt";
3284 0 : imageExts(PB)=".pb";
3285 0 : imageExts(FORWARDGRID)=".forward";
3286 0 : imageExts(BACKWARDGRID)=".backward";
3287 0 : imageExts(IMAGEPBCOR)=".image.pbcor";
3288 :
3289 0 : itsParentPsf = itsPsf;
3290 0 : itsParentModel=itsModel;
3291 0 : itsParentResidual=itsResidual;
3292 0 : itsParentWeight=itsWeight;
3293 0 : itsParentImage=itsImage;
3294 0 : itsParentSumWt=itsSumWt;
3295 0 : itsParentMask=itsMask;
3296 0 : itsParentImagePBcor=itsImagePBcor;
3297 :
3298 : // cout << "parent shape : " << itsParentImageShape << " shape : " << itsImageShape << endl;
3299 0 : itsParentImageShape = itsImageShape;
3300 0 : itsParentCoordSys = itsCoordSys;
3301 :
3302 0 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 ) { itsImageShape=IPosition(4,0,0,0,0); }
3303 :
3304 0 : itsOpened=0;
3305 :
3306 0 : itsPSFSideLobeLevel=0.0;
3307 0 : itsReadLock=nullptr;
3308 0 : itsDataPolRep=StokesImageUtil::UNKNOWN; //Should throw an exception if it is not initialized properly
3309 0 : }
3310 :
3311 :
3312 0 : void SIImageStore::regridToModelImage( ImageInterface<Float> &inputimage, Int term )
3313 : {
3314 : try
3315 : {
3316 :
3317 : //output coords
3318 0 : IPosition outshape = itsImageShape;
3319 0 : CoordinateSystem outcsys = itsCoordSys;
3320 0 : DirectionCoordinate outDirCsys = outcsys.directionCoordinate();
3321 0 : SpectralCoordinate outSpecCsys = outcsys.spectralCoordinate();
3322 :
3323 : // do regrid
3324 0 : IPosition axes(3,0, 1, 2);
3325 0 : IPosition inshape = inputimage.shape();
3326 0 : CoordinateSystem incsys = inputimage.coordinates();
3327 0 : DirectionCoordinate inDirCsys = incsys.directionCoordinate();
3328 0 : SpectralCoordinate inSpecCsys = incsys.spectralCoordinate();
3329 :
3330 0 : Vector<Int> dirAxes = CoordinateUtil::findDirectionAxes(incsys);
3331 0 : axes(0) = dirAxes(0);
3332 0 : axes(1) = dirAxes(1);
3333 0 : axes(2) = CoordinateUtil::findSpectralAxis(incsys);
3334 :
3335 : try {
3336 0 : ImageRegrid<Float> imregrid;
3337 0 : imregrid.regrid( *(model(term)), Interpolate2D::LINEAR, axes, inputimage );
3338 0 : } catch (AipsError &x) {
3339 0 : throw(AipsError("ImageRegrid error : "+ x.getMesg()));
3340 0 : }
3341 :
3342 0 : }catch(AipsError &x)
3343 : {
3344 0 : throw(AipsError("Error in regridding input model image to target coordsys : " + x.getMesg()));
3345 0 : }
3346 0 : }
3347 :
3348 : //
3349 : //---------------------------------------------------------------
3350 : //
3351 0 : void SIImageStore::makePersistent(String& fileName)
3352 : {
3353 0 : LogIO logIO(LogOrigin("SIImageStore", "makePersistent"));
3354 0 : ofstream outFile; outFile.open(fileName.c_str(),std::ofstream::out);
3355 0 : if (!outFile) logIO << "Failed to open filed \"" << fileName << "\"" << LogIO::EXCEPTION;
3356 : // String itsImageName;
3357 0 : outFile << "itsImageNameBase: " << itsImageName << endl;
3358 :
3359 : //IPosition itsImageShape;
3360 0 : outFile << "itsImageShape: " << itsImageShape.nelements() << " ";
3361 0 : for (uInt i=0;i<itsImageShape.nelements(); i++) outFile << itsImageShape(i) << " "; outFile << endl;
3362 :
3363 : // Don't know what to do with this. Looks like this gets
3364 : // filled-in from one of the images. So load this from one of the
3365 : // images if they exist?
3366 : //CoordinateSystem itsCoordSys;
3367 :
3368 : // Int itsNFacets;
3369 0 : outFile << "itsNFacets: " << itsNFacets << endl;
3370 0 : outFile << "itsUseWeight: " << itsUseWeight << endl;
3371 :
3372 :
3373 : // Misc Information to go into the header.
3374 : // Record itsMiscInfo;
3375 0 : itsMiscInfo.print(outFile);
3376 :
3377 : // std::shared_ptr<ImageInterface<Float> > itsMask, itsPsf, itsModel, itsResidual, itsWeight, itsImage, itsSumWt;
3378 : // std::shared_ptr<ImageInterface<Complex> > itsForwardGrid, itsBackwardGrid;
3379 :
3380 0 : Vector<Bool> ImageExists(MAX_IMAGE_IDS);
3381 0 : if ( ! itsMask ) ImageExists(MASK)=False;
3382 0 : if ( ! itsPsf ) ImageExists(PSF)=False;
3383 0 : if ( ! itsModel ) ImageExists(MODEL)=False;
3384 0 : if ( ! itsResidual ) ImageExists(RESIDUAL)=False;
3385 0 : if ( ! itsWeight ) ImageExists(WEIGHT)=False;
3386 0 : if ( ! itsImage ) ImageExists(IMAGE)=False;
3387 0 : if ( ! itsSumWt ) ImageExists(SUMWT)=False;
3388 0 : if ( ! itsGridWt ) ImageExists(GRIDWT)=False;
3389 0 : if ( ! itsPB ) ImageExists(PB)=False;
3390 :
3391 0 : if ( ! itsForwardGrid ) ImageExists(FORWARDGRID)=False;
3392 0 : if ( ! itsBackwardGrid ) ImageExists(BACKWARDGRID)=False;
3393 :
3394 0 : outFile << "ImagesExist: " << ImageExists << endl;
3395 0 : }
3396 : //
3397 : //---------------------------------------------------------------
3398 : //
3399 0 : void SIImageStore::recreate(String& fileName)
3400 : {
3401 0 : LogIO logIO(LogOrigin("SIImageStore", "recreate"));
3402 0 : ifstream inFile; inFile.open(fileName.c_str(),std::ofstream::out);
3403 0 : if (!inFile) logIO << "Failed to open filed \"" << fileName << "\"" << LogIO::EXCEPTION;
3404 :
3405 0 : String token;
3406 0 : inFile >> token; if (token == "itsImageNameBase:") inFile >> itsImageName;
3407 :
3408 0 : inFile >> token;
3409 0 : if (token=="itsImageShape:")
3410 : {
3411 : Int n;
3412 0 : inFile >> n;
3413 0 : itsImageShape.resize(n);
3414 0 : for (Int i=0;i<n; i++) inFile >> itsImageShape(i);
3415 : }
3416 :
3417 : // Int itsNFacets;
3418 0 : inFile >> token; if (token=="itsNFacets:") inFile >> itsNFacets;
3419 0 : inFile >> token; if (token=="itsUseWeight:") inFile >> itsUseWeight;
3420 :
3421 0 : Bool coordSysLoaded=False;
3422 0 : String itsName;
3423 : try
3424 : {
3425 0 : itsName=itsImageName+imageExts(MASK);casa::openImage(itsName, itsMask);
3426 0 : if (coordSysLoaded==False) {itsCoordSys=itsMask->coordinates(); itsMiscInfo=itsImage->miscInfo();coordSysLoaded=True;}
3427 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3428 : try
3429 : {
3430 0 : itsName=itsImageName+imageExts(MODEL);casa::openImage(itsName, itsModel);
3431 0 : if (coordSysLoaded==False) {itsCoordSys=itsModel->coordinates(); itsMiscInfo=itsModel->miscInfo();coordSysLoaded=True;}
3432 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3433 : try
3434 : {
3435 0 : itsName=itsImageName+imageExts(RESIDUAL);casa::openImage(itsName, itsResidual);
3436 0 : if (coordSysLoaded==False) {itsCoordSys=itsResidual->coordinates(); itsMiscInfo=itsResidual->miscInfo();coordSysLoaded=True;}
3437 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3438 : try
3439 : {
3440 0 : itsName=itsImageName+imageExts(WEIGHT);casa::openImage(itsName, itsWeight);
3441 0 : if (coordSysLoaded==False) {itsCoordSys=itsWeight->coordinates(); itsMiscInfo=itsWeight->miscInfo();coordSysLoaded=True;}
3442 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3443 : try
3444 : {
3445 0 : itsName=itsImageName+imageExts(IMAGE);casa::openImage(itsName, itsImage);
3446 0 : if (coordSysLoaded==False) {itsCoordSys=itsImage->coordinates(); itsMiscInfo=itsImage->miscInfo();coordSysLoaded=True;}
3447 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3448 : try
3449 : {
3450 0 : itsName=itsImageName+imageExts(SUMWT);casa::openImage(itsName, itsSumWt);
3451 0 : if (coordSysLoaded==False) {itsCoordSys=itsSumWt->coordinates(); itsMiscInfo=itsSumWt->miscInfo();coordSysLoaded=True;}
3452 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3453 : try
3454 : {
3455 0 : itsName=itsImageName+imageExts(PSF);casa::openImage(itsName, itsPsf);
3456 0 : if (coordSysLoaded==False) {itsCoordSys=itsPsf->coordinates(); itsMiscInfo=itsPsf->miscInfo();coordSysLoaded=True;}
3457 0 : } catch (AipsIO& x) {logIO << "\"" << itsName << "\" not found." << LogIO::WARN;};
3458 : try
3459 : {
3460 0 : casa::openImage(itsImageName+imageExts(FORWARDGRID), itsForwardGrid);
3461 0 : casa::openImage(itsImageName+imageExts(BACKWARDGRID), itsBackwardGrid);
3462 : }
3463 0 : catch (AipsError& x)
3464 : {
3465 0 : logIO << "Did not find forward and/or backward grid. Just say'n..." << LogIO::POST;
3466 0 : }
3467 :
3468 0 : }
3469 : //////////////
3470 0 : Bool SIImageStore::intersectComplexImage(const String& ComplexImageName){
3471 0 : Vector<Int> whichStokes(0);
3472 0 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
3473 0 : whichStokes, itsDataPolRep);
3474 :
3475 :
3476 : //cerr <<"itsDataPolRep " << itsDataPolRep << endl;
3477 :
3478 0 : CountedPtr<PagedImage<Complex> > compliantImage =nullptr;
3479 : {
3480 0 : PagedImage<Complex>inputImage(ComplexImageName);
3481 0 : IPosition theShape=itsImageShape;
3482 0 : theShape(0)=inputImage.shape()(0);
3483 0 : theShape(1)=inputImage.shape()(1);
3484 0 : CoordinateSystem inpcsys=inputImage.coordinates();
3485 0 : Vector<Double> refpix=cimageCoord.referencePixel();
3486 0 : refpix(0)+=(theShape(0)-itsImageShape(0))/2.0;
3487 0 : refpix(1)+=(theShape(1)-itsImageShape(1))/2.0;
3488 0 : cimageCoord.setReferencePixel(refpix);
3489 0 : String tmpImage=File::newUniqueName(".", "TempImage").absoluteName();
3490 0 : compliantImage=new PagedImage<Complex>(theShape, cimageCoord, tmpImage);
3491 0 : compliantImage->set(0.0);
3492 0 : IPosition iblc(theShape.nelements(),0);
3493 0 : IPosition itrc=theShape-1;
3494 : //cerr << "blc " << iblc << " trc " << itrc << " shape " << theShape << endl;
3495 0 : LCBox lbox(iblc, itrc, theShape);
3496 0 : ImageRegion imagreg(WCBox(lbox, cimageCoord));
3497 :
3498 :
3499 0 : SubImage<Complex> subim(inputImage, imagreg, false);
3500 : //cerr << "shapes " << inputImage.shape() << " sub " << subim.shape() << " compl " << compliantImage->shape() << endl;
3501 0 : compliantImage->copyData(subim);
3502 :
3503 0 : }
3504 0 : TableUtil::deleteTable(ComplexImageName);
3505 0 : compliantImage->rename(ComplexImageName);
3506 0 : return True;
3507 :
3508 0 : }
3509 :
3510 :
3511 : //////////////////////////////////////////////////////////////////////////////////////////////////////
3512 : //////////////////////////////////////////////////////////////////////////////////////////////////////
3513 :
3514 : } //# NAMESPACE CASA - END
3515 :
|
