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 17 : SynthesisNormalizer::SynthesisNormalizer() :
62 17 : itsImages(std::shared_ptr<SIImageStore>()),
63 17 : itsPartImages(Vector<std::shared_ptr<SIImageStore> >()),
64 17 : itsImageName(""),
65 17 : itsPartImageNames(Vector<String>(0)),
66 17 : itsPBLimit(0.2),
67 17 : itsMapperType("default"),
68 17 : itsNTaylorTerms(1),
69 17 : itsNFacets(1),
70 51 : itsIsSingleDish(False)
71 : {
72 17 : itsFacetImageStores.resize(0);
73 17 : itsPBLimit = 0.35;
74 17 : }
75 :
76 17 : SynthesisNormalizer::~SynthesisNormalizer()
77 : {
78 34 : LogIO os( LogOrigin("SynthesisNormalizer","destructor",WHERE) );
79 17 : os << LogIO::DEBUG1 << "SynthesisNormalizer destroyed" << LogIO::POST;
80 17 : SynthesisUtilMethods::getResource("End SynthesisNormalizer");
81 :
82 17 : }
83 :
84 :
85 17 : void SynthesisNormalizer::setupNormalizer(Record normpars)
86 : {
87 34 : LogIO os( LogOrigin("SynthesisNormalizer","setupNormalizer",WHERE) );
88 :
89 : try {
90 17 : if ( normpars.isDefined("psfcutoff") ) { // A single string
91 17 : normpars.get( RecordFieldId("psfcutoff") , itsPsfcutoff );
92 : } else {
93 0 : itsPsfcutoff=0.35;
94 : }
95 :
96 17 : if ( normpars.isDefined("imagename") ) { // A single string
97 17 : itsImageName = normpars.asString( RecordFieldId("imagename"));
98 : } else {
99 0 : throw AipsError("imagename not specified");
100 : }
101 :
102 17 : if ( normpars.isDefined("partimagenames") ) { // A vector of strings
103 0 : normpars.get( RecordFieldId("partimagenames") , itsPartImageNames );
104 : } else {
105 17 : itsPartImageNames.resize(0);
106 : }
107 :
108 17 : if ( normpars.isDefined("pblimit") ) {
109 17 : normpars.get( RecordFieldId("pblimit") , itsPBLimit );
110 : } else {
111 0 : itsPBLimit = 0.2;
112 : }
113 :
114 17 : if ( normpars.isDefined("normtype") ) { // A single string
115 17 : itsNormType = normpars.asString( RecordFieldId("normtype"));
116 : } else {
117 0 : itsNormType = "flatnoise"; // flatnoise, flatsky
118 : }
119 :
120 : { // Debug comments
121 : // cout << "Chosen normtype : " << itsNormType << endl;
122 :
123 : // For multi-term choices.
124 : // 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 :
133 17 : if ( normpars.isDefined("deconvolver") ) { // A single string
134 17 : String dec = normpars.asString( RecordFieldId("deconvolver"));
135 17 : if (dec == "mtmfs") { itsMapperType="multiterm"; }
136 17 : else itsMapperType="default";
137 17 : } else {
138 0 : itsMapperType = "default";
139 : }
140 :
141 17 : if ( normpars.isDefined("nterms") ) { // A single int
142 17 : itsNTaylorTerms = normpars.asuInt( RecordFieldId("nterms")); }
143 : else {
144 0 : itsNTaylorTerms = 1;
145 : }
146 :
147 17 : if ( normpars.isDefined("facets") ) { // A single int
148 17 : itsNFacets = normpars.asuInt( RecordFieldId("facets"));
149 : } else {
150 0 : itsNFacets = 1;
151 : }
152 :
153 17 : if ( normpars.isDefined("restoringbeam") ) {
154 17 : if (normpars.dataType("restoringbeam") == TpString) {
155 0 : itsUseBeam = normpars.asString( RecordFieldId("restoringbeam") );
156 : } else {
157 17 : itsUseBeam = "";
158 : }
159 : } // else: keep itsUseBeam unchanged
160 :
161 17 : if ( normpars.isDefined("makesingledishnormalizer") ) { // A single bool
162 17 : itsIsSingleDish =
163 17 : normpars.asBool( RecordFieldId("makesingledishnormalizer") );
164 : } else {
165 0 : itsIsSingleDish = False;
166 : }
167 :
168 : }
169 0 : catch (AipsError &x) {
170 0 : throw AipsError(
171 0 : "Error in reading gather/scatter parameters: " + x.getMesg()
172 0 : );
173 0 : }
174 :
175 17 : } // end of setupNormalizer
176 :
177 :
178 34 : void SynthesisNormalizer::gatherImages(Bool dopsf, Bool doresidual, Bool dodensity)
179 : {
180 : // cout << " partimagenames :" << itsPartImageNames << endl;
181 :
182 34 : Bool needToGatherImages = setupImagesOnDisk();
183 : //if(dopsf)
184 : // cerr << "GATHER " << itsImages->psf()->shape() << endl;
185 :
186 34 : if( needToGatherImages )
187 : {
188 0 : LogIO os( LogOrigin("SynthesisNormalizer", "gatherImages",WHERE) );
189 :
190 0 : AlwaysAssert( itsPartImages.nelements()>0 , AipsError );
191 : //Bool doresidual = !dopsf;
192 0 : Bool doweight = dopsf ; //|| ( doresidual && ! itsImages->hasSensitivity() );
193 : // Bool doweight = dopsf || ( doresidual && itsImages->getUseWeightImage(*(itsPartImages[0]->residual())) );
194 :
195 0 : os << "Gather "<< (doresidual?"residual":"") << ( (dopsf&&doresidual)?",":"")
196 0 : << (dopsf?"psf":"") << ( (dopsf&&doweight)?",":"")
197 : << (doweight?"weight":"")<< " images : " << itsPartImageNames
198 0 : << " onto :" << itsImageName << LogIO::POST;
199 :
200 : // Add intelligence to modify all only the first time, but later, only residual;
201 0 : itsImages->resetImages( /*psf*/dopsf, /*residual*/doresidual, /*weight*/doweight );
202 :
203 0 : for( uInt part=0;part<itsPartImages.nelements();part++)
204 : {
205 0 : itsImages->addImages( itsPartImages[part], /*psf*/dopsf, /*residual*/doresidual, /*weight*/doweight, /*griddedwt*/dodensity );
206 0 : itsPartImages[part]->releaseLocks();
207 : }
208 :
209 0 : }// end of image gathering.
210 :
211 : // Normalize by the weight image.
212 : // divideResidualByWeight();
213 34 : itsImages->releaseLocks();
214 :
215 34 : }// end of gatherImages
216 :
217 17 : void SynthesisNormalizer::gatherPB()
218 : {
219 17 : if( itsPartImageNames.nelements()>0 )
220 : {
221 :
222 : try
223 : {
224 0 : itsPartImages[0] = makeImageStore( itsPartImageNames[0] );
225 : }
226 0 : catch(AipsError &x)
227 : {
228 0 : throw(AipsError("Cannot construct ImageStore for "+itsPartImageNames[0]+". "+x.what()));
229 0 : }
230 :
231 : try{
232 0 : LatticeExpr<Float> thepb( *(itsPartImages[0]->pb()) );
233 0 : LatticeLocker lock1(*(itsImages->pb()), FileLocker::Write);
234 0 : itsImages->pb()->copyData(thepb);
235 :
236 0 : }
237 0 : catch(AipsError &x)
238 : {
239 0 : throw(AipsError("Cannot copy the PB : "+x.getMesg()));
240 0 : }
241 :
242 : }
243 17 : }
244 :
245 17 : void SynthesisNormalizer::scatterModel()
246 : {
247 :
248 34 : LogIO os( LogOrigin("SynthesisNormalizer", "scatterModel",WHERE) );
249 :
250 17 : setupImagesOnDisk(); // To open up and initialize itsPartImages.
251 :
252 : // os << "In ScatterModel : " << itsPartImages.nelements() << " for " << itsPartImageNames << LogIO::POST;
253 :
254 17 : if( itsPartImages.nelements() > 0 ) // && itsImages->doesImageExist(modelName) )
255 : {
256 0 : os << "Send the model from : " << itsImageName << " to all nodes :" << itsPartImageNames << LogIO::POST;
257 :
258 : // Make the list of model images. This list is of length >1 only for multi-term runs.
259 : // Vector<String> modelNames( itsImages->getNTaylorTerms() );
260 : // if( modelNames.nelements() ==1 ) modelNames[0] = itsImages->getName()+".model";
261 : // if( modelNames.nelements() > 1 )
262 : // {
263 : // for( uInt nt=0;nt<itsImages->getNTaylorTerms();nt++)
264 : // modelNames[nt] = itsImages->getName()+".model.tt" + String::toString(nt);
265 : // }
266 :
267 :
268 0 : Vector<String> modelNames( itsImages->getNTaylorTerms() );
269 0 : if( itsImages->getType()=="default" ) modelNames[0] = itsImages->getName()+".model";
270 0 : if( itsImages->getType()=="multiterm" )
271 : {
272 0 : for( uInt nt=0;nt<itsImages->getNTaylorTerms();nt++)
273 0 : modelNames[nt] = itsImages->getName()+".model.tt" + String::toString(nt);
274 : }
275 :
276 :
277 :
278 0 : for( uInt part=0;part<itsPartImages.nelements();part++)
279 : {
280 0 : itsPartImages[part]->setModelImage( modelNames );
281 0 : itsPartImages[part]->releaseLocks();
282 : }
283 0 : itsImages->releaseLocks();
284 0 : }
285 :
286 17 : }// end of scatterModel
287 :
288 : //----------------
289 :
290 0 : void SynthesisNormalizer::gatherWeightDensity(){
291 :
292 0 : Bool needToGatherImages = setupImagesOnDisk();
293 : // if(dopsf)
294 : // cerr << "GATHER " << itsImages->psf()->shape() << endl;
295 :
296 0 : if (needToGatherImages) {
297 0 : LogIO os(LogOrigin("SynthesisNormalizer", "gatherWeightDensity", WHERE));
298 :
299 0 : itsImages->gridwt()->set(0.0);
300 0 : Record iminfo(itsImages->gridwt()->miscInfo());
301 0 : Vector<Float> d2out;
302 0 : Vector<Float> f2out;
303 0 : if(iminfo.isDefined("d2")){
304 0 : iminfo.get("d2", d2out);
305 0 : iminfo.get("f2", f2out);
306 : }
307 0 : Vector<Float> f2part(f2out.nelements(), 0.0);
308 0 : bool recalcF2 = (d2out.nelements() >0 && d2out[0] == 1.0); // need to recalculate F2 for normalized Briggs
309 0 : if(recalcF2)
310 0 : f2out.set(0.0);
311 0 : Int nxim = itsImages->gridwt()->shape()(0);
312 0 : Int nyim = itsImages->gridwt()->shape()(1);
313 0 : for (uInt part = 0; part < itsPartImages.nelements(); part++) {
314 0 : itsImages->addImages(itsPartImages[part], false,
315 : /*residual*/ false, /*weight*/ false,
316 : /*griddedwt*/ true);
317 0 : Record partiminfo(itsPartImages[part]->gridwt()->miscInfo());
318 0 : mergeWeightDensityInfo(iminfo, partiminfo);
319 0 : if (recalcF2) {
320 0 : Vector<Float> f2in;
321 0 : (itsPartImages[part]->gridwt()->miscInfo()).get("f2", f2in);
322 0 : Array<Float> arr;
323 0 : if (d2out.nelements() > 1) {
324 0 : IPosition blc(5, 0, 0, 0, 0, 0);
325 0 : IPosition shp(5, nxim, nyim, 1, 1, 1);
326 0 : for (uInt k = 0; k < d2out.nelements(); ++k) {
327 0 : blc[4] = k;
328 0 : itsPartImages[part]->gridwt()->getSlice(arr, blc, shp, True);
329 0 : f2in[k] = f2in[k] * sum(arr * arr);
330 : }
331 0 : }
332 : else{ //not multifield weighting
333 0 : itsPartImages[part]->gridwt()->get(arr, True);
334 0 : f2in[0] = f2in[0] * sum(arr * arr);
335 : }
336 0 : f2out += f2in;
337 0 : }
338 0 : itsPartImages[part]->releaseLocks();
339 0 : }
340 0 : if(recalcF2){
341 0 : Array<Float> arr;
342 0 : if (d2out.nelements() > 1) {
343 0 : IPosition blc(5, 0, 0, 0, 0, 0);
344 0 : IPosition shp(5, nxim, nyim, 1, 1, 1);
345 0 : for (uInt k = 0; k < d2out.nelements(); ++k) {
346 0 : blc[4] = k;
347 0 : itsImages->gridwt()->getSlice(arr, blc, shp, True);
348 0 : Float sumarr2 = sum(arr * arr);
349 0 : if(sumarr2> 0)
350 0 : f2out[k] = f2out[k] / sumarr2;
351 : }
352 0 : }
353 : else{
354 0 : itsImages->gridwt()->get(arr, True);
355 0 : Float sumarr2 = sum(arr * arr);
356 0 : if (sumarr2 > 0)
357 0 : f2out[0] = f2out[0] / sumarr2;
358 : }
359 0 : iminfo.define("f2", f2out);
360 :
361 0 : }
362 0 : itsImages->gridwt()->setMiscInfo(iminfo);
363 :
364 0 : } // end of image gathering.
365 :
366 : // Normalize by the weight image.
367 : // divideResidualByWeight();
368 0 : itsImages->releaseLocks();
369 :
370 0 : }//end of gatherweightdensity
371 :
372 : //-------------------------
373 :
374 :
375 0 : string SynthesisNormalizer::scatterWeightDensity()
376 : {
377 :
378 0 : LogIO os( LogOrigin("SynthesisNormalizer", "scatterWeightDensity",WHERE) );
379 :
380 0 : setupImagesOnDisk(); // To open up and initialize itsPartImages.
381 0 : string weightname = "";
382 : // os << "In ScatterModel : " << itsPartImages.nelements() << " for " <<
383 : // itsPartImageNames << LogIO::POST;
384 :
385 0 : if( itsPartImages.nelements() > 0 )
386 : {
387 0 : os << "Send the gridded weight from : " << itsImageName << " to all nodes :" << itsPartImageNames << LogIO::POST;
388 :
389 0 : for( uInt part=0;part<itsPartImages.nelements();part++)
390 : {
391 0 : itsPartImages[part]->setWeightDensity( itsImages );
392 0 : itsPartImages[part]->releaseLocks();
393 : }
394 0 : weightname=itsImages->gridwt()->name();
395 0 : itsImages->releaseLocks();
396 : }
397 0 : return weightname;
398 0 : } // end of gatherImages
399 :
400 17 : void SynthesisNormalizer::divideResidualByWeight()
401 : {
402 34 : LogIO os( LogOrigin("SynthesisNormalizer", "divideResidualByWeight",WHERE) );
403 :
404 17 : if( itsNFacets==1) {
405 17 : itsImages->divideResidualByWeight( itsPBLimit, itsNormType );
406 : }
407 : else {
408 0 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
409 0 : { itsFacetImageStores[facet]->divideResidualByWeight( itsPBLimit , itsNormType ); }
410 : }
411 17 : itsImages->releaseLocks();
412 :
413 17 : }
414 :
415 0 : void SynthesisNormalizer::divideResidualByWeightSD()
416 : {
417 0 : LogIO os( LogOrigin("SynthesisNormalizer", "divideResidualByWeightSD",WHERE) );
418 :
419 0 : if( itsNFacets==1) {
420 0 : itsImages->divideResidualByWeightSD( itsPBLimit );
421 : }
422 : else {
423 0 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
424 0 : { itsFacetImageStores[facet]->divideResidualByWeightSD( itsPBLimit ); }
425 : }
426 0 : itsImages->releaseLocks();
427 :
428 0 : }
429 0 : void SynthesisNormalizer::makePSFBeamset(){
430 0 : LogIO os(LogOrigin("SynthesisNormalizer", "dividePSFByWeight", WHERE));
431 : try{
432 0 : if(!itsImages){
433 0 : itsImages = makeImageStore( itsImageName, false );
434 : }
435 :
436 :
437 : }
438 0 : catch(AipsError& x){
439 :
440 0 : throw(AipsError("Programmers error no psf is made on disk and trying to fit"));
441 0 : }
442 0 : itsImages->makeImageBeamSet(itsPsfcutoff);
443 0 : itsImages->releaseLocks();
444 0 : }
445 :
446 :
447 17 : void SynthesisNormalizer::divideWeightBySumWt() {
448 :
449 34 : LogIO os(LogOrigin("SynthesisNormalizer", "divideWeightBySumWt", WHERE));
450 17 : itsImages->divideWeightBySumwt();
451 :
452 17 : }
453 :
454 :
455 :
456 17 : void SynthesisNormalizer::dividePSFByWeight()
457 : {
458 34 : LogIO os( LogOrigin("SynthesisNormalizer", "dividePSFByWeight",WHERE) );
459 : {
460 :
461 17 : LatticeLocker lock1 (*(itsImages->psf()), FileLocker::Write);
462 :
463 17 : if( itsNFacets==1) {
464 17 : itsImages->dividePSFByWeight(itsPBLimit);
465 : }
466 : else {
467 : // Since PSFs are normed just by their max, this sequence is OK.
468 0 : setPsfFromOneFacet();
469 0 : itsImages->dividePSFByWeight(itsPBLimit);
470 : }
471 17 : }// release of lock1 otherwise releaselock below will cause it to core
472 : //dump as the psf pointer goes dangling
473 : // Check PSF quality by fitting beams
474 : {
475 17 : itsImages->calcSensitivity();
476 :
477 17 : itsImages->makeImageBeamSet(itsPsfcutoff);
478 17 : Bool verbose(False);
479 17 : if (itsUseBeam=="common") verbose=True;
480 17 : itsImages->printBeamSet(verbose);
481 : }
482 17 : itsImages->releaseLocks();
483 :
484 17 : }
485 :
486 17 : void SynthesisNormalizer::normalizePrimaryBeam()
487 : {
488 34 : LogIO os( LogOrigin("SynthesisNormalizer", "normalizePrimaryBeam",WHERE) );
489 :
490 17 : if( itsImages==NULL ) { itsImages = makeImageStore( itsImageName ); }
491 :
492 17 : gatherPB();
493 :
494 : // Irrespective of facets.
495 17 : itsImages->normalizePrimaryBeam(itsPBLimit);
496 17 : }
497 :
498 :
499 17 : void SynthesisNormalizer::divideModelByWeight()
500 : {
501 34 : LogIO os( LogOrigin("SynthesisNormalizer", "divideModelByWeight",WHERE) );
502 17 : if( ! itsImages )
503 : {
504 : //os << LogIO::WARN << "No imagestore yet. Trying to construct, so that the starting model can be divided by wt if needed...." << LogIO::POST;
505 :
506 : try
507 : {
508 0 : itsImages = makeImageStore( itsImageName );
509 : }
510 0 : catch(AipsError &x)
511 : {
512 0 : throw(AipsError("Cannot construct ImageStore for "+itsImageName+". "+x.what()));
513 0 : }
514 : // return;
515 : }
516 17 : if( itsNFacets==1) {
517 17 : itsImages->divideModelByWeight( itsPBLimit, itsNormType );
518 : }
519 : else {
520 0 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
521 0 : { itsFacetImageStores[facet]->divideModelByWeight( itsPBLimit, itsNormType ); }
522 : }
523 17 : }
524 :
525 17 : void SynthesisNormalizer::multiplyModelByWeight()
526 : {
527 : // LogIO os( LogOrigin("SynthesisNormalizer", "multiplyModelByWeight",WHERE) );
528 17 : if( itsNFacets==1) {
529 17 : itsImages->multiplyModelByWeight( itsPBLimit , itsNormType );
530 : }
531 : else {
532 0 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; facet++ )
533 0 : { itsFacetImageStores[facet]->multiplyModelByWeight( itsPBLimit , itsNormType); }
534 : }
535 17 : }
536 :
537 :
538 0 : std::shared_ptr<SIImageStore> SynthesisNormalizer::getImageStore()
539 : {
540 0 : LogIO os( LogOrigin("SynthesisNormalizer", "getImageStore", WHERE) );
541 0 : return itsImages;
542 0 : }
543 :
544 0 : void SynthesisNormalizer::setImageStore( SIImageStore* imstore )
545 : {
546 0 : LogIO os( LogOrigin("SynthesisNormalizer", "setImageStore", WHERE) );
547 0 : itsImages.reset( imstore );
548 0 : }
549 :
550 0 : void SynthesisNormalizer::setImageStore( std::shared_ptr<SIImageStore>& imstore )
551 : {
552 0 : LogIO os( LogOrigin("SynthesisNormalizer", "setImageStore", WHERE) );
553 0 : itsImages= imstore ;
554 :
555 0 : }
556 :
557 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
558 : //// Internal Functions start here. These are not visible to the tool layer.
559 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
560 :
561 :
562 51 : Bool SynthesisNormalizer::setupImagesOnDisk() {
563 102 : LogIO os( LogOrigin("SynthesisNormalizer","setupImagesOnDisk",WHERE) );
564 :
565 51 : Bool needToGatherImages = false;
566 :
567 51 : String err("");
568 :
569 : // Check if full images exist, and open them if possible.
570 51 : Bool foundFullImage = false;
571 : try {
572 51 : itsImages = makeImageStore( itsImageName );
573 51 : foundFullImage = true;
574 : }
575 0 : catch (AipsError &x) {
576 : //throw( AipsError("Error in constructing a Deconvolver : "+x.getMesg()) );
577 0 : err = err += String(x.getMesg()) + "\n";
578 0 : foundFullImage = false;
579 0 : }
580 :
581 : os << LogIO::DEBUG2
582 : << " Found full images : " << foundFullImage
583 51 : << LogIO::POST;
584 :
585 : // Check if part images exist
586 51 : Bool foundPartImages = itsPartImageNames.nelements() > 0 ? true : false ;
587 51 : itsPartImages.resize( itsPartImageNames.nelements() );
588 :
589 51 : for ( uInt part=0; part < itsPartImageNames.nelements() ; part++ ) {
590 : try {
591 0 : itsPartImages[part] = makeImageStore ( itsPartImageNames[part] );
592 0 : foundPartImages |= true;
593 : }
594 0 : catch (AipsError &x) {
595 : // throw AipsError(
596 : // "Error in constructing a Deconvolver : " + x.getMesg()
597 : // );
598 0 : err = err += String(x.getMesg()) + "\n";
599 0 : foundPartImages = false;
600 0 : }
601 : }
602 :
603 : os << LogIO::DEBUG2
604 : << " Found part images : " << foundPartImages
605 51 : << LogIO::POST;
606 :
607 51 : if ( not foundPartImages ) {
608 51 : if ( foundFullImage and itsPartImageNames.nelements() > 0 ) {
609 : // Pick the coordsys, etc from fullImage,
610 : // and construct new/fresh partial images.
611 : os << LogIO::DEBUG2
612 : << "Found full image, but no partial images."
613 : << " Making partImStores for : " << itsPartImageNames
614 0 : << LogIO::POST;
615 :
616 0 : String imopen = itsImages->getName() + ".residual"
617 0 : + ( (itsMapperType == "multiterm") ? ".tt0" : "" );
618 0 : Directory imdir( imopen );
619 0 : if ( not imdir.exists() ) {
620 0 : imopen = itsImages->getName() + ".psf"
621 0 : + ( (itsMapperType == "multiterm") ? ".tt0" : "" );
622 0 : Directory imdir2( imopen );
623 0 : if ( not imdir2.exists() )
624 0 : throw AipsError(
625 : "Cannot find partial image psf or residual for "
626 0 : + itsImages->getName() + err
627 0 : );
628 0 : }
629 :
630 0 : PagedImage<Float> temppart(imopen);
631 :
632 : Bool useweightimage =
633 0 : itsImages->getUseWeightImage( *(itsImages->sumwt()) );
634 0 : for ( uInt part=0; part<itsPartImageNames.nelements(); part++ ) {
635 0 : itsPartImages[part] = makeImageStore(
636 0 : itsPartImageNames[part], temppart, useweightimage);
637 : }
638 0 : foundPartImages = True;
639 0 : } else {
640 51 : itsPartImages.resize(0);
641 51 : foundPartImages = False;
642 : }
643 : } else { // Check that all have the same shape.
644 0 : AlwaysAssert( itsPartImages.nelements() > 0 , AipsError );
645 0 : IPosition tempshape = itsPartImages[0]->getShape();
646 0 : for ( uInt part=1; part<itsPartImages.nelements(); part++ ) {
647 0 : if ( tempshape != itsPartImages[part]->getShape() ) {
648 0 : throw AipsError(
649 0 : "Shapes of partial images to be combined, do not match" + err
650 0 : );
651 : }
652 : }
653 0 : }
654 :
655 : // Make sure all images exist and are consistent with each other.
656 : // At the end, itsImages should be valid
657 51 : if ( foundPartImages ) { // Partial Images exist. Check that 'full' exists,
658 : // and do the gather.
659 0 : if ( foundFullImage ) { // Full image exists.
660 : // Just check that shapes match with parts.
661 : os << LogIO::DEBUG2
662 : << "Partial and Full images exist."
663 0 : << "Checking if part images have the same shape as the full image : ";
664 :
665 0 : IPosition fullshape = itsImages->getShape();
666 0 : for ( uInt part=0; part < itsPartImages.nelements() ; part++ ) {
667 0 : IPosition partshape = itsPartImages[part]->getShape();
668 0 : if ( partshape != fullshape ) {
669 : os << LogIO::DEBUG2
670 0 : << "NO" << LogIO::POST;
671 0 : throw AipsError(
672 : "Shapes of the partial and full images on disk do not match."
673 0 : " Cannot gather" + err
674 0 : );
675 : }
676 0 : }
677 0 : os << LogIO::DEBUG2 << "Yes" << LogIO::POST;
678 :
679 0 : } else { // Full image does not exist.
680 : // Need to make it, using the shape and coords of part[0]
681 : os << LogIO::DEBUG2
682 : << "Only partial images exist. Need to make full images"
683 0 : << LogIO::POST ;
684 :
685 0 : AlwaysAssert( itsPartImages.nelements() > 0, AipsError );
686 :
687 : // Find an image to open and pick csys,shape from.
688 0 : String imopen = itsPartImageNames[0] + ".residual"
689 0 : + ( (itsMapperType == "multiterm") ? ".tt0" : "" );
690 0 : Directory imdir( imopen );
691 0 : if ( not imdir.exists() ) {
692 0 : imopen = itsPartImageNames[0] + ".psf"
693 0 : + ( (itsMapperType == "multiterm") ? ".tt0" : "" );
694 0 : Directory imdir2( imopen );
695 0 : if ( not imdir2.exists() ) {
696 0 : imopen = itsPartImageNames[0] + ".gridwt";
697 0 : Directory imdir3( imopen );
698 0 : if ( not imdir3.exists() )
699 0 : throw AipsError(
700 : "Cannot find partial image psf or residual or gridwt for "
701 0 : + itsPartImageNames[0] + err
702 0 : );
703 0 : }
704 0 : }
705 :
706 0 : PagedImage<Float> temppart( imopen );
707 :
708 : Bool useweightimage =
709 0 : itsPartImages[0]->getUseWeightImage( *(itsPartImages[0]->sumwt()) );
710 0 : itsImages = makeImageStore (itsImageName, temppart, useweightimage);
711 0 : foundFullImage = true;
712 0 : }
713 :
714 : // By now, all partial images and the full images exist on disk,
715 : // and have the same shape.
716 0 : needToGatherImages = true;
717 : } else { // No partial images supplied. Operating only with full images.
718 :
719 51 : if ( foundFullImage ) {
720 : os << LogIO::DEBUG2
721 51 : << "Full images exist : " << itsImageName
722 51 : << LogIO::POST;
723 : }
724 : else { // No full image on disk either. Error.
725 0 : throw AipsError(
726 0 : "No images named " + itsImageName
727 0 : + " found on disk. No partial images found either." + err
728 0 : );
729 : }
730 : }
731 :
732 : // Remove ?
733 51 : itsImages->psf();
734 :
735 51 : itsImages->validate();
736 :
737 : // Set up facet Imstores..... if needed
738 51 : if ( itsNFacets > 1 ) {
739 : // First, make sure that full images have been allocated
740 : // before trying to make references.....
741 : // if ( not itsImages->checkValidity(
742 : // true/*psf*/,
743 : // true/*res*/,
744 : // true/*wgt*/,
745 : // true/*model*/,
746 : // false/*image*/,
747 : // false/*mask*/,
748 : // true/*sumwt*/ ) ) {
749 : // throw AipsError(
750 : // "Internal Error : Invalid ImageStore for "
751 : // + itsImages->getName())
752 : // );
753 : // }
754 :
755 : // Array<Float> ttt;
756 : // (itsImages->sumwt())->get(ttt);
757 : // cout << "SUMWT full : " << ttt << endl;
758 :
759 0 : itsFacetImageStores.resize( itsNFacets*itsNFacets );
760 0 : for ( uInt facet=0; facet<itsNFacets*itsNFacets; ++facet ) {
761 0 : itsFacetImageStores[facet] =
762 0 : itsImages->getSubImageStore(facet, itsNFacets);
763 : // Array<Float> qqq;
764 : // itsFacetImageStores[facet]->sumwt()->get(qqq);
765 : // cout << "SUMWTs for " << facet << " : " << qqq << endl;
766 : }
767 : }
768 :
769 : os << LogIO::DEBUG2
770 : << "Need to Gather ? " << needToGatherImages
771 51 : << LogIO::POST;
772 :
773 51 : itsImages->releaseLocks();
774 :
775 51 : return needToGatherImages;
776 51 : } // end of setupImagesOnDisk
777 :
778 :
779 51 : std::shared_ptr<SIImageStore> SynthesisNormalizer::makeImageStore(
780 : const String &imagename,
781 : const bool useweightimage)
782 : {
783 : // The constructors use ignoresumwt so use the !useweightimage
784 51 : if ( itsMapperType == "multiterm" ) {
785 : return std::shared_ptr<SIImageStore>(
786 : new SIImageStoreMultiTerm(
787 : imagename,
788 : itsNTaylorTerms,
789 : /* ignorefacets= */ true,
790 0 : /* ignoresumwt= */ not useweightimage
791 0 : )
792 0 : );
793 : } else {
794 : return std::shared_ptr<SIImageStore>(
795 : new SIImageStore(
796 : imagename,
797 : /* ignorefacets= */ true,
798 51 : /* noRequireSumwt= */ not useweightimage,
799 51 : /* makeSingleDishStore= */ itsIsSingleDish
800 51 : )
801 51 : );
802 : }
803 : itsImages->releaseLocks();
804 : }
805 :
806 :
807 : /**
808 : * build a new ImageStore, whether SIImageStore or SIImageStoreMultiTerm, borrowing
809 : * image information from one partial image.
810 : *
811 : * @param imagename image name for the new SIStorageManager
812 : * @param part partial image from which miscinfo, etc. data will be borrowed
813 : * @param useweightimage useweight option for the new SIStorageManager
814 : *
815 : * @return A new SIImageStore object for the image name given.
816 : */
817 0 : std::shared_ptr<SIImageStore> SynthesisNormalizer::makeImageStore(const String &imagename,
818 : const PagedImage<Float> &part,
819 : Bool useweightimage)
820 : {
821 : // borrow shape, coord, imageinfo and miscinfo
822 0 : auto shape = part.shape();
823 0 : auto csys = part.coordinates();
824 0 : auto objectname = part.imageInfo().objectName();
825 0 : auto miscinfo = part.miscInfo();
826 0 : if (itsMapperType == "multiterm") {
827 : std::shared_ptr<SIImageStore> multiTermStore =
828 0 : std::make_shared<SIImageStoreMultiTerm>(
829 0 : imagename, csys, shape, objectname, miscinfo, itsNFacets, false,
830 0 : itsNTaylorTerms, useweightimage);
831 0 : return multiTermStore;
832 0 : } else {
833 : return std::make_shared<SIImageStore>(imagename, csys, shape, objectname,
834 0 : miscinfo, false, useweightimage);
835 : }
836 0 : }
837 :
838 : //
839 : //---------------------------------------------------------------
840 : //
841 0 : void SynthesisNormalizer::setPsfFromOneFacet()
842 : {
843 : // Copy the PSF from one facet to the center of the full image, to use for the minor cycle
844 : //
845 : // This code segment will work for single and multi-term
846 : // - for single term, the index will always be 0, and SIImageStore's access functions know this.
847 : // - for multi-term, the index will be the taylor index and SIImageStoreMultiTerm knows this.
848 :
849 : /// itsImages
850 : /// itsFacetImageStores
851 :
852 : {
853 0 : IPosition shape=(itsImages->psf(0))->shape();
854 0 : IPosition blc(4, 0, 0, 0, 0);
855 0 : IPosition trc=shape-1;
856 0 : for(uInt tix=0; tix<2 * itsImages->getNTaylorTerms() - 1; tix++)
857 : {
858 0 : TempImage<Float> onepsf((itsFacetImageStores[0]->psf(tix))->shape(),
859 0 : (itsFacetImageStores[0]->psf(tix))->coordinates());
860 0 : onepsf.copyData(*(itsFacetImageStores[0]->psf(tix)));
861 : //now set the original to 0 as we have a copy of one facet psf
862 0 : (itsImages->psf(tix))->set(0.0);
863 0 : blc[0]=(shape[0]-(onepsf.shape()[0]))/2;
864 0 : trc[0]=onepsf.shape()[0]+blc[0]-1;
865 0 : blc[1]=(shape[1]-(onepsf.shape()[1]))/2;
866 0 : trc[1]=onepsf.shape()[1]+blc[1]-1;
867 0 : Slicer sl(blc, trc, Slicer::endIsLast);
868 0 : SubImage<Float> sub(*(itsImages->psf(tix)), sl, true);
869 0 : sub.copyData(onepsf);
870 0 : }
871 0 : }
872 :
873 : // cout << "In setPsfFromOneFacet : sumwt : " << itsImages->sumwt()->get() << endl;
874 :
875 0 : }
876 :
877 0 : void SynthesisNormalizer::mergeWeightDensityInfo(Record& outinfo, const Record& partinfo){
878 0 : uInt mapsize=0;
879 0 : if(outinfo.isDefined("multimapsize"))
880 0 : outinfo.get("multimapsize", mapsize);
881 : uInt partsize;
882 0 : if(!partinfo.isDefined("multimapsize"))
883 0 : return;
884 0 : partinfo.get("multimapsize", partsize);
885 0 : std::vector<String> key2add;
886 0 : std::vector<Int> val2add;
887 0 : for (uInt k = 0; k < partsize; ++k){
888 0 : String key;
889 0 : partinfo.get("key" + String::toString(k), key);
890 0 : Bool hasKey = False;
891 0 : for (uInt j = 0; j < mapsize; ++j) {
892 0 : String imkey;
893 0 : outinfo.get("key" + String::toString(j), imkey);
894 0 : hasKey = hasKey || (imkey == key);
895 0 : }
896 0 : if (!hasKey) {
897 0 : key2add.push_back(key);
898 : Int val;
899 0 : partinfo.get("val" + String::toString(k), val);
900 0 : val2add.push_back(val);
901 : }
902 0 : }
903 0 : if (key2add.size() > 0) {
904 0 : for (uInt k = mapsize; k < (mapsize + key2add.size()); ++k) {
905 0 : String key = "key" + String::toString(k);
906 0 : outinfo.define(key, key2add[k - mapsize]);
907 0 : String val = "val" + String::toString(k);
908 0 : outinfo.define(val, val2add[k - mapsize]);
909 0 : }
910 0 : mapsize += key2add.size();
911 0 : outinfo.define("multimapsize", mapsize);
912 : }
913 0 : }
914 :
915 : } //# NAMESPACE CASA - END
916 :
|
