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 0 : SIImageStore::SIImageStore()
111 : {
112 :
113 0 : itsPsf.reset( );
114 0 : itsModel.reset( );
115 0 : itsResidual.reset( );
116 0 : itsWeight.reset( );
117 0 : itsImage.reset( );
118 0 : itsMask.reset( );
119 0 : itsGridWt.reset( );
120 0 : itsPB.reset( );
121 0 : itsImagePBcor.reset();
122 0 : itsTempWorkIm.reset();
123 :
124 0 : itsSumWt.reset( );
125 0 : itsOverWrite = False;
126 0 : itsUseWeight = False;
127 0 : itsPBScaleFactor = 1.0;
128 :
129 0 : itsNFacets = 1;
130 0 : itsFacetId = 0;
131 0 : itsNChanChunks = 1;
132 0 : itsChanId = 0;
133 0 : itsNPolChunks = 1;
134 0 : itsPolId = 0;
135 :
136 0 : itsImageShape = IPosition(4,0,0,0,0);
137 0 : itsImageName = String("");
138 0 : itsCoordSys = CoordinateSystem();
139 0 : itsMiscInfo = Record();
140 :
141 0 : itsIsSingleDishStore = False;
142 :
143 0 : init();
144 : // validate();
145 :
146 0 : }
147 :
148 : // Used from SynthesisNormalizer::makeImageStore()
149 73 : 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 73 : )
160 : // TODO : Add parameter to indicate weight image shape.
161 : {
162 146 : LogIO os( LogOrigin("SIImageStore","Open new Images",WHERE) );
163 :
164 :
165 73 : itsPsf.reset( );
166 73 : itsModel.reset( );
167 73 : itsResidual.reset( );
168 73 : itsWeight.reset( );
169 73 : itsImage.reset( );
170 73 : itsMask.reset( );
171 73 : itsGridWt.reset( );
172 73 : itsPB.reset( );
173 73 : itsImagePBcor.reset( );
174 73 : itsTempWorkIm.reset();
175 :
176 73 : itsSumWt.reset( );
177 73 : itsOverWrite = False; // Hard Coding this. See CAS-6937. overwrite;
178 73 : itsUseWeight = useweightimage;
179 73 : itsPBScaleFactor = 1.0;
180 :
181 73 : itsNFacets = 1;
182 73 : itsFacetId = 0;
183 73 : itsNChanChunks = 1;
184 73 : itsChanId = 0;
185 73 : itsNPolChunks = 1;
186 73 : itsPolId = 0;
187 :
188 73 : itsImageName = imagename;
189 73 : itsCoordSys = imcoordsys;
190 73 : itsImageShape = imshape;
191 73 : itsObjectName = objectname;
192 73 : itsMiscInfo = miscinfo;
193 :
194 73 : itsIsSingleDishStore = issingledishstore;
195 :
196 73 : init();
197 :
198 73 : validate();
199 73 : }
200 :
201 : // Used from SynthesisNormalizer::makeImageStore()
202 : // This constructor creates an Image Store from images on disk
203 838 : SIImageStore::SIImageStore(
204 : const String &imagename,
205 : const Bool ignorefacets,
206 : const Bool noRequireSumwt,
207 838 : const Bool makeSingleDishStore)
208 : {
209 1676 : LogIO os( LogOrigin("SIImageStore","Open existing Images",WHERE) );
210 :
211 : { // Initialize some members
212 838 : itsImageName = imagename;
213 :
214 838 : itsPsf.reset( );
215 838 : itsModel.reset( );
216 838 : itsResidual.reset( );
217 838 : itsWeight.reset( );
218 838 : itsImage.reset( );
219 838 : itsMask.reset( );
220 838 : itsGridWt.reset( );
221 838 : itsPB.reset( );
222 838 : itsImagePBcor.reset( );
223 838 : itsTempWorkIm.reset();
224 838 : itsMiscInfo = Record();
225 :
226 838 : itsSumWt.reset( );
227 838 : itsNFacets = 1;
228 838 : itsFacetId = 0;
229 838 : itsNChanChunks = 1;
230 838 : itsChanId = 0;
231 838 : itsNPolChunks = 1;
232 838 : itsPolId = 0;
233 :
234 838 : itsOverWrite = False;
235 :
236 838 : itsIsSingleDishStore = makeSingleDishStore;
237 : }
238 :
239 : // Need to do this init now so that imageExts is initialized
240 838 : 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 838 : 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 838 : std::shared_ptr<ImageInterface<Float> > imptr;
257 :
258 838 : auto haveImage = False;
259 838 : for (auto imageId : imageIds) {
260 838 : const auto imageName = imageFullName(imageId);
261 838 : if (doesImageExist(imageName)) {
262 838 : if (imageId == SIImageStore::GRIDWT) {
263 0 : constexpr auto preserveOldComment = True;
264 : // How can this be right ?
265 : }
266 838 : buildImage(imptr, imageName);
267 838 : haveImage = True;
268 838 : break;
269 : }
270 838 : }
271 :
272 838 : if (haveImage) {
273 838 : itsObjectName = imptr->imageInfo().objectName();
274 838 : itsImageShape = imptr->shape();
275 838 : itsCoordSys = imptr->coordinates();
276 838 : 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 838 : }
291 :
292 : { // Handle special case: psf or residual exist
293 : // Should we update things here when itsIsSingleDishStore = True ?
294 838 : if ( doesImageExist( imageFullName(RESIDUAL) )
295 838 : or doesImageExist( imageFullName(PSF) ) ) {
296 838 : if ( doesImageExist( imageFullName(SUMWT)) ) {
297 838 : std::shared_ptr<ImageInterface<Float> > imptr;
298 838 : buildImage(imptr, imageFullName(SUMWT) );
299 838 : itsNFacets = imptr->shape()[0];
300 838 : itsFacetId = 0;
301 838 : itsUseWeight = getUseWeightImage( *imptr );
302 838 : 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 838 : itsCoordSys = imptr->coordinates();
306 838 : itsMiscInfo =imptr->miscInfo();
307 1676 : if ( itsUseWeight
308 838 : 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 838 : } else {
315 0 : 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 0 : << LogIO::POST;
324 0 : std::shared_ptr<ImageInterface<Float> > imptr;
325 0 : if ( doesImageExist( imageFullName(RESIDUAL) ) ) {
326 0 : buildImage(imptr, imageFullName(RESIDUAL) );
327 : }
328 : else {
329 0 : buildImage(imptr, imageFullName(PSF) );
330 : }
331 :
332 0 : itsNFacets = 1;
333 0 : itsFacetId = 0;
334 0 : itsUseWeight = False;
335 0 : itsPBScaleFactor = 1.0;
336 0 : itsCoordSys = imptr->coordinates();
337 0 : itsMiscInfo = imptr->miscInfo();
338 0 : }
339 : }
340 : } // if psf or residual exist
341 : } // Handle special case
342 :
343 838 : if (ignorefacets == True) itsNFacets = 1;
344 :
345 : // Why again ?
346 838 : init();
347 :
348 838 : validate();
349 838 : }
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 273 : 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 273 : const Bool useweightimage)
374 : {
375 :
376 :
377 273 : itsPsf=psfim;
378 273 : itsModel=modelim;
379 273 : itsResidual=residim;
380 : /* if(residim){
381 : LatticeLocker lock1(*itsResidual, FileLocker::Read);
382 : cerr << "RESIDMAX " << max(itsResidual->get()) << " " << max(residim->get()) << endl;
383 : }*/
384 273 : itsWeight=weightim;
385 273 : itsImage=restoredim;
386 273 : itsMask=maskim;
387 :
388 273 : itsSumWt=sumwtim;
389 :
390 273 : itsGridWt=gridwtim;
391 273 : itsPB=pbim;
392 273 : itsImagePBcor=restoredpbcorim;
393 273 : itsTempWorkIm=tempworkimage;
394 :
395 273 : itsPBScaleFactor=1.0;// No need to set properly here as it will be computed in makePSF.
396 :
397 273 : itsObjectName = objectname;
398 273 : itsMiscInfo = miscinfo;
399 :
400 273 : itsNFacets = nfacets;
401 273 : itsFacetId = facet;
402 273 : itsNChanChunks = nchanchunks;
403 273 : itsChanId = chan;
404 273 : itsNPolChunks = npolchunks;
405 273 : itsPolId = pol;
406 :
407 273 : itsOverWrite=False;
408 273 : itsUseWeight=useweightimage;
409 :
410 273 : itsParentImageShape = imshape;
411 273 : itsImageShape = imshape;
412 :
413 273 : itsParentCoordSys = csys;
414 273 : itsCoordSys = csys; // Hopefully this doesn't change for a reference image
415 273 : itsImageName = imagename;
416 :
417 273 : itsIsSingleDishStore = False;
418 :
419 : //-----------------------
420 273 : init(); // Connect parent pointers to the images.
421 : //-----------------------
422 :
423 : // Set these to null, to be set later upon first access.
424 273 : itsPsf.reset( );
425 273 : itsModel.reset( );
426 273 : itsResidual.reset( );
427 273 : itsWeight.reset( );
428 273 : itsImage.reset( );
429 273 : itsMask.reset( );
430 273 : itsSumWt.reset( );
431 273 : itsPB.reset( );
432 :
433 273 : validate();
434 273 : }
435 :
436 2222 : void SIImageStore::validate()
437 : {
438 : // There are two valid states. Check for at least one of them.
439 2222 : Bool inValidState = False;
440 :
441 2222 : stringstream oss;
442 : { // Initialize error message
443 : oss
444 2222 : << "shape:" << itsImageShape
445 2222 : << " parentimageshape:" << itsParentImageShape
446 2222 : << " nfacets:" << itsNFacets << "x" << itsNFacets
447 2222 : << " facetid:" << itsFacetId
448 2222 : << " nchanchunks:" << itsNChanChunks << " chanid:" << itsChanId
449 2222 : << " npolchunks:" << itsNPolChunks << " polid:" << itsPolId
450 2222 : << " coord-dim:" << itsCoordSys.nPixelAxes()
451 2222 : << " psf/res:" << (hasPsf() or hasResidual());
452 2222 : if ( hasSumWt() ) oss << " sumwtshape : " << sumwt()->shape();
453 2222 : oss << endl;
454 : }
455 :
456 : try {
457 :
458 2222 : if ( itsCoordSys.nPixelAxes() != 4 ) inValidState = False;
459 :
460 : // (1) Regular imagestore
461 2222 : if (
462 2222 : itsNFacets == 1 and itsFacetId == 0
463 2222 : and itsNChanChunks == 1 and itsChanId == 0
464 2222 : and itsNPolChunks == 1 and itsPolId == 0 ) {
465 2222 : Bool sumWtShapeOK = hasSumWt() and sumwt()->shape()[0] == 1;
466 2222 : if ( itsImageShape.isEqual(itsParentImageShape)
467 2222 : and ( sumWtShapeOK or not hasSumWt() )
468 4444 : and itsParentImageShape.product() > 0 ) inValidState = True;
469 2222 : }
470 : // (2) Reference Sub Imagestore
471 0 : else if (
472 0 : ( itsNFacets > 1 and itsFacetId >= 0 )
473 0 : or ( itsNChanChunks > 1 and itsChanId >= 0 )
474 0 : or ( itsNPolChunks > 1 and itsPolId >= 0 ) ) {
475 : // If shape is still unset, even when the first image has been made....
476 : Bool imgShapeOK1 =
477 0 : ( itsImageShape.product() > 0 and ( hasPsf() or hasResidual() ) );
478 : Bool imgShapeOK2 =
479 0 : ( itsImageShape.isEqual(IPosition(4,0,0,0,0)) and
480 0 : ( not hasPsf() and not hasResidual() )
481 0 : );
482 0 : Bool imgShapeOK = imgShapeOK1 or imgShapeOK2;
483 : Bool sumWtShapeOK =
484 0 : hasSumWt() and sumwt()->shape()[0] == 1; // One facet only.
485 :
486 0 : if ( imgShapeOK and ( itsParentImageShape.product() > 0 )
487 0 : and ( itsFacetId < itsNFacets*itsNFacets )
488 : // and (itsChanId<=itsNChanChunks) // chanchunks can be larger...
489 0 : and ( itsPolId < itsNPolChunks )
490 0 : 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 2222 : if ( not inValidState ) {
500 0 : throw AipsError(
501 0 : "Internal Error : Invalid ImageStore state : " + oss.str()
502 0 : );
503 : }
504 :
505 2222 : }
506 :
507 :
508 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
509 :
510 273 : 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 273 : 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 5073 : 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 5073 : std::shared_ptr<ImageInterface<Float> > imPtr;
528 5073 : IPosition useShape( itsParentImageShape );
529 :
530 5073 : if( dosumwt ) // change shape to sumwt image shape.
531 : {
532 838 : useShape[0] = nfacetsperside;
533 838 : 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 5073 : Bool dbg=False;
553 5073 : if( doesImageExist( imagenamefull ) )
554 : {
555 : ///// not used since itsOverWrite is hardcoded to FALSE (CAS-6937)
556 4565 : 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 4565 : if(dbg) cout << "Trying to open existing image : "<< imagenamefull << endl;
589 : try{
590 4565 : buildImage( imPtr, imagenamefull ) ;
591 :
592 4565 : 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 3800 : if( useShape.product()>0 && ! useShape.isEqual( imPtr->shape() ) )
598 : {
599 0 : ostringstream oo1,oo2;
600 0 : oo1 << useShape; oo2 << imPtr->shape();
601 0 : throw( AipsError( "There is a shape mismatch between existing images ("+oo2.str()+") and current parameters ("+oo1.str()+"). If you are attempting to restart a run with a new image shape, please change imagename and supply the old model or mask as inputs (via the startmodel or mask parameters) so that they can be regridded to the new shape before continuing." ) );
602 0 : }
603 :
604 :
605 :
606 3800 : 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 0 : 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 0 : if ( ! itsParentCoordSys.near( imPtr->coordinates() ) )
627 : {
628 : //cout << " Still different..." << endl;
629 0 : throw( AipsError( "There is a coordinate system mismatch between existing images on disk and current parameters ("+itsParentCoordSys.errorMessage()+"). If you are attempting to restart a run, please change imagename and supply the old model or mask as inputs (via the startmodel or mask parameters) so that they can be regridded to the new coordinate system before continuing. " ) );
630 : }
631 : //}
632 : }
633 : }// not dosumwt
634 : }
635 0 : catch (AipsError &x){
636 0 : throw( AipsError("Cannot open existing image : "+imagenamefull+" : " + x.getMesg() ) );
637 0 : }
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 508 : if(dbg) cout << "Trying to open new image named : " << imagenamefull << endl;
649 : try{
650 508 : buildImage(imPtr, useShape, itsParentCoordSys, imagenamefull) ;
651 : // initialize to zeros...
652 : {
653 508 : LatticeLocker lock1(*imPtr, FileLocker::Write);
654 508 : imPtr->set(0.0);
655 508 : imPtr->flush();
656 508 : imPtr->unlock();
657 508 : }
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 10146 : return imPtr;
698 5073 : }
699 :
700 508 : void SIImageStore::buildImage(std::shared_ptr<ImageInterface<Float> > &imptr, IPosition shape,
701 : CoordinateSystem csys, const String name)
702 : {
703 1016 : LogIO os( LogOrigin("SIImageStore", "Open non-existing image", WHERE) );
704 508 : os <<"Opening image, name: " << name << LogIO::DEBUG1;
705 :
706 508 : itsOpened++;
707 : //For some reason cannot open a new paged image with UserNoread directly
708 : {
709 508 : PagedImage<Float> godot(shape, csys, name);
710 508 : godot.unlock();
711 508 : }
712 508 : 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 508 : imptr.reset( new PagedImage<Float> (name, locktype) );
718 : {
719 508 : LatticeLocker lock1(*imptr, FileLocker::Write);
720 508 : initMetaInfo(imptr, name);
721 508 : imptr->unlock();
722 508 : }
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 508 : }
737 :
738 6241 : void SIImageStore::buildImage(std::shared_ptr<ImageInterface<Float> > &imptr, const String name)
739 : {
740 12482 : LogIO os(LogOrigin("SIImageStore", "Open existing Images", WHERE));
741 6241 : os <<"Opening image, name: " << name << LogIO::DEBUG1;
742 :
743 6241 : itsOpened++;
744 6241 : if ( Table::isReadable(name) ) {
745 6241 : 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 6241 : Table table(name, locktype);
761 6241 : String type = table.tableInfo().type();
762 6241 : if ( type != TableInfo::type(TableInfo::PAGEDIMAGE) ) {
763 :
764 0 : imptr.reset( new PagedImage<Float>( table ) );
765 0 : imptr->unlock();
766 0 : return;
767 : }
768 6241 : }
769 :
770 6241 : LatticeBase* latt =ImageOpener::openImage(name);
771 6241 : if (not latt) {
772 0 : throw AipsError("Error in opening Image : "+name);
773 : }
774 6241 : DataType dtype = latt->dataType();
775 6241 : if (dtype == TpFloat) {
776 6241 : 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 6241 : }
805 :
806 : /**
807 : * Sets ImageInfo and MiscInfo on an image
808 : *
809 : * @param imptr image to initialize
810 : */
811 508 : 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 508 : LatticeLocker lock1(*imptr, FileLocker::Write);
817 508 : if (not itsObjectName.empty()) {
818 508 : ImageInfo info = imptr->imageInfo();
819 508 : info.setObjectName(itsObjectName);
820 508 : imptr->setImageInfo(info);
821 508 : imptr->setMiscInfo(itsMiscInfo);
822 508 : } 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 508 : imptr->unlock();
830 508 : }
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 635 : 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 635 : Int nx_facets=Int(sqrt(Double(nfacets)));
853 1270 : LogIO os( LogOrigin("SynthesisImager","makeFacet") );
854 : // Make the output image
855 635 : Slicer imageSlicer;
856 :
857 : // Add checks for all dimensions..
858 635 : 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 635 : IPosition imshp=image.shape();
863 635 : IPosition blc(imshp.nelements(), 0);
864 635 : IPosition trc=imshp-1;
865 635 : IPosition inc(imshp.nelements(), 1);
866 :
867 : /// Facets
868 635 : Int facetx = facet % nx_facets;
869 635 : Int facety = (facet - facetx) / nx_facets;
870 635 : Int sizex = imshp(0) / nx_facets;
871 635 : Int sizey = imshp(1) / nx_facets;
872 635 : blc(1) = facety * sizey;
873 635 : trc(1) = blc(1) + sizey - 1;
874 635 : blc(0) = facetx * sizex;
875 635 : trc(0) = blc(0) + sizex - 1;
876 :
877 : /// Pol Chunks
878 635 : Int sizepol = imshp(2) / npolchunks;
879 635 : blc(2) = pol * sizepol;
880 635 : trc(2) = blc(2) + sizepol - 1;
881 :
882 : /// Chan Chunks
883 635 : Int sizechan = imshp(3) / nchanchunks;
884 635 : blc(3) = chan * sizechan;
885 635 : trc(3) = blc(3) + sizechan - 1;
886 :
887 635 : LCBox::verify(blc, trc, inc, imshp);
888 635 : Slicer imslice(blc, trc, inc, Slicer::endIsLast);
889 :
890 : // Now create the sub image
891 635 : std::shared_ptr<ImageInterface<Float> > referenceImage( new SubImage<Float>(image, imslice, True) );
892 : {
893 635 : LatticeLocker lock1(*referenceImage, FileLocker::Write);
894 635 : referenceImage->setMiscInfo(image.miscInfo());
895 635 : referenceImage->setUnits(image.units());
896 :
897 635 : }
898 :
899 : // cout << "Made Ref subimage of shape : " << referenceImage->shape() << endl;
900 :
901 635 : return referenceImage;
902 :
903 635 : }
904 :
905 :
906 :
907 : //////////////////////////////////////////////////////////////////////////////////////////////////////
908 :
909 2022 : SIImageStore::~SIImageStore()
910 : {
911 2022 : }
912 :
913 : //////////////////////////////////////////////////////////////////////////////////////////////////////
914 : //////////////////////////////////////////////////////////////////////////////////////////////////////
915 :
916 3000 : Bool SIImageStore::releaseLocks()
917 : {
918 : //LogIO os( LogOrigin("SIImageStore","releaseLocks",WHERE) );
919 :
920 : // String fname( itsImageName+String(".info") );
921 : // makePersistent( fname );
922 :
923 3000 : if( itsPsf ) releaseImage( itsPsf );
924 3000 : if( itsModel ) { releaseImage( itsModel ); }
925 3000 : if( itsResidual ) releaseImage( itsResidual );
926 3000 : if( itsImage ) releaseImage( itsImage );
927 3000 : if( itsWeight ) releaseImage( itsWeight );
928 3000 : if( itsMask ) releaseImage( itsMask );
929 3000 : if( itsSumWt ) releaseImage( itsSumWt );
930 3000 : if( itsGridWt ) releaseImage( itsGridWt );
931 3000 : if( itsPB ) releaseImage( itsPB );
932 3000 : if( itsImagePBcor ) releaseImage( itsImagePBcor );
933 :
934 3000 : 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 4181 : 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 4181 : im->flush();
952 : //os << LogIO::WARN << "clear cache" << LogIO::POST;
953 4181 : im->clearCache();
954 : //os << LogIO::WARN << "unlock" << LogIO::POST;
955 4181 : im->unlock();
956 : //os << LogIO::WARN << "tempClose" << LogIO::POST;
957 : //im->tempClose();
958 : //os << LogIO::WARN << "NULL" << LogIO::POST;
959 4181 : im.reset(); // This was added to allow modification by modules independently
960 4181 : }
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 73 : Long SIImageStore::estimateRAM(){
977 : //right now this is estimated at 2MB for the 2 complex lattices;
978 73 : return Long(2000);
979 : }
980 0 : void SIImageStore::setModelImage( const Vector<String> &modelnames)
981 : {
982 0 : LogIO os( LogOrigin("SIImageStore","setModelImage",WHERE) );
983 :
984 0 : 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 0 : if( modelnames.nelements()==1 ) { setModelImageOne( modelnames[0] ); }
990 0 : }
991 :
992 :
993 :
994 0 : void SIImageStore::setModelImageOne( const String &modelname , Int nterm)
995 : {
996 0 : LogIO os( LogOrigin("SIImageStore","setModelImageOne",WHERE) );
997 :
998 0 : if(modelname==String("")) return;
999 :
1000 0 : Bool multiterm=False;
1001 0 : if(nterm>-1)multiterm=True;
1002 0 : if(nterm==-1)nterm=0;
1003 :
1004 0 : Directory immodel( modelname ); // +String(".model") );
1005 0 : 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 0 : std::shared_ptr<ImageInterface<Float> > newmodel;
1015 0 : buildImage(newmodel, modelname);
1016 : // in master
1017 : //std::shared_ptr<PagedImage<Float> > newmodel( new PagedImage<Float>( modelname ) ); //+String(".model") ) );
1018 :
1019 0 : Bool hasMask = newmodel->isMasked(); /// || newmodel->hasPixelMask() ;
1020 :
1021 0 : if( hasMask )
1022 : {
1023 :
1024 0 : 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 0 : TempImage<Float> maskmodel( newmodel->shape(), newmodel->coordinates() );
1036 0 : IPosition inshape = newmodel->shape();
1037 0 : for(Int pol=0; pol<inshape[2]; pol++)
1038 : {
1039 0 : for(Int chan=0; chan<inshape[3]; chan++)
1040 : {
1041 0 : IPosition pos(4,0,0,pol,chan);
1042 : std::shared_ptr<ImageInterface<Float> > subim=makeSubImage(0,1,
1043 0 : chan, inshape[3],
1044 0 : pol, inshape[2],
1045 0 : (*newmodel) );
1046 :
1047 : std::shared_ptr<ImageInterface<Float> > submaskmodel=makeSubImage(0,1,
1048 0 : chan, inshape[3],
1049 0 : pol, inshape[2],
1050 0 : maskmodel );
1051 :
1052 0 : ArrayLattice<Bool> pixmask( subim->getMask() );
1053 0 : LatticeExpr<Float> masked( (*subim) * iif(pixmask,1.0,0.0) );
1054 0 : submaskmodel->copyData( masked );
1055 0 : }
1056 : }
1057 :
1058 :
1059 :
1060 : // Check shapes, coordsys with those of other images. If different, try to re-grid here.
1061 0 : 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 0 : os << "Copying input model " << modelname << " to " << itsImageName << ".model" << ((multiterm)?".tt"+String::toString(nterm) :"") << LogIO::POST;
1069 0 : model(nterm)->copyData( LatticeExpr<Float> (maskmodel) );
1070 : }
1071 :
1072 :
1073 0 : } 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 0 : if( (newmodel->shape() != model(nterm)->shape()) || (! itsCoordSys.near(newmodel->coordinates() )) )
1083 : {
1084 0 : os << "Regridding input model " << modelname << " to target coordinate system for " << itsImageName << ".model" << ((multiterm)?".tt"+String::toString(nterm) :"") << LogIO::POST;
1085 0 : regridToModelImage( *newmodel, nterm );
1086 : }
1087 : else
1088 : {
1089 0 : os << "Copying input model " << modelname << " to " << itsImageName << ".model" << ((multiterm)?".tt"+String::toString(nterm) :"") << LogIO::POST;
1090 0 : model(nterm)->copyData( LatticeExpr<Float> (*newmodel) );
1091 : }
1092 : }//nomask
1093 0 : }
1094 :
1095 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1096 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1097 :
1098 1607 : IPosition SIImageStore::getShape()
1099 : {
1100 1607 : return itsImageShape;
1101 : }
1102 :
1103 759 : String SIImageStore::getName()
1104 : {
1105 759 : return itsImageName;
1106 : }
1107 :
1108 3571 : String SIImageStore::imageFullName(IMAGE_IDS imageId)
1109 : {
1110 3571 : return itsImageName + imageExts(imageId);
1111 : }
1112 :
1113 1111 : uInt SIImageStore::getNTaylorTerms(Bool /*dopsf*/)
1114 : {
1115 1111 : 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 20836 : 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 20836 : Bool sw=False;
1135 20836 : if( label.contains(imageExts(SUMWT)) ) sw=True;
1136 :
1137 20836 : if( !ptr )
1138 : {
1139 : //cout << itsNFacets << " " << itsNChanChunks << " " << itsNPolChunks << endl;
1140 5346 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
1141 : {
1142 0 : if( ! parentptr )
1143 : {
1144 : //cout << "Making parent : " << itsImageName+label << " sw : " << sw << endl;
1145 0 : parentptr = openImage(itsImageName+label , itsOverWrite, sw, itsNFacets );
1146 0 : 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 0 : ptr = makeSubImage( itsFacetId, itsNFacets*itsNFacets,
1152 : itsChanId, itsNChanChunks,
1153 : itsPolId, itsNPolChunks,
1154 0 : *parentptr );
1155 0 : if( ! sw )
1156 : {
1157 0 : itsImageShape = ptr->shape(); // over and over again.... FIX.
1158 0 : itsCoordSys = ptr->coordinates();
1159 0 : itsMiscInfo = ptr->miscInfo();
1160 : }
1161 :
1162 : //cout << "accessImage : " << label << " : sumwt : " << sw << " : shape : " << itsImageShape << endl;
1163 :
1164 : }
1165 : else
1166 : { //no facets of chanchunks
1167 5346 : if(!parentptr){
1168 : ///coordsys for psf can be different ...shape should be the same.
1169 5073 : ptr = openImage(itsImageName+label , itsOverWrite, sw, 1, !(label.contains(imageExts(PSF)))); }
1170 : else{
1171 273 : 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 20836 : }
1184 :
1185 :
1186 2546 : std::shared_ptr<ImageInterface<Float> > SIImageStore::psf(uInt /*nterm*/)
1187 : {
1188 2546 : accessImage( itsPsf, itsParentPsf, imageExts(PSF) );
1189 :
1190 2546 : return itsPsf;
1191 : }
1192 8115 : std::shared_ptr<ImageInterface<Float> > SIImageStore::residual(uInt /*nterm*/)
1193 : {
1194 8115 : accessImage( itsResidual, itsParentResidual, imageExts(RESIDUAL) );
1195 : // Record mi = itsResidual->miscInfo(); ostringstream oss;mi.print(oss);cout<<"MiscInfo(res) : " << oss.str() << endl;
1196 8115 : return itsResidual;
1197 : }
1198 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::weight(uInt /*nterm*/)
1199 : {
1200 0 : accessImage( itsWeight, itsParentWeight, imageExts(WEIGHT) );
1201 0 : return itsWeight;
1202 : }
1203 :
1204 2187 : std::shared_ptr<ImageInterface<Float> > SIImageStore::sumwt(uInt /*nterm*/)
1205 : {
1206 :
1207 2187 : accessImage( itsSumWt, itsParentSumWt, imageExts(SUMWT) );
1208 :
1209 2187 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
1210 0 : { itsUseWeight = getUseWeightImage( *itsParentSumWt );}
1211 2187 : setUseWeightImage( *itsSumWt , itsUseWeight); // Sets a flag in the SumWt image.
1212 :
1213 2187 : return itsSumWt;
1214 : }
1215 :
1216 3016 : std::shared_ptr<ImageInterface<Float> > SIImageStore::model(uInt /*nterm*/)
1217 : {
1218 3016 : accessImage( itsModel, itsParentModel, imageExts(MODEL) );
1219 3016 : LatticeLocker lock1(*itsModel, FileLocker::Write);
1220 : // Set up header info the first time.
1221 3016 : itsModel->setUnits("Jy/pixel");
1222 :
1223 6032 : return itsModel;
1224 3016 : }
1225 438 : std::shared_ptr<ImageInterface<Float> > SIImageStore::image(uInt /*nterm*/)
1226 : {
1227 438 : accessImage( itsImage, itsParentImage, imageExts(IMAGE) );
1228 438 : LatticeLocker lock1(*itsImage, FileLocker::Write);
1229 438 : itsImage->setUnits(Unit("Jy/beam"));
1230 876 : return itsImage;
1231 438 : }
1232 :
1233 3360 : std::shared_ptr<ImageInterface<Float> > SIImageStore::mask(uInt /*nterm*/)
1234 : {
1235 3360 : accessImage( itsMask, itsParentMask, imageExts(MASK) );
1236 3360 : return itsMask;
1237 : }
1238 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::gridwt(IPosition newshape)
1239 :
1240 : {
1241 0 : if(newshape.empty()){
1242 0 : accessImage( itsGridWt, itsParentGridWt, imageExts(GRIDWT) );
1243 : }
1244 : else{
1245 0 : if(!itsGridWt || (itsGridWt && (itsGridWt->shape() != newshape))){
1246 0 : itsGridWt.reset(); // deassign previous one hopefully it'll close it
1247 0 : CoordinateSystem newcoordsys=itsCoordSys;
1248 0 : if(newshape.nelements() > 4){
1249 0 : Matrix<Double> pc(1,1); pc = 1.0;
1250 0 : LinearCoordinate zc(Vector<String>(1, "FIELD_ORDER"), Vector<String>(1, ""), Vector<Double>(1, 0.0), Vector<Double>(1,1.0), pc, Vector<Double>(1,0.0));
1251 0 : newcoordsys.addCoordinate(zc);
1252 0 : }
1253 0 : itsGridWt.reset(new PagedImage<Float>(newshape, newcoordsys, itsImageName+ imageExts(GRIDWT)));
1254 0 : initMetaInfo(itsGridWt, itsImageName+ imageExts(GRIDWT));
1255 :
1256 0 : }
1257 : }
1258 : /// change the coordinate system here, to uv.
1259 0 : return itsGridWt;
1260 : }
1261 :
1262 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::tempworkimage(uInt /*term*/){
1263 0 : if(itsTempWorkIm) return itsTempWorkIm;
1264 0 : itsTempWorkIm.reset(new PagedImage<Float>(itsImageShape, itsCoordSys, itsImageName+ ".work.temp"));
1265 0 : static_cast<PagedImage<Float>* > (itsTempWorkIm.get())->set(0.0);
1266 0 : itsTempWorkIm->flush();
1267 0 : static_cast<PagedImage<Float>* > (itsTempWorkIm.get())->table().markForDelete();
1268 0 : return itsTempWorkIm;
1269 : }
1270 :
1271 1174 : std::shared_ptr<ImageInterface<Float> > SIImageStore::pb(uInt /*nterm*/)
1272 : {
1273 1174 : accessImage( itsPB, itsParentPB, imageExts(PB) );
1274 1174 : return itsPB;
1275 : }
1276 0 : std::shared_ptr<ImageInterface<Float> > SIImageStore::imagepbcor(uInt /*nterm*/)
1277 : {
1278 0 : accessImage( itsImagePBcor, itsParentImagePBcor, imageExts(IMAGEPBCOR) );
1279 0 : LatticeLocker lock1(*itsImagePBcor, FileLocker::Write);
1280 0 : itsImagePBcor->setUnits(Unit("Jy/beam"));
1281 0 : return itsImagePBcor;
1282 0 : }
1283 :
1284 1092 : std::shared_ptr<ImageInterface<Complex> > SIImageStore::forwardGrid(uInt /*nterm*/){
1285 1092 : if( itsForwardGrid ) // && (itsForwardGrid->shape() == itsImageShape))
1286 : {
1287 : // cout << "Forward grid has shape : " << itsForwardGrid->shape() << endl;
1288 1022 : return itsForwardGrid;
1289 : }
1290 70 : Vector<Int> whichStokes(0);
1291 70 : IPosition cimageShape;
1292 70 : cimageShape=itsImageShape;
1293 70 : 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 70 : if(freqframe != MFrequency::LSRK && freqframe!=MFrequency::Undefined && freqframe!=MFrequency::REST)
1296 0 : { itsCoordSys.setSpectralConversion("LSRK"); }
1297 70 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
1298 70 : whichStokes, itsDataPolRep);
1299 70 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
1300 70 : cimageShape(2)=whichStokes.nelements();
1301 : //cout << "Making forward grid of shape : " << cimageShape << " for imshape : " << itsImageShape << endl;
1302 70 : itsForwardGrid.reset( new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()) );
1303 : //if(image())
1304 70 : if(hasRestored()){
1305 0 : LatticeLocker lock1(*itsForwardGrid, FileLocker::Write);
1306 0 : itsForwardGrid->setImageInfo((image())->imageInfo());
1307 :
1308 0 : }
1309 70 : return itsForwardGrid;
1310 70 : }
1311 :
1312 838 : std::shared_ptr<ImageInterface<Complex> > SIImageStore::backwardGrid(uInt /*nterm*/){
1313 838 : if( itsBackwardGrid ) //&& (itsBackwardGrid->shape() == itsImageShape))
1314 : {
1315 : // cout << "Backward grid has shape : " << itsBackwardGrid->shape() << endl;
1316 765 : return itsBackwardGrid;
1317 : }
1318 73 : Vector<Int> whichStokes(0);
1319 73 : IPosition cimageShape;
1320 73 : cimageShape=itsImageShape;
1321 73 : 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 73 : if(freqframe != MFrequency::LSRK && freqframe!=MFrequency::Undefined && freqframe!=MFrequency::REST )
1324 0 : { itsCoordSys.setSpectralConversion("LSRK"); }
1325 73 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
1326 73 : whichStokes, itsDataPolRep);
1327 73 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
1328 73 : cimageShape(2)=whichStokes.nelements();
1329 : //cout << "Making backward grid of shape : " << cimageShape << " for imshape : " << itsImageShape << endl;
1330 73 : itsBackwardGrid.reset( new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()) );
1331 : //if(image())
1332 73 : if(hasRestored()){
1333 0 : LatticeLocker lock1(*itsBackwardGrid, FileLocker::Write);
1334 0 : itsBackwardGrid->setImageInfo((image())->imageInfo());
1335 0 : }
1336 73 : return itsBackwardGrid;
1337 73 : }
1338 143 : Double SIImageStore::memoryBeforeLattice(){
1339 : //Calculate how much memory to use per temporary images before disking
1340 143 : return 1.0; /// in MB
1341 : }
1342 143 : IPosition SIImageStore::tileShape(){
1343 : //Need to have settable stuff here or algorith to determine this
1344 143 : 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 9486 : Bool SIImageStore::doesImageExist(String imagename)
1349 : {
1350 18972 : LogIO os( LogOrigin("SIImageStore","doesImageExist",WHERE) );
1351 9486 : Directory image( imagename );
1352 18972 : return image.exists();
1353 9486 : }
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 0 : 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 0 : LatticeExprNode le( sqrt(max( *(weight(0)) )) );
1449 0 : return le.getFloat();
1450 :
1451 0 : }
1452 :
1453 0 : 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 0 : CountedPtr<ImageInterface<Float> > subim=makeSubImage(0, 1,
1459 0 : chan, itsImageShape[3],
1460 0 : pol, itsImageShape[2],
1461 0 : *weight(0));
1462 :
1463 0 : return sqrt(max(subim->get()));
1464 0 : }
1465 :
1466 0 : void SIImageStore::makePBFromWeight(const Float pblimit)
1467 : {
1468 0 : LogIO os( LogOrigin("SIImageStore","makePBFromWeight",WHERE) );
1469 :
1470 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1471 : {
1472 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1473 : {
1474 :
1475 0 : itsPBScaleFactor = getPbMax(pol,chan);
1476 :
1477 0 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1478 : else {
1479 :
1480 0 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1481 0 : chan, itsImageShape[3],
1482 0 : pol, itsImageShape[2],
1483 0 : *weight() );
1484 0 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
1485 0 : chan, itsImageShape[3],
1486 0 : pol, itsImageShape[2],
1487 0 : *pb() );
1488 :
1489 :
1490 0 : LatticeExpr<Float> normed( sqrt(abs(*wtsubim)) / itsPBScaleFactor );
1491 0 : LatticeExpr<Float> limited( iif( normed > fabs(pblimit) , normed, 0.0 ) );
1492 0 : pbsubim->copyData( limited );
1493 0 : }// if not zero
1494 : }
1495 : }
1496 :
1497 0 : 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 0 : LatticeExpr<Bool> pbmask( iif( *pb() > fabs(pblimit) , True , False ) );
1506 : //MSK//
1507 0 : createMask( pbmask, pb() );
1508 0 : pb()->pixelMask().unlock();
1509 0 : }
1510 :
1511 : }
1512 0 : weight()->unlock();
1513 0 : pb()->unlock();
1514 0 : }
1515 :
1516 73 : void SIImageStore::makePBImage(const Float pblimit)
1517 : {
1518 146 : LogIO os( LogOrigin("SIImageStore","makePBImage",WHERE) );
1519 :
1520 73 : if( hasPB() )
1521 : {
1522 :
1523 146 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1524 : {
1525 146 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1526 : {
1527 :
1528 73 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
1529 73 : chan, itsImageShape[3],
1530 73 : pol, itsImageShape[2],
1531 146 : *pb() );
1532 :
1533 : // Norm by the max
1534 73 : LatticeExprNode elmax= max( *pbsubim );
1535 73 : 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 73 : if(fmax>1.0)
1539 : {
1540 2 : LatticeExpr<Float> normed( (*pbsubim) / fmax );
1541 2 : LatticeExpr<Float> limited( iif( normed > fabs(pblimit) , normed, 0.0 ) );
1542 1 : pbsubim->copyData( limited );
1543 1 : }
1544 : else
1545 : {
1546 144 : LatticeExpr<Float> limited( iif((*pbsubim) > fabs(pblimit) , (*pbsubim), 0.0 ) );
1547 72 : pbsubim->copyData( limited );
1548 72 : }
1549 73 : }
1550 : }
1551 :
1552 73 : if((pb()->getDefaultMask()==""))// && pblimit >= 0.0)
1553 : {
1554 : // removeMask( pb() );
1555 : //MSK//
1556 146 : LatticeExpr<Bool> pbmask( iif( *pb() > fabs(pblimit) , True , False ) );
1557 : //MSK//
1558 73 : createMask( pbmask, pb() );
1559 73 : pb()->pixelMask().unlock();
1560 73 : }
1561 :
1562 : }// if hasPB
1563 73 : pb()->unlock();
1564 :
1565 73 : }
1566 :
1567 73 : 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 73 : LatticeLocker lock1(*outimage, FileLocker::Write);
1574 73 : ImageRegion outreg = outimage->makeMask("mask0",False,True);
1575 73 : LCRegion& outmask=outreg.asMask();
1576 73 : outmask.copyData(lemask);
1577 73 : outimage->defineRegion("mask0",outreg, RegionHandler::Masks, True);
1578 :
1579 73 : outimage->setDefaultMask("mask0");
1580 :
1581 73 : outimage->unlock();
1582 73 : outimage->tempClose();
1583 :
1584 : // outimage->table().unmarkForDelete();
1585 73 : }
1586 0 : catch (const AipsError& x) {
1587 0 : throw(AipsError("Error in creating internal T/F mask : " + x.getMesg() ));
1588 0 : }
1589 :
1590 73 : return True;
1591 : }
1592 :
1593 73 : 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 73 : if( (inimage->getDefaultMask()).matches("mask0") ) // input mask exists.
1601 0 : {LatticeLocker lock1(*outimage, FileLocker::Write);
1602 0 : 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 0 : ImageRegion outreg=outimage->makeMask("mask0", False, True);
1611 0 : LCRegion& outmask=outreg.asMask();
1612 0 : outmask.copyData(inimage->getRegion("mask0").asLCRegion());
1613 0 : outimage->defineRegion("mask0",outreg, RegionHandler::Masks,True);
1614 0 : outimage->setDefaultMask("mask0");
1615 0 : }
1616 : }
1617 0 : catch (const AipsError& x) {
1618 0 : throw(AipsError("Error in copying internal T/F mask : " + x.getMesg() ));
1619 0 : }
1620 :
1621 73 : return True;
1622 : }
1623 :
1624 0 : 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 0 : LatticeLocker lock1(*im, FileLocker::Write);
1634 0 : if (im-> getDefaultMask() != String(""))
1635 : {
1636 0 : String strung=im->getDefaultMask();
1637 0 : im->setDefaultMask("");
1638 0 : im->removeRegion(strung);
1639 0 : }
1640 0 : if( im->hasRegion("mask0") )
1641 : {
1642 0 : im->removeRegion("mask0");
1643 : }
1644 0 : }
1645 0 : catch (const AipsError& x)
1646 : {
1647 0 : throw(AipsError("Error in deleting internal T/F mask : " + x.getMesg() ));
1648 0 : }
1649 0 : }
1650 0 : void SIImageStore:: rescaleResolution(Int chan,
1651 : ImageInterface<Float>& image,
1652 : const GaussianBeam& newbeam,
1653 : const GaussianBeam& oldbeam){
1654 :
1655 0 : LogIO os( LogOrigin("SIImageStore","rescaleResolution",WHERE) );
1656 0 : GaussianBeam toBeUsed(Quantity(0.0, "arcsec"),Quantity(0.0, "arcsec"),
1657 0 : Quantity(0.0, "deg")) ;
1658 : try {
1659 0 : 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 0 : if( samesize )
1669 : {
1670 0 : os << LogIO::NORMAL2 << "Input and output beam sizes are the same for Channel : " << chan << ". Not convolving residuals." << LogIO::POST;
1671 : }
1672 : else
1673 : {
1674 0 : Double pixwidth=sqrt(image.coordinates().increment()(0)*image.coordinates().increment()(0)+image.coordinates().increment()(1)*image.coordinates().increment()(1));
1675 :
1676 0 : if(toBeUsed.getMinor(image.coordinates().worldAxisUnits()[0]) > pixwidth)
1677 : {
1678 : //cerr << "old beam area " << oldbeam.getArea("rad2") << " new beam " << newbeam.getArea("rad2") << endl;
1679 0 : StokesImageUtil::Convolve(image, toBeUsed, True);
1680 0 : image.copyData(LatticeExpr<Float>(image*newbeam.getArea("rad2")/ oldbeam.getArea("rad2")));
1681 : }
1682 : }
1683 : }
1684 0 : catch (const AipsError& x) {
1685 : //throw(AipsError("Cannot convolve to new beam: may be smaller than old beam : " + x.getMesg() ));
1686 0 : os << LogIO::WARN << "Cannot convolve to new beam for Channel : " << chan << " : " << x.getMesg() << LogIO::POST;
1687 0 : }
1688 :
1689 :
1690 0 : }
1691 :
1692 73 : void SIImageStore::divideWeightBySumwt()
1693 : {
1694 146 : LogIO os(LogOrigin("SIImageStore", "divideWeightBySumWt", WHERE));
1695 :
1696 :
1697 73 : if( itsUseWeight )
1698 : {
1699 :
1700 0 : divideImageByWeightVal( *weight() );
1701 : }
1702 :
1703 73 : }
1704 :
1705 :
1706 73 : void SIImageStore::dividePSFByWeight(const Float /* pblimit*/)
1707 : {
1708 146 : LogIO os( LogOrigin("SIImageStore","dividePSFByWeight",WHERE) );
1709 :
1710 73 : LatticeLocker lock1 (*(psf()), FileLocker::Write);
1711 73 : 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 73 : (psf())->unlock();
1721 :
1722 73 : }
1723 :
1724 73 : void SIImageStore::normalizePrimaryBeam(const Float pblimit)
1725 : {
1726 146 : LogIO os( LogOrigin("SIImageStore","normalizePrimaryBeam",WHERE) );
1727 :
1728 : // cout << "In dividePSFByWeight : itsUseWeight : " << itsUseWeight << endl;
1729 73 : if( itsUseWeight )
1730 : {
1731 :
1732 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1733 : {
1734 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1735 : {
1736 0 : os << LogIO::NORMAL1 << "Scale factor for [C" +String::toString(chan) + ":P" + String::toString(pol) + "] to keep the model image w.r.to a PB of max=1 is " << getPbMax(pol,chan) << LogIO::POST;
1737 : }//chan
1738 : }//pol
1739 :
1740 0 : makePBFromWeight(pblimit);
1741 :
1742 : }//if itsUseWeight
1743 73 : else { makePBImage(pblimit); } // OR... just check that it exists already.
1744 73 : pb()->unlock();
1745 :
1746 73 : }
1747 :
1748 : // Make another for the PSF too.
1749 346 : void SIImageStore::divideResidualByWeight(Float pblimit, String normtype) {
1750 :
1751 692 : LogIO os(LogOrigin("SIImageStore", "divideResidualByWeight", WHERE));
1752 346 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1753 :
1754 0 : auto logTemplate = [&](Int const &chan, Int const &pol, string const &normalizer, string const &result) {
1755 0 : os << LogIO::NORMAL1
1756 0 : << "[C" + String::toString(chan) + ":P" + String::toString(pol) + "] "
1757 0 : << "Dividing " << itsImageName + String(".residual") << " by "
1758 : << "[ " << normalizer << " ] "
1759 0 : << "to get " << result << "." << LogIO::POST;
1760 0 : };
1761 :
1762 : // Normalize by the sumwt, per plane.
1763 346 : Bool didNorm = divideImageByWeightVal(*residual());
1764 :
1765 346 : if (itsUseWeight) {
1766 0 : for (Int pol = 0; pol < itsImageShape[2]; pol++) {
1767 0 : for (Int chan = 0; chan < itsImageShape[3]; chan++) {
1768 0 : itsPBScaleFactor = getPbMax(pol, chan);
1769 :
1770 0 : if (itsPBScaleFactor <= 0) {
1771 : os << LogIO::NORMAL1
1772 : << "Skipping normalization for C:" << chan << " P:" << pol
1773 0 : << " because pb max is zero " << LogIO::POST;
1774 :
1775 : } else {
1776 0 : CountedPtr<ImageInterface<Float> > wtsubim = makeSubImage(0, 1,
1777 0 : chan, itsImageShape[3],
1778 0 : pol, itsImageShape[2],
1779 0 : *weight());
1780 0 : CountedPtr<ImageInterface<Float> > ressubim = makeSubImage(0, 1,
1781 0 : chan, itsImageShape[3],
1782 0 : pol, itsImageShape[2],
1783 0 : *residual());
1784 0 : LatticeExpr<Float> ratio;
1785 0 : Float scalepb = 1.0;
1786 :
1787 0 : if (normtype == "flatnoise") {
1788 0 : logTemplate(chan, pol,
1789 0 : "sqrt(weightimage) * " + String::toString(itsPBScaleFactor),
1790 : "flat noise with unit pb peak");
1791 :
1792 0 : LatticeExpr<Float> deno = itsPBScaleFactor * sqrt(abs(LatticeExpr<Float>(*(wtsubim))));
1793 0 : scalepb = fabs(pblimit) * itsPBScaleFactor * itsPBScaleFactor;
1794 0 : ratio = iif(deno > scalepb, (*(ressubim) / deno), 0.0);
1795 :
1796 0 : } 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 0 : 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 0 : } // 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 346 : 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 346 : if ((residual()->getDefaultMask() == "") && hasPB() && pblimit >=0.0) {
1840 0 : copyMask(pb(), residual());
1841 : }
1842 :
1843 346 : if ((pblimit < 0.0) && (residual()->getDefaultMask()).matches("mask0")) {
1844 0 : removeMask(residual());
1845 : }
1846 :
1847 346 : residual()->unlock();
1848 346 : }
1849 :
1850 0 : void SIImageStore::divideResidualByWeightSD(Float pblimit) {
1851 :
1852 0 : LogIO os(LogOrigin("SIImageStore", "divideResidualByWeightSD", WHERE));
1853 0 : LatticeLocker lock1(*(residual()), FileLocker::Write);
1854 :
1855 0 : if (itsUseWeight) {
1856 0 : LatticeExpr<Float> deno = LatticeExpr<Float>(*weight());
1857 0 : LatticeExpr<Float> ratio = iif(deno > 0.0, *(residual()) / deno, 0.0);
1858 0 : residual()->copyData(ratio);
1859 0 : }
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 0 : if ((residual()->getDefaultMask() == "") && hasPB() && pblimit >=0.0) {
1869 0 : copyMask(pb(), residual());
1870 : }
1871 :
1872 0 : if ((pblimit < 0.0) && (residual()->getDefaultMask()).matches("mask0")) {
1873 0 : removeMask(residual());
1874 : }
1875 :
1876 0 : residual()->unlock();
1877 0 : }
1878 :
1879 346 : void SIImageStore::divideModelByWeight(Float pblimit, const String normtype)
1880 : {
1881 692 : LogIO os( LogOrigin("SIImageStore","divideModelByWeight",WHERE) );
1882 :
1883 :
1884 692 : if(itsUseWeight // only when needed
1885 346 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
1886 : {
1887 :
1888 0 : 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 0 : else if( normtype=="flatnoise"){
1897 :
1898 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
1899 : {
1900 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
1901 : {
1902 :
1903 0 : itsPBScaleFactor = getPbMax(pol,chan);
1904 : // cout << " pbscale : " << itsPBScaleFactor << endl;
1905 0 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
1906 : else {
1907 :
1908 0 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
1909 0 : chan, itsImageShape[3],
1910 0 : pol, itsImageShape[2],
1911 0 : *weight() );
1912 0 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
1913 0 : chan, itsImageShape[3],
1914 0 : pol, itsImageShape[2],
1915 0 : *model() );
1916 0 : os << LogIO::NORMAL1 ;
1917 0 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
1918 0 : os << "Dividing " << itsImageName+String(".model") ;
1919 0 : os << " by [ sqrt(weight) / " << itsPBScaleFactor ;
1920 0 : os <<" ] to get to flat sky model before prediction" << LogIO::POST;
1921 :
1922 :
1923 0 : LatticeExpr<Float> deno( sqrt( abs(*(wtsubim)) ) / itsPBScaleFactor );
1924 :
1925 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1926 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1927 0 : 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 0 : 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 0 : 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 0 : }// 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 346 : }
2011 :
2012 346 : void SIImageStore::multiplyModelByWeight(Float pblimit, const String normtype)
2013 : {
2014 692 : LogIO os( LogOrigin("SIImageStore","multiplyModelByWeight",WHERE) );
2015 :
2016 692 : if(itsUseWeight // only when needed
2017 346 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
2018 : {
2019 0 : 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 0 : else if( normtype=="flatnoise"){
2024 :
2025 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2026 : {
2027 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
2028 : {
2029 :
2030 0 : itsPBScaleFactor = getPbMax(pol,chan);
2031 : // cout << " pbscale : " << itsPBScaleFactor << endl;
2032 0 : if(itsPBScaleFactor<=0){os << LogIO::NORMAL1 << "Skipping normalization for C:" << chan << " P:" << pol << " because pb max is zero " << LogIO::POST;}
2033 : else {
2034 :
2035 0 : CountedPtr<ImageInterface<Float> > wtsubim=makeSubImage(0,1,
2036 0 : chan, itsImageShape[3],
2037 0 : pol, itsImageShape[2],
2038 0 : *weight() );
2039 0 : CountedPtr<ImageInterface<Float> > modsubim=makeSubImage(0,1,
2040 0 : chan, itsImageShape[3],
2041 0 : pol, itsImageShape[2],
2042 0 : *model() );
2043 :
2044 :
2045 :
2046 0 : os << LogIO::NORMAL1 ;
2047 0 : os << "[C" +String::toString(chan) + ":P" + String::toString(pol) + "] ";
2048 0 : os << "Multiplying " << itsImageName+String(".model") ;
2049 0 : os << " by [ sqrt(weight) / " << itsPBScaleFactor;
2050 0 : os << " ] to take model back to flat noise with unit pb peak." << LogIO::POST;
2051 :
2052 0 : LatticeExpr<Float> deno( sqrt( abs(*(wtsubim)) ) / itsPBScaleFactor );
2053 :
2054 0 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
2055 0 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
2056 0 : 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 0 : modsubim->copyData(ratio);
2061 :
2062 0 : }// 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 346 : }
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 143 : void SIImageStore::makeImageBeamSet(Float psfcutoff, const Bool forcefit)
2133 : {
2134 143 : clock_t begin = clock();
2135 :
2136 286 : LogIO os( LogOrigin("SIImageStore","getPSFGaussian",WHERE) );
2137 : // For all chans/pols, call getPSFGaussian() and put it into ImageBeamSet(chan,pol).
2138 143 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2139 143 : uInt nx = itsImageShape[0];
2140 143 : uInt ny = itsImageShape[1];
2141 143 : uInt npol = itsImageShape[2];
2142 143 : uInt nchan = itsImageShape[3];
2143 143 : ImageInfo ii = psf()->imageInfo();
2144 143 : ImageBeamSet iibeamset=ii.getBeamSet();
2145 143 : if(iibeamset.nchan()==nchan && iibeamset.nstokes()==npol && forcefit==False){
2146 70 : itsPSFBeams=iibeamset;
2147 70 : itsRestoredBeams=iibeamset;
2148 70 : return;
2149 : }
2150 73 : itsPSFBeams.resize( nchan, npol );
2151 73 : itsRestoredBeams.resize(nchan, npol);
2152 : // cout << "makeImBeamSet : imshape : " << itsImageShape << endl;
2153 :
2154 73 : String blankpsfs="";
2155 146 : for( uInt chanid=0; chanid<nchan;chanid++) {
2156 146 : for( uInt polid=0; polid<npol; polid++ ) {
2157 73 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
2158 :
2159 73 : IPosition substart(4,0,0,polid,chanid);
2160 73 : IPosition substop(4,nx-1,ny-1,polid,chanid);
2161 73 : Slicer psfslice(substart, substop,Slicer::endIsLast);
2162 146 : SubImage<Float> subPsf( *psf() , psfslice, True );
2163 73 : GaussianBeam beam;
2164 :
2165 73 : Bool tryfit=True;
2166 :
2167 73 : LatticeExprNode le( max(subPsf) );
2168 73 : tryfit = le.getFloat()>0.0;
2169 73 : if(tryfit)
2170 : {
2171 : try
2172 : {
2173 73 : StokesImageUtil::FitGaussianPSF( subPsf, beam,psfcutoff );
2174 73 : itsPSFBeams.setBeam( chanid, polid, beam );
2175 73 : 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 0 : blankpsfs += "[C" +String::toString(chanid) + ":P" + String::toString(polid) + "] ";
2188 : }
2189 :
2190 73 : }// end of pol loop
2191 : }// end of chan loop
2192 :
2193 73 : if( blankpsfs.length() >0 )
2194 0 : 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 73 : 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 73 : if(defaultbeam.getArea("rad2")==0.0){
2207 0 : Quantity majax(1e-06,"arcsec"),minax(1e-06,"arcsec"),pa(0.0,"deg");
2208 0 : defaultbeam.setMajorMinor(majax,minax);
2209 0 : defaultbeam.setPA(pa);
2210 0 : }
2211 146 : for( uInt chanid=0; chanid<nchan;chanid++) {
2212 146 : for( uInt polid=0; polid<npol; polid++ ) {
2213 73 : if( (itsPSFBeams.getBeam(chanid, polid)).isNull() )
2214 0 : { itsPSFBeams.setBeam( chanid, polid, defaultbeam );
2215 0 : 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 73 : ii.setBeams( itsPSFBeams );
2243 : {
2244 73 : LatticeLocker lock1(*(psf()), FileLocker::Write);
2245 73 : psf()->setImageInfo(ii);
2246 73 : psf()->unlock();
2247 73 : }
2248 73 : clock_t end = clock();
2249 73 : os << LogIO::NORMAL << "Time to fit Gaussian to PSF " << double(end - begin) / CLOCKS_PER_SEC << LogIO::POST;
2250 283 : }// 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 0 : void SIImageStore::setBeamSet(const ImageBeamSet& bs){
2275 :
2276 0 : itsPSFBeams=bs;
2277 0 : }
2278 :
2279 70 : ImageBeamSet SIImageStore::getBeamSet(Float psfcutoff)
2280 : {
2281 70 : IPosition beamshp = itsPSFBeams.shape();
2282 70 : AlwaysAssert( beamshp.nelements()==2 , AipsError );
2283 70 : if( beamshp[0]==0 || beamshp[1]==0 ) {makeImageBeamSet(psfcutoff);}
2284 140 : return itsPSFBeams;
2285 70 : }
2286 :
2287 73 : void SIImageStore::printBeamSet(Bool verbose)
2288 : {
2289 146 : LogIO os( LogOrigin("SIImageStore","printBeamSet",WHERE) );
2290 73 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2291 73 : if( itsImageShape[3] == 1 && itsImageShape[2]==1 )
2292 : {
2293 73 : GaussianBeam beam = itsPSFBeams.getBeam();
2294 73 : os << "Beam : " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2295 73 : }
2296 : else
2297 : {
2298 : // per CAS-11415, verbose=True when restoringbeam='common'
2299 0 : if( verbose )
2300 : {
2301 0 : ostringstream oss;
2302 0 : oss<<"Beam";
2303 0 : Int nchan = itsImageShape[3];
2304 0 : for( Int chanid=0; chanid<nchan;chanid++) {
2305 0 : Int polid=0;
2306 : // for( Int polid=0; polid<npol; polid++ ) {
2307 0 : GaussianBeam beam = itsPSFBeams.getBeam( chanid, polid );
2308 0 : oss << " [C" << chanid << "]: " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg";
2309 0 : }//for chanid
2310 0 : os << oss.str() << LogIO::POST;
2311 0 : }
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 0 : 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 73 : }
2340 :
2341 : /////////////////////////////// Restore all planes.
2342 :
2343 73 : void SIImageStore::restore(GaussianBeam& rbeam, String& usebeam, uInt term, Float psfcutoff)
2344 : {
2345 146 : LogIO os( LogOrigin("SIImageStore","restore",WHERE) );
2346 : // << ". Optionally, PB-correct too." << LogIO::POST;
2347 :
2348 73 : AlwaysAssert( itsImageShape.nelements() == 4, AipsError );
2349 73 : Int nx = itsImageShape[0];
2350 73 : Int ny = itsImageShape[1];
2351 73 : Int npol = itsImageShape[2];
2352 73 : 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 73 : IPosition beamset = itsPSFBeams.shape();
2364 73 : AlwaysAssert( beamset.nelements()==2 , AipsError );
2365 73 : if( beamset[0] != nchan || beamset[1] != npol )
2366 : {
2367 :
2368 : // Get PSF Beams....
2369 3 : ImageInfo ii = psf()->imageInfo();
2370 3 : itsPSFBeams = ii.getBeamSet();
2371 3 : itsRestoredBeams=itsPSFBeams;
2372 3 : IPosition beamset2 = itsPSFBeams.shape();
2373 :
2374 3 : AlwaysAssert( beamset2.nelements()==2 , AipsError );
2375 3 : 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 3 : }
2382 :
2383 : // toggle for printing common beam to the log
2384 73 : 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 73 : if( rbeam.isNull() && usebeam=="common") {
2390 0 : os << "Getting common beam" << LogIO::POST;
2391 0 : rbeam = CasaImageBeamSet(itsPSFBeams).getCommonBeam();
2392 0 : os << "Common Beam : " << rbeam.getMajor(Unit("arcsec")) << " arcsec, " << rbeam.getMinor(Unit("arcsec"))<< " arcsec, " << rbeam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2393 0 : printcommonbeam=True;
2394 : }
2395 73 : 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 0 : itsRestoredBeams=ImageBeamSet(rbeam);
2404 0 : GaussianBeam beam = itsRestoredBeams.getBeam();
2405 :
2406 : //if commonbeam has not printed in the log
2407 0 : if (!printcommonbeam) {
2408 0 : os << "Common Beam : " << beam.getMajor(Unit("arcsec")) << " arcsec, " << beam.getMinor(Unit("arcsec"))<< " arcsec, " << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
2409 0 : printcommonbeam=True; // to avoid duplicate the info is printed to th elog
2410 : }
2411 0 : }// 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 73 : Bool emptymodel = isModelEmpty();
2418 73 : 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 73 : LatticeLocker lock1(*(image(term)), FileLocker::Write);
2422 73 : LatticeLocker lock2(*(residual(term)), FileLocker::Write);
2423 73 : LatticeLocker lock3(*(model(term)), FileLocker::Read);
2424 : //// Use beamset to restore
2425 146 : for( Int chanid=0; chanid<nchan;chanid++) {
2426 146 : for( Int polid=0; polid<npol; polid++ ) {
2427 :
2428 73 : IPosition substart(4,0,0,polid,chanid);
2429 73 : IPosition substop(4,nx-1,ny-1,polid,chanid);
2430 73 : Slicer imslice(substart, substop,Slicer::endIsLast);
2431 146 : SubImage<Float> subRestored( *image(term) , imslice, True );
2432 146 : SubImage<Float> subModel( *model(term) , imslice, False );
2433 146 : SubImage<Float> subResidual( *residual(term) , imslice, True );
2434 :
2435 73 : 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 73 : if(!printcommonbeam) {
2439 73 : 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 73 : subRestored.set(0.0);
2446 73 : if( !emptymodel ) {
2447 : // Copy model into it
2448 70 : subRestored.copyData( LatticeExpr<Float>( subModel ) );
2449 : // Smooth model by beam
2450 70 : StokesImageUtil::Convolve( subRestored, beam);
2451 : }
2452 : // Add residual image
2453 73 : 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 0 : TempImage<Float> tmpSubResidualCopy( IPosition(4,nx,ny,1,1), subResidual.coordinates());
2457 0 : tmpSubResidualCopy.copyData( subResidual );
2458 :
2459 0 : rescaleResolution(chanid, tmpSubResidualCopy, beam, itsPSFBeams.getBeam(chanid, polid));
2460 0 : subRestored.copyData( LatticeExpr<Float>( subRestored + tmpSubResidualCopy ) );
2461 0 : }
2462 : else// if no need to rescale residuals, just add the residuals.
2463 : {
2464 :
2465 :
2466 73 : 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 73 : }// end of pol loop
2477 : }// end of chan loop
2478 :
2479 : try
2480 : {
2481 : //MSK//
2482 73 : if( hasPB() )
2483 : {
2484 73 : if( (image(term)->getDefaultMask()).matches("mask0") ) removeMask( image(term) );
2485 73 : copyMask(residual(term),image(term));
2486 : }
2487 :
2488 : // if(hasPB()){copyMask(residual(term),image(term));}
2489 73 : ImageInfo iminf = image(term)->imageInfo();
2490 : //iminf.setBeams( itsRestoredBeams);
2491 :
2492 73 : os << "Beam Set : Single beam : " << itsRestoredBeams.hasSingleBeam() << " Multi-beam : " << itsRestoredBeams.hasMultiBeam() << LogIO::DEBUG2;
2493 :
2494 73 : iminf.removeRestoringBeam();
2495 :
2496 73 : if( itsRestoredBeams.hasSingleBeam() )
2497 73 : { iminf.setRestoringBeam( itsRestoredBeams.getBeam() );}
2498 : else
2499 0 : {iminf.setBeams( itsRestoredBeams);}
2500 :
2501 73 : image(term)->setImageInfo(iminf);
2502 :
2503 73 : }
2504 0 : catch(AipsError &x)
2505 : {
2506 0 : throw( AipsError("Restoration Error : " + x.getMesg() ) );
2507 0 : }
2508 :
2509 73 : }// 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 0 : void SIImageStore::pbcor(uInt term)
2752 : {
2753 :
2754 0 : LogIO os( LogOrigin("SIImageStore","pbcor",WHERE) );
2755 :
2756 0 : 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 0 : LatticeLocker lock1(*(image(term)), FileLocker::Write);
2763 0 : LatticeLocker lock2(*(pb(term)), FileLocker::Write);
2764 0 : LatticeLocker lock3(*(imagepbcor(term)), FileLocker::Write);
2765 :
2766 0 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2767 : {
2768 0 : for(Int chan=0; chan<itsImageShape[3]; chan++)
2769 : {
2770 :
2771 0 : CountedPtr<ImageInterface<Float> > restoredsubim=makeSubImage(0,1,
2772 0 : chan, itsImageShape[3],
2773 0 : pol, itsImageShape[2],
2774 0 : *image(term) );
2775 0 : CountedPtr<ImageInterface<Float> > pbsubim=makeSubImage(0,1,
2776 0 : chan, itsImageShape[3],
2777 0 : pol, itsImageShape[2],
2778 0 : *pb() );
2779 :
2780 0 : CountedPtr<ImageInterface<Float> > pbcorsubim=makeSubImage(0,1,
2781 0 : chan, itsImageShape[3],
2782 0 : pol, itsImageShape[2],
2783 0 : *imagepbcor(term) );
2784 :
2785 :
2786 0 : LatticeExprNode pbmax( max( *pbsubim ) );
2787 0 : Float pbmaxval = pbmax.getFloat();
2788 :
2789 0 : 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 0 : LatticeExpr<Float> thepbcor( iif( *(pbsubim) > 0.0 , (*(restoredsubim))/(*(pbsubim)) , 0.0 ) );
2798 0 : pbcorsubim->copyData( thepbcor );
2799 :
2800 0 : }// if not zero
2801 0 : }//chan
2802 : }//pol
2803 :
2804 : // Copy over the PB mask.
2805 0 : if((imagepbcor(term)->getDefaultMask()=="") && hasPB())
2806 0 : {copyMask(pb(),imagepbcor(term));}
2807 :
2808 : // Set restoring beam info
2809 0 : ImageInfo iminf = image(term)->imageInfo();
2810 : //iminf.setBeams( itsRestoredBeams );
2811 0 : imagepbcor(term)->setImageInfo(iminf);
2812 :
2813 0 : }// 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 838 : Bool SIImageStore::getUseWeightImage(ImageInterface<Float>& target)
2839 : {
2840 838 : Record miscinfo = target.miscInfo();
2841 : Bool useweightimage;
2842 838 : if( miscinfo.isDefined("useweightimage") && miscinfo.dataType("useweightimage")==TpBool )
2843 838 : { miscinfo.get( "useweightimage", useweightimage ); }
2844 0 : else { useweightimage = False; }
2845 :
2846 838 : return useweightimage;
2847 838 : }
2848 : /*
2849 : Bool SIImageStore::getUseWeightImage()
2850 : {
2851 : if( ! itsParentSumWt )
2852 : return False;
2853 : else
2854 : return getUseWeightImage( *itsParentSumWt );
2855 : }
2856 : */
2857 2187 : void SIImageStore::setUseWeightImage(ImageInterface<Float>& target, Bool useweightimage)
2858 : {
2859 2187 : Record miscinfo = target.miscInfo();
2860 2187 : miscinfo.define("useweightimage", useweightimage);
2861 2187 : LatticeLocker lock1(target, FileLocker::Write);
2862 2187 : target.setMiscInfo( miscinfo );
2863 2187 : }
2864 :
2865 :
2866 :
2867 346 : Bool SIImageStore::divideImageByWeightVal( ImageInterface<Float>& target )
2868 : {
2869 :
2870 346 : Array<Float> lsumwt;
2871 346 : sumwt()->get( lsumwt , False ); // For MT, this will always pick the zeroth sumwt, which it should.
2872 346 : LatticeLocker lock2(target, FileLocker::Write);
2873 :
2874 346 : IPosition imshape = target.shape();
2875 :
2876 : //cerr << " SumWt : " << lsumwt << " sumwtshape : " << lsumwt.shape() << " image shape : " << imshape << endl;
2877 :
2878 346 : AlwaysAssert( lsumwt.shape()[2] == imshape[2] , AipsError ); // polplanes
2879 346 : AlwaysAssert( lsumwt.shape()[3] == imshape[3] , AipsError ); // chanplanes
2880 :
2881 346 : Bool div=False; // flag to signal if division actually happened, or weights are all 1.0
2882 :
2883 692 : for(Int pol=0; pol<lsumwt.shape()[2]; pol++)
2884 : {
2885 692 : for(Int chan=0; chan<lsumwt.shape()[3]; chan++)
2886 : {
2887 346 : IPosition pos(4,0,0,pol,chan);
2888 346 : 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 346 : chan, lsumwt.shape()[3],
2893 346 : pol, lsumwt.shape()[2],
2894 692 : target );
2895 346 : if ( lsumwt(pos) > 1e-07 ) {
2896 692 : LatticeExpr<Float> le( (*subim)/lsumwt(pos) );
2897 346 : subim->copyData( le );
2898 346 : }
2899 : else {
2900 0 : subim->set(0.0);
2901 : }
2902 346 : div=True;
2903 346 : }
2904 346 : }
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 346 : return div;
2913 346 : }
2914 :
2915 73 : void SIImageStore::normPSF(Int term)
2916 : {
2917 :
2918 146 : for(Int pol=0; pol<itsImageShape[2]; pol++)
2919 : {
2920 146 : 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 146 : chan, itsImageShape[3],
2926 73 : pol, itsImageShape[2],
2927 146 : (*psf(term)) );
2928 :
2929 : std::shared_ptr<ImageInterface<Float> > subim0=makeSubImage(0,1,
2930 146 : chan, itsImageShape[3],
2931 73 : pol, itsImageShape[2],
2932 146 : (*psf(0)) );
2933 :
2934 73 : LatticeLocker lock1(*(subim), FileLocker::Write);
2935 :
2936 73 : LatticeExprNode themax( max(*(subim0)) );
2937 73 : Float maxim = themax.getFloat();
2938 :
2939 73 : if ( maxim > 1e-07 )
2940 : {
2941 146 : LatticeExpr<Float> normed( (*(subim)) / maxim );
2942 73 : subim->copyData( normed );
2943 73 : }
2944 : else
2945 : {
2946 0 : subim->set(0.0);
2947 : }
2948 73 : }//chan
2949 : }//pol
2950 :
2951 73 : }
2952 :
2953 73 : void SIImageStore::calcSensitivity()
2954 : {
2955 146 : LogIO os( LogOrigin("SIImageStore","calcSensitivity",WHERE) );
2956 :
2957 73 : Array<Float> lsumwt;
2958 73 : sumwt()->get( lsumwt , False ); // For MT, this will always pick the zeroth sumwt, which it should.
2959 :
2960 73 : IPosition shp( lsumwt.shape() );
2961 : //cout << "Sumwt shape : " << shp << " : " << lsumwt << endl;
2962 : //AlwaysAssert( shp.nelements()==4 , AipsError );
2963 :
2964 73 : os << "[" << itsImageName << "] Theoretical sensitivity (Jy/bm):" ;
2965 :
2966 73 : IPosition it(4,0,0,0,0);
2967 146 : for ( it[0]=0; it[0]<shp[0]; it[0]++)
2968 146 : for ( it[1]=0; it[1]<shp[1]; it[1]++)
2969 146 : for ( it[2]=0; it[2]<shp[2]; it[2]++)
2970 146 : for ( it[3]=0; it[3]<shp[3]; it[3]++)
2971 : {
2972 73 : if( shp[0]>1 ){os << "f"<< it[0]+(it[1]*shp[0]) << ":" ;}
2973 73 : if( shp[3]>1 ) { os << "c"<< it[3] << ":"; }
2974 73 : if( shp[2]>1 ) { os << "p"<< it[2]<< ":" ; }
2975 73 : if( lsumwt( it ) > 1e-07 )
2976 : {
2977 73 : os << 1.0/(sqrt(lsumwt(it))) << " " ;
2978 : }
2979 : else
2980 : {
2981 0 : os << "none" << " ";
2982 : }
2983 : }
2984 :
2985 73 : os << LogIO::POST;
2986 :
2987 : // Array<Float> sens = 1.0/sqrtsumwt;
2988 :
2989 :
2990 73 : }
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 959 : Float SIImageStore::getPeakResidual()
3003 : {
3004 1918 : LogIO os( LogOrigin("SIImageStore","getPeakResidual",WHERE) );
3005 959 : LatticeLocker lock2 (*(residual()), FileLocker::Read);
3006 1918 : ArrayLattice<Bool> pixelmask(residual()->getMask());
3007 1918 : LatticeExpr<Float> resd(iif(pixelmask,abs(*residual()),0));
3008 : //LatticeExprNode pres( max(abs( *residual() ) ));
3009 959 : LatticeExprNode pres( max(resd) );
3010 959 : Float maxresidual = pres.getFloat();
3011 :
3012 959 : return maxresidual;
3013 959 : }
3014 :
3015 616 : Float SIImageStore::getPeakResidualWithinMask()
3016 : {
3017 1232 : LogIO os( LogOrigin("SIImageStore","getPeakResidualWithinMask",WHERE) );
3018 : //Float minresmask, maxresmask, minres, maxres;
3019 : //findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
3020 :
3021 1232 : ArrayLattice<Bool> pixelmask(residual()->getMask());
3022 : //findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
3023 1232 : LatticeExpr<Float> resd(iif(((*mask()) > 0) && pixelmask ,abs((*residual())),0));
3024 616 : LatticeExprNode pres( max(resd) );
3025 616 : Float maxresidual = pres.getFloat();
3026 : //Float maxresidual=max( abs(maxresmask), abs(minresmask) );
3027 616 : return maxresidual;
3028 616 : }
3029 :
3030 : // Calculate the total model flux
3031 616 : Float SIImageStore::getModelFlux(uInt term)
3032 : {
3033 : // LogIO os( LogOrigin("SIImageStore","getModelFlux",WHERE) );
3034 616 : LatticeLocker lock2 (*(model(term)), FileLocker::Read);
3035 1232 : LatticeExprNode mflux( sum( *model(term) ) );
3036 616 : Float modelflux = mflux.getFloat();
3037 : // Float modelflux = sum( model(term)->get() );
3038 :
3039 616 : return modelflux;
3040 616 : }
3041 :
3042 : // Check for non-zero model (this is different from getting model flux, for derived SIIMMT)
3043 419 : Bool SIImageStore::isModelEmpty()
3044 : {
3045 419 : if( !itsModel && (! hasModel()) ) return True;
3046 343 : else return ( fabs( getModelFlux(0) ) < 1e-08 );
3047 : }
3048 :
3049 :
3050 0 : void SIImageStore::setPSFSidelobeLevel(const Float level){
3051 :
3052 0 : itsPSFSideLobeLevel=level;
3053 0 : }
3054 : // Calculate the PSF sidelobe level...
3055 826 : Float SIImageStore::getPSFSidelobeLevel()
3056 : {
3057 1652 : LogIO os( LogOrigin("SIImageStore","getPSFSidelobeLevel",WHERE) );
3058 : //cerr << "*****PSF sidelobe "<< itsPSFSideLobeLevel << endl;
3059 : /// Calculate only once, store and return for all subsequent calls.
3060 826 : if( itsPSFSideLobeLevel == 0.0 )
3061 : {
3062 :
3063 70 : ImageBeamSet thebeams = getBeamSet();
3064 70 : LatticeLocker lock2 (*(psf()), FileLocker::Read);
3065 :
3066 : //------------------------------------------------------------
3067 70 : IPosition oneplaneshape( itsImageShape );
3068 70 : AlwaysAssert( oneplaneshape.nelements()==4, AipsError );
3069 70 : oneplaneshape[2]=1; oneplaneshape[3]=1;
3070 70 : TempImage<Float> psfbeam( oneplaneshape, itsCoordSys );
3071 :
3072 : // In a loop through channels, subtract out or mask out the main lobe
3073 70 : Float allmin=0.0, allmax=0.0;
3074 140 : for(Int pol=0; pol<itsImageShape[2]; pol++)
3075 : {
3076 140 : for(Int chan=0; chan<itsImageShape[3]; chan++)
3077 : {
3078 : std::shared_ptr<ImageInterface<Float> > onepsf=makeSubImage(0,1,
3079 140 : chan, itsImageShape[3],
3080 70 : pol, itsImageShape[2],
3081 140 : (*psf()) );
3082 :
3083 :
3084 70 : GaussianBeam beam = thebeams.getBeam( chan, pol );
3085 70 : Vector<Float> abeam(3); // Holds bmaj, bmin, pa in asec, asec, deg
3086 70 : abeam[0] = beam.getMajor().get("arcsec").getValue() * C::arcsec;
3087 70 : abeam[1] = beam.getMinor().get("arcsec").getValue() * C::arcsec;
3088 70 : abeam[2] = (beam.getPA().get("deg").getValue() + 90.0)* C::degree;
3089 :
3090 : //cout << "Beam : " << abeam << endl;
3091 :
3092 70 : StokesImageUtil::MakeGaussianPSF( psfbeam, abeam, False);
3093 :
3094 : // storeImg( String("psfbeam.im"), psfbeam );
3095 :
3096 : //Subtract from PSF plane
3097 140 : 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 70 : LatticeExprNode minval_le( min( *onepsf ) );
3105 70 : LatticeExprNode maxval_le( max( delobed ) );
3106 :
3107 70 : Float minval = minval_le.getFloat();
3108 70 : Float maxval = maxval_le.getFloat();
3109 :
3110 70 : if( minval < allmin ) allmin = minval;
3111 70 : if( maxval > allmax ) allmax = maxval;
3112 :
3113 : //cout << "Chan : " << chan << " minval : " << minval << " maxval : " << maxval << endl;
3114 :
3115 70 : }//chan
3116 : }//pol
3117 :
3118 : //------------------------------------------------------------
3119 :
3120 70 : itsPSFSideLobeLevel = max( fabs(allmin), fabs(allmax) );
3121 :
3122 : //os << "PSF min : " << allmin << " max : " << allmax << " psfsidelobelevel : " << itsPSFSideLobeLevel << LogIO::POST;
3123 :
3124 70 : }// if changed.
3125 :
3126 : // LatticeExprNode psfside( min( *psf() ) );
3127 : // itsPSFSideLobeLevel = fabs( psfside.getFloat() );
3128 :
3129 : //cout << "PSF sidelobe level : " << itsPSFSideLobeLevel << endl;
3130 826 : return itsPSFSideLobeLevel;
3131 826 : }
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 0 : Array<Double> SIImageStore::calcRobustRMS(Array<Double>& mdns, const Float pbmasklevel, const Bool fastcalc)
3148 : {
3149 0 : LogIO os( LogOrigin("SIImageStore","calcRobustRMS",WHERE) );
3150 0 : Record* regionPtr=0;
3151 0 : String LELmask("");
3152 0 : LatticeLocker lockres (*(residual()), FileLocker::Read);
3153 0 : ArrayLattice<Bool> pbmasklat(residual()->shape());
3154 0 : pbmasklat.set(False);
3155 0 : LatticeExpr<Bool> pbmask(pbmasklat);
3156 0 : if (hasPB()) {
3157 : // set bool mask: False = masked
3158 0 : pbmask = LatticeExpr<Bool> (iif(*pb() > pbmasklevel, True, False));
3159 0 : os << LogIO::DEBUG1 << "pbmask to be attached===> nfalse(pbmask.getMask())="<<nfalse(pbmask.getMask())<<endl;
3160 : }
3161 :
3162 :
3163 0 : Record thestats;
3164 0 : 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 0 : auto tempRes = std::shared_ptr<TempImage<Float>>(new TempImage<Float>(residual()->shape(), residual()->coordinates(), SDMaskHandler::memoryToUse()));
3168 0 : tempRes->copyData(*residual());
3169 0 : tempRes->attachMask(pbmask);
3170 : //thestats = SDMaskHandler::calcImageStatistics(*residual(), LELmask, regionPtr, True);
3171 0 : thestats = SDMaskHandler::calcImageStatistics(*tempRes, LELmask, regionPtr, True);
3172 0 : }
3173 : else { // older way to calculate
3174 : // use the new statistic calculation algorithm
3175 0 : Vector<Bool> dummyvec;
3176 : // TT: 2018.08.01 using revised version (the older version of this is renameed to calcRobustImageStatisticsOld)
3177 0 : thestats = SDMaskHandler::calcRobustImageStatistics(*residual(), *mask(), pbmask, LELmask, regionPtr, True, dummyvec);
3178 0 : }
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 0 : Array<Double>rmss, mads;
3195 : //thestats.get(RecordFieldId("max"), maxs);
3196 0 : thestats.get(RecordFieldId("rms"), rmss);
3197 0 : thestats.get(RecordFieldId("medabsdevmed"), mads);
3198 0 : thestats.get(RecordFieldId("median"), mdns);
3199 :
3200 : //os << LogIO::DEBUG1 << "Max : " << maxs << LogIO::POST;
3201 0 : os << LogIO::DEBUG1 << "RMS : " << rmss << LogIO::POST;
3202 0 : 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 0 : return mads*1.4826;
3208 0 : }
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 686 : Float SIImageStore::getMaskSum()
3254 : {
3255 1372 : LogIO os( LogOrigin("SIImageStore","getMaskSum",WHERE) );
3256 686 : LatticeLocker lock2 (*(mask()), FileLocker::Read);
3257 1372 : LatticeExprNode msum( sum( *mask() ) );
3258 686 : Float masksum = msum.getFloat();
3259 :
3260 : // Float masksum = sum( mask()->get() );
3261 :
3262 686 : return masksum;
3263 686 : }
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 2022 : void SIImageStore::init()
3334 : {
3335 2022 : imageExts.resize(MAX_IMAGE_IDS);
3336 :
3337 2022 : imageExts(MASK) = ".mask";
3338 2022 : imageExts(PSF) = ".psf";
3339 2022 : imageExts(MODEL) = ".model";
3340 2022 : if (not itsIsSingleDishStore) {
3341 2022 : imageExts(RESIDUAL) = ".residual";
3342 2022 : imageExts(IMAGE) = ".image";
3343 : }
3344 : else {
3345 : // The initial residual image IS the single-dish image
3346 0 : imageExts(RESIDUAL) = ".image";
3347 : // Make sure we have no duplicates in the vector
3348 : // Not sure what should be done here
3349 0 : imageExts(IMAGE) = ".wrongly-deconvolved-single-dish-image";
3350 : }
3351 2022 : imageExts(WEIGHT) = ".weight";
3352 2022 : imageExts(SUMWT) = ".sumwt";
3353 2022 : imageExts(GRIDWT) = ".gridwt";
3354 2022 : imageExts(PB) = ".pb";
3355 2022 : imageExts(FORWARDGRID) = ".forward";
3356 2022 : imageExts(BACKWARDGRID) = ".backward";
3357 2022 : imageExts(IMAGEPBCOR) = ".image.pbcor";
3358 :
3359 2022 : itsParentPsf = itsPsf;
3360 2022 : itsParentModel = itsModel;
3361 2022 : itsParentResidual = itsResidual;
3362 2022 : itsParentWeight = itsWeight;
3363 2022 : itsParentImage = itsImage;
3364 2022 : itsParentSumWt = itsSumWt;
3365 2022 : itsParentMask = itsMask;
3366 2022 : itsParentImagePBcor = itsImagePBcor;
3367 :
3368 : // cout << "parent shape : " << itsParentImageShape
3369 : // << " shape : " << itsImageShape << endl;
3370 2022 : itsParentImageShape = itsImageShape;
3371 2022 : itsParentCoordSys = itsCoordSys;
3372 :
3373 2022 : if ( itsNFacets>1 or itsNChanChunks>1 or itsNPolChunks>1 ) {
3374 0 : itsImageShape = IPosition(4,0,0,0,0);
3375 : }
3376 :
3377 2022 : itsOpened = 0;
3378 :
3379 2022 : itsPSFSideLobeLevel = 0.0;
3380 2022 : itsReadLock = nullptr;
3381 2022 : itsDataPolRep = StokesImageUtil::UNKNOWN; // Should throw an exception if
3382 : // it is not initialized properly
3383 2022 : }
3384 :
3385 :
3386 0 : void SIImageStore::regridToModelImage( ImageInterface<Float> &inputimage, Int term )
3387 : {
3388 : try
3389 : {
3390 :
3391 : //output coords
3392 0 : IPosition outshape = itsImageShape;
3393 0 : CoordinateSystem outcsys = itsCoordSys;
3394 0 : DirectionCoordinate outDirCsys = outcsys.directionCoordinate();
3395 0 : SpectralCoordinate outSpecCsys = outcsys.spectralCoordinate();
3396 :
3397 : // do regrid
3398 0 : IPosition axes(3,0, 1, 2);
3399 0 : IPosition inshape = inputimage.shape();
3400 0 : CoordinateSystem incsys = inputimage.coordinates();
3401 0 : DirectionCoordinate inDirCsys = incsys.directionCoordinate();
3402 0 : SpectralCoordinate inSpecCsys = incsys.spectralCoordinate();
3403 :
3404 0 : Vector<Int> dirAxes = CoordinateUtil::findDirectionAxes(incsys);
3405 0 : axes(0) = dirAxes(0);
3406 0 : axes(1) = dirAxes(1);
3407 0 : axes(2) = CoordinateUtil::findSpectralAxis(incsys);
3408 :
3409 : try {
3410 0 : ImageRegrid<Float> imregrid;
3411 0 : imregrid.regrid( *(model(term)), Interpolate2D::LINEAR, axes, inputimage );
3412 0 : } catch (AipsError &x) {
3413 0 : throw(AipsError("ImageRegrid error : "+ x.getMesg()));
3414 0 : }
3415 :
3416 0 : }catch(AipsError &x)
3417 : {
3418 0 : throw(AipsError("Error in regridding input model image to target coordsys : " + x.getMesg()));
3419 0 : }
3420 0 : }
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 0 : Bool SIImageStore::intersectComplexImage(const String& ComplexImageName){
3545 0 : Vector<Int> whichStokes(0);
3546 0 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
3547 0 : whichStokes, itsDataPolRep);
3548 :
3549 :
3550 : //cerr <<"itsDataPolRep " << itsDataPolRep << endl;
3551 :
3552 0 : CountedPtr<PagedImage<Complex> > compliantImage =nullptr;
3553 : {
3554 0 : PagedImage<Complex>inputImage(ComplexImageName);
3555 0 : IPosition theShape=itsImageShape;
3556 0 : theShape(0)=inputImage.shape()(0);
3557 0 : theShape(1)=inputImage.shape()(1);
3558 0 : CoordinateSystem inpcsys=inputImage.coordinates();
3559 0 : Vector<Double> refpix=cimageCoord.referencePixel();
3560 0 : refpix(0)+=(theShape(0)-itsImageShape(0))/2.0;
3561 0 : refpix(1)+=(theShape(1)-itsImageShape(1))/2.0;
3562 0 : cimageCoord.setReferencePixel(refpix);
3563 0 : String tmpImage=File::newUniqueName(".", "TempImage").absoluteName();
3564 0 : compliantImage=new PagedImage<Complex>(theShape, cimageCoord, tmpImage);
3565 0 : compliantImage->set(0.0);
3566 0 : IPosition iblc(theShape.nelements(),0);
3567 0 : IPosition itrc=theShape-1;
3568 : //cerr << "blc " << iblc << " trc " << itrc << " shape " << theShape << endl;
3569 0 : LCBox lbox(iblc, itrc, theShape);
3570 0 : ImageRegion imagreg(WCBox(lbox, cimageCoord));
3571 :
3572 :
3573 0 : SubImage<Complex> subim(inputImage, imagreg, false);
3574 : //cerr << "shapes " << inputImage.shape() << " sub " << subim.shape() << " compl " << compliantImage->shape() << endl;
3575 0 : compliantImage->copyData(subim);
3576 :
3577 0 : }
3578 0 : TableUtil::deleteTable(ComplexImageName);
3579 0 : compliantImage->rename(ComplexImageName);
3580 0 : return True;
3581 :
3582 0 : }
3583 :
3584 :
3585 : //////////////////////////////////////////////////////////////////////////////////////////////////////
3586 : //////////////////////////////////////////////////////////////////////////////////////////////////////
3587 :
3588 : } //# NAMESPACE CASA - END
3589 :
|
