Line data Source code
1 : //# MFMSCleanImageSkyModel.cc: Implementation of MFMSCleanImageSkyModel 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: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <synthesis/MeasurementComponents/MFMSCleanImageSkyModel.h>
30 : #include <casacore/images/Images/PagedImage.h>
31 : #include <casacore/casa/OS/File.h>
32 : #include <casacore/lattices/LEL/LatticeExpr.h>
33 : #include <casacore/lattices/LEL/LatticeExprNode.h>
34 : #include <casacore/images/Images/SubImage.h>
35 : #include <casacore/casa/Arrays/IPosition.h>
36 : #include <casacore/lattices/LRegions/LCBox.h>
37 : #include <synthesis/MeasurementEquations/ImageMSCleaner.h>
38 : #include <synthesis/MeasurementEquations/SkyEquation.h>
39 : #include <casacore/casa/Exceptions/Error.h>
40 : #include <casacore/casa/BasicSL/String.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 : #include <casacore/casa/Quanta/Quantum.h>
43 :
44 : #include <sstream>
45 : #include <casacore/casa/Logging/LogMessage.h>
46 : #include <casacore/casa/Logging/LogSink.h>
47 : #include <casacore/casa/Logging/LogIO.h>
48 :
49 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
50 :
51 : using namespace casacore;
52 : namespace casa {
53 :
54 :
55 : // Some constructors
56 :
57 0 : MFMSCleanImageSkyModel::MFMSCleanImageSkyModel()
58 : : method_p(NSCALES), nscales_p(4), progress_p(0),
59 0 : stopLargeNegatives_p(2), stopPointMode_p(-1)
60 : {
61 :
62 0 : donePSF_p=false;
63 0 : modified_p=true;
64 0 : getScales();
65 :
66 0 : };
67 :
68 0 : MFMSCleanImageSkyModel::MFMSCleanImageSkyModel(const Int nscales,
69 : const Int sln,
70 : const Int spm,
71 0 : const Float inbias)
72 : : method_p(NSCALES), nscales_p(nscales), progress_p(0),
73 0 : stopLargeNegatives_p(sln), stopPointMode_p(spm),smallScaleBias_p(inbias)
74 : {
75 :
76 0 : donePSF_p=false;
77 0 : modified_p=true;
78 0 : getScales();
79 :
80 0 : };
81 :
82 0 : MFMSCleanImageSkyModel::MFMSCleanImageSkyModel(const Vector<Float>& userScaleSizes,
83 : const Int sln,
84 : const Int spm,
85 0 : const Float inbias)
86 : : method_p(USERVECTOR), userScaleSizes_p(userScaleSizes), progress_p(0),
87 0 : stopLargeNegatives_p(sln), stopPointMode_p(spm), smallScaleBias_p(inbias)
88 : {
89 :
90 0 : donePSF_p=false;
91 0 : modified_p=true;
92 0 : getScales();
93 :
94 0 : };
95 :
96 0 : MFMSCleanImageSkyModel::~MFMSCleanImageSkyModel()
97 : {
98 0 : if (progress_p) {
99 0 : delete progress_p;
100 : }
101 :
102 0 : if(componentList_p) delete componentList_p; componentList_p=0;
103 0 : for (Int thismodel=0;thismodel<nmodels_p;thismodel++) {
104 0 : if(residualImage_p[thismodel]) delete residualImage_p[thismodel];
105 0 : residualImage_p[thismodel]=0;
106 0 : if(cimage_p[thismodel]) delete cimage_p[thismodel];
107 0 : cimage_p[thismodel]=0;
108 0 : for(Int numXFR=0;numXFR<maxNumXFR_p;numXFR++) {
109 0 : if(cxfr_p[thismodel*maxNumXFR_p+numXFR])
110 0 : delete cxfr_p[thismodel*maxNumXFR_p+numXFR];
111 0 : cxfr_p[thismodel*maxNumXFR_p+numXFR]=0;
112 : }
113 0 : if(gS_p[thismodel]) delete gS_p[thismodel]; gS_p[thismodel]=0;
114 0 : if(psf_p[thismodel]) delete psf_p[thismodel]; psf_p[thismodel]=0;
115 0 : if(ggS_p[thismodel]) delete ggS_p[thismodel]; ggS_p[thismodel]=0;
116 0 : if(fluxScale_p[thismodel]) delete fluxScale_p[thismodel];
117 0 : fluxScale_p[thismodel]=0;
118 0 : if(work_p[thismodel]) delete work_p[thismodel]; work_p[thismodel]=0;
119 0 : if(deltaimage_p[thismodel]) delete deltaimage_p[thismodel];
120 0 : deltaimage_p[thismodel]=0;
121 0 : if(weight_p[thismodel]) delete weight_p[thismodel]; weight_p[thismodel]=0;
122 0 : if(beam_p[thismodel]) delete beam_p[thismodel]; beam_p[thismodel]=0;
123 : }
124 :
125 0 : };
126 :
127 :
128 : // Clean solver
129 0 : Bool MFMSCleanImageSkyModel::solve(SkyEquation& se) {
130 :
131 0 : LogIO os(LogOrigin("MFMSCleanImageSkyModel","solve"));
132 : /*backed out for now
133 : if(modified_p){
134 : makeNewtonRaphsonStep(se, false);
135 : }
136 :
137 : if(numberIterations() < 1){
138 : return true;
139 : }
140 : */
141 : //Make the PSFs, one per field
142 0 : if(!donePSF_p){
143 : os << LogIO::NORMAL // Loglevel PROGRESS
144 0 : << "Making approximate PSFs" << LogIO::POST;
145 0 : makeApproxPSFs(se);
146 : }
147 :
148 0 : Int converged=true;
149 0 : Int converging=0;
150 :
151 : // Validate PSFs for each field
152 0 : Vector<Float> psfmax(numberOfModels()); psfmax=0.0;
153 0 : Vector<Float> psfmin(numberOfModels()); psfmin=1.0;
154 0 : Vector<Float> resmax(numberOfModels());
155 0 : Vector<Float> resmin(numberOfModels());
156 :
157 0 : Float maxSidelobe=0.0;
158 : Int model;
159 0 : for (model=0;model<numberOfModels();model++) {
160 0 : if(isSolveable(model)) {
161 0 : IPosition blc(PSF(model).shape().nelements(), 0);
162 0 : IPosition trc(PSF(model).shape()-1);
163 0 : blc(2) = 0; trc(2) = 0;
164 0 : blc(3) = 0; trc(3) = 0;
165 :
166 0 : SubImage<Float> subPSF;
167 0 : Int k =0;
168 0 : Int numchan= PSF(model).shape()(3);
169 : //PSF of first non zero plane
170 0 : while(psfmax(model) < 0.1 && k< numchan){
171 0 : blc(3)= k;
172 0 : trc(3)=k;
173 0 : LCBox onePlane(blc, trc, PSF(model).shape());
174 :
175 0 : subPSF=SubImage<Float> ( PSF(model), onePlane, true);
176 : {
177 0 : LatticeExprNode node = max(subPSF);
178 0 : psfmax(model) = node.getFloat();
179 : }
180 0 : ++k;
181 : }
182 : {
183 0 : LatticeExprNode node = min(subPSF);
184 0 : psfmin(model) = node.getFloat();
185 : }
186 : os << LogIO::NORMAL // Loglevel INFO
187 : << "Model " << model+1 << ": max, min PSF = "
188 0 : << psfmax(model) << ", " << psfmin(model) << endl;
189 0 : if(abs(psfmin(model))>maxSidelobe) maxSidelobe=abs(psfmin(model));
190 : }
191 : }
192 0 : os << LogIO::POST;
193 :
194 0 : Float absmax=threshold();
195 0 : Block< Matrix<Int> > iterations(numberOfModels());
196 0 : Int maxIterations=0;
197 0 : Int oldMaxIterations=-1;
198 :
199 : // Loop over major cycles
200 : /*if (displayProgress_p) {
201 : if(progress_p) delete progress_p;
202 : progress_p=0;
203 : progress_p = new LatticeCleanProgress( pgplotter_p );
204 : }*/
205 :
206 0 : Block<CountedPtr<ImageMSCleaner > > cleaner(numberOfModels());
207 0 : cleaner=0;
208 :
209 0 : Int cycle=0;
210 0 : Bool stop=false;
211 0 : Bool lastCycleWriteModel=false;
212 :
213 0 : while(absmax>=threshold()&&maxIterations<numberIterations()&&!stop) {
214 :
215 : os << LogIO::NORMAL // Loglevel PROGRESS
216 0 : << "*** Starting major cycle " << cycle+1 << LogIO::POST;
217 0 : cycle++;
218 :
219 :
220 : // Make the residual images. We do an incremental update
221 : // for cycles after the first one. If we have only one
222 : // model then we use convolutions to speed the processing
223 0 : Bool incremental(cycle>1);
224 0 : if(modified_p){
225 0 : if (!incremental||(itsSubAlgorithm == "full")) {
226 :
227 : os << LogIO::NORMAL1 // Loglevel INFO
228 : << "Using visibility-subtraction for residual calculation"
229 0 : << LogIO::POST;
230 0 : makeNewtonRaphsonStep(se, false, (numberIterations()<1)?true:False);
231 :
232 : }
233 : else {
234 : os << LogIO::NORMAL1 // Loglevel INFO
235 : << "Using XFR-based shortcut for residual calculation"
236 0 : << LogIO::POST;
237 0 : makeNewtonRaphsonStep(se, true);
238 : }
239 :
240 : }
241 0 : if(numberIterations() < 1){
242 : // Why waste the time to set up
243 0 : return true;
244 : }
245 :
246 0 : absmax=maxField(resmax, resmin);
247 :
248 0 : for (model=0;model<numberOfModels();model++) {
249 : os << LogIO::NORMAL // Loglevel INFO
250 : << "Model " << model+1 << ": max, min residuals (over all pixels) = "
251 0 : << resmax(model) << ", " << resmin(model) << endl;
252 : }
253 0 : os << LogIO::POST;
254 :
255 : // Can we stop?
256 0 : if(absmax<threshold()) {
257 : os << LogIO::NORMAL // Loglevel PROGRESS?
258 : << "Reached stopping peak point source residual = "
259 0 : << absmax << LogIO::POST;
260 0 : stop=true;
261 0 : if(cycle >1)
262 0 : lastCycleWriteModel=true;
263 : }
264 : else {
265 : // Calculate the threshold for this cycle. Add a safety factor
266 0 : Float fudge = cycleFactor_p*maxSidelobe;
267 0 : if (fudge > 0.8) fudge = 0.8; // painfully slow!
268 0 : Quantity fThreshold(fudge*100, "%");
269 0 : Quantity aThreshold(threshold(), "Jy");
270 : os << LogIO::NORMAL // Loglevel INFO
271 0 : << "Maximum point source residual = " << absmax << LogIO::POST;
272 : os << LogIO::NORMAL // Loglevel INFO
273 : << "Cleaning scale flux down to maximum of " << fThreshold.getValue("%")
274 0 : << " % and " << aThreshold.getValue("Jy") << " Jy" << LogIO::POST;
275 :
276 0 : for (model=0;model<numberOfModels();model++) {
277 :
278 0 : Int nchan=image(model).shape()(3);
279 0 : Int npol=image(model).shape()(2);
280 :
281 0 : IPosition blcDirty(image(model).shape().nelements(), 0);
282 0 : IPosition trcDirty(image(model).shape()-1);
283 :
284 0 : if(cycle==1) {
285 0 : iterations[model].resize(nchan,npol);
286 0 : iterations[model]=0;
287 : }
288 :
289 : // Initialize delta image
290 0 : deltaImage(model).set(0.0);
291 :
292 : // Only process solveable models
293 0 : if(isSolveable(model)) {
294 :
295 : os << LogIO::NORMAL // Loglevel PROGRESS
296 0 : << "Processing model " << model+1 << LogIO::POST;
297 :
298 : // If mask exists, use it;
299 : // If not, use the fluxScale image to figure out
300 0 : Bool doMask = false;
301 0 : Bool mustDeleteMask = false;
302 0 : ImageInterface<Float> *maskPointer = 0;
303 :
304 0 : if (hasMask(model)) {
305 0 : doMask = true;
306 0 : maskPointer = &mask(model);
307 0 : } else if (doFluxScale(model)) {
308 0 : doMask = true;
309 0 : mustDeleteMask = true;
310 :
311 0 : maskPointer = new TempImage<Float> ( fluxScale(model).shape(),
312 0 : fluxScale(model).coordinates());
313 0 : maskPointer->copyData( (LatticeExpr<Float>)
314 0 : (iif( (fluxScale(model) > 0.0), 1.0, 0.0)) );
315 : }
316 :
317 : // Now clean each channel and each pol
318 0 : for (Int chan=0; chan<nchan; chan++) {
319 0 : if(psfmax(model)>0.0) {
320 : // We could keep a cleaner per channel but for the moment
321 : // we simply make a new one for each channel
322 0 : if(nchan>1) {
323 : os << LogIO::NORMAL // Loglevel PROGRESS
324 0 : <<"Processing channel "<<chan<<" of 0 to "<<nchan-1<<LogIO::POST;
325 0 : if(!cleaner[model].null()) cleaner[model]=0;
326 : }
327 :
328 0 : blcDirty(3) = chan;
329 0 : trcDirty(3) = chan;
330 0 : blcDirty(2) = 0; trcDirty(2) = 0;
331 0 : LCBox firstPlane(blcDirty, trcDirty, image(model).shape());
332 0 : SubImage<Float> subPSF( PSF(model), firstPlane);
333 :
334 0 : for (Int pol=0; pol<npol; pol++) {
335 0 : blcDirty(2) = pol; trcDirty(2) = pol;
336 : // The PSF should be the same for each polarization so we
337 : // can use the existing cleaner (unlike the spectral case)
338 0 : if(npol>1) {
339 : os << LogIO::NORMAL // Loglevel PROGRESS
340 : <<"Processing polarization "<<pol+1<<" of "<< npol
341 0 : <<LogIO::POST;
342 : }
343 0 : LCBox onePlane(blcDirty, trcDirty, image(model).shape());
344 :
345 0 : SubImage<Float> subImage( image(model), onePlane, true);
346 0 : SubImage<Float> subResid( residual(model), onePlane);
347 0 : SubImage<Float> subDeltaImage( deltaImage(model), onePlane, true);
348 0 : SubImage<Float> subMask;
349 0 : Bool skipThisPlane=false;
350 0 : if (doMask) {
351 0 : subMask = SubImage<Float> ( *maskPointer, onePlane, true);
352 0 : if(max(subMask).getFloat() <= 0.0){
353 0 : skipThisPlane=true;
354 : }
355 : }
356 0 : if(!skipThisPlane){
357 0 : if(!cleaner[model].null()) {
358 : os << LogIO::NORMAL2 // Loglevel PROGRESS
359 : << "Updating multiscale cleaner with new residual images"
360 0 : << LogIO::POST;
361 0 : cleaner[model]->update(subResid);
362 : }
363 : else {
364 : os << LogIO::NORMAL2 // Loglevel PROGRESS
365 : << "Creating multiscale cleaner with psf and residual images"
366 0 : << LogIO::POST;
367 0 : cleaner[model]=new ImageMSCleaner(subPSF, subResid);
368 : //setScales(*(cleaner[model]));
369 0 : cleaner[model]->setscales(userScaleSizes_p);
370 0 : cleaner[model]->setSmallScaleBias(smallScaleBias_p);
371 0 : if (doMask) {
372 0 : cleaner[model]->setMask(subMask);
373 : }
374 : }
375 0 : subDeltaImage.set(0.0);
376 :
377 0 : cleaner[model]->setcontrol(CleanEnums::MULTISCALE, numberIterations(), gain(),
378 : aThreshold, fThreshold);
379 :
380 0 : if (cycleSpeedup_p > 1) {
381 : os << LogIO::NORMAL // Loglevel INFO
382 0 : << "cycleSpeedup is " << cycleSpeedup_p << LogIO::POST;
383 0 : cleaner[model]->speedup(cycleSpeedup_p);
384 : }
385 0 : cleaner[model]->startingIteration( iterations[model](chan,pol) );
386 0 : if (cycle <= stopLargeNegatives_p) {
387 0 : cleaner[model]->stopAtLargeScaleNegative();
388 : }
389 0 : cleaner[model]->stopPointMode(stopPointMode_p);
390 0 : cleaner[model]->ignoreCenterBox(true);
391 0 : converging=cleaner[model]->clean(subDeltaImage, "fullmsclean", numberIterations(), gain(), aThreshold, fThreshold, displayProgress_p );
392 : //diverging
393 0 : if(converging==-3)
394 0 : stop=true;
395 : //reduce scales on main field only
396 0 : if(converging==-2 && model==0){
397 0 : if(userScaleSizes_p.nelements() > 1){
398 0 : userScaleSizes_p.resize(userScaleSizes_p.nelements()-1, true);
399 0 : cleaner[model]->setscales(userScaleSizes_p);
400 : }
401 : else
402 0 : stop=true;
403 : }
404 0 : iterations[model](chan,pol)=cleaner[model]->numberIterations();
405 0 : maxIterations=(iterations[model](chan,pol)>maxIterations) ?
406 0 : iterations[model](chan,pol) : maxIterations;
407 : os << LogIO::NORMAL // Loglevel INFO
408 0 : << "Clean used " << iterations[model](chan,pol)
409 : << " iterations"
410 0 : << LogIO::POST;
411 0 : modified_p=true;
412 :
413 0 : subImage.copyData( LatticeExpr<Float>( subImage + subDeltaImage));
414 :
415 0 : if (cleaner[model]->queryStopPointMode()) {
416 0 : stop = true;
417 : os << LogIO::NORMAL // Loglevel INFO
418 : << "MSClean terminating because we hit "
419 : << stopPointMode_p
420 0 : << " consecutive compact sources" << LogIO::POST;
421 : }
422 : //if (doMask) {
423 : // delete subMaskPointer;
424 : //}
425 : }
426 : }
427 : }
428 : }
429 0 : if (mustDeleteMask) {
430 0 : delete maskPointer;
431 : }
432 : } // if solveable
433 : } // for model
434 0 : if(maxIterations==oldMaxIterations){
435 : os << LogIO::NORMAL // Loglevel INFO
436 : << "MSClean terminating because components search stopped "
437 0 : << LogIO::POST;
438 0 : stop=true;
439 : }
440 : else{
441 0 : oldMaxIterations=maxIterations;
442 : }
443 : }
444 : }
445 :
446 0 : if(modified_p || lastCycleWriteModel) {
447 : os << LogIO::NORMAL // Loglevel PROGRESS?
448 0 : << "Finalizing residual images for all fields" << LogIO::POST;
449 0 : makeNewtonRaphsonStep(se, false, true);
450 0 : Float finalabsmax=maxField(resmax, resmin);
451 0 : setNumberIterations(numberIterations()-maxIterations);
452 0 : converged=(finalabsmax < 1.05 * threshold());
453 0 : setThreshold(finalabsmax);
454 : os << LogIO::NORMAL // Loglevel INFO
455 0 : << "Final maximum residual = " << finalabsmax << LogIO::POST;
456 :
457 0 : for (model=0;model<numberOfModels();model++) {
458 : os << LogIO::NORMAL // Loglevel INFO
459 : << "Model " << model+1 << ": max, min residuals (over all pixels) = "
460 0 : << resmax(model) << ", " << resmin(model) << endl;
461 : }
462 0 : os << LogIO::POST;
463 : }
464 : else {
465 0 : os << "Residual images for all fields are up-to-date" << LogIO::POST;
466 : }
467 :
468 : //if(cleaner) delete cleaner; cleaner=0;
469 :
470 0 : return(converged);
471 : };
472 :
473 :
474 : void
475 0 : MFMSCleanImageSkyModel::getScales()
476 : {
477 0 : LogIO os(LogOrigin("MFCleanImageSkyModel", "getScales"));
478 0 : if (method_p == USERVECTOR) {
479 0 : if (userScaleSizes_p.nelements() <= 0) {
480 : os << LogIO::SEVERE
481 : << "Need at least one scale for method uservector"
482 0 : << LogIO::POST;
483 : }
484 0 : os << "Creating scales from uservector method: " << LogIO::POST;
485 : }
486 : else {
487 0 : if (nscales_p <= 0) nscales_p = 1;
488 0 : Vector<Float> scaleSizes(nscales_p);
489 : os << "Creating " << nscales_p <<
490 0 : " scales from powerlaw nscales method" << LogIO::POST;
491 0 : scaleSizes(0) = 0.0;
492 0 : Float scaleInc = 2.0;
493 0 : for(uInt scale = 1; scale < nscales_p; ++scale)
494 0 : scaleSizes[scale] = scaleInc * pow(10.0, (Float(scale) - 2.0)/2.0);
495 :
496 : //store the scales as user setscales..in case we need to reduce scales
497 0 : userScaleSizes_p.resize();
498 0 : userScaleSizes_p=scaleSizes;
499 0 : method_p=USERVECTOR;
500 : }
501 :
502 0 : for(uInt scale = 0; scale < userScaleSizes_p.nelements(); ++scale)
503 : os << LogIO::NORMAL <<
504 0 : "scale " << scale+1 << " = " << userScaleSizes_p(scale)
505 : << " pixels"
506 0 : << LogIO::POST;
507 0 : };
508 :
509 : /*
510 : // Inlined here and not in the .h because the .h doesn't #include LatticeCleaner.
511 : inline void MFMSCleanImageSkyModel::setScales(LatticeCleaner<Float>& cleaner)
512 : {
513 : cleaner.setscales(userScaleSizes_p);
514 : }
515 :
516 : */
517 : } //#End casa namespace
|
