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: 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 :
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 1796 : SynthesisNormalizer::SynthesisNormalizer() :
62 1796 : itsImages(std::shared_ptr<SIImageStore>()),
63 1796 : itsPartImages(Vector<std::shared_ptr<SIImageStore> >()),
64 1796 : itsImageName(""),
65 1796 : itsPartImageNames(Vector<String>(0)),
66 1796 : itsPBLimit(0.2),
67 1796 : itsMapperType("default"),
68 1796 : itsNTaylorTerms(1),
69 5388 : itsNFacets(1)
70 : {
71 1796 : itsFacetImageStores.resize(0);
72 1796 : itsPBLimit = 0.35;
73 1796 : }
74 :
75 1796 : SynthesisNormalizer::~SynthesisNormalizer()
76 : {
77 3592 : LogIO os( LogOrigin("SynthesisNormalizer","destructor",WHERE) );
78 1796 : os << LogIO::DEBUG1 << "SynthesisNormalizer destroyed" << LogIO::POST;
79 1796 : SynthesisUtilMethods::getResource("End SynthesisNormalizer");
80 :
81 1796 : }
82 :
83 :
84 1796 : void SynthesisNormalizer::setupNormalizer(Record normpars)
85 : {
86 3592 : LogIO os( LogOrigin("SynthesisNormalizer","setupNormalizer",WHERE) );
87 :
88 : try
89 : {
90 1796 : if( normpars.isDefined("psfcutoff") ) // A single string
91 : {
92 1796 : normpars.get( RecordFieldId("psfcutoff") , itsPsfcutoff );
93 : }else
94 : {
95 0 : itsPsfcutoff=0.35;
96 : }
97 :
98 :
99 :
100 1796 : if( normpars.isDefined("imagename") ) // A single string
101 1796 : { itsImageName = normpars.asString( RecordFieldId("imagename")); }
102 : else
103 0 : {throw( AipsError("imagename not specified")); }
104 :
105 1796 : if( normpars.isDefined("partimagenames") ) // A vector of strings
106 0 : { normpars.get( RecordFieldId("partimagenames") , itsPartImageNames ); }
107 : else
108 1796 : { itsPartImageNames.resize(0); }
109 :
110 1796 : if( normpars.isDefined("pblimit") )
111 : {
112 1796 : normpars.get( RecordFieldId("pblimit") , itsPBLimit );
113 : }
114 : else
115 0 : { itsPBLimit = 0.2; }
116 :
117 1796 : if( normpars.isDefined("normtype") ) // A single string
118 1796 : { 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 1796 : if( normpars.isDefined("deconvolver") ) // A single string
133 1796 : { String dec = normpars.asString( RecordFieldId("deconvolver"));
134 1796 : if (dec == "mtmfs") { itsMapperType="multiterm"; }
135 1670 : else itsMapperType="default";
136 1796 : }
137 : else
138 0 : { itsMapperType = "default";}
139 :
140 1796 : if( normpars.isDefined("nterms") ) // A single int
141 1796 : { itsNTaylorTerms = normpars.asuInt( RecordFieldId("nterms")); }
142 : else
143 0 : { itsNTaylorTerms = 1;}
144 :
145 1796 : if( normpars.isDefined("facets") ) // A single int
146 1796 : { itsNFacets = normpars.asuInt( RecordFieldId("facets")); }
147 : else
148 0 : { itsNFacets = 1;}
149 :
150 1796 : if( normpars.isDefined("restoringbeam") )
151 : {
152 1796 : if (normpars.dataType("restoringbeam")==TpString) {
153 46 : itsUseBeam = normpars.asString( RecordFieldId("restoringbeam") ); }
154 : else
155 1750 : { itsUseBeam = "";}
156 : }
157 : }
158 0 : catch(AipsError &x)
159 : {
160 0 : throw( AipsError("Error in reading gather/scatter parameters: "+x.getMesg()) );
161 0 : }
162 :
163 1796 : }//end of setupParSync
164 :
165 :
166 1461 : void SynthesisNormalizer::gatherImages(Bool dopsf, Bool doresidual, Bool dodensity)
167 : {
168 : // cout << " partimagenames :" << itsPartImageNames << endl;
169 :
170 1461 : Bool needToGatherImages = setupImagesOnDisk();
171 :
172 1461 : 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 0 : }// end of image gathering.
196 :
197 : // Normalize by the weight image.
198 : // divideResidualByWeight();
199 1461 : itsImages->releaseLocks();
200 :
201 1461 : }// end of gatherImages
202 :
203 703 : void SynthesisNormalizer::gatherPB()
204 : {
205 703 : 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 0 : }
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 0 : }
223 0 : catch(AipsError &x)
224 : {
225 0 : throw(AipsError("Cannot copy the PB : "+x.getMesg()));
226 0 : }
227 :
228 : }
229 703 : }
230 :
231 929 : void SynthesisNormalizer::scatterModel()
232 : {
233 :
234 1858 : LogIO os( LogOrigin("SynthesisNormalizer", "scatterModel",WHERE) );
235 :
236 929 : setupImagesOnDisk(); // To open up and initialize itsPartImages.
237 :
238 : // os << "In ScatterModel : " << itsPartImages.nelements() << " for " << itsPartImageNames << LogIO::POST;
239 :
240 929 : 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 0 : }
271 :
272 929 : }// 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 1462 : void SynthesisNormalizer::divideResidualByWeight()
299 : {
300 2924 : LogIO os( LogOrigin("SynthesisNormalizer", "divideResidualByWeight",WHERE) );
301 :
302 1462 : if( itsNFacets==1) {
303 1454 : 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 1462 : itsImages->releaseLocks();
310 :
311 1462 : }
312 :
313 144 : void SynthesisNormalizer::divideResidualByWeightSD()
314 : {
315 288 : LogIO os( LogOrigin("SynthesisNormalizer", "divideResidualByWeightSD",WHERE) );
316 :
317 144 : if( itsNFacets==1) {
318 144 : itsImages->divideResidualByWeightSD( itsPBLimit );
319 : }
320 : else {
321 0 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
322 0 : { itsFacetImageStores[facet]->divideResidualByWeightSD( itsPBLimit ); }
323 : }
324 144 : itsImages->releaseLocks();
325 :
326 144 : }
327 7 : void SynthesisNormalizer::makePSFBeamset(){
328 14 : LogIO os(LogOrigin("SynthesisNormalizer", "dividePSFByWeight", WHERE));
329 : try{
330 7 : if(!itsImages){
331 7 : itsImages = makeImageStore( itsImageName, false );
332 : }
333 :
334 : }
335 0 : catch(AipsError& x){
336 :
337 0 : throw(AipsError("Programmers error no psf is made on disk and trying to fit"));
338 0 : }
339 7 : itsImages->makeImageBeamSet(itsPsfcutoff);
340 7 : itsImages->releaseLocks();
341 7 : }
342 696 : void SynthesisNormalizer::dividePSFByWeight()
343 : {
344 1392 : LogIO os( LogOrigin("SynthesisNormalizer", "dividePSFByWeight",WHERE) );
345 : {
346 :
347 696 : LatticeLocker lock1 (*(itsImages->psf()), FileLocker::Write);
348 :
349 696 : if( itsNFacets==1) {
350 692 : itsImages->dividePSFByWeight(itsPBLimit);
351 : }
352 : else {
353 : // Since PSFs are normed just by their max, this sequence is OK.
354 4 : setPsfFromOneFacet();
355 4 : itsImages->dividePSFByWeight(itsPBLimit);
356 : }
357 696 : }// release of lock1 otherwise releaselock below will cause it to core
358 : //dump as the psf pointer goes dangling
359 : // Check PSF quality by fitting beams
360 : {
361 696 : itsImages->calcSensitivity();
362 :
363 696 : itsImages->makeImageBeamSet(itsPsfcutoff);
364 696 : Bool verbose(False);
365 696 : if (itsUseBeam=="common") verbose=True;
366 696 : itsImages->printBeamSet(verbose);
367 : }
368 696 : itsImages->releaseLocks();
369 :
370 696 : }
371 :
372 703 : void SynthesisNormalizer::normalizePrimaryBeam()
373 : {
374 1406 : LogIO os( LogOrigin("SynthesisNormalizer", "normalizePrimaryBeam",WHERE) );
375 :
376 703 : if( itsImages==NULL ) { itsImages = makeImageStore( itsImageName ); }
377 :
378 703 : gatherPB();
379 :
380 : // Irrespective of facets.
381 703 : itsImages->normalizePrimaryBeam(itsPBLimit);
382 703 : }
383 :
384 :
385 1000 : void SynthesisNormalizer::divideModelByWeight()
386 : {
387 2000 : LogIO os( LogOrigin("SynthesisNormalizer", "divideModelByWeight",WHERE) );
388 1000 : if( ! itsImages )
389 : {
390 : //os << LogIO::WARN << "No imagestore yet. Trying to construct, so that the starting model can be divided by wt if needed...." << LogIO::POST;
391 :
392 : try
393 : {
394 22 : itsImages = makeImageStore( itsImageName );
395 : }
396 0 : catch(AipsError &x)
397 : {
398 0 : throw(AipsError("Cannot construct ImageStore for "+itsImageName+". "+x.what()));
399 0 : }
400 : // return;
401 : }
402 1000 : if( itsNFacets==1) {
403 992 : itsImages->divideModelByWeight( itsPBLimit, itsNormType );
404 : }
405 : else {
406 64 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
407 56 : { itsFacetImageStores[facet]->divideModelByWeight( itsPBLimit, itsNormType ); }
408 : }
409 1000 : }
410 :
411 929 : void SynthesisNormalizer::multiplyModelByWeight()
412 : {
413 : // LogIO os( LogOrigin("SynthesisNormalizer", "multiplyModelByWeight",WHERE) );
414 929 : if( itsNFacets==1) {
415 921 : itsImages->multiplyModelByWeight( itsPBLimit , itsNormType );
416 : }
417 : else {
418 64 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
419 56 : { itsFacetImageStores[facet]->multiplyModelByWeight( itsPBLimit , itsNormType); }
420 : }
421 929 : }
422 :
423 :
424 0 : std::shared_ptr<SIImageStore> SynthesisNormalizer::getImageStore()
425 : {
426 0 : LogIO os( LogOrigin("SynthesisNormalizer", "getImageStore", WHERE) );
427 0 : return itsImages;
428 0 : }
429 :
430 0 : void SynthesisNormalizer::setImageStore( SIImageStore* imstore )
431 : {
432 0 : LogIO os( LogOrigin("SynthesisNormalizer", "setImageStore", WHERE) );
433 0 : itsImages.reset( imstore );
434 0 : }
435 :
436 912 : void SynthesisNormalizer::setImageStore( std::shared_ptr<SIImageStore>& imstore )
437 : {
438 1824 : LogIO os( LogOrigin("SynthesisNormalizer", "setImageStore", WHERE) );
439 912 : itsImages= imstore ;
440 :
441 912 : }
442 :
443 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
444 : //// Internal Functions start here. These are not visible to the tool layer.
445 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
446 :
447 :
448 2390 : Bool SynthesisNormalizer::setupImagesOnDisk()
449 : {
450 4780 : LogIO os( LogOrigin("SynthesisNormalizer","setupImagesOnDisk",WHERE) );
451 :
452 2390 : Bool needToGatherImages=false;
453 :
454 2390 : String err("");
455 :
456 : // Check if full images exist, and open them if possible.
457 2390 : Bool foundFullImage=false;
458 : try
459 : {
460 2390 : itsImages = makeImageStore( itsImageName );
461 2390 : foundFullImage = true;
462 : }
463 0 : catch(AipsError &x)
464 : {
465 : //throw( AipsError("Error in constructing a Deconvolver : "+x.getMesg()) );
466 0 : err = err += String(x.getMesg()) + "\n";
467 0 : foundFullImage = false;
468 0 : }
469 :
470 2390 : os << LogIO::DEBUG2 << " Found full images : " << foundFullImage << LogIO::POST;
471 :
472 : // Check if part images exist
473 2390 : Bool foundPartImages = itsPartImageNames.nelements()>0 ? true : false ;
474 2390 : itsPartImages.resize( itsPartImageNames.nelements() );
475 :
476 2390 : for ( uInt part=0; part < itsPartImageNames.nelements() ; part++ )
477 : {
478 : try
479 : {
480 0 : itsPartImages[part] = makeImageStore ( itsPartImageNames[part] );
481 0 : foundPartImages |= true;
482 : }
483 0 : catch(AipsError &x)
484 : {
485 : //throw( AipsError("Error in constructing a Deconvolver : "+x.getMesg()) );
486 0 : err = err += String(x.getMesg()) + "\n";
487 0 : foundPartImages = false;
488 0 : }
489 : }
490 :
491 2390 : os << LogIO::DEBUG2 << " Found part images : " << foundPartImages << LogIO::POST;
492 :
493 2390 : if( foundPartImages == false)
494 : {
495 2390 : if( foundFullImage == true && itsPartImageNames.nelements()>0 )
496 : {
497 : // Pick the coordsys, etc from fullImage, and construct new/fresh partial images.
498 0 : os << LogIO::DEBUG2 << "Found full image, but no partial images. Make partImStores for : " << itsPartImageNames << LogIO::POST;
499 :
500 0 : String imopen = itsImages->getName()+".residual"+((itsMapperType=="multiterm")?".tt0":"");
501 0 : Directory imdir( imopen );
502 0 : if( ! imdir.exists() )
503 : {
504 0 : imopen = itsImages->getName()+".psf"+((itsMapperType=="multiterm")?".tt0":"");
505 0 : Directory imdir2( imopen );
506 0 : if( ! imdir2.exists() )
507 0 : throw(AipsError("Cannot find partial image psf or residual for " +itsImages->getName() +err));
508 0 : }
509 :
510 0 : PagedImage<Float> temppart(imopen);
511 :
512 0 : Bool useweightimage = itsImages->getUseWeightImage( *(itsImages->sumwt()) );
513 0 : for( uInt part=0; part<itsPartImageNames.nelements(); part++ )
514 : {
515 0 : itsPartImages[part] = makeImageStore (itsPartImageNames[part], temppart,
516 0 : useweightimage);
517 : }
518 0 : foundPartImages = True;
519 0 : }
520 : else
521 : {
522 2390 : itsPartImages.resize(0);
523 2390 : foundPartImages = False;
524 : }
525 : }
526 : else // Check that all have the same shape.
527 : {
528 0 : AlwaysAssert( itsPartImages.nelements() > 0 , AipsError );
529 0 : IPosition tempshape = itsPartImages[0]->getShape();
530 0 : for( uInt part=1; part<itsPartImages.nelements(); part++ )
531 : {
532 0 : if( tempshape != itsPartImages[part]->getShape() )
533 : {
534 0 : throw( AipsError("Shapes of partial images to be combined, do not match" + err) );
535 : }
536 : }
537 0 : }
538 :
539 :
540 : // Make sure all images exist and are consistent with each other. At the end, itsImages should be valid
541 2390 : if( foundPartImages == true ) // Partial Images exist. Check that 'full' exists, and do the gather.
542 : {
543 0 : if ( foundFullImage == true ) // Full image exists. Just check that shapes match with parts.
544 : {
545 0 : os << LogIO::DEBUG2 << "Partial and Full images exist. Checking if part images have the same shape as the full image : ";
546 0 : IPosition fullshape = itsImages->getShape();
547 :
548 0 : for ( uInt part=0; part < itsPartImages.nelements() ; part++ )
549 : {
550 0 : IPosition partshape = itsPartImages[part]->getShape();
551 0 : if( partshape != fullshape )
552 : {
553 0 : os << LogIO::DEBUG2<< "NO" << LogIO::POST;
554 0 : throw( AipsError("Shapes of the partial and full images on disk do not match. Cannot gather" + err) );
555 : }
556 0 : }
557 0 : os << LogIO::DEBUG2 << "Yes" << LogIO::POST;
558 :
559 0 : }
560 : else // Full image does not exist. Need to make it, using the shape and coords of part[0]
561 : {
562 0 : os << LogIO::DEBUG2 << "Only partial images exist. Need to make full images" << LogIO::POST;
563 :
564 0 : AlwaysAssert( itsPartImages.nelements() > 0, AipsError );
565 :
566 : // Find an image to open and pick csys,shape from.
567 0 : String imopen = itsPartImageNames[0]+".residual"+((itsMapperType=="multiterm")?".tt0":"");
568 0 : Directory imdir( imopen );
569 0 : if( ! imdir.exists() )
570 : {
571 0 : imopen = itsPartImageNames[0]+".psf"+((itsMapperType=="multiterm")?".tt0":"");
572 0 : Directory imdir2( imopen );
573 0 : if( ! imdir2.exists() )
574 : {
575 0 : imopen = itsPartImageNames[0]+".gridwt";
576 0 : Directory imdir3( imopen );
577 0 : if( ! imdir3.exists() )
578 0 : throw(AipsError("Cannot find partial image psf or residual or gridwt for " + itsPartImageNames[0]+err));
579 0 : }
580 :
581 0 : }
582 :
583 0 : PagedImage<Float> temppart( imopen );
584 :
585 0 : Bool useweightimage = itsPartImages[0]->getUseWeightImage( *(itsPartImages[0]->sumwt()) );
586 0 : itsImages = makeImageStore (itsImageName, temppart, useweightimage);
587 0 : foundFullImage = true;
588 0 : }
589 :
590 : // By now, all partial images and the full images exist on disk, and have the same shape.
591 0 : needToGatherImages=true;
592 :
593 : }
594 : else // No partial images supplied. Operating only with full images.
595 : {
596 2390 : if ( foundFullImage == true )
597 : {
598 2390 : os << LogIO::DEBUG2 << "Full images exist : " << itsImageName << LogIO::POST;
599 : }
600 : else // No full image on disk either. Error.
601 : {
602 0 : throw( AipsError("No images named " + itsImageName + " found on disk. No partial images found either."+err) );
603 : }
604 : }
605 :
606 : // Remove ?
607 2390 : itsImages->psf();
608 2390 : itsImages->validate();
609 :
610 : // Set up facet Imstores..... if needed
611 2390 : if( itsNFacets>1 )
612 : {
613 :
614 : // First, make sure that full images have been allocated before trying to make references.....
615 : // if( ! itsImages->checkValidity(true/*psf*/, true/*res*/,true/*wgt*/,true/*model*/,false/*image*/,false/*mask*/,true/*sumwt*/ ) )
616 : // { throw(AipsError("Internal Error : Invalid ImageStore for " + itsImages->getName())); }
617 :
618 : // Array<Float> ttt;
619 : // (itsImages->sumwt())->get(ttt);
620 : // cout << "SUMWT full : " << ttt << endl;
621 :
622 20 : itsFacetImageStores.resize( itsNFacets*itsNFacets );
623 172 : for( uInt facet=0; facet<itsNFacets*itsNFacets; ++facet )
624 : {
625 152 : itsFacetImageStores[facet] = itsImages->getSubImageStore(facet, itsNFacets);
626 :
627 : //Array<Float> qqq;
628 : //itsFacetImageStores[facet]->sumwt()->get(qqq);
629 : // cout << "SUMWTs for " << facet << " : " << qqq << endl;
630 :
631 : }
632 : }
633 2390 : os << LogIO::DEBUG2 << "Need to Gather ? " << needToGatherImages << LogIO::POST;
634 2390 : itsImages->releaseLocks();
635 2390 : return needToGatherImages;
636 2390 : }// end of setupImagesOnDisk
637 :
638 :
639 2708 : std::shared_ptr<SIImageStore> SynthesisNormalizer::makeImageStore(const String &imagename , const bool useweightimage)
640 : {
641 : //The constructors use ignoresumwt so use the !useweightimage
642 2708 : if( itsMapperType == "multiterm" )
643 458 : { return std::shared_ptr<SIImageStore>(new SIImageStoreMultiTerm( imagename, itsNTaylorTerms, true, !useweightimage )); }
644 : else
645 2250 : { return std::shared_ptr<SIImageStore>(new SIImageStore( imagename, true /*ignorefacets*/, !useweightimage )); }
646 : itsImages->releaseLocks();
647 : }
648 :
649 :
650 : /**
651 : * build a new ImageStore, whether SIImageStore or SIImageStoreMultiTerm, borrowing
652 : * image information from one partial image.
653 : *
654 : * @param imagename image name for the new SIStorageManager
655 : * @param part partial image from which miscinfo, etc. data will be borrowed
656 : * @param useweightimage useweight option for the new SIStorageManager
657 : *
658 : * @return A new SIImageStore object for the image name given.
659 : */
660 0 : std::shared_ptr<SIImageStore> SynthesisNormalizer::makeImageStore(const String &imagename,
661 : const PagedImage<Float> &part,
662 : Bool useweightimage)
663 : {
664 : // borrow shape, coord, imageinfo and miscinfo
665 0 : auto shape = part.shape();
666 0 : auto csys = part.coordinates();
667 0 : auto objectname = part.imageInfo().objectName();
668 0 : auto miscinfo = part.miscInfo();
669 0 : if( itsMapperType == "multiterm" )
670 : {
671 : std::shared_ptr<SIImageStore> multiTermStore =
672 0 : std::make_shared<SIImageStoreMultiTerm>(imagename, csys, shape, objectname,
673 0 : miscinfo, itsNFacets, false, itsNTaylorTerms, useweightimage );
674 0 : return multiTermStore;
675 0 : }
676 : else
677 : {
678 : return std::make_shared<SIImageStore>(imagename, csys, shape, objectname, miscinfo,
679 0 : false, useweightimage);
680 : }
681 0 : }
682 :
683 : //
684 : //---------------------------------------------------------------
685 : //
686 4 : void SynthesisNormalizer::setPsfFromOneFacet()
687 : {
688 : // Copy the PSF from one facet to the center of the full image, to use for the minor cycle
689 : //
690 : // This code segment will work for single and multi-term
691 : // - for single term, the index will always be 0, and SIImageStore's access functions know this.
692 : // - for multi-term, the index will be the taylor index and SIImageStoreMultiTerm knows this.
693 :
694 : /// itsImages
695 : /// itsFacetImageStores
696 :
697 : {
698 4 : IPosition shape=(itsImages->psf(0))->shape();
699 4 : IPosition blc(4, 0, 0, 0, 0);
700 4 : IPosition trc=shape-1;
701 10 : for(uInt tix=0; tix<2 * itsImages->getNTaylorTerms() - 1; tix++)
702 : {
703 12 : TempImage<Float> onepsf((itsFacetImageStores[0]->psf(tix))->shape(),
704 18 : (itsFacetImageStores[0]->psf(tix))->coordinates());
705 6 : onepsf.copyData(*(itsFacetImageStores[0]->psf(tix)));
706 : //now set the original to 0 as we have a copy of one facet psf
707 6 : (itsImages->psf(tix))->set(0.0);
708 6 : blc[0]=(shape[0]-(onepsf.shape()[0]))/2;
709 6 : trc[0]=onepsf.shape()[0]+blc[0]-1;
710 6 : blc[1]=(shape[1]-(onepsf.shape()[1]))/2;
711 6 : trc[1]=onepsf.shape()[1]+blc[1]-1;
712 6 : Slicer sl(blc, trc, Slicer::endIsLast);
713 12 : SubImage<Float> sub(*(itsImages->psf(tix)), sl, true);
714 6 : sub.copyData(onepsf);
715 6 : }
716 4 : }
717 :
718 : // cout << "In setPsfFromOneFacet : sumwt : " << itsImages->sumwt()->get() << endl;
719 :
720 4 : }
721 :
722 :
723 :
724 : } //# NAMESPACE CASA - END
725 :
|
