Line data Source code
1 : //# WBCleanImageSkyModel.cc: Implementation of WBCleanImageSkyModel class
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library 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 Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 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: WBCleanImageSkyModel.cc 13615 2010-12-20 14:04:00 UrvashiRV$
27 : //# v2.6 : Added psf-patch support to reduce memory footprint.
28 :
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 : #include <synthesis/MeasurementComponents/WBCleanImageSkyModel.h>
31 : #include <synthesis/MeasurementEquations/CubeSkyEquation.h>
32 : #include <casacore/casa/OS/File.h>
33 : #include <synthesis/MeasurementEquations/SkyEquation.h>
34 : #include <synthesis/TransformMachines/StokesImageUtil.h>
35 : #include <synthesis/MeasurementEquations/LatticeModel.h>
36 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
37 : #include <casacore/casa/Exceptions/Error.h>
38 : #include <casacore/casa/BasicSL/String.h>
39 : #include <casacore/casa/Utilities/Assert.h>
40 :
41 : #include <casacore/images/Images/PagedImage.h>
42 : #include <casacore/images/Images/SubImage.h>
43 : #include <casacore/images/Regions/ImageRegion.h>
44 : #include <casacore/images/Regions/RegionManager.h>
45 :
46 : #include <casacore/images/Regions/WCBox.h>
47 :
48 : #include <casacore/measures/Measures/Quality.h>
49 : #include <casacore/coordinates/Coordinates/QualityCoordinate.h>
50 : #include <casacore/images/Images/ImageUtilities.h>
51 :
52 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
53 :
54 : #include <msvis/MSVis/VisSet.h>
55 : #include <msvis/MSVis/VisSetUtil.h>
56 :
57 : #include <casacore/ms/MeasurementSets/MSColumns.h>
58 :
59 : #include <sstream>
60 :
61 : #include <casacore/casa/Logging/LogMessage.h>
62 : #include <casacore/casa/Logging/LogSink.h>
63 :
64 : #include <casacore/casa/OS/HostInfo.h>
65 :
66 :
67 : using namespace casacore;
68 : namespace casa { //# NAMESPACE CASA - BEGIN
69 : #define TMR(a) "[User: " << a.user() << "] [System: " << a.system() << "] [Real: " << a.real() << "]"
70 :
71 : #define MIN(a,b) ((a)<=(b) ? (a) : (b))
72 : #define MAX(a,b) ((a)>=(b) ? (a) : (b))
73 :
74 : /*************************************
75 : * Constructor
76 : *************************************/
77 0 : WBCleanImageSkyModel::WBCleanImageSkyModel()
78 : {
79 0 : initVars();
80 0 : nscales_p=4;
81 0 : ntaylor_p=2;
82 0 : refFrequency_p=1.42e+09;
83 0 : scaleSizes_p.resize(0);
84 0 : };
85 0 : WBCleanImageSkyModel::WBCleanImageSkyModel(const Int ntaylor,const Int nscales,const Double reffreq)
86 : {
87 0 : initVars();
88 0 : nscales_p=nscales;
89 0 : ntaylor_p=ntaylor;
90 0 : refFrequency_p=reffreq;
91 0 : scaleSizes_p.resize(0);
92 0 : };
93 0 : WBCleanImageSkyModel::WBCleanImageSkyModel(const Int ntaylor,const Vector<Float>& userScaleSizes,const Double reffreq)
94 : {
95 0 : initVars();
96 0 : if(userScaleSizes.nelements()==0) os << "No scales set !" << LogIO::WARN;
97 0 : if(userScaleSizes[0]!=0.0)
98 : {
99 0 : nscales_p=userScaleSizes.nelements()+1;
100 0 : scaleSizes_p.resize(nscales_p);
101 0 : for(Int i=1;i<nscales_p;i++) scaleSizes_p[i] = userScaleSizes[i-1];
102 0 : scaleSizes_p[0]=0.0;
103 : }
104 : else
105 : {
106 0 : nscales_p=userScaleSizes.nelements();
107 0 : scaleSizes_p.resize(nscales_p);
108 0 : for(Int i=0;i<nscales_p;i++) scaleSizes_p[i] = userScaleSizes[i];
109 : }
110 0 : ntaylor_p=ntaylor;
111 0 : refFrequency_p=reffreq;
112 0 : };
113 :
114 0 : void WBCleanImageSkyModel::initVars()
115 : {
116 0 : adbg=0;
117 0 : ddbg=1; // output per iteration
118 0 : tdbg=0;
119 :
120 0 : modified_p=true;
121 0 : memoryMB_p = Double(HostInfo::memoryTotal(true)/1024)/(2.0);
122 0 : donePSF_p=false;
123 0 : doneMTMCinit_p=false;
124 :
125 0 : numbermajorcycles_p=0;
126 0 : nfields_p=1;
127 0 : lc_p.resize(0);
128 :
129 0 : os = LogIO(LogOrigin("WBCleanImageSkyModel","solve"));
130 :
131 0 : setAlgorithm("MSMFS");
132 :
133 0 : }
134 :
135 : /*************************************
136 : * Destructor
137 : *************************************/
138 0 : WBCleanImageSkyModel::~WBCleanImageSkyModel()
139 : {
140 0 : lc_p.resize(0);
141 : //cout << "WBCleanImageSkyModel destructor " << endl;
142 :
143 : /*
144 : if(nmodels_p > numberOfModels())
145 : {
146 : resizeWorkArrays(numberOfModels());
147 : nmodels_p = numberOfModels();
148 : }
149 : */
150 :
151 0 : };
152 :
153 : /*************************************
154 : * Solver
155 : *************************************/
156 0 : Bool WBCleanImageSkyModel::solve(SkyEquation& se)
157 : {
158 0 : os << "MSMFS algorithm (v2.6) with " << ntaylor_p << " Taylor coefficients and Reference Frequency of " << refFrequency_p << " Hz" << LogIO::POST;
159 0 : Int stopflag=0;
160 0 : Int nchan=0,npol=0;
161 :
162 : /* Gather shape information */
163 0 : nmodels_p = numberOfModels();
164 :
165 0 : if(nmodels_p % ntaylor_p != 0)
166 : {
167 0 : os << "Incorrect number of input models " << LogIO::EXCEPTION;
168 0 : os << "NModels = N_fields x N_taylor" << LogIO::EXCEPTION;
169 0 : AlwaysAssert((nmodels_p % ntaylor_p == 0), AipsError);
170 : }
171 :
172 : /* Check supplied bandwidth-ratio and print warnings if needed */
173 0 : checkParameters();
174 :
175 : /* Calc the number of fields */
176 0 : nfields_p = nmodels_p/ntaylor_p;
177 : ///// NOTE : Code is written with loops on 'fields' to support a future implementation
178 : ///// Disable it for now, since it has not been tested.
179 : //AlwaysAssert(nfields_p==1, AipsError);
180 :
181 : //cout << "Number of fields : " << nfields_p << endl;
182 :
183 : /* Current restriction : one pol and one chan-plane */
184 0 : for(Int model=0;model<nmodels_p;model++)
185 : {
186 0 : nx = image(model).shape()(0);
187 0 : ny = image(model).shape()(1);
188 0 : npol=image(model).shape()(2);
189 0 : nchan=image(model).shape()(3);
190 0 : if(nchan > 1) os << "Cannot process more than one output channel" << LogIO::EXCEPTION;
191 0 : if(npol > 1) os << "Cannot process more than one output polarization" << LogIO::EXCEPTION;
192 0 : AlwaysAssert((nchan==1), AipsError);
193 0 : AlwaysAssert((npol==1), AipsError);
194 : }
195 :
196 : //------------------------------- For 'active' ---------------------------------------------------------------///
197 : /* Calculate the initial residual image for all models. */
198 0 : if( se.isNewFTM() == false )
199 : {
200 0 : if(!donePSF_p)
201 : {
202 0 : os << "Calculating initial residual images..." << LogIO::POST;
203 0 : solveResiduals(se,(numberIterations()<1)?true:False);
204 : }
205 : else
206 : {
207 : /*
208 : if(numbermajorcycles_p>0)
209 : {
210 : os << "RE-Calculating residual images because previous residuals have been modified in-place during restoration to be 'coefficient residuals'." << LogIO::POST;
211 : solveResiduals(se,(numberIterations()<1)?true:False);
212 : }
213 : */
214 : }
215 : }
216 : //------------------------------- For 'active' ---------------------------------------------------------------///
217 :
218 : /* Create the Point Spread Functions */
219 0 : if(!donePSF_p)
220 : {
221 : /* Resize the work arrays to calculate extra PSFs */
222 0 : Int original_nmodels = nmodels_p;
223 0 : nmodels_p = original_nmodels + nfields_p * (ntaylor_p - 1);
224 0 : resizeWorkArrays(nmodels_p);
225 :
226 : try
227 : {
228 : /* Make the 2N-1 PSFs */
229 0 : os << "Calculating spectral PSFs..." << LogIO::POST;
230 0 : makeSpectralPSFs(se, numberIterations()<0?true:False);
231 : }
232 0 : catch(AipsError &x)
233 : {
234 : /* Resize the work arrays to normal size - the destructors use 'nmodels_p' on other lists */
235 0 : nmodels_p = original_nmodels;
236 0 : resizeWorkArrays(nmodels_p);
237 0 : os << "Could not make PSFs. Please check image co-ordinate system : " << x.getMesg() << LogIO::EXCEPTION;
238 0 : }
239 :
240 0 : if(adbg) cout << "Shape of lc_p : " << lc_p.nelements() << endl;
241 : /* Initialize MTMC, allocate memory, and send in all 2N-1 psfs */
242 0 : if(lc_p.nelements()==0 && numberIterations()>=0)
243 : {
244 0 : lc_p.resize(nfields_p);
245 0 : Bool state=true;
246 : /* Initialize the MultiTermMatrixCleaners */
247 0 : for(Int thismodel=0;thismodel<nfields_p;thismodel++)
248 : {
249 0 : lc_p[thismodel].setscales(scaleSizes_p);
250 0 : lc_p[thismodel].setntaylorterms(ntaylor_p);
251 0 : nx = image(thismodel).shape()(0);
252 0 : ny = image(thismodel).shape()(1);
253 0 : state &= lc_p[thismodel].initialise(nx,ny); // allocates memory once....
254 : }
255 0 : if( !state ) // initialise will return false if there is any internal inconsistency with settings so far.
256 : {
257 0 : lc_p.resize(0);
258 0 : nmodels_p = original_nmodels;
259 0 : resizeWorkArrays(nmodels_p);
260 0 : os << LogIO::SEVERE << "Could not initialize MS-MFS minor cycle" << LogIO::EXCEPTION;
261 0 : return false; // redundant
262 : }
263 :
264 : /* Send all 2N-1 PSFs into the MultiTermLatticeCleaner */
265 0 : for(Int thismodel=0;thismodel<nfields_p;thismodel++)
266 : {
267 0 : for (Int order=0;order<2*ntaylor_p-1;order++)
268 : {
269 : // This should be doing a reference only. Make sure of this.
270 0 : Matrix<Float> tempMat;
271 0 : Array<Float> tempArr;
272 0 : (PSF(getModelIndex(thismodel,order))).get(tempArr,true);
273 0 : tempMat.reference(tempArr);
274 :
275 0 : lc_p[thismodel].setpsf( order , tempMat );
276 0 : }
277 : }
278 0 : doneMTMCinit_p = true;
279 : }
280 :
281 : /* Resize the work arrays to normal size - for residual comps, etc. */
282 0 : nmodels_p = original_nmodels;
283 0 : resizeWorkArrays(nmodels_p);
284 : ///// TODO : Make sure this releases memory as it is supposed to.
285 :
286 0 : donePSF_p=true;
287 : }
288 :
289 :
290 : //------------------------------ For new FTMs --------------------------------------------------------------//
291 : ///// Consider doing lc_p.setmodel and getmodel instead of dividing/multiplying the avgPB in NewMTFT::initializeToVis
292 0 : if( se.isNewFTM() == true )
293 : {
294 : /* Calculate the initial residual image for all models. */
295 0 : if( numbermajorcycles_p==0 )
296 : {
297 0 : os << "Calculating initial residual images..." << LogIO::POST;
298 0 : solveResiduals(se,(numberIterations()<1)?true:False);
299 : }
300 :
301 : }
302 : //------------------------------- For new FTMs -------------------------------------------------------------//
303 :
304 : /* Return if niter=0 */
305 : /* Check if this is an interactive-clean run, or if niter=0 */
306 0 : if(adbg) cout << "NumberIterations - before any cycles: " << numberIterations() << endl;
307 0 : if(numberIterations() < 1)
308 : {
309 0 : return true;
310 : }
311 :
312 : /* Check that lc_p's have been initialized by now.. */
313 0 : if(doneMTMCinit_p == false)
314 : {
315 0 : os << LogIO::SEVERE << "MultiTermMatrixCleaners are un-initialized, perhaps because of a previous im.clean(niter=-1) call. Please close imager, re-open it, and run with niter>=0" << LogIO::POST;
316 : }
317 :
318 : /* Set up the Mask image */
319 0 : for(Int thismodel=0;thismodel<nfields_p;thismodel++)
320 : {
321 0 : if(hasMask(getModelIndex(thismodel,0)))
322 : {
323 0 : if(adbg) os << "Sending in the mask for lc_p : " << thismodel << LogIO::POST;
324 :
325 0 : Matrix<Float> tempMat;
326 0 : Array<Float> tempArr;
327 0 : (mask(getModelIndex(thismodel,0))).get(tempArr,true);
328 0 : tempMat.reference(tempArr);
329 :
330 0 : lc_p[thismodel].setmask( tempMat );
331 0 : }
332 : }
333 :
334 : /******************* START MAJOR CYCLE LOOP *****************/
335 : /* Logic for number of iterations in the minor-cycle vs major cycle vs total,
336 : how they relate to convergence thresholds, and how they change with
337 : interactive vs non-interactive use (user specified total niter, vs niters per
338 : major cycle...) is a bit of a mess. Needs cleanup. Right now, hard-coded for
339 : nfields=1 */
340 0 : previous_maxresidual_p=1e+10;
341 0 : Int index=0;
342 0 : Float fractionOfPsf=0.1;
343 0 : Vector<Int> iterationcount(nfields_p); iterationcount=0;
344 0 : Vector<Int> moreiterations(nfields_p); moreiterations=0;
345 0 : Vector<Int> thiscycleniter(nfields_p); thiscycleniter=0;
346 0 : for(Int itercountmaj=0;itercountmaj<1000;itercountmaj++)
347 : {
348 0 : numbermajorcycles_p ++;
349 0 : os << "**** Major Cycle " << numbermajorcycles_p << LogIO::POST;
350 0 : thiscycleniter=0;
351 : /* Compute stopping threshold for this major cycle */
352 0 : stopflag = static_cast<casacore::Int>(computeFluxLimit(fractionOfPsf));
353 : /* If the peak residual is already less than the user-threshold, stop */
354 0 : if(stopflag==1) break;
355 : /* If we detect divergence across major cycles, stop */
356 0 : if(stopflag==-1) break;
357 :
358 0 : for(Int thismodel=0;thismodel<nfields_p;thismodel++) // For now, nfields_p=1 always.
359 : {
360 :
361 0 : if( !isSolveable(getModelIndex(thismodel,0)) )
362 : {
363 : // This field is not to be cleaned in this set of minor-cycle iterations
364 0 : continue;
365 : }
366 :
367 : /* Number of iterations left to do */
368 0 : moreiterations[thismodel] = numberIterations() - iterationcount[thismodel];
369 :
370 : /* If all iterations are done for this run, stop - for all fields*/
371 0 : if(moreiterations[thismodel] <=0 ) {stopflag=-1;break;}
372 :
373 : /* Send in the current model and residual */
374 0 : for (Int order=0;order<ntaylor_p;order++)
375 : {
376 0 : index = getModelIndex(thismodel,order);
377 :
378 0 : Matrix<Float> tempMat;
379 0 : Array<Float> tempArr;
380 :
381 0 : (residual(index)).get(tempArr,true);
382 0 : tempMat.reference(tempArr);
383 0 : lc_p[thismodel].setresidual(order,tempMat);
384 :
385 0 : (image(index)).get(tempArr,true);
386 0 : tempMat.reference(tempArr);
387 0 : lc_p[thismodel].setmodel(order,tempMat);
388 0 : }
389 :
390 : /* Deconvolve */
391 : /* Return value is the number of minor-cycle iterations done */
392 0 : os << "Starting Minor Cycle iterations for field : " << thismodel << LogIO::POST;
393 0 : thiscycleniter[thismodel] = lc_p[thismodel].mtclean(moreiterations[thismodel], fractionOfPsf, gain(), threshold());
394 :
395 : /* A signal for a singular Hessian : stop */
396 0 : if(thiscycleniter[thismodel]==-2) { stopflag=-2; break;}
397 :
398 : /* A signal for convergence with no iterations */
399 : /* One field may have converged, while others go on... so 'continue' */
400 : // if(thiscycleniter[thismodel]==0) { stopflag=1; break; }
401 0 : if(thiscycleniter[thismodel]==0) { stopflag=1; continue; }
402 :
403 : ///* A signal for divergence. Save current model and stop (force a major cycle). */
404 : //if(thiscycleniter==-1) { stopflag=1; }
405 :
406 : /* Increment the minor-cycle iteration counter */
407 0 : iterationcount[thismodel] += thiscycleniter[thismodel];
408 :
409 : /* Get out the updated model */
410 0 : for (Int order=0;order<ntaylor_p;order++)
411 : {
412 0 : index = getModelIndex(thismodel,order);
413 0 : Matrix<Float> tempMod;
414 0 : lc_p[thismodel].getmodel(order,tempMod);
415 0 : image(index).put(tempMod);
416 0 : }
417 :
418 : }// end of model loop
419 :
420 0 : if(adbg) cout << "end of major cycle : " << itercountmaj << " -> iterationcount : " << iterationcount << " thiscycleniter : " << thiscycleniter << " stopflag : " << stopflag << endl;
421 :
422 : /* Exit without further ado if MTMC cannot invert matrices */
423 0 : if(stopflag == -2)
424 : {
425 0 : os << "Cannot invert Multi-Term Hessian matrix. Please check the reference-frequency and ensure that the number of frequency-channels in the selected data >= nterms" << LogIO::WARN << LogIO::POST;
426 0 : break;
427 : }
428 :
429 : /* If reached 1000 major cycles - assume something is wrong */
430 0 : if(itercountmaj==999) os << " Reached the allowed maximum of 1000 major cycles " << LogIO::POST;
431 :
432 : /* Do the prediction and residual computation for all models. */
433 : /* Exit if we're already at the user-specified flux threshold, or reached numberIterations() */
434 0 : if( abs(stopflag)==1 || max(iterationcount) >= numberIterations() || itercountmaj==999)
435 : {
436 0 : os << "Calculating final residual images..." << LogIO::POST;
437 : /* Call 'solveResiduals' with modelToMS = true to write the model to the MS */
438 0 : solveResiduals(se,true);
439 0 : break;
440 : }
441 : else
442 : {
443 0 : os << "Calculating new residual images..." << LogIO::POST;
444 0 : solveResiduals(se);
445 : }
446 :
447 : /*
448 : if( se.isNewFTM() == true )
449 : {
450 : // The model will get PB-corrected in-place before prediction.
451 : // This needs to be reset to what the preceeding minor cycle got,
452 : // for incremental additions in the next cycle.
453 : saveCurrentModels();
454 : }
455 : */
456 :
457 : }
458 : /******************* END MAJOR CYCLE LOOP *****************/
459 :
460 : /* Compute and write alpha,beta results to disk */
461 : ///writeResultsToDisk();
462 : // calculateCoeffResiduals();
463 :
464 : /* stopflag=1 -> stopped because the user-threshold was reached */
465 : /* stopflag=-1 -> stopped because number of iterations was reached */
466 : /* stopflag=-2 -> stopped because of singular Hessian */
467 : /* stopflag=0 -> reached maximum number of major-cycle iterations !! */
468 0 : if(stopflag>0) return(true);
469 0 : else return(false);
470 0 : } // END OF SOLVE
471 :
472 :
473 0 : void WBCleanImageSkyModel::saveCurrentModels()
474 : {
475 :
476 0 : for(Int thismodel=0;thismodel<nfields_p;thismodel++)
477 : {
478 : /* Get out the updated model */
479 0 : for (Int order=0;order<ntaylor_p;order++)
480 : {
481 0 : Int index = getModelIndex(thismodel,order);
482 0 : Matrix<Float> tempMod;
483 0 : lc_p[thismodel].getmodel(order,tempMod);
484 0 : image(index).put(tempMod);
485 0 : }
486 : }// end of model loop
487 :
488 0 : }// end of saveCurrentModels
489 :
490 :
491 : /* Calculate stopping threshold for the current set of minor-cycle iterations */
492 : /* Note : The stopping threshold is computed from the residual image.
493 : However, MTMC checks on the peak of (residual * psf * scale_0).
494 : These values can differ slightly.
495 : TODO : Perhaps move this function back into MTMC.
496 : (It was brought out to conform to the other devonvolvers..)
497 : */
498 0 : Float WBCleanImageSkyModel::computeFluxLimit(Float &fractionOfPsf)
499 : {
500 0 : LogIO os(LogOrigin("WBCleanImageSkyModel", "computeFluxLimit", WHERE));
501 0 : Float maxres=0.0,maxval=0.0,minval=0.0;
502 0 : Vector<Float> fmaxres(nfields_p); fmaxres=0.0;
503 :
504 : /* Measure the peak residual across all fields */
505 0 : for(Int field=0;field<nfields_p;field++)
506 : {
507 0 : Int index = getModelIndex(field,0);
508 0 : Array<Float> tempArr;
509 0 : (residual(index)).get(tempArr,true);
510 0 : IPosition maxpos(tempArr.shape()),minpos(tempArr.shape());
511 0 : minMax(minval,maxval, minpos, maxpos,tempArr);
512 0 : fmaxres[field]=maxval;
513 0 : }
514 :
515 : /* Pick the peak residual across all fields */
516 0 : maxres = max(fmaxres);
517 0 : if(adbg) cout << "Peak Residual across fields (over all pixels) : " << maxres << endl;
518 :
519 : /* If we've already converged, return */
520 0 : if(maxres < threshold())
521 : {
522 0 : os << "Peak residual : " << maxres << " is lower than the user-specified stopping threshold : " << threshold() << LogIO::POST;
523 0 : return 1;
524 : }
525 :
526 : /* If we detect divergence across major cycles, warn or STOP */
527 0 : if(fabs( ((maxres - previous_maxresidual_p)/previous_maxresidual_p ) > 10.0 ))
528 : {
529 0 : os << "Peak residual : " << maxres << " has increased by more than a factor of 10 across major cycles. Could be diverging. Stopping" << LogIO::POST;
530 0 : return -1;
531 : }
532 0 : if(fabs( ((maxres - previous_maxresidual_p)/previous_maxresidual_p ) > 2.0 ))
533 : {
534 0 : os << "Peak residual : " << maxres << " has increased across major cycles. Could be diverging, but continuing..." << LogIO::POST;
535 : }
536 0 : previous_maxresidual_p = maxres;
537 :
538 : /* Find PSF sidelobe level */
539 0 : Float maxpsfside=0.0;
540 0 : for(uInt field=0;static_cast<Int>(field)<nfields_p;field++)
541 : {
542 0 : Int index = getModelIndex(field,0);
543 : /* abs(min of the PSF) will be treated as the max sidelobe level */
544 0 : Array<Float> tempArr;
545 0 : (PSF(index)).get(tempArr,true);
546 0 : IPosition maxpos(tempArr.shape()),minpos(tempArr.shape());
547 0 : minMax(minval,maxval, minpos, maxpos,tempArr);
548 0 : maxpsfside = max(maxpsfside, abs(minval));
549 0 : }
550 :
551 0 : fractionOfPsf = min(cycleMaxPsfFraction_p, cycleFactor_p * maxpsfside);
552 0 : if(fractionOfPsf>0.8) fractionOfPsf=0.8;
553 : // cyclethreshold = max(threshold(), fractionOfPsf * maxres);
554 :
555 : // if(adbg)
556 : {
557 0 : os << "Peak Residual (all pixels) : " << maxres << " User Threshold : " << threshold() << " Max PSF Sidelobe : " << abs(minval) << " User maxPsfFraction : " << cycleMaxPsfFraction_p << " User cyclefactor : " << cycleFactor_p << " fractionOfPsf = min(maxPsfFraction, PSFsidelobe x cyclefactor) : " << fractionOfPsf << LogIO::POST;
558 : // os << "Stopping threshold for this major cycle min(user threshold , fractionOfPsf x Max Residual) : " << cyclethreshold << endl;
559 : }
560 :
561 0 : return 0;
562 0 : }
563 :
564 : /***********************************************************************/
565 0 : Bool WBCleanImageSkyModel::calculateAlphaBeta(const Vector<String> &restoredNames,
566 : const Vector<String> &residualNames)
567 : {
568 0 : LogIO os(LogOrigin("WBCleanImageSkyModel", "calculateAlphaBeta", WHERE));
569 : //UNUSED: Int index=0;
570 0 : Bool writeerror=true;
571 :
572 :
573 0 : for(Int field=0;field<nfields_p;field++)
574 : {
575 0 : Int baseindex = getModelIndex(field,0);
576 0 : String alphaname,betaname, alphaerrorname;
577 0 : if( ( (restoredNames[baseindex]).substr( (restoredNames[baseindex]).length()-3 , 3 ) ).matches("tt0") )
578 : {
579 0 : alphaname = (restoredNames[baseindex]).substr(0,(restoredNames[baseindex]).length()-3) + "alpha";
580 0 : betaname = (restoredNames[baseindex]).substr(0,(restoredNames[baseindex]).length()-3) + "beta";
581 0 : if(writeerror) alphaerrorname = alphaname+".error";
582 : }
583 : else
584 : {
585 0 : alphaname = (restoredNames[baseindex]) + String(".alpha");
586 0 : betaname = (restoredNames[baseindex]) + String(".beta");
587 : }
588 :
589 : /* Create empty alpha image, and alpha error image */
590 0 : PagedImage<Float> imalpha(image(baseindex).shape(),image(baseindex).coordinates(),alphaname);
591 0 : imalpha.set(0.0);
592 :
593 : /* Open restored images */
594 0 : PagedImage<Float> imtaylor0(restoredNames[getModelIndex(field,0)]);
595 0 : PagedImage<Float> imtaylor1(restoredNames[getModelIndex(field,1)]);
596 :
597 : /* Open first residual image */
598 0 : PagedImage<Float> residual0(residualNames[getModelIndex(field,0)]);
599 :
600 : /* Create a mask - make this adapt to the signal-to-noise */
601 0 : LatticeExprNode leMaxRes=max(residual0);
602 0 : Float maxres = leMaxRes.getFloat();
603 : // Threshold is either 10% of the peak residual (psf sidelobe level) or
604 : // user threshold, if deconvolution has gone that deep.
605 0 : Float specthreshold = MAX( threshold()*5 , maxres/5.0 );
606 0 : os << "Calculating spectral parameters for Intensity > MAX(threshold*5,peakresidual/5) = " << specthreshold << " Jy/beam" << LogIO::POST;
607 0 : LatticeExpr<Float> mask1(iif((imtaylor0)>(specthreshold),1.0,0.0));
608 0 : LatticeExpr<Float> mask0(iif((imtaylor0)>(specthreshold),0.0,1.0));
609 :
610 : /////// Calculate alpha
611 0 : LatticeExpr<Float> alphacalc( ((imtaylor1)*mask1)/((imtaylor0)+(mask0)) );
612 0 : imalpha.copyData(alphacalc);
613 :
614 : // Set the restoring beam for alpha
615 0 : ImageInfo ii = imalpha.imageInfo();
616 0 : ii.setRestoringBeam( (imtaylor0.imageInfo()).restoringBeam() );
617 0 : imalpha.setImageInfo(ii);
618 : //imalpha.setUnits(Unit("Spectral Index"));
619 0 : imalpha.table().unmarkForDelete();
620 :
621 : // Make a mask for the alpha image
622 0 : LatticeExpr<Bool> lemask(iif((imtaylor0 > specthreshold) , true, false));
623 :
624 0 : createMask(lemask, imalpha);
625 0 : os << "Written Spectral Index Image : " << alphaname << LogIO::POST;
626 :
627 :
628 : ////// Calculate error on alpha
629 0 : if(writeerror)
630 : {
631 0 : PagedImage<Float> imalphaerror(image(baseindex).shape(),image(baseindex).coordinates(),alphaerrorname);
632 0 : imalphaerror.set(0.0);
633 0 : PagedImage<Float> residual1(residualNames[getModelIndex(field,1)]);
634 :
635 0 : LatticeExpr<Float> alphacalcerror( abs(alphacalc) * sqrt( ( (residual0*mask1)/(imtaylor0+mask0) )*( (residual0*mask1)/(imtaylor0+mask0) ) + ( (residual1*mask1)/(imtaylor1+mask0) )*( (residual1*mask1)/(imtaylor1+mask0) ) ) );
636 0 : imalphaerror.copyData(alphacalcerror);
637 0 : imalphaerror.setImageInfo(ii);
638 0 : createMask(lemask, imalphaerror);
639 0 : imalphaerror.table().unmarkForDelete();
640 0 : os << "Written Spectral Index Error Image : " << alphaerrorname << LogIO::POST;
641 :
642 : // mergeDataError( imalpha, imalphaerror, alphaerrorname+".new" );
643 :
644 0 : }
645 :
646 : ////// Calculate beta
647 0 : if(ntaylor_p>2)
648 : {
649 0 : PagedImage<Float> imbeta(image(baseindex).shape(),image(baseindex).coordinates(),betaname);
650 0 : imbeta.set(0.0);
651 0 : PagedImage<Float> imtaylor2(restoredNames[getModelIndex(field,2)]);
652 :
653 0 : LatticeExpr<Float> betacalc( ((imtaylor2)*mask1)/((imtaylor0)+(mask0))-0.5*(imalpha)*(imalpha-1.0) );
654 0 : imbeta.copyData(betacalc);
655 0 : imbeta.setImageInfo(ii);
656 : //imbeta.setUnits(Unit("Spectral Curvature"));
657 0 : createMask(lemask, imbeta);
658 0 : imbeta.table().unmarkForDelete();
659 :
660 0 : os << "Written Spectral Curvature Image : " << betaname << LogIO::POST;
661 0 : }
662 :
663 0 : }// field loop
664 :
665 0 : return 0;
666 :
667 0 : }
668 : /***********************************************************************/
669 :
670 0 : Bool WBCleanImageSkyModel::createMask(LatticeExpr<Bool> &lemask, ImageInterface<Float> &outimage)
671 : {
672 0 : ImageRegion outreg = outimage.makeMask("mask0",false,true);
673 0 : LCRegion& outmask=outreg.asMask();
674 0 : outmask.copyData(lemask);
675 0 : outimage.defineRegion("mask0",outreg, RegionHandler::Masks, true);
676 0 : outimage.setDefaultMask("mask0");
677 0 : return true;
678 0 : }
679 :
680 :
681 : /***********************************************************************/
682 0 : Bool WBCleanImageSkyModel::calculateCoeffResiduals()
683 : {
684 0 : for(Int field=0;field<nfields_p;field++)
685 : {
686 0 : Int baseindex = getModelIndex(field,0);
687 :
688 :
689 : // Send in the final residuals
690 0 : for (Int order=0;order<ntaylor_p;order++)
691 : {
692 0 : Int index = getModelIndex(field,order);
693 0 : Matrix<Float> tempMat;
694 0 : Array<Float> tempArr;
695 0 : (residual(index)).get(tempArr,true);
696 0 : tempMat.reference(tempArr);
697 0 : lc_p[baseindex].setresidual(order,tempMat);
698 0 : }
699 :
700 : // Compute principal solution in-place.
701 0 : lc_p[baseindex].computeprincipalsolution();
702 :
703 : // Get the new residuals
704 0 : for (Int order=0;order<ntaylor_p;order++)
705 : {
706 0 : Int index = getModelIndex(field,order);
707 0 : Matrix<Float> tempMod;
708 0 : lc_p[baseindex].getresidual(order,tempMod);
709 0 : residual(index).put(tempMod);
710 0 : }
711 :
712 :
713 : /*
714 :
715 : //Apply Inverse Hessian to the residuals
716 : IPosition gip(4,image(baseindex).shape()[0],image(baseindex).shape()[1],1,1);
717 : Matrix<Double> invhessian;
718 : lc_p[field].getinvhessian(invhessian);
719 : //cout << "Inverse Hessian : " << invhessian << endl;
720 :
721 : Int tindex;
722 : LatticeExprNode len_p;
723 : PtrBlock<TempLattice<Float>* > coeffresiduals(ntaylor_p); //,smoothresiduals(ntaylor_p);
724 : for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
725 : {
726 : coeffresiduals[taylor1] = new TempLattice<Float>(gip,memoryMB_p);
727 : }
728 :
729 : //Apply the inverse Hessian to the residuals
730 : for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
731 : {
732 : len_p = LatticeExprNode(0.0);
733 : for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
734 : {
735 : tindex = getModelIndex(field,taylor2);
736 : len_p = len_p + LatticeExprNode((Float)(invhessian)(taylor1,taylor2)*(residual(tindex)));
737 : }
738 : (*coeffresiduals[taylor1]).copyData(LatticeExpr<Float>(len_p));
739 : }
740 :
741 : //Fill in the residual images with these coefficient residuals
742 : for(Int taylor=0;taylor<ntaylor_p;taylor++)
743 : {
744 : tindex = getModelIndex(field,taylor);
745 : (residual(tindex)).copyData(LatticeExpr<Float>(*coeffresiduals[taylor]));
746 : }
747 :
748 : for(uInt i=0;i<coeffresiduals.nelements();i++)
749 : {
750 : if(coeffresiduals[i]) delete coeffresiduals[i];
751 : }
752 : */
753 :
754 :
755 : }//end of field loop
756 0 : os << "Converting final residuals to 'coefficient residuals', for restoration" << LogIO::POST;
757 0 : return true;
758 : }//end of calculateCoeffResiduals
759 :
760 : /***********************************************************************/
761 : ///// Write alpha and beta to disk. Calculate from smoothed model + residuals.
762 0 : Int WBCleanImageSkyModel::writeResultsToDisk()
763 : {
764 0 : if(ntaylor_p<=1) return 0;
765 :
766 0 : os << "Output Taylor-coefficient images : " << imageNames << LogIO::POST;
767 : // if(ntaylor_p==2) os << "Calculating Spectral Index" << LogIO::POST;
768 : // if(ntaylor_p>2) os << "Calculating Spectral Index and Curvature" << LogIO::POST;
769 :
770 :
771 0 : GaussianBeam beam;
772 0 : Int index=0;
773 :
774 0 : for(Int field=0;field<nfields_p;field++)
775 : {
776 0 : PtrBlock<TempLattice<Float>* > smoothed;
777 0 : if(ntaylor_p>2) smoothed.resize(3);
778 0 : else if(ntaylor_p==2) smoothed.resize(2);
779 :
780 0 : Int baseindex = getModelIndex(field,0);
781 0 : String alphaname,betaname;
782 0 : if( ( (imageNames[baseindex]).substr( (imageNames[baseindex]).length()-3 , 3 ) ).matches("tt0") )
783 : {
784 0 : alphaname = (imageNames[baseindex]).substr(0,(imageNames[baseindex]).length()-3) + "alpha";
785 0 : betaname = (imageNames[baseindex]).substr(0,(imageNames[baseindex]).length()-3) + "beta";
786 : }
787 : else
788 : {
789 0 : alphaname = (imageNames[baseindex]) + String(".alpha");
790 0 : betaname = (imageNames[baseindex]) + String(".beta");
791 : }
792 :
793 0 : StokesImageUtil::FitGaussianPSF(PSF(baseindex), beam);
794 :
795 : /* Create empty alpha image */
796 0 : PagedImage<Float> imalpha(image(baseindex).shape(),image(baseindex).coordinates(),alphaname);
797 0 : imalpha.set(0.0);
798 :
799 : /* Apply Inverse Hessian to the residuals */
800 0 : IPosition gip(4,image(baseindex).shape()[0],image(baseindex).shape()[1],1,1);
801 0 : Matrix<Double> invhessian;
802 0 : lc_p[field].getinvhessian(invhessian);
803 : //cout << "Inverse Hessian : " << invhessian << endl;
804 :
805 : Int tindex;
806 0 : LatticeExprNode len_p;
807 0 : PtrBlock<TempLattice<Float>* > coeffresiduals(ntaylor_p); //,smoothresiduals(ntaylor_p);
808 0 : for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
809 : {
810 0 : coeffresiduals[taylor1] = new TempLattice<Float>(gip,memoryMB_p);
811 : }
812 :
813 : /* Apply the inverse Hessian to the residuals */
814 0 : for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
815 : {
816 0 : len_p = LatticeExprNode(0.0);
817 0 : for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
818 : {
819 0 : tindex = getModelIndex(field,taylor2);
820 0 : len_p = len_p + LatticeExprNode((Float)(invhessian)(taylor1,taylor2)*(residual(tindex)));
821 : }
822 0 : (*coeffresiduals[taylor1]).copyData(LatticeExpr<Float>(len_p));
823 : }
824 :
825 : /* Fill in the residual images with these coefficient residuals */
826 0 : for(Int taylor=0;taylor<ntaylor_p;taylor++)
827 : {
828 0 : tindex = getModelIndex(field,taylor);
829 0 : (residual(taylor)).copyData(LatticeExpr<Float>(*coeffresiduals[taylor]));
830 : }
831 :
832 :
833 : /* Smooth the model images and add the above coefficient residuals */
834 0 : for(uInt i=0;i<smoothed.nelements();i++)
835 : {
836 0 : smoothed[i] = new TempLattice<Float>(gip,memoryMB_p);
837 :
838 0 : index = getModelIndex(field,i);
839 0 : LatticeExpr<Float> cop(image(index));
840 0 : imalpha.copyData(cop);
841 0 : StokesImageUtil::Convolve(imalpha, beam);
842 : //cout << "Clean Beam from WBC : " << bmaj << " , " << bmin << " , " << bpa << endl;
843 : //cout << "SkyModel internally-recorded beam for index " << index << " : " << beam(index) << endl;
844 : //LatticeExpr<Float> le(imalpha);
845 0 : LatticeExpr<Float> le(imalpha+( *coeffresiduals[i] ));
846 0 : (*smoothed[i]).copyData(le);
847 0 : }
848 :
849 :
850 : /* Create a mask - make this adapt to the signal-to-noise */
851 0 : LatticeExprNode leMaxRes=max(residual( getModelIndex(field,0) ));
852 0 : Float maxres = leMaxRes.getFloat();
853 : // Threshold is either 10% of the peak residual (psf sidelobe level) or
854 : // user threshold, if deconvolution has gone that deep.
855 0 : Float specthreshold = MAX( threshold()*5 , maxres/5.0 );
856 0 : os << "Calculating spectral parameters for Intensity > MAX(threshold*5,peakresidual/5) = " << specthreshold << " Jy/beam" << LogIO::POST;
857 0 : LatticeExpr<Float> mask1(iif((*smoothed[0])>(specthreshold),1.0,0.0));
858 0 : LatticeExpr<Float> mask0(iif((*smoothed[0])>(specthreshold),0.0,1.0));
859 :
860 : /* Calculate alpha and beta */
861 0 : LatticeExpr<Float> alphacalc( ((*smoothed[1])*mask1)/((*smoothed[0])+(mask0)) );
862 0 : imalpha.copyData(alphacalc);
863 :
864 0 : ImageInfo ii = imalpha.imageInfo();
865 0 : ii.setRestoringBeam(beam);
866 :
867 0 : imalpha.setImageInfo(ii);
868 : //imalpha.setUnits(Unit("Spectral Index"));
869 0 : imalpha.table().unmarkForDelete();
870 0 : os << "Written Spectral Index Image : " << alphaname << LogIO::POST;
871 :
872 0 : if(ntaylor_p>2)
873 : {
874 0 : PagedImage<Float> imbeta(image(baseindex).shape(),image(baseindex).coordinates(),betaname);
875 0 : imbeta.set(0.0);
876 :
877 0 : LatticeExpr<Float> betacalc( ((*smoothed[2])*mask1)/((*smoothed[0])+(mask0))-0.5*(imalpha)*(imalpha-1.0) );
878 0 : imbeta.copyData(betacalc);
879 :
880 0 : imbeta.setImageInfo(ii);
881 : //imbeta.setUnits(Unit("Spectral Curvature"));
882 0 : imbeta.table().unmarkForDelete();
883 0 : os << "Written Spectral Curvature Image : " << betaname << LogIO::POST;
884 0 : }
885 :
886 : /* Print out debugging info for center pixel */
887 : /*
888 : IPosition cgip(4,512,512,0,0);
889 : IPosition cgip2(4,490,542,0,0);
890 : for(Int i=0;i<ntaylor_p;i++)
891 : {
892 : cout << "Extended : " << endl;
893 : cout << "Original residual : " << i << " : " << residual( getModelIndex(model,i) ).getAt(cgip) << endl;
894 : cout << "Smoothed residual : " << i << " : " << (*smoothresiduals[i]).getAt(cgip) << endl;
895 : cout << "Coeff residual : " << i << " : " << (*coeffresiduals[i]).getAt(cgip) << endl;
896 : cout << "Point : " << endl;
897 : cout << "Original residual : " << i << " : " << residual( getModelIndex(model,i) ).getAt(cgip2) << endl;
898 : cout << "Smoothed residual : " << i << " : " << (*smoothresiduals[i]).getAt(cgip2) << endl;
899 : cout << "Coeff residual : " << i << " : " << (*coeffresiduals[i]).getAt(cgip2) << endl;
900 : }
901 : */
902 :
903 :
904 : /* Clean up temp arrays */
905 0 : for(uInt i=0;i<smoothed.nelements();i++) if(smoothed[i]) delete smoothed[i];
906 0 : for(uInt i=0;i<coeffresiduals.nelements();i++)
907 : {
908 0 : if(coeffresiduals[i]) delete coeffresiduals[i];
909 : }
910 :
911 0 : }// field loop
912 0 : return 0;
913 0 : }
914 :
915 : /***********************************************************************/
916 : /*************************************
917 : * Make Residuals and compute the current peak
918 : *************************************/
919 0 : Bool WBCleanImageSkyModel::solveResiduals(SkyEquation& se, Bool modelToMS)
920 : {
921 0 : blankOverlappingModels();
922 0 : makeNewtonRaphsonStep(se,false,modelToMS);
923 0 : restoreOverlappingModels();
924 :
925 0 : return true;
926 : }
927 : /***********************************************************************/
928 :
929 : /*************************************
930 : * Make Residuals
931 : *************************************/
932 : // Currently - all normalized by ggS from taylor0.
933 : // SAME WITH MAKE_SPECTRAL_PSFS
934 :
935 0 : Bool WBCleanImageSkyModel::makeNewtonRaphsonStep(SkyEquation& se, Bool incremental, Bool modelToMS)
936 : {
937 0 : LogIO os(LogOrigin("WBCleanImageSkyModel", "makeNewtonRaphsonStep"));
938 0 : se.gradientsChiSquared(incremental, modelToMS);
939 :
940 0 : Int index=0,baseindex=0;
941 0 : LatticeExpr<Float> le;
942 0 : for(Int thismodel=0;thismodel<nfields_p;thismodel++)
943 : {
944 0 : baseindex = getModelIndex(thismodel,0);
945 0 : for(Int taylor=0;taylor<ntaylor_p;taylor++)
946 : {
947 : /* Normalize by the Taylor 0 weight image */
948 0 : index = getModelIndex(thismodel,taylor);
949 : //cout << "Normalizing image " << index << " with " << baseindex << endl;
950 : //cout << "Shapes : " << ggS(index).shape() << gS(baseindex).shape() << endl;
951 :
952 : // UUU
953 : // LatticeExpr<Float> le(iif(ggS(baseindex)>(0.0), -gS(index)/ggS(baseindex), 0.0));
954 :
955 : ///cout << "WBC : makeNRs : isnormalized : " << isImageNormalized() << endl;
956 :
957 0 : if (isImageNormalized()) le = LatticeExpr<Float>(gS(index));
958 0 : else le = LatticeExpr<Float>(iif(ggS(baseindex)>(0.0),
959 0 : -gS(index)/ggS(baseindex), 0.0));
960 0 : residual(index).copyData(le);
961 :
962 : //LatticeExprNode LEN = max( residual(index) );
963 : //LatticeExprNode len2 = max( ggS(baseindex) );
964 : //cout << "Max Residual : " << LEN.getFloat() << " Max ggS : " << len2.getFloat() << endl;
965 :
966 :
967 : //storeAsImg(String("Weight.")+String::toString(thismodel)+String(".")+String::toString(taylor),ggS(index));
968 : //storeAsImg(String("TstResidual.")+String::toString(thismodel)+String(".")+String::toString(taylor),residual(index));
969 : }
970 : }
971 0 : modified_p=false;
972 0 : return true;
973 0 : }
974 :
975 : /*************************************
976 : * Make Approx PSFs.
977 : *************************************/
978 : // The normalization ignores that done in makeSimplePSFs in the Sky Eqn
979 : // and recalculates it from gS and ggS.
980 0 : Int WBCleanImageSkyModel::makeSpectralPSFs(SkyEquation& se, Bool writeToDisk)
981 : {
982 0 : LogIO os(LogOrigin("WBCleanImageSkyModel", "makeSpectralPSFs"));
983 0 : if(!donePSF_p)
984 : {
985 : //make sure the psf images are made
986 0 : for (Int thismodel=0;thismodel<nmodels_p;thismodel++)
987 : {
988 0 : PSF(thismodel);
989 : }
990 : }
991 :
992 : // All 2N-1 psfs will get made here....
993 0 : se.makeApproxPSF(psf_p);
994 :
995 : // Normalize
996 0 : Float normfactor=1.0;
997 0 : Int index=0,baseindex=0;
998 0 : LatticeExpr<Float> le;
999 0 : for (Int thismodel=0;thismodel<nfields_p;thismodel++)
1000 : {
1001 0 : normfactor=1.0;
1002 0 : baseindex = getModelIndex(thismodel,0);
1003 0 : for(Int taylor=0;taylor<2*ntaylor_p-1;taylor++)
1004 : {
1005 : /* Normalize by the Taylor 0 weight image */
1006 0 : index = getModelIndex(thismodel,taylor);
1007 : //cout << "Normalizing PSF " << index << " with " << baseindex << endl;
1008 : //cout << "Shapes : " << ggS(index).shape() << gS(baseindex).shape() << endl;
1009 : //le = LatticeExpr<Float>(iif(ggS(baseindex)>(0.0), gS(index)/ggS(baseindex), 0.0));
1010 0 : if (isImageNormalized()) le = LatticeExpr<Float>(gS(index));
1011 0 : else le = LatticeExpr<Float>(iif(ggS(baseindex)>(0.0),
1012 0 : gS(index)/ggS(baseindex), 0.0));
1013 0 : PSF(index).copyData(le);
1014 0 : if(taylor==0)
1015 : {
1016 0 : LatticeExprNode maxPSF=max(PSF(index));
1017 0 : normfactor = maxPSF.getFloat();
1018 0 : if(adbg) os << "Normalize PSFs for field " << thismodel << " by " << normfactor << LogIO::POST;
1019 0 : }
1020 0 : LatticeExpr<Float> lenorm(PSF(index)/normfactor);
1021 0 : PSF(index).copyData(lenorm);
1022 0 : LatticeExprNode maxPSF2=max(PSF(index));
1023 0 : Float maxpsf=maxPSF2.getFloat();
1024 0 : if(adbg) os << "Psf for Model " << thismodel << " and Taylor " << taylor << " has peak " << maxpsf << LogIO::POST;
1025 :
1026 0 : if(writeToDisk)
1027 : {
1028 : //need unique name as multiple processes can be running in the
1029 : //same workingdir
1030 0 : String tmpName=image_p[thismodel]->name();
1031 0 : tmpName.del(String(".model.tt0"), 0);
1032 0 : storeAsImg(tmpName+String(".TempPsf.")+String::toString(taylor),PSF(index));
1033 0 : }
1034 0 : }
1035 :
1036 : // index = getModelIndex(thismodel,0);
1037 : //beam(thismodel)=0.0;
1038 0 : if(!StokesImageUtil::FitGaussianPSF(PSF(baseindex),beam(thismodel)))
1039 0 : os << "Beam fit failed: using default" << LogIO::POST;
1040 : }
1041 0 : return 0;
1042 0 : }
1043 :
1044 :
1045 :
1046 : /***********************************************************************/
1047 :
1048 : /*************************************
1049 : * Store a Templattice as an image
1050 : *************************************/
1051 :
1052 0 : Int WBCleanImageSkyModel::storeAsImg(String fileName, ImageInterface<Float>& theImg)
1053 : {
1054 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), fileName);
1055 0 : LatticeExpr<Float> le(theImg);
1056 0 : tmp.copyData(le);
1057 0 : return 0;
1058 0 : }
1059 : /*
1060 : Int WBCleanImageSkyModel::storeTLAsImg(String fileName, TempLattice<Float> &TL, ImageInterface<Float>& theImg)
1061 : {
1062 : PagedImage<Float> tmp(TL.shape(), theImg.coordinates(), fileName);
1063 : LatticeExpr<Float> le(TL);
1064 : tmp.copyData(le);
1065 : return 0;
1066 : }
1067 :
1068 : Int WBCleanImageSkyModel::storeTLAsImg(String fileName, TempLattice<Complex> &TL, ImageInterface<Float>& theImg)
1069 : {
1070 : PagedImage<Complex> tmp(TL.shape(), theImg.coordinates(), fileName);
1071 : LatticeExpr<Complex> le(TL);
1072 : tmp.copyData(le);
1073 : return 0;
1074 : }
1075 : */
1076 :
1077 : /**************************
1078 : Resize Work Arrays to calculate extra PSFs.
1079 : (Resizing back to a shorter list must free memory correctly).
1080 : *************************/
1081 0 : Bool WBCleanImageSkyModel::resizeWorkArrays(Int length)
1082 : {
1083 0 : Int originallength = gS_p.nelements();
1084 :
1085 0 : if(length < originallength) // Clean up extra arrays
1086 : {
1087 0 : for(Int i = length; i < originallength; ++i)
1088 : {
1089 0 : if(psf_p[i]){delete psf_p[i]; psf_p[i]=0;}
1090 0 : if(image_p[i]){delete image_p[i]; image_p[i]=0;}
1091 0 : if(cimage_p[i]){delete cimage_p[i]; cimage_p[i]=0;}
1092 0 : if(gS_p[i]){delete gS_p[i]; gS_p[i]=0;}
1093 0 : if(ggS_p[i]){delete ggS_p[i]; ggS_p[i]=0;}
1094 0 : if(work_p[i]){delete work_p[i]; work_p[i]=0;}
1095 0 : if(fluxScale_p[i]){delete fluxScale_p[i]; fluxScale_p[i]=0;}
1096 : }
1097 :
1098 : }
1099 :
1100 0 : psf_p.resize(length,true);
1101 0 : image_p.resize(length,true);
1102 0 : cimage_p.resize(length,true);
1103 0 : gS_p.resize(length,true);
1104 0 : ggS_p.resize(length,true);
1105 0 : work_p.resize(length,true);
1106 0 : fluxScale_p.resize(length,true);
1107 :
1108 0 : if(length > originallength) // Add extra arrays // TODO : This part can go - I think !!!
1109 : {
1110 0 : for(Int i = originallength; i < length; ++i)
1111 : {
1112 0 : psf_p[i]=0;gS_p[i]=0;ggS_p[i]=0;cimage_p[i]=0;work_p[i]=0;fluxScale_p[i]=0;
1113 :
1114 0 : Int ind = getFieldIndex(i);
1115 : TempImage<Float>* imptr =
1116 0 : new TempImage<Float> (IPosition(image_p[ind]->ndim(),
1117 0 : image_p[ind]->shape()(0),
1118 0 : image_p[ind]->shape()(1),
1119 0 : image_p[ind]->shape()(2),
1120 0 : image_p[ind]->shape()(3)),
1121 0 : image_p[ind]->coordinates(), memoryMB_p);
1122 0 : AlwaysAssert(imptr, AipsError);
1123 0 : image_p[i] = imptr;
1124 : }
1125 : }
1126 :
1127 : ///cout << "Memory held by ImageSkyModel : " << 6*length << " float + " << 1*length << " complex images " << endl;
1128 :
1129 0 : return true;
1130 : }
1131 :
1132 : /************************************************************************************
1133 : Check some input parameters and print warnings for the user
1134 : fbw = (fmax-fmin)/fref.
1135 : if(fbw < 0.1) and nterms>2
1136 : => lever-arm may be insufficient for more than alpha.
1137 : => polynomial fit will work but alpha interpretation may not be ok.
1138 : if(ref < fmin or ref > fmax)
1139 : => polynomial fit will work, but alpha interpretation will not be right.
1140 : if(nchan==1) or fbw = 0, then ask to use only nterms=1, or choose more than one chan
1141 :
1142 : ***********************************************************************************/
1143 0 : Bool WBCleanImageSkyModel::checkParameters()
1144 : {
1145 : /* Check ntaylor_p, nrefFrequency_p with min and max freq from the image-coords */
1146 :
1147 0 : for(uInt i=0; i<image(0).coordinates().nCoordinates(); i++)
1148 : {
1149 0 : if( image(0).coordinates().type(i) == Coordinate::SPECTRAL )
1150 : {
1151 0 : SpectralCoordinate speccoord(image(0).coordinates().spectralCoordinate(i));
1152 0 : Double startfreq=0.0,startpixel=-0.5;
1153 0 : Double endfreq=0.0,endpixel=+0.5;
1154 0 : speccoord.toWorld(startfreq,startpixel);
1155 0 : speccoord.toWorld(endfreq,endpixel);
1156 0 : Float fbw = (endfreq - startfreq)/refFrequency_p;
1157 : //os << "Freq range of the mfs channel : " << startfreq << " -> " << endfreq << endl;
1158 : //cout << "Fractional bandwidth : " << fbw << endl;
1159 :
1160 0 : os << "Fractional Bandwidth : " << fbw*100 << " %." << LogIO::POST;
1161 :
1162 : /*
1163 : if(fbw < 0.1 && ntaylor_p == 2 )
1164 : os << "Fractional Bandwidth is " << fbw*100 << " %. Please check that the flux variation across the chosen frequency range (" << startfreq << " Hz to " << endfreq << " Hz) is at least twice the single-channel noise-level. If not, please use nterms=1." << LogIO::WARN << LogIO::POST;
1165 :
1166 : if(fbw < 0.1 && ntaylor_p > 2)
1167 : os << "Fractional Bandwidth is " << fbw*100 << " %. Please check that (a) the flux variation across the chosen frequency range (" << startfreq << " Hz to " << endfreq << " Hz) is at least twice the single-channel noise-level, and (b) a " << ntaylor_p << "-term Taylor-polynomial fit across this frequency range is appropriate. " << LogIO::WARN << LogIO::POST;
1168 : */
1169 0 : if(refFrequency_p < startfreq || refFrequency_p > endfreq)
1170 0 : os << "A Reference frequency of " << refFrequency_p << "Hz is outside the frequency range of the selected data (" << startfreq << " Hz to " << endfreq << " Hz). A power-law interpretation of the resulting Taylor-coefficients may not be accurate." << LogIO::POST;
1171 :
1172 0 : }
1173 : }
1174 :
1175 :
1176 0 : return true;
1177 : }
1178 :
1179 : // Copied from MFCleanImageSkyModel
1180 0 : void WBCleanImageSkyModel::blankOverlappingModels(){
1181 0 : if(nfields_p == 1) return;
1182 :
1183 0 : for(Int taylor=0;taylor<ntaylor_p;taylor++)
1184 : {
1185 0 : for (Int field=0;field<nfields_p-1; ++field)
1186 : {
1187 0 : Int model=getModelIndex(field,taylor);
1188 0 : CoordinateSystem cs0=image(model).coordinates();
1189 0 : IPosition iblc0(image(model).shape().nelements(),0);
1190 :
1191 0 : IPosition itrc0(image(model).shape());
1192 0 : itrc0=itrc0-Int(1);
1193 0 : LCBox lbox0(iblc0, itrc0, image(model).shape());
1194 :
1195 0 : ImageRegion imagreg0(WCBox(lbox0, cs0));
1196 0 : for (Int nextfield=field+1; nextfield < nfields_p; ++nextfield)
1197 : {
1198 0 : Int nextmodel = getModelIndex(nextfield, taylor);
1199 0 : CoordinateSystem cs=image(nextmodel).coordinates();
1200 0 : IPosition iblc(image(nextmodel).shape().nelements(),0);
1201 :
1202 0 : IPosition itrc(image(nextmodel).shape());
1203 0 : itrc=itrc-Int(1);
1204 :
1205 0 : LCBox lbox(iblc, itrc, image(nextmodel).shape());
1206 :
1207 0 : ImageRegion imagreg(WCBox(lbox, cs));
1208 : try{
1209 0 : LatticeRegion latReg=imagreg.toLatticeRegion(image(model).coordinates(), image(model).shape());
1210 0 : ArrayLattice<Bool> pixmask(latReg.get());
1211 0 : SubImage<Float> partToMask(image(model), imagreg, true);
1212 0 : LatticeExpr<Float> myexpr(iif(pixmask, 0.0, partToMask) );
1213 0 : partToMask.copyData(myexpr);
1214 :
1215 0 : }
1216 0 : catch(...){
1217 : //no overlap you think ?
1218 : //cout << "Did i fail " << endl;
1219 0 : continue;
1220 0 : }
1221 :
1222 :
1223 0 : }
1224 :
1225 :
1226 :
1227 0 : }// for field
1228 : }// for taylor
1229 : }
1230 :
1231 : // Copied from MFCleanImageSkyModel
1232 0 : void WBCleanImageSkyModel::restoreOverlappingModels(){
1233 0 : LogIO os(LogOrigin("WBImageSkyModel","restoreOverlappingModels"));
1234 0 : if(nfields_p == 1) return;
1235 :
1236 0 : for(Int taylor=0;taylor<ntaylor_p;taylor++)
1237 : {
1238 :
1239 0 : for (Int field=0;field<nfields_p-1; ++field)
1240 : {
1241 0 : Int model = getModelIndex(field,taylor);
1242 0 : CoordinateSystem cs0=image(model).coordinates();
1243 0 : IPosition iblc0(image(model).shape().nelements(),0);
1244 0 : IPosition itrc0(image(model).shape());
1245 0 : itrc0=itrc0-Int(1);
1246 0 : LCBox lbox0(iblc0, itrc0, image(model).shape());
1247 0 : ImageRegion imagreg0(WCBox(lbox0, cs0));
1248 :
1249 0 : for (Int nextfield=field+1; nextfield < nfields_p; ++nextfield)
1250 : {
1251 0 : Int nextmodel = getModelIndex(nextfield,taylor);
1252 0 : CoordinateSystem cs=image(nextmodel).coordinates();
1253 0 : IPosition iblc(image(nextmodel).shape().nelements(),0);
1254 :
1255 0 : IPosition itrc(image(nextmodel).shape());
1256 0 : itrc=itrc-Int(1);
1257 :
1258 0 : LCBox lbox(iblc, itrc, image(nextmodel).shape());
1259 :
1260 0 : ImageRegion imagreg(WCBox(lbox, cs));
1261 : try{
1262 0 : LatticeRegion latReg0=imagreg0.toLatticeRegion(image(nextmodel).coordinates(), image(nextmodel).shape());
1263 0 : LatticeRegion latReg=imagreg.toLatticeRegion(image(model).coordinates(), image(model).shape());
1264 0 : ArrayLattice<Bool> pixmask(latReg.get());
1265 :
1266 0 : SubImage<Float> partToMerge(image(nextmodel), imagreg0, true);
1267 0 : SubImage<Float> partToUnmask(image(model), imagreg, true);
1268 :
1269 0 : LatticeExpr<Float> myexpr0(iif(pixmask,partToMerge,partToUnmask));
1270 0 : partToUnmask.copyData(myexpr0);
1271 :
1272 0 : }
1273 0 : catch(AipsError x){
1274 : /*
1275 : os << LogIO::WARN
1276 : << "no overlap or failure of copying the clean components"
1277 : << x.getMesg()
1278 : << LogIO::POST;
1279 : */
1280 0 : continue;
1281 0 : }
1282 0 : }// for field - inner
1283 0 : }// for field - outer
1284 : }// for taylor
1285 0 : }
1286 :
1287 0 : Bool WBCleanImageSkyModel::mergeDataError(ImageInterface<Float> &data, ImageInterface<Float> &error, const String &outImg)
1288 : ///Bool WBCleanImageSkyModel::mergeDataError(const String &dataImg, const String &errorImg, const String &outImg)
1289 : {
1290 0 : LogIO os(LogOrigin("WBImageSkyModel",__FUNCTION__));
1291 :
1292 : // open the data and the error image
1293 : // ImageInterface<Float> *data = new PagedImage<Float>(dataImg, TableLock::AutoNoReadLocking);
1294 : // ImageInterface<Float> *error = new PagedImage<Float>(errorImg, TableLock::AutoNoReadLocking);
1295 :
1296 : // create the tiled shape for the output image
1297 0 : IPosition newShape=IPosition(data.shape());
1298 0 : newShape.append(IPosition(1, 2));
1299 0 : TiledShape tiledShape(newShape);
1300 :
1301 : // create the coordinate system for the output image
1302 0 : CoordinateSystem newCSys = data.coordinates();
1303 0 : Vector<Int> quality(2);
1304 0 : quality(0) = Quality::DATA;
1305 0 : quality(1) = Quality::ERROR;
1306 0 : QualityCoordinate qualAxis(quality);
1307 0 : newCSys.addCoordinate(qualAxis);
1308 :
1309 0 : Array<Float> outData=Array<Float>(newShape, 0.0f);
1310 0 : Array<Bool> outMask;
1311 :
1312 : // get all the data values
1313 0 : Array<Float> inData;
1314 0 : Slicer inSlicer(IPosition((data.shape()).size(), 0), IPosition(data.shape()));
1315 0 : data.doGetSlice(inData, inSlicer);
1316 :
1317 : // define in the output array
1318 : // the slicers for data and error
1319 : Int qualCooPos, qualIndex;
1320 0 : qualCooPos = newCSys.findCoordinate(Coordinate::QUALITY);
1321 0 : (newCSys.qualityCoordinate(qualCooPos)).toPixel(qualIndex, Quality::DATA);
1322 0 : IPosition outStart(newShape.size(), 0);
1323 0 : outStart(newShape.size()-1)=qualIndex;
1324 0 : IPosition outLength(newShape);
1325 0 : outLength(newShape.size()-1)=1;
1326 0 : Slicer outDataSlice(outStart, outLength);
1327 0 : (newCSys.qualityCoordinate(qualCooPos)).toPixel(qualIndex, Quality::ERROR);
1328 0 : outStart(newShape.size()-1)=qualIndex;
1329 0 : Slicer outErrorSlice(outStart, outLength);
1330 :
1331 : // add the data values to the output array
1332 0 : outData(outDataSlice) = inData.addDegenerate(1);
1333 :
1334 : // get all the error values
1335 0 : error.doGetSlice(inData, inSlicer);
1336 :
1337 :
1338 : // add the error values to the output array
1339 0 : outData(outErrorSlice) = inData.addDegenerate(1);
1340 :
1341 : // check whether a mask is necessary
1342 0 : if (data.hasPixelMask() || error.hasPixelMask()){
1343 0 : Array<Bool> inMask;
1344 :
1345 0 : outMask=Array<Bool>(newShape, true);
1346 :
1347 : // make the mask for the data values
1348 0 : if (data.hasPixelMask()){
1349 0 : inMask = (data.pixelMask()).get();
1350 : }
1351 : else{
1352 0 : inMask = Array<Bool>(data.shape(), true);
1353 : }
1354 :
1355 : // add the data mask to the output
1356 0 : outMask(outDataSlice) = inMask.addDegenerate(1);
1357 :
1358 : // make the mask for the error values
1359 0 : if (error.hasPixelMask()){
1360 0 : inMask = (error.pixelMask()).get();
1361 : }
1362 : else{
1363 0 : inMask = Array<Bool>(error.shape(), true);
1364 : }
1365 :
1366 : // add the data mask to the output
1367 0 : outMask(outErrorSlice) = inMask.addDegenerate(1);
1368 0 : }
1369 :
1370 :
1371 : // write out the combined image
1372 0 : ImageUtilities::writeImage(tiledShape, newCSys, outImg, outData, os, outMask);
1373 :
1374 : // delete data;
1375 : // delete error;
1376 :
1377 0 : return true;
1378 0 : }
1379 :
1380 : } //# NAMESPACE CASA - END
1381 :
|
