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