Line data Source code
1 : //# SynthesisNormalizer.cc: Implementation of Gather/Scatter for parallel major cycles
2 : //# Copyright (C) 1997-2008
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Exceptions/Error.h>
29 : #include <iostream>
30 : #include <sstream>
31 :
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 :
36 : #include <casacore/casa/Logging.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : #include <casacore/casa/Logging/LogMessage.h>
39 : #include <casacore/casa/Logging/LogSink.h>
40 : #include <casacore/casa/Logging/LogMessage.h>
41 :
42 : #include <casacore/casa/OS/DirectoryIterator.h>
43 : #include <casacore/casa/OS/File.h>
44 : #include <casacore/casa/OS/Path.h>
45 :
46 : #include <casacore/casa/OS/HostInfo.h>
47 : #include <casacore/lattices/Lattices/LatticeLocker.h>
48 : #include <casacore/images/Images/TempImage.h>
49 : #include <casacore/images/Images/SubImage.h>
50 : #include <casacore/images/Regions/ImageRegion.h>
51 :
52 : #include <synthesis/ImagerObjects/SynthesisNormalizer.h>
53 :
54 : #include <sys/types.h>
55 : #include <unistd.h>
56 : using namespace std;
57 :
58 : using namespace casacore;
59 : namespace casa { //# NAMESPACE CASA - BEGIN
60 :
61 1723 : SynthesisNormalizer::SynthesisNormalizer() :
62 : itsImages(std::shared_ptr<SIImageStore>()),
63 : itsPartImages(Vector<std::shared_ptr<SIImageStore> >()),
64 : itsImageName(""),
65 : itsPartImageNames(Vector<String>(0)),
66 : itsPBLimit(0.2),
67 : itsMapperType("default"),
68 : itsNTaylorTerms(1),
69 1723 : itsNFacets(1)
70 : {
71 1723 : itsFacetImageStores.resize(0);
72 1723 : itsPBLimit = 0.35;
73 1723 : }
74 :
75 1723 : SynthesisNormalizer::~SynthesisNormalizer()
76 : {
77 3446 : LogIO os( LogOrigin("SynthesisNormalizer","destructor",WHERE) );
78 1723 : os << LogIO::DEBUG1 << "SynthesisNormalizer destroyed" << LogIO::POST;
79 1723 : SynthesisUtilMethods::getResource("End SynthesisNormalizer");
80 :
81 1723 : }
82 :
83 :
84 1723 : void SynthesisNormalizer::setupNormalizer(Record normpars)
85 : {
86 5169 : LogIO os( LogOrigin("SynthesisNormalizer","setupNormalizer",WHERE) );
87 :
88 : try
89 : {
90 1723 : if( normpars.isDefined("psfcutoff") ) // A single string
91 : {
92 1723 : normpars.get( RecordFieldId("psfcutoff") , itsPsfcutoff );
93 : }else
94 : {
95 0 : itsPsfcutoff=0.35;
96 : }
97 :
98 :
99 :
100 1723 : if( normpars.isDefined("imagename") ) // A single string
101 1723 : { itsImageName = normpars.asString( RecordFieldId("imagename")); }
102 : else
103 0 : {throw( AipsError("imagename not specified")); }
104 :
105 1723 : if( normpars.isDefined("partimagenames") ) // A vector of strings
106 0 : { normpars.get( RecordFieldId("partimagenames") , itsPartImageNames ); }
107 : else
108 1723 : { itsPartImageNames.resize(0); }
109 :
110 1723 : if( normpars.isDefined("pblimit") )
111 : {
112 1723 : normpars.get( RecordFieldId("pblimit") , itsPBLimit );
113 : }
114 : else
115 0 : { itsPBLimit = 0.2; }
116 :
117 1723 : if( normpars.isDefined("normtype") ) // A single string
118 1723 : { itsNormType = normpars.asString( RecordFieldId("normtype")); }
119 : else
120 0 : { itsNormType = "flatnoise";} // flatnoise, flatsky
121 :
122 : // cout << "Chosen normtype : " << itsNormType << endl;
123 :
124 : // For multi-term choices. Try to eliminate, after making imstores hold aux descriptive info.
125 : /*
126 : if( normpars.isDefined("mtype") ) // A single string
127 : { itsMapperType = normpars.asString( RecordFieldId("mtype")); }
128 : else
129 : { itsMapperType = "default";}
130 : */
131 :
132 1723 : if( normpars.isDefined("deconvolver") ) // A single string
133 3446 : { String dec = normpars.asString( RecordFieldId("deconvolver"));
134 1723 : if (dec == "mtmfs") { itsMapperType="multiterm"; }
135 1604 : else itsMapperType="default";
136 : }
137 : else
138 0 : { itsMapperType = "default";}
139 :
140 1723 : if( normpars.isDefined("nterms") ) // A single int
141 1723 : { itsNTaylorTerms = normpars.asuInt( RecordFieldId("nterms")); }
142 : else
143 0 : { itsNTaylorTerms = 1;}
144 :
145 1723 : if( normpars.isDefined("facets") ) // A single int
146 1723 : { itsNFacets = normpars.asuInt( RecordFieldId("facets")); }
147 : else
148 0 : { itsNFacets = 1;}
149 :
150 1723 : if( normpars.isDefined("restoringbeam") )
151 : {
152 1723 : if (normpars.dataType("restoringbeam")==TpString) {
153 50 : itsUseBeam = normpars.asString( RecordFieldId("restoringbeam") ); }
154 : else
155 1673 : { itsUseBeam = "";}
156 : }
157 : }
158 0 : catch(AipsError &x)
159 : {
160 0 : throw( AipsError("Error in reading gather/scatter parameters: "+x.getMesg()) );
161 : }
162 :
163 1723 : }//end of setupParSync
164 :
165 :
166 1420 : void SynthesisNormalizer::gatherImages(Bool dopsf, Bool doresidual, Bool dodensity)
167 : {
168 : // cout << " partimagenames :" << itsPartImageNames << endl;
169 :
170 1420 : Bool needToGatherImages = setupImagesOnDisk();
171 :
172 1420 : if( needToGatherImages )
173 : {
174 0 : LogIO os( LogOrigin("SynthesisNormalizer", "gatherImages",WHERE) );
175 :
176 0 : AlwaysAssert( itsPartImages.nelements()>0 , AipsError );
177 : //Bool doresidual = !dopsf;
178 0 : Bool doweight = dopsf ; //|| ( doresidual && ! itsImages->hasSensitivity() );
179 : // Bool doweight = dopsf || ( doresidual && itsImages->getUseWeightImage(*(itsPartImages[0]->residual())) );
180 :
181 0 : os << "Gather "<< (doresidual?"residual":"") << ( (dopsf&&doresidual)?",":"")
182 0 : << (dopsf?"psf":"") << ( (dopsf&&doweight)?",":"")
183 : << (doweight?"weight":"")<< " images : " << itsPartImageNames
184 0 : << " onto :" << itsImageName << LogIO::POST;
185 :
186 : // Add intelligence to modify all only the first time, but later, only residual;
187 0 : itsImages->resetImages( /*psf*/dopsf, /*residual*/doresidual, /*weight*/doweight );
188 :
189 0 : for( uInt part=0;part<itsPartImages.nelements();part++)
190 : {
191 0 : itsImages->addImages( itsPartImages[part], /*psf*/dopsf, /*residual*/doresidual, /*weight*/doweight, /*griddedwt*/dodensity );
192 0 : itsPartImages[part]->releaseLocks();
193 : }
194 :
195 : }// end of image gathering.
196 :
197 : // Normalize by the weight image.
198 : // divideResidualByWeight();
199 1420 : itsImages->releaseLocks();
200 :
201 1420 : }// end of gatherImages
202 :
203 676 : void SynthesisNormalizer::gatherPB()
204 : {
205 676 : if( itsPartImageNames.nelements()>0 )
206 : {
207 :
208 : try
209 : {
210 0 : itsPartImages[0] = makeImageStore( itsPartImageNames[0] );
211 : }
212 0 : catch(AipsError &x)
213 : {
214 0 : throw(AipsError("Cannot construct ImageStore for "+itsPartImageNames[0]+". "+x.what()));
215 : }
216 :
217 : try{
218 0 : LatticeExpr<Float> thepb( *(itsPartImages[0]->pb()) );
219 0 : LatticeLocker lock1(*(itsImages->pb()), FileLocker::Write);
220 0 : itsImages->pb()->copyData(thepb);
221 :
222 : }
223 0 : catch(AipsError &x)
224 : {
225 0 : throw(AipsError("Cannot copy the PB : "+x.getMesg()));
226 : }
227 :
228 : }
229 676 : }
230 :
231 914 : void SynthesisNormalizer::scatterModel()
232 : {
233 :
234 2742 : LogIO os( LogOrigin("SynthesisNormalizer", "scatterModel",WHERE) );
235 :
236 914 : setupImagesOnDisk(); // To open up and initialize itsPartImages.
237 :
238 : // os << "In ScatterModel : " << itsPartImages.nelements() << " for " << itsPartImageNames << LogIO::POST;
239 :
240 914 : if( itsPartImages.nelements() > 0 ) // && itsImages->doesImageExist(modelName) )
241 : {
242 0 : os << "Send the model from : " << itsImageName << " to all nodes :" << itsPartImageNames << LogIO::POST;
243 :
244 : // Make the list of model images. This list is of length >1 only for multi-term runs.
245 : // Vector<String> modelNames( itsImages->getNTaylorTerms() );
246 : // if( modelNames.nelements() ==1 ) modelNames[0] = itsImages->getName()+".model";
247 : // if( modelNames.nelements() > 1 )
248 : // {
249 : // for( uInt nt=0;nt<itsImages->getNTaylorTerms();nt++)
250 : // modelNames[nt] = itsImages->getName()+".model.tt" + String::toString(nt);
251 : // }
252 :
253 :
254 0 : Vector<String> modelNames( itsImages->getNTaylorTerms() );
255 0 : if( itsImages->getType()=="default" ) modelNames[0] = itsImages->getName()+".model";
256 0 : if( itsImages->getType()=="multiterm" )
257 : {
258 0 : for( uInt nt=0;nt<itsImages->getNTaylorTerms();nt++)
259 0 : modelNames[nt] = itsImages->getName()+".model.tt" + String::toString(nt);
260 : }
261 :
262 :
263 :
264 0 : for( uInt part=0;part<itsPartImages.nelements();part++)
265 : {
266 0 : itsPartImages[part]->setModelImage( modelNames );
267 0 : itsPartImages[part]->releaseLocks();
268 : }
269 0 : itsImages->releaseLocks();
270 : }
271 :
272 914 : }// end of scatterModel
273 :
274 0 : void SynthesisNormalizer::scatterWeightDensity()
275 : {
276 :
277 0 : LogIO os( LogOrigin("SynthesisNormalizer", "scatterWeightDensity",WHERE) );
278 :
279 0 : setupImagesOnDisk(); // To open up and initialize itsPartImages.
280 :
281 : // os << "In ScatterModel : " << itsPartImages.nelements() << " for " << itsPartImageNames << LogIO::POST;
282 :
283 0 : if( itsPartImages.nelements() > 0 )
284 : {
285 0 : os << "Send the gridded weight from : " << itsImageName << " to all nodes :" << itsPartImageNames << LogIO::POST;
286 :
287 0 : for( uInt part=0;part<itsPartImages.nelements();part++)
288 : {
289 0 : itsPartImages[part]->setWeightDensity( itsImages );
290 0 : itsPartImages[part]->releaseLocks();
291 : }
292 0 : itsImages->releaseLocks();
293 : }
294 0 : }// end of gatherImages
295 :
296 :
297 :
298 1428 : void SynthesisNormalizer::divideResidualByWeight()
299 : {
300 4284 : LogIO os( LogOrigin("SynthesisNormalizer", "divideResidualByWeight",WHERE) );
301 :
302 1428 : if( itsNFacets==1) {
303 1420 : itsImages->divideResidualByWeight( itsPBLimit, itsNormType );
304 : }
305 : else {
306 64 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
307 56 : { itsFacetImageStores[facet]->divideResidualByWeight( itsPBLimit , itsNormType ); }
308 : }
309 1428 : itsImages->releaseLocks();
310 :
311 1428 : }
312 :
313 129 : void SynthesisNormalizer::divideResidualByWeightSD()
314 : {
315 387 : LogIO os( LogOrigin("SynthesisNormalizer", "divideResidualByWeightSD",WHERE) );
316 :
317 129 : if( itsNFacets==1) {
318 129 : itsImages->divideResidualByWeightSD( itsPBLimit );
319 : }
320 : else {
321 0 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
322 0 : { itsFacetImageStores[facet]->divideResidualByWeightSD( itsPBLimit ); }
323 : }
324 129 : itsImages->releaseLocks();
325 :
326 129 : }
327 :
328 676 : void SynthesisNormalizer::dividePSFByWeight()
329 : {
330 2028 : LogIO os( LogOrigin("SynthesisNormalizer", "dividePSFByWeight",WHERE) );
331 : {
332 :
333 1352 : LatticeLocker lock1 (*(itsImages->psf()), FileLocker::Write);
334 :
335 676 : if( itsNFacets==1) {
336 672 : itsImages->dividePSFByWeight(itsPBLimit);
337 : }
338 : else {
339 : // Since PSFs are normed just by their max, this sequence is OK.
340 4 : setPsfFromOneFacet();
341 4 : itsImages->dividePSFByWeight(itsPBLimit);
342 : }
343 : }// release of lock1 otherwise releaselock below will cause it to core
344 : //dump as the psf pointer goes dangling
345 : // Check PSF quality by fitting beams
346 : {
347 676 : itsImages->calcSensitivity();
348 :
349 676 : itsImages->makeImageBeamSet(itsPsfcutoff);
350 676 : Bool verbose(False);
351 676 : if (itsUseBeam=="common") verbose=True;
352 676 : itsImages->printBeamSet(verbose);
353 : }
354 676 : itsImages->releaseLocks();
355 :
356 676 : }
357 :
358 676 : void SynthesisNormalizer::normalizePrimaryBeam()
359 : {
360 2028 : LogIO os( LogOrigin("SynthesisNormalizer", "normalizePrimaryBeam",WHERE) );
361 :
362 676 : if( itsImages==NULL ) { itsImages = makeImageStore( itsImageName ); }
363 :
364 676 : gatherPB();
365 :
366 : // Irrespective of facets.
367 676 : itsImages->normalizePrimaryBeam(itsPBLimit);
368 676 : }
369 :
370 :
371 980 : void SynthesisNormalizer::divideModelByWeight()
372 : {
373 2940 : LogIO os( LogOrigin("SynthesisNormalizer", "divideModelByWeight",WHERE) );
374 980 : if( ! itsImages )
375 : {
376 : //os << LogIO::WARN << "No imagestore yet. Trying to construct, so that the starting model can be divided by wt if needed...." << LogIO::POST;
377 :
378 : try
379 : {
380 22 : itsImages = makeImageStore( itsImageName );
381 : }
382 0 : catch(AipsError &x)
383 : {
384 0 : throw(AipsError("Cannot construct ImageStore for "+itsImageName+". "+x.what()));
385 : }
386 : // return;
387 : }
388 980 : if( itsNFacets==1) {
389 972 : itsImages->divideModelByWeight( itsPBLimit, itsNormType );
390 : }
391 : else {
392 64 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
393 56 : { itsFacetImageStores[facet]->divideModelByWeight( itsPBLimit, itsNormType ); }
394 : }
395 980 : }
396 :
397 914 : void SynthesisNormalizer::multiplyModelByWeight()
398 : {
399 : // LogIO os( LogOrigin("SynthesisNormalizer", "multiplyModelByWeight",WHERE) );
400 914 : if( itsNFacets==1) {
401 906 : itsImages->multiplyModelByWeight( itsPBLimit , itsNormType );
402 : }
403 : else {
404 64 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
405 56 : { itsFacetImageStores[facet]->multiplyModelByWeight( itsPBLimit , itsNormType); }
406 : }
407 914 : }
408 :
409 :
410 0 : std::shared_ptr<SIImageStore> SynthesisNormalizer::getImageStore()
411 : {
412 0 : LogIO os( LogOrigin("SynthesisNormalizer", "getImageStore", WHERE) );
413 0 : return itsImages;
414 : }
415 :
416 0 : void SynthesisNormalizer::setImageStore( SIImageStore* imstore )
417 : {
418 0 : LogIO os( LogOrigin("SynthesisNormalizer", "setImageStore", WHERE) );
419 0 : itsImages.reset( imstore );
420 0 : }
421 :
422 879 : void SynthesisNormalizer::setImageStore( std::shared_ptr<SIImageStore>& imstore )
423 : {
424 2637 : LogIO os( LogOrigin("SynthesisNormalizer", "setImageStore", WHERE) );
425 879 : itsImages= imstore ;
426 :
427 879 : }
428 :
429 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
430 : //// Internal Functions start here. These are not visible to the tool layer.
431 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
432 :
433 :
434 2334 : Bool SynthesisNormalizer::setupImagesOnDisk()
435 : {
436 7002 : LogIO os( LogOrigin("SynthesisNormalizer","setupImagesOnDisk",WHERE) );
437 :
438 2334 : Bool needToGatherImages=false;
439 :
440 2334 : String err("");
441 :
442 : // Check if full images exist, and open them if possible.
443 2334 : Bool foundFullImage=false;
444 : try
445 : {
446 2334 : itsImages = makeImageStore( itsImageName );
447 2334 : foundFullImage = true;
448 : }
449 0 : catch(AipsError &x)
450 : {
451 : //throw( AipsError("Error in constructing a Deconvolver : "+x.getMesg()) );
452 0 : err = err += String(x.getMesg()) + "\n";
453 0 : foundFullImage = false;
454 : }
455 :
456 2334 : os << LogIO::DEBUG2 << " Found full images : " << foundFullImage << LogIO::POST;
457 :
458 : // Check if part images exist
459 2334 : Bool foundPartImages = itsPartImageNames.nelements()>0 ? true : false ;
460 2334 : itsPartImages.resize( itsPartImageNames.nelements() );
461 :
462 2334 : for ( uInt part=0; part < itsPartImageNames.nelements() ; part++ )
463 : {
464 : try
465 : {
466 0 : itsPartImages[part] = makeImageStore ( itsPartImageNames[part] );
467 0 : foundPartImages |= true;
468 : }
469 0 : catch(AipsError &x)
470 : {
471 : //throw( AipsError("Error in constructing a Deconvolver : "+x.getMesg()) );
472 0 : err = err += String(x.getMesg()) + "\n";
473 0 : foundPartImages = false;
474 : }
475 : }
476 :
477 2334 : os << LogIO::DEBUG2 << " Found part images : " << foundPartImages << LogIO::POST;
478 :
479 2334 : if( foundPartImages == false)
480 : {
481 2334 : if( foundFullImage == true && itsPartImageNames.nelements()>0 )
482 : {
483 : // Pick the coordsys, etc from fullImage, and construct new/fresh partial images.
484 0 : os << LogIO::DEBUG2 << "Found full image, but no partial images. Make partImStores for : " << itsPartImageNames << LogIO::POST;
485 :
486 0 : String imopen = itsImages->getName()+".residual"+((itsMapperType=="multiterm")?".tt0":"");
487 0 : Directory imdir( imopen );
488 0 : if( ! imdir.exists() )
489 : {
490 0 : imopen = itsImages->getName()+".psf"+((itsMapperType=="multiterm")?".tt0":"");
491 0 : Directory imdir2( imopen );
492 0 : if( ! imdir2.exists() )
493 0 : throw(AipsError("Cannot find partial image psf or residual for " +itsImages->getName() +err));
494 : }
495 :
496 0 : PagedImage<Float> temppart(imopen);
497 :
498 0 : Bool useweightimage = itsImages->getUseWeightImage( *(itsImages->sumwt()) );
499 0 : for( uInt part=0; part<itsPartImageNames.nelements(); part++ )
500 : {
501 0 : itsPartImages[part] = makeImageStore (itsPartImageNames[part], temppart,
502 0 : useweightimage);
503 : }
504 0 : foundPartImages = True;
505 : }
506 : else
507 : {
508 2334 : itsPartImages.resize(0);
509 2334 : foundPartImages = False;
510 : }
511 : }
512 : else // Check that all have the same shape.
513 : {
514 0 : AlwaysAssert( itsPartImages.nelements() > 0 , AipsError );
515 0 : IPosition tempshape = itsPartImages[0]->getShape();
516 0 : for( uInt part=1; part<itsPartImages.nelements(); part++ )
517 : {
518 0 : if( tempshape != itsPartImages[part]->getShape() )
519 : {
520 0 : throw( AipsError("Shapes of partial images to be combined, do not match" + err) );
521 : }
522 : }
523 : }
524 :
525 :
526 : // Make sure all images exist and are consistent with each other. At the end, itsImages should be valid
527 2334 : if( foundPartImages == true ) // Partial Images exist. Check that 'full' exists, and do the gather.
528 : {
529 0 : if ( foundFullImage == true ) // Full image exists. Just check that shapes match with parts.
530 : {
531 0 : os << LogIO::DEBUG2 << "Partial and Full images exist. Checking if part images have the same shape as the full image : ";
532 0 : IPosition fullshape = itsImages->getShape();
533 :
534 0 : for ( uInt part=0; part < itsPartImages.nelements() ; part++ )
535 : {
536 0 : IPosition partshape = itsPartImages[part]->getShape();
537 0 : if( partshape != fullshape )
538 : {
539 0 : os << LogIO::DEBUG2<< "NO" << LogIO::POST;
540 0 : throw( AipsError("Shapes of the partial and full images on disk do not match. Cannot gather" + err) );
541 : }
542 : }
543 0 : os << LogIO::DEBUG2 << "Yes" << LogIO::POST;
544 :
545 : }
546 : else // Full image does not exist. Need to make it, using the shape and coords of part[0]
547 : {
548 0 : os << LogIO::DEBUG2 << "Only partial images exist. Need to make full images" << LogIO::POST;
549 :
550 0 : AlwaysAssert( itsPartImages.nelements() > 0, AipsError );
551 :
552 : // Find an image to open and pick csys,shape from.
553 0 : String imopen = itsPartImageNames[0]+".residual"+((itsMapperType=="multiterm")?".tt0":"");
554 0 : Directory imdir( imopen );
555 0 : if( ! imdir.exists() )
556 : {
557 0 : imopen = itsPartImageNames[0]+".psf"+((itsMapperType=="multiterm")?".tt0":"");
558 0 : Directory imdir2( imopen );
559 0 : if( ! imdir2.exists() )
560 : {
561 0 : imopen = itsPartImageNames[0]+".gridwt";
562 0 : Directory imdir3( imopen );
563 0 : if( ! imdir3.exists() )
564 0 : throw(AipsError("Cannot find partial image psf or residual or gridwt for " + itsPartImageNames[0]+err));
565 : }
566 :
567 : }
568 :
569 0 : PagedImage<Float> temppart( imopen );
570 :
571 0 : Bool useweightimage = itsPartImages[0]->getUseWeightImage( *(itsPartImages[0]->sumwt()) );
572 0 : itsImages = makeImageStore (itsImageName, temppart, useweightimage);
573 0 : foundFullImage = true;
574 : }
575 :
576 : // By now, all partial images and the full images exist on disk, and have the same shape.
577 0 : needToGatherImages=true;
578 :
579 : }
580 : else // No partial images supplied. Operating only with full images.
581 : {
582 2334 : if ( foundFullImage == true )
583 : {
584 2334 : os << LogIO::DEBUG2 << "Full images exist : " << itsImageName << LogIO::POST;
585 : }
586 : else // No full image on disk either. Error.
587 : {
588 0 : throw( AipsError("No images named " + itsImageName + " found on disk. No partial images found either."+err) );
589 : }
590 : }
591 :
592 : // Remove ?
593 2334 : itsImages->psf();
594 2334 : itsImages->validate();
595 :
596 : // Set up facet Imstores..... if needed
597 2334 : if( itsNFacets>1 )
598 : {
599 :
600 : // First, make sure that full images have been allocated before trying to make references.....
601 : // if( ! itsImages->checkValidity(true/*psf*/, true/*res*/,true/*wgt*/,true/*model*/,false/*image*/,false/*mask*/,true/*sumwt*/ ) )
602 : // { throw(AipsError("Internal Error : Invalid ImageStore for " + itsImages->getName())); }
603 :
604 : // Array<Float> ttt;
605 : // (itsImages->sumwt())->get(ttt);
606 : // cout << "SUMWT full : " << ttt << endl;
607 :
608 20 : itsFacetImageStores.resize( itsNFacets*itsNFacets );
609 172 : for( uInt facet=0; facet<itsNFacets*itsNFacets; ++facet )
610 : {
611 152 : itsFacetImageStores[facet] = itsImages->getSubImageStore(facet, itsNFacets);
612 :
613 : //Array<Float> qqq;
614 : //itsFacetImageStores[facet]->sumwt()->get(qqq);
615 : // cout << "SUMWTs for " << facet << " : " << qqq << endl;
616 :
617 : }
618 : }
619 2334 : os << LogIO::DEBUG2 << "Need to Gather ? " << needToGatherImages << LogIO::POST;
620 2334 : itsImages->releaseLocks();
621 4668 : return needToGatherImages;
622 : }// end of setupImagesOnDisk
623 :
624 :
625 2636 : std::shared_ptr<SIImageStore> SynthesisNormalizer::makeImageStore(const String &imagename )
626 : {
627 2636 : if( itsMapperType == "multiterm" )
628 449 : { return std::shared_ptr<SIImageStore>(new SIImageStoreMultiTerm( imagename, itsNTaylorTerms, true )); }
629 : else
630 2187 : { return std::shared_ptr<SIImageStore>(new SIImageStore( imagename, true )); }
631 : itsImages->releaseLocks();
632 : }
633 :
634 :
635 : /**
636 : * build a new ImageStore, whether SIImageStore or SIImageStoreMultiTerm, borrowing
637 : * image information from one partial image.
638 : *
639 : * @param imagename image name for the new SIStorageManager
640 : * @param part partial image from which miscinfo, etc. data will be borrowed
641 : * @param useweightimage useweight option for the new SIStorageManager
642 : *
643 : * @return A new SIImageStore object for the image name given.
644 : */
645 0 : std::shared_ptr<SIImageStore> SynthesisNormalizer::makeImageStore(const String &imagename,
646 : const PagedImage<Float> &part,
647 : Bool useweightimage)
648 : {
649 : // borrow shape, coord, imageinfo and miscinfo
650 0 : auto shape = part.shape();
651 0 : auto csys = part.coordinates();
652 0 : auto objectname = part.imageInfo().objectName();
653 0 : auto miscinfo = part.miscInfo();
654 0 : if( itsMapperType == "multiterm" )
655 : {
656 : std::shared_ptr<SIImageStore> multiTermStore =
657 0 : std::make_shared<SIImageStoreMultiTerm>(imagename, csys, shape, objectname,
658 0 : miscinfo, itsNFacets, false, itsNTaylorTerms, useweightimage );
659 0 : return multiTermStore;
660 : }
661 : else
662 : {
663 : return std::make_shared<SIImageStore>(imagename, csys, shape, objectname, miscinfo,
664 0 : false, useweightimage);
665 : }
666 : }
667 :
668 : //
669 : //---------------------------------------------------------------
670 : //
671 4 : void SynthesisNormalizer::setPsfFromOneFacet()
672 : {
673 : // Copy the PSF from one facet to the center of the full image, to use for the minor cycle
674 : //
675 : // This code segment will work for single and multi-term
676 : // - for single term, the index will always be 0, and SIImageStore's access functions know this.
677 : // - for multi-term, the index will be the taylor index and SIImageStoreMultiTerm knows this.
678 :
679 : /// itsImages
680 : /// itsFacetImageStores
681 :
682 : {
683 8 : IPosition shape=(itsImages->psf(0))->shape();
684 8 : IPosition blc(4, 0, 0, 0, 0);
685 8 : IPosition trc=shape-1;
686 10 : for(uInt tix=0; tix<2 * itsImages->getNTaylorTerms() - 1; tix++)
687 : {
688 12 : TempImage<Float> onepsf((itsFacetImageStores[0]->psf(tix))->shape(),
689 24 : (itsFacetImageStores[0]->psf(tix))->coordinates());
690 6 : onepsf.copyData(*(itsFacetImageStores[0]->psf(tix)));
691 : //now set the original to 0 as we have a copy of one facet psf
692 6 : (itsImages->psf(tix))->set(0.0);
693 6 : blc[0]=(shape[0]-(onepsf.shape()[0]))/2;
694 6 : trc[0]=onepsf.shape()[0]+blc[0]-1;
695 6 : blc[1]=(shape[1]-(onepsf.shape()[1]))/2;
696 6 : trc[1]=onepsf.shape()[1]+blc[1]-1;
697 12 : Slicer sl(blc, trc, Slicer::endIsLast);
698 18 : SubImage<Float> sub(*(itsImages->psf(tix)), sl, true);
699 6 : sub.copyData(onepsf);
700 : }
701 : }
702 :
703 : // cout << "In setPsfFromOneFacet : sumwt : " << itsImages->sumwt()->get() << endl;
704 :
705 4 : }
706 :
707 :
708 :
709 : } //# NAMESPACE CASA - END
710 :
|
