Line data Source code
1 : //# MFCleanImageSkyModel.cc: Implementation of MFCleanImageSkyModel class
2 : //# Copyrighogt (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$
27 :
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <casacore/casa/Arrays/ArrayLogical.h>
30 : #include <casacore/casa/Arrays/Matrix.h>
31 : #include <synthesis/MeasurementComponents/MFCleanImageSkyModel.h>
32 : #include <casacore/images/Images/PagedImage.h>
33 : #include <casacore/images/Images/SubImage.h>
34 : #include <casacore/images/Regions/ImageRegion.h>
35 : #include <casacore/images/Regions/WCBox.h>
36 : #include <casacore/casa/OS/File.h>
37 : #include <casacore/lattices/LEL/LatticeExpr.h>
38 : #include <casacore/lattices/LEL/LatticeExprNode.h>
39 : #include <casacore/lattices/Lattices/LatticeStepper.h>
40 : #include <casacore/lattices/Lattices/LatticeIterator.h>
41 : #include <casacore/lattices/LRegions/RegionType.h>
42 : #include <casacore/lattices/Lattices/SubLattice.h>
43 : #include <casacore/lattices/LRegions/LCBox.h>
44 :
45 : #include <synthesis/MeasurementEquations/SkyEquation.h>
46 : #include <casacore/casa/Exceptions/Error.h>
47 : #include <casacore/casa/BasicSL/String.h>
48 : #include <casacore/casa/Utilities/Assert.h>
49 : #include <casacore/casa/Utilities/CountedPtr.h>
50 :
51 : #include <sstream>
52 :
53 : #include <casacore/casa/Logging/LogMessage.h>
54 : #include <casacore/casa/Logging/LogSink.h>
55 : #include <casacore/casa/Logging/LogIO.h>
56 : #include <casacore/casa/BasicSL/Constants.h>
57 :
58 : #include <msvis/MSVis/StokesVector.h>
59 : #include <synthesis/MeasurementEquations/ConvolutionEquation.h>
60 : #include <synthesis/MeasurementEquations/ClarkCleanModel.h>
61 :
62 :
63 :
64 : using namespace casacore;
65 : namespace casa {
66 :
67 : #define NEED_UNDERSCORES
68 : #if defined(NEED_UNDERSCORES)
69 : #define hclean hclean_
70 : #endif
71 : extern "C" {
72 : void hclean(Float*, Float*, Float*, int*, Float*, int*, int*, int*,
73 : int*, int*, int*, int*, int*, int*, int*, Float*, Float*,
74 : Float*, void *, void *);
75 : };
76 :
77 :
78 : // This is defined in HogbomCleanImageSkyModel.cc
79 : void HogbomCleanImageSkyModelstopnow (Int *yes);
80 :
81 : // This is defined in HogbomCleanImageSkyModel.cc
82 : void HogbomCleanImageSkyModelmsgput(Int* npol, Int* pol, Int* iter, Int* px, Int* py, Float* fMaxVal);
83 :
84 :
85 0 : Int MFCleanImageSkyModel::add(ImageInterface<Float>& image,
86 : const Int maxNumXfr)
87 : {
88 0 : return CleanImageSkyModel::add(image, maxNumXfr);
89 : };
90 :
91 0 : Bool MFCleanImageSkyModel::addMask(Int image,
92 : ImageInterface<Float>& mask)
93 : {
94 0 : return CleanImageSkyModel::addMask(image, mask);
95 : };
96 :
97 0 : Bool MFCleanImageSkyModel::addResidual(Int image,
98 : ImageInterface<Float>& residual)
99 : {
100 0 : return ImageSkyModel::addResidual(image, residual);
101 : };
102 :
103 : // Clean solver
104 0 : Bool MFCleanImageSkyModel::solve(SkyEquation& se) {
105 :
106 :
107 0 : LogIO os(LogOrigin("MFCleanImageSkyModel","solve"));
108 0 : Bool converged=true;
109 : //Make the PSFs, one per field
110 : /*back out for now
111 : if(modified_p){
112 : blankOverlappingModels();
113 : makeNewtonRaphsonStep(se, false);
114 : }
115 : if(numberIterations() < 1){
116 : // Why waste the time to set up
117 : return true;
118 : }
119 : */
120 :
121 0 : if(!donePSF_p) {
122 0 : makeApproxPSFs(se);
123 : }
124 :
125 : // Validate PSFs for each field
126 0 : Vector<Float> psfmax(numberOfModels()); psfmax=0.0;
127 0 : Vector<Float> psfmaxouter(numberOfModels()); psfmaxouter=0.0;
128 0 : Vector<Float> psfmin(numberOfModels()); psfmin=1.0;
129 0 : Vector<Float> resmax(numberOfModels());
130 0 : Vector<Float> resmin(numberOfModels());
131 0 : Int psfpatch=51;
132 0 : Float maxSidelobe=0.0;
133 : Int model;
134 0 : os << LogIO::NORMAL1; // Loglevel INFO
135 0 : for (model=0;model<numberOfModels();model++) {
136 0 : if(isSolveable(model)) {
137 :
138 0 : IPosition blc(PSF(model).shape().nelements(), 0);
139 0 : IPosition trc(PSF(model).shape()-1);
140 0 : blc(2) = 0; trc(2) = 0;
141 0 : blc(3) = 0; trc(3) = 0;
142 :
143 0 : Float maxMaxPSF=max(PSF(model)).getFloat();
144 0 : SubLattice<Float> subPSF;
145 0 : Int k =0;
146 0 : Int numchan= PSF(model).shape()(3);
147 : //PSF of first non zero plane
148 0 : while(psfmax(model) < (0.8*maxMaxPSF) && k< numchan){
149 0 : blc(3)= k;
150 0 : trc(3)=k;
151 0 : LCBox onePlane(blc, trc, PSF(model).shape());
152 0 : subPSF=SubLattice<Float> ( PSF(model), onePlane, true);
153 : {
154 0 : LatticeExprNode node = max(subPSF);
155 0 : psfmax(model) = node.getFloat();
156 0 : }
157 0 : ++k;
158 0 : }
159 : {
160 0 : LatticeExprNode node = min(subPSF);
161 0 : psfmin(model) = node.getFloat();
162 0 : }
163 : // 4 pixels: pretty arbitrary, but only look for sidelobes
164 : // outside the inner (2n+1) * (2n+1) square
165 0 : Int ncent=4;
166 : {//locate peak size
167 0 : CoordinateSystem cs=PSF(model).coordinates();
168 0 : Vector<String> unitas=cs.worldAxisUnits();
169 0 : unitas(0)="arcsec";
170 0 : unitas(1)="arcsec";
171 0 : cs.setWorldAxisUnits(unitas);
172 : //Get the minimum increment in arcsec
173 0 : Double incr=abs(min(cs.increment()(0), cs.increment()(1)));
174 0 : if(incr > 0.0){
175 0 : GaussianBeam beamModel=beam(model)(0,0);
176 0 : ncent=max(ncent, Int(ceil(beamModel.getMajor("arcsec")/incr)));
177 0 : ncent=max(ncent, Int(ceil(beamModel.getMinor("arcsec")/incr)));
178 0 : }
179 0 : }
180 0 : psfpatch=3*ncent+1;
181 :
182 0 : psfmaxouter(model) = maxOuter(subPSF, ncent);
183 : os << "Model " << model << ": max, min, maxOuter PSF = "
184 0 : << psfmax(model) << ", " << psfmin(model) << ", " <<
185 0 : psfmaxouter(model) << endl;
186 0 : if(abs(psfmin(model))>maxSidelobe) maxSidelobe=abs(psfmin(model));
187 0 : if(psfmaxouter(model) > maxSidelobe) maxSidelobe= psfmaxouter(model );
188 : // os << " : PSFPatch = " << psfpatch << LogIO::POST;
189 0 : }
190 : }
191 0 : os << LogIO::POST;
192 :
193 0 : Float absmax=threshold();
194 0 : Float oldabsmax=C::flt_max;
195 0 : Float cycleThreshold=0.0;
196 0 : Block< Vector<Int> > iterations(numberOfModels());
197 0 : Int maxIterations=0;
198 0 : Int oldMaxIterations=0;
199 :
200 : // Loop over major cycles
201 0 : Int cycle=0;
202 0 : Bool stop=false;
203 0 : Bool diverging=false;
204 :
205 0 : if (displayProgress_p) {
206 0 : progress_p = new ClarkCleanProgress( pgplotter_p );
207 : } else {
208 0 : progress_p = 0;
209 : }
210 :
211 : //Making the mask matrices
212 0 : Int nmodels=numberOfModels();
213 0 : Block< Vector<Int> > xbeg(nmodels), ybeg(nmodels), xend(nmodels),
214 0 : yend(nmodels);
215 0 : Vector<Bool> isCubeMask(nmodels);
216 0 : isCubeMask.set(false);
217 : Int nx, ny;
218 0 : Int npol=image(0).shape()(2);
219 0 : Int nchan=image(0).shape()(3);
220 : //Int polloop=1;
221 0 : Vector<String> stokesID(npol, "");
222 0 : if (!doPolJoint_p) {
223 : //polloop=npol;
224 0 : CoordinateSystem cs=image(0).coordinates();
225 0 : Int stokesindex=cs.findCoordinate(Coordinate::STOKES);
226 0 : StokesCoordinate stcoord=cs.stokesCoordinate(stokesindex);
227 0 : for (Int jj=0; jj < npol ; ++jj){
228 0 : stokesID(jj)=Stokes::name(Stokes::type(stcoord.stokes()(jj)));
229 : }
230 0 : }
231 :
232 : //merge the overlapping part of the masks if any
233 0 : mergeOverlappingMasks();
234 :
235 :
236 0 : PtrBlock<Matrix<Float> *> lmaskCube(nmodels*nchan);
237 :
238 0 : for (model=0;model<numberOfModels();model++) {
239 :
240 0 : nx=image(model).shape()(0);
241 0 : ny=image(model).shape()(1);
242 0 : npol=image(model).shape()(2);
243 0 : nchan=image(model).shape()(3);
244 0 : xbeg[model].resize(nchan);
245 0 : xend[model].resize(nchan);
246 0 : ybeg[model].resize(nchan);
247 0 : yend[model].resize(nchan);
248 0 : for (Int chan=0; chan < nchan; ++chan){
249 0 : xbeg[model][chan]=0;
250 0 : ybeg[model][chan]=0;
251 0 : xend[model][chan]=nx-1;
252 0 : yend[model][chan]=ny-1;
253 0 : lmaskCube[chan*nmodels+model]=0;
254 : }
255 0 : if(hasMask(model)){
256 0 : Int mx=mask(model).shape()(0);
257 0 : Int my=mask(model).shape()(1);
258 0 : Int mpol=mask(model).shape()(2);
259 0 : if((mx != nx) || (my != ny) || (mpol != npol)){
260 0 : throw(AipsError("Mask image shape is not the same as dirty image"));
261 : }
262 0 : LatticeStepper mls(mask(model).shape(),
263 0 : IPosition(4, nx, ny, 1, 1),
264 0 : IPosition(4, 0, 1, 2, 3));
265 0 : RO_LatticeIterator<Float>maskli(mask(model), mls);
266 0 : maskli.reset();
267 0 : if(nchan > 1){
268 0 : if(nchan==mask(model).shape()(3)){
269 0 : isCubeMask[model]=true;
270 0 : for (Int chan=0; chan < nchan; ++chan){
271 0 : if(lmaskCube[chan*nmodels+model])
272 0 : delete lmaskCube[chan*nmodels+model];
273 0 : lmaskCube[chan*nmodels+model]=0;
274 0 : lmaskCube[chan*nmodels+model]=makeMaskMatrix(nx,ny, maskli,
275 0 : xbeg[model][chan],
276 0 : xend[model][chan],
277 0 : ybeg[model][chan],
278 0 : yend[model][chan]);
279 0 : maskli++;
280 : }
281 : }
282 : else{
283 : //Using first plane mask
284 0 : if(lmaskCube[model])
285 0 : delete lmaskCube[model];
286 0 : lmaskCube[model]=0;
287 0 : lmaskCube[model]=makeMaskMatrix(nx,ny, maskli,
288 0 : xbeg[model][0],
289 0 : xend[model][0],
290 0 : ybeg[model][0],
291 0 : yend[model][0]);
292 0 : for (Int chan=1; chan < nchan; ++chan){
293 0 : xbeg[model][chan]=xbeg[model][0];
294 0 : ybeg[model][chan]=ybeg[model][0];
295 0 : xend[model][chan]=xend[model][0];
296 0 : yend[model][chan]=yend[model][0];
297 : }
298 :
299 : }
300 : }
301 : else{
302 0 : if(lmaskCube[model])
303 0 : delete lmaskCube[model];
304 0 : lmaskCube[model]=0;
305 0 : lmaskCube[model]=makeMaskMatrix(nx,ny, maskli,
306 0 : xbeg[model][0],
307 0 : xend[model][0],
308 0 : ybeg[model][0],
309 0 : yend[model][0]);
310 :
311 :
312 : }
313 0 : }
314 :
315 :
316 :
317 : }
318 :
319 :
320 0 : Vector<Float> maxggS(nmodels,0.0);
321 0 : Bool lastCycleWriteModel=false;
322 :
323 0 : while(absmax>=threshold()&&maxIterations<numberIterations()&&!stop && !diverging) {
324 : os << LogIO::NORMAL
325 : << "*** Starting major cycle " << cycle
326 0 : << LogIO::POST; // Loglevel PROGRESS
327 0 : cycle++;
328 :
329 : // Make the residual images. We do an incremental update
330 : // for cycles after the first one. If we have only one
331 : // model then we use convolutions to speed the processing
332 0 : os << LogIO::NORMAL2 << "Making residual images for all fields" << LogIO::POST; // Loglevel PROGRESS
333 0 : if(modified_p){
334 0 : blankOverlappingModels();
335 0 : makeNewtonRaphsonStep(se, false, (numberIterations()<1)?true:False);
336 : //makeNewtonRaphsonStep(se, (cycle>1));
337 : }
338 :
339 0 : if(numberIterations() < 1){
340 : // Why waste the time to set up
341 0 : return true;
342 : }
343 :
344 :
345 : if(0) {
346 : ostringstream name;
347 : name << "Residual" << cycle;
348 : PagedImage<Float> thisCycleImage(residual(0).shape(),
349 : residual(0).coordinates(),
350 : name.str());
351 : thisCycleImage.copyData(residual(0));
352 : }
353 :
354 0 : if(cycle==1) {
355 0 : for (Int model=0;model<numberOfModels();model++) {
356 0 : LatticeExprNode LEN = max(LatticeExpr<Float>(ggS(model)));
357 0 : maxggS[model] = LEN.getFloat();
358 :
359 : if(0) {
360 : ostringstream name;
361 : name << "ggS";
362 : PagedImage<Float> thisImage(residual(0).shape(),
363 : residual(0).coordinates(),
364 : name.str());
365 : thisImage.copyData(ggS(0));
366 : }
367 0 : }
368 : os << LogIO::NORMAL1 // Loglevel INFO
369 0 : << "Maximum sensitivity = " << 1.0/sqrt(max(maxggS))
370 0 : << " Jy/beam" << LogIO::POST;
371 : }
372 :
373 : if(0) {
374 : ostringstream name;
375 : name << "ggSxResidual" << cycle;
376 : PagedImage<Float> thisCycleImage(residual(0).shape(),
377 : residual(0).coordinates(),
378 : name.str());
379 : thisCycleImage.copyData(LatticeExpr<Float>(sqrt(ggS(0)/maxggS[0])*residual(0)));
380 : }
381 :
382 0 : absmax=maxField(resmax, resmin);
383 0 : if(cycle >1){
384 : //check if its 5% above previous value
385 0 : if(absmax < 1.000005*oldabsmax)
386 0 : oldabsmax=absmax;
387 : else{
388 0 : diverging=true;
389 : os << LogIO::WARN
390 0 : << "Clean not converging " << LogIO::POST;
391 : }
392 : }
393 :
394 0 : os << LogIO::NORMAL1; // Loglevel INFO
395 0 : for (model=0;model<numberOfModels();model++) {
396 : os << "Model " << model << ": max, min (weighted) residuals (over all pixels) = "
397 0 : << resmax(model) << ", " << resmin(model) << endl;
398 : }
399 0 : os << LogIO::POST;
400 :
401 : // Can we stop?
402 0 : if(absmax<threshold()) {
403 : os << LogIO::NORMAL // Loglevel PROGRESS
404 0 : << "Reached stopping peak residual = " << absmax << LogIO::POST;
405 0 : stop=true;
406 0 : if(cycle >1)
407 0 : lastCycleWriteModel=true;
408 : }
409 0 : else if(!diverging) {
410 :
411 :
412 : // Calculate the threshold for this cycle. Add a safety factor
413 :
414 : // fractionOfPsf controls how deep the cleaning should go.
415 : // There arse two user-controls.
416 : // cycleFactor_p : scale factor for the PSF sidelobe level.
417 : // 1 : clean down to the psf sidelobe level
418 : // <1 : go deeper
419 : // >1 : shallower : stop sooner.
420 : // Default : 1.5
421 : // cycleMaxPsfFraction_p : scale factor as a fraction of the PSF peak
422 : // must be 0.0 < xx < 1.0 (obviously)
423 : // Default : 0.8
424 0 : Float fractionOfPsf = min(cycleMaxPsfFraction_p, cycleFactor_p * maxSidelobe);
425 0 : if (fractionOfPsf > (Float)0.8)
426 : {
427 : // os << LogIO::WARN << "PSF fraction for threshold computation is too high : " << fractionOfPsf << ". Forcing to 0.8 to ensure that the threshold is smaller than the peak residual !" << LogIO::POST;
428 0 : os << LogIO::NORMAL << "Current values of max-PSF-fraction, cycle-factor and max PSF sidelobe level result in a stopping threshold more than 80% of the peak residual. Forcing the maximum threshold to 80% of the peak." << LogIO::POST;
429 0 : fractionOfPsf = 0.8; // painfully slow!
430 : }
431 :
432 0 : os << LogIO::NORMAL << "The minor-cycle threshold is MAX[ 0.95 x " << threshold() << " , peak residual x " << fractionOfPsf << " ] " << LogIO::POST;
433 :
434 0 : cycleThreshold=max(0.95*threshold(), fractionOfPsf * absmax);
435 : os << LogIO::NORMAL << "Maximum residual = " << absmax // Loglevel INFO
436 0 : << ", cleaning down to " << cycleThreshold << LogIO::POST;
437 :
438 0 : for (model=0;model<numberOfModels();model++) {
439 :
440 0 : nx=image(model).shape()(0);
441 0 : ny=image(model).shape()(1);
442 0 : npol=image(model).shape()(2);
443 0 : nchan=image(model).shape()(3);
444 :
445 0 : Int npolcube=npol;
446 0 : if(!doPolJoint_p){
447 0 : npolcube=1;
448 : }
449 0 : if(cycle==1) {
450 0 : if(doPolJoint_p)
451 0 : iterations[model].resize(nchan);
452 : else
453 0 : iterations[model].resize(nchan*npol);
454 0 : iterations[model].resize(nchan*npol);
455 0 : iterations[model]=0;
456 : }
457 :
458 : // Initialize delta image
459 0 : deltaImage(model).set(0.0);
460 :
461 : // Only process solveable models
462 0 : if(isSolveable(model)) {
463 :
464 :
465 0 : if((abs(resmax(model))>cycleThreshold)||
466 0 : (abs(resmin(model))>cycleThreshold)) {
467 :
468 : os << LogIO::NORMAL
469 : << "Processing model " << model
470 0 : << LogIO::POST; // Loglevel PROGRESS
471 :
472 0 : IPosition onePlane(4, nx, ny, 1, 1);
473 0 : IPosition oneCube(4, nx, ny, npolcube, 1);
474 :
475 : //AlwaysAssert((npol==1)||(npol==2)||(npol==4), AipsError);
476 :
477 : // Loop over all channels. We only want the PSF for the first
478 : // polarization so we iterate over polarization LAST
479 :
480 0 : LatticeStepper psfls(PSF(model).shape(), onePlane,
481 0 : IPosition(4,0,1,3,2));
482 0 : RO_LatticeIterator<Float> psfli(PSF(model),psfls);
483 :
484 0 : LatticeStepper fsls(ggS(model).shape(), oneCube,
485 0 : IPosition(4,0,1,3,2));
486 0 : RO_LatticeIterator<Float> ggSli(ggS(model),fsls);
487 :
488 0 : LatticeStepper ls(image(model).shape(), oneCube,
489 0 : IPosition(4, 0, 1, 3, 2));
490 0 : LatticeIterator<Float> imageStepli(residual(model), ls);
491 0 : LatticeIterator<Float> imageli(image(model), ls);
492 0 : LatticeIterator<Float> deltaimageli(deltaImage(model), ls);
493 :
494 : // If we are using a mask: assign the first ones
495 0 : Matrix<Float> maskArray;
496 0 : Bool domask=false;
497 0 : if(hasMask(model)) {
498 0 : maskArray.resize();
499 : //maskArray.resize(lmaskCube[model]->shape());
500 0 : maskArray.reference(*(lmaskCube[model]));
501 0 : domask=true;
502 : }
503 :
504 : // Now clean each channel
505 0 : Int chan=0;
506 0 : Int ipol=0;
507 0 : for (imageStepli.reset(),imageli.reset(),psfli.reset(),
508 0 : deltaimageli.reset(),ggSli.reset();
509 0 : !imageStepli.atEnd();
510 0 : imageStepli++,imageli++,psfli++,deltaimageli++,
511 0 : ggSli++,chan++) {
512 0 : if(!doPolJoint_p){
513 0 : if(chan==0){
514 : os << LogIO::NORMAL // Loglevel PROGRESS
515 0 : << "Doing stokes "<< stokesID(ipol) << " image" <<LogIO::POST;
516 : }
517 0 : if(chan==nchan){
518 0 : chan=0;
519 0 : psfli.reset();
520 0 : ++ipol;
521 : os << LogIO::NORMAL2 // Loglevel PROGRESS
522 0 : << "Doing stokes "<< stokesID(ipol) << " image" <<LogIO::POST;
523 : }
524 : }
525 0 : if(hasMask(model) && isCubeMask[model] && chan >0) {
526 0 : maskArray.reference(*(lmaskCube[chan*nmodels+model]));
527 : }
528 :
529 0 : if(psfmax(model)>0.0) {
530 :
531 0 : if(nchan>1) {
532 : os << LogIO::NORMAL // Loglevel PROGRESS
533 0 : << "Processing channel "<<chan<<" of "<<nchan<<LogIO::POST;
534 : }
535 :
536 : // Renormalize by the weights
537 0 : Matrix<Float> wmask(nx, ny, 0.0);
538 : //Trouble..
539 0 : Array<Float> residu(imageStepli.cursor());
540 0 : Cube<Float> resid(residu.nonDegenerate(3));
541 0 : Cube<Float> weight(sqrt(ggSli.cursor().nonDegenerate(3)/maxggS[model]));
542 :
543 :
544 0 : for (Int iy=0;iy<ny;iy++) {
545 0 : for (Int ix=0;ix<nx;ix++) {
546 0 : for (Int pol=0;pol<npolcube;pol++) {
547 : //resid(ix,iy,pol)*=weight(ix,iy,pol);
548 0 : if(weight(ix,iy,pol)>0.01) {
549 0 : wmask(ix,iy)=1.0;
550 : }
551 : }
552 : }
553 : }
554 0 : if (domask) {
555 0 : wmask*=maskArray;
556 : }
557 : Bool delete_its;
558 0 : Float* limageStep_data=resid.getStorage(delete_its);
559 :
560 0 : if (itsSubAlgorithm == "hogbom") {
561 : Bool delete_itdi, delete_itp, delete_itm;
562 : Float *ldeltaimage_data;
563 0 : ldeltaimage_data=deltaimageli.rwCursor().getStorage(delete_itdi);
564 0 : const Float *lpsf_data(0), *lmask_data(0);
565 0 : lmask_data=wmask.getStorage(delete_itm);
566 0 : lpsf_data=psfli.cursor().getStorage(delete_itp);
567 0 : Int niter=numberIterations();
568 0 : Int starting_iteration = iterations[model](chan*npolcube+ipol);
569 : Int ending_iteration;
570 :
571 0 : Float g=gain();
572 0 : Float thres=cycleThreshold;
573 0 : Int fxbeg=xbeg[model][chan]+1;
574 0 : Int fxend=xend[model][chan]+1;
575 0 : Int fybeg=ybeg[model][chan]+1;
576 0 : Int fyend=yend[model][chan]+1;
577 :
578 0 : Int domaskI=1;
579 0 : hclean(ldeltaimage_data, limageStep_data, (Float*)lpsf_data,
580 : &domaskI, (Float*)lmask_data,
581 : &nx, &ny, &npolcube,
582 : &fxbeg, &fxend, &fybeg, &fyend, &niter,
583 : &starting_iteration, &ending_iteration,
584 : &g, &thres, &cycleSpeedup_p,
585 : (void*) &HogbomCleanImageSkyModelmsgput,
586 : (void*) &HogbomCleanImageSkyModelstopnow);
587 :
588 : os << LogIO::NORMAL // Loglevel PROGRESS
589 0 : << "Finished Hogbom clean inner cycle " << LogIO::POST;
590 :
591 0 : imageStepli.rwCursor().putStorage (limageStep_data, delete_its);
592 0 : deltaimageli.rwCursor().putStorage (ldeltaimage_data, delete_itdi);
593 :
594 0 : Cube<Float> step(deltaimageli.rwCursor().nonDegenerate(3));
595 : ///Need to go
596 : /*
597 : for (Int iy=0;iy<ny;iy++) {
598 : for (Int ix=0;ix<nx;ix++) {
599 : for (Int pol=0;pol<npol;pol++) {
600 : if(weight(ix,iy,pol)>0.0) {
601 : step(ix,iy,pol)/=weight(ix,iy,pol);
602 : }
603 : }
604 : }
605 : }
606 : deltaimageli.rwCursor()=step.addDegenerate(1);
607 : */
608 : /////
609 : //imageli.rwCursor()+=deltaimageli.cursor();
610 0 : psfli.cursor().freeStorage (lpsf_data, delete_itp);
611 0 : if (domask) {
612 0 : maskArray.freeStorage (lmask_data, delete_itm);
613 : }
614 0 : iterations[model](chan*npolcube+ipol) = ending_iteration;
615 0 : maxIterations=(iterations[model](chan*npolcube+ipol)>maxIterations) ?
616 0 : iterations[model](chan*npolcube+ipol) : maxIterations;
617 0 : modified_p=true;
618 0 : } else { // clark is the default for now
619 :
620 : // Now make a convolution equation for this
621 : // residual image and psf and then clean it
622 0 : ConvolutionEquation eqn(psfli.matrixCursor(),
623 0 : residu);
624 0 : deltaimageli.rwCursor().set(Float(0.0));
625 0 : ClarkCleanModel cleaner(deltaimageli.rwCursor());
626 0 : cleaner.setMask(wmask);
627 0 : cleaner.setGain(gain());
628 0 : cleaner.setNumberIterations(numberIterations());
629 0 : cleaner.setInitialNumberIterations(iterations[model](chan*npolcube+ipol));
630 0 : cleaner.setThreshold(cycleThreshold);
631 0 : cleaner.setPsfPatchSize(IPosition(2,psfpatch));
632 0 : cleaner.setMaxNumberMajorCycles(10);
633 0 : cleaner.setHistLength(1024);
634 0 : cleaner.setMaxNumPix(32*1024);
635 0 : cleaner.setChoose(false);
636 0 : cleaner.setCycleSpeedup(cycleSpeedup_p);
637 0 : cleaner.setSpeedup(0.0);
638 0 : if ( displayProgress_p ) {
639 0 : cleaner.setProgress( *progress_p );
640 : }
641 0 : cleaner.singleSolve(eqn, residu);
642 0 : iterations[model](chan*npolcube+ipol)=cleaner.numberIterations();
643 0 : maxIterations=(iterations[model](chan*npolcube+ipol)>maxIterations) ?
644 0 : iterations[model](chan*npolcube+ipol) : maxIterations;
645 0 : modified_p=true;
646 :
647 0 : cleaner.getModel(deltaimageli.rwCursor());
648 : //imageli.rwCursor()+=deltaimageli.cursor();
649 : //eqn.residual(imageStepli.rwCursor(), cleaner);
650 :
651 : os << LogIO::NORMAL // Loglevel PROGRESS
652 0 : <<"Finished Clark clean inner cycle " << LogIO::POST;
653 0 : }
654 0 : if(maxIterations==0) {
655 0 : stop=true;
656 : }
657 : else{
658 0 : stop=false;
659 : }
660 : // Upped to NORMAL per CAS-2017.
661 : // cerr << model << " " << chan << " " << npolcube << " " << ipol
662 : // << " " << iterations.nelements()
663 : // << " " << iterations[0].shape()
664 : // << endl;
665 0 : if (iterations[model](chan*npolcube+ipol) > 0) {
666 : os << LogIO::NORMAL << "Clean used " // Loglevel PROGRESS
667 0 : << iterations[model](chan*npolcube+ipol) << " iterations"
668 : << " to approach a threshhold of " << cycleThreshold
669 0 : << LogIO::POST;
670 : }
671 0 : }
672 : }
673 0 : }
674 : else {
675 : os << LogIO::NORMAL // Loglevel INFO
676 : << "No need to clean model " << model
677 : << ": peak residual below threshhold"
678 0 : << LogIO::POST;
679 : }
680 : }
681 : }
682 :
683 0 : if(maxIterations != oldMaxIterations){
684 0 : oldMaxIterations=maxIterations;
685 0 : for(Int model=0; model < numberOfModels(); ++model){
686 0 : image(model).copyData( LatticeExpr<Float>((image(model))+(deltaImage(model))));
687 0 : os << LogIO::NORMAL << LatticeExprNode(sum(deltaImage(model))).getFloat() // Loglevel INFO
688 : << " Jy <- cleaned in this cycle for model "
689 : << model
690 0 : << " (Total flux : " << LatticeExprNode(sum(image(model))).getFloat() << "Jy)"
691 0 : << LogIO::POST;
692 :
693 : }
694 :
695 : }
696 : else {
697 : os << LogIO::NORMAL
698 : << "No more iterations left in this major cycle - stopping now"
699 0 : << LogIO::POST;
700 0 : stop=true;
701 0 : converged=true;
702 : }
703 : }
704 : if(0) {
705 : ostringstream name;
706 : name << "Delta" << cycle;
707 : PagedImage<Float> thisCycleDeltaImage(residual(1).shape(),
708 : residual(1).coordinates(),
709 : name.str());
710 : thisCycleDeltaImage.copyData(deltaImage(1));
711 : }
712 : }
713 0 : if (progress_p) delete progress_p;
714 :
715 0 : for(model=0; model < nmodels; ++model){
716 0 : image(model).clearCache();
717 0 : PSF(model).clearCache();
718 0 : ggS(model).clearCache();
719 0 : residual(model).clearCache();
720 0 : deltaImage(model).clearCache();
721 0 : for(Int chan=0; chan < nchan; ++chan){
722 0 : if(lmaskCube[chan*nmodels+model] !=0)
723 0 : delete lmaskCube[chan*nmodels+model];
724 : }
725 : }
726 :
727 :
728 0 : setNumberIterations(maxIterations);
729 0 : if(modified_p || lastCycleWriteModel) {
730 : os << LogIO::NORMAL2 // Loglevel PROGRESS
731 0 : << "Finalizing residual images for all fields" << LogIO::POST;
732 0 : blankOverlappingModels();
733 0 : makeNewtonRaphsonStep(se, false, true); //committing model to MS
734 0 : restoreOverlappingModels();
735 0 : Float finalabsmax=maxField(resmax, resmin);
736 0 : converged=(finalabsmax < 1.05 * threshold());
737 0 : setThreshold(finalabsmax);
738 0 : os << LogIO::NORMAL << "Final maximum residual = " << finalabsmax << LogIO::POST; // Loglevel INFO
739 :
740 0 : os << LogIO::NORMAL; // Loglevel INFO
741 0 : for (model=0;model<numberOfModels();model++) {
742 : os << "Model " << model << ": max, min residuals = "
743 0 : << resmax(model) << ", " << resmin(model) << " clean flux " << LatticeExprNode(sum(image(model))).getFloat() << endl;
744 : }
745 0 : }
746 : else {
747 : os << LogIO::NORMAL2 // Loglevel PROGRESS
748 0 : << "Residual images for all fields are up-to-date" << LogIO::POST;
749 : }
750 :
751 0 : os << LogIO::POST;
752 0 : if(stop)
753 0 : converged=true;
754 :
755 0 : return(converged);
756 0 : };
757 :
758 :
759 0 : void MFCleanImageSkyModel::blankOverlappingModels(){
760 0 : if(numberOfModels() == 1)
761 0 : return;
762 : /////////////////////
763 : /*
764 : for (Int model=0;model<(numberOfModels()); ++model) {
765 :
766 :
767 : LatticeExprNode maxIm=max(deltaImage(model));
768 : cout << "MAX model " << model << " " << maxIm.getFloat() << endl;
769 : }
770 : */
771 : //////////////
772 0 : for (Int model=0;model<(numberOfModels()-1); ++model) {
773 0 : CoordinateSystem cs0=image(model).coordinates();
774 0 : IPosition iblc0(image(model).shape().nelements(),0);
775 :
776 0 : IPosition itrc0(image(model).shape());
777 0 : itrc0=itrc0-Int(1);
778 0 : LCBox lbox0(iblc0, itrc0, image(model).shape());
779 :
780 0 : ImageRegion imagreg0(WCBox(lbox0, cs0));
781 0 : for (Int nextmodel=model+1; nextmodel < numberOfModels(); ++nextmodel){
782 0 : CoordinateSystem cs=image(nextmodel).coordinates();
783 0 : IPosition iblc(image(nextmodel).shape().nelements(),0);
784 :
785 0 : IPosition itrc(image(nextmodel).shape());
786 0 : itrc=itrc-Int(1);
787 :
788 0 : LCBox lbox(iblc, itrc, image(nextmodel).shape());
789 :
790 0 : ImageRegion imagreg(WCBox(lbox, cs));
791 :
792 :
793 : try{
794 0 : LatticeRegion latReg=imagreg.toLatticeRegion(image(model).coordinates(), image(model).shape());
795 :
796 0 : SubImage<Float> partToMask(image(model), imagreg, true);
797 0 : ArrayLattice<Bool> pixmask(latReg.get());
798 0 : LatticeExpr<Float> myexpr(iif(pixmask, 0.0, partToMask) );
799 0 : partToMask.copyData(myexpr);
800 :
801 0 : }
802 0 : catch(...){
803 : //no overlap you think ?
804 : /*
805 : os << LogIO::WARN
806 : << "no overlap or failure of copying the clean components"
807 : << x.getMesg()
808 : << LogIO::POST;
809 : */
810 0 : continue;
811 0 : }
812 :
813 :
814 0 : }
815 :
816 :
817 :
818 0 : }
819 : }
820 :
821 : // test method - reversing blanking of overlapped regions
822 0 : void MFCleanImageSkyModel::restoreOverlappingModels(){
823 0 : LogIO os(LogOrigin("MFCleanImageSkyModel","restoreOverlappingModels"));
824 0 : if(numberOfModels() == 1)
825 0 : return;
826 0 : for (Int model=0;model<(numberOfModels()-1); ++model) {
827 0 : CoordinateSystem cs0=image(model).coordinates();
828 0 : IPosition iblc0(image(model).shape().nelements(),0);
829 0 : IPosition itrc0(image(model).shape());
830 0 : itrc0=itrc0-Int(1);
831 0 : LCBox lbox0(iblc0, itrc0, image(model).shape());
832 0 : ImageRegion imagreg0(WCBox(lbox0, cs0));
833 :
834 0 : for (Int nextmodel=model+1; nextmodel < numberOfModels(); ++nextmodel){
835 0 : CoordinateSystem cs=image(nextmodel).coordinates();
836 0 : IPosition iblc(image(nextmodel).shape().nelements(),0);
837 :
838 0 : IPosition itrc(image(nextmodel).shape());
839 0 : itrc=itrc-Int(1);
840 :
841 0 : LCBox lbox(iblc, itrc, image(nextmodel).shape());
842 :
843 0 : ImageRegion imagreg(WCBox(lbox, cs));
844 : try{
845 0 : LatticeRegion latReg0=imagreg0.toLatticeRegion(image(nextmodel).coordinates(), image(nextmodel).shape());
846 0 : LatticeRegion latReg=imagreg.toLatticeRegion(image(model).coordinates(), image(model).shape());
847 :
848 :
849 0 : SubImage<Float> partToMerge(image(nextmodel), imagreg0, true);
850 0 : SubImage<Float> partToUnmask(image(model), imagreg, true);
851 :
852 0 : ArrayLattice<Bool> pixmask(latReg.get());
853 0 : LatticeExpr<Float> myexpr0(iif(pixmask,partToMerge,partToUnmask));
854 0 : partToUnmask.copyData(myexpr0);
855 :
856 0 : }
857 0 : catch(AipsError x){
858 : /*
859 : os << LogIO::WARN
860 : << "no overlap or failure of copying the clean components"
861 : << x.getMesg()
862 : << LogIO::POST;
863 : */
864 0 : continue;
865 0 : }
866 0 : }
867 0 : }
868 0 : }
869 :
870 0 : void MFCleanImageSkyModel::mergeOverlappingMasks(){
871 :
872 : using casacore::operator<;
873 : using casacore::max;
874 :
875 0 : LogIO os(LogOrigin("MFCleanImageSkyModel","restoreOverlappingModels"));
876 0 : if(numberOfModels() == 1)
877 0 : return;
878 0 : for (Int model=0;model<(numberOfModels()-1); ++model) {
879 0 : if(hasMask(model)){
880 0 : CoordinateSystem cs0=mask(model).coordinates();
881 0 : IPosition iblc0(mask(model).shape().nelements(),0);
882 :
883 0 : IPosition itrc0(mask(model).shape());
884 0 : itrc0=itrc0-Int(1);
885 0 : LCBox lbox0(iblc0, itrc0, mask(model).shape());
886 :
887 0 : ImageRegion imagreg0(WCBox(lbox0, cs0));
888 0 : for (Int nextmodel=model+1; nextmodel < numberOfModels(); ++nextmodel){
889 0 : if(hasMask(nextmodel)){
890 0 : CoordinateSystem cs=mask(nextmodel).coordinates();
891 0 : IPosition iblc(mask(nextmodel).shape().nelements(),0);
892 :
893 0 : IPosition itrc(mask(nextmodel).shape());
894 0 : itrc=itrc-Int(1);
895 :
896 0 : LCBox lbox(iblc, itrc, image(nextmodel).shape());
897 :
898 0 : ImageRegion imagreg(WCBox(lbox, cs));
899 :
900 : try{
901 :
902 : LatticeRegion latRef1 = imagreg0.toLatticeRegion(
903 0 : (mask(nextmodel)).coordinates(),
904 0 : (mask(nextmodel)).shape() );
905 : LatticeRegion latRef2 = imagreg.toLatticeRegion(
906 0 : (mask(model)).coordinates(),
907 0 : (mask(model)).shape() );
908 :
909 0 : SubImage<Float> partToMerge(mask(nextmodel), imagreg0, true);
910 0 : SubImage<Float> partToMask(mask(model), imagreg, true);
911 0 : LatticeExpr<Float> myexpr0(max(partToMerge,partToMask));
912 0 : partToMask.copyData(myexpr0);
913 0 : partToMerge.copyData(myexpr0);
914 :
915 0 : }
916 0 : catch(AipsError x){
917 : //most probably no overlap
918 : /*
919 : os << LogIO::WARN
920 : << "no overlap or failure of copying the clean components"
921 : << x.getMesg()
922 : << LogIO::POST;
923 : */
924 :
925 0 : continue;
926 0 : }
927 0 : }
928 : }// for model
929 0 : }//hasmask(model)
930 : }
931 0 : }
932 :
933 :
934 :
935 0 : Float MFCleanImageSkyModel::maxOuter(Lattice<Float> & lat, const uInt nCenter )
936 : {
937 0 : Float myMax=0.0;
938 0 : Float myMin=0.0;
939 :
940 : /*
941 : TempLattice<Float> mask(lat.shape());
942 : mask.set(1.0);
943 :
944 : IPosition pos(4,0,0,0,0 );
945 : uInt nxc = lat.shape()(0)/2;
946 : uInt nyc = lat.shape()(1)/2;
947 : for (uInt ix = -nCenter; ix < nCenter; ix++) {
948 : for (uInt iy = -nCenter; iy < nCenter; iy++) {
949 : pos(0) = nxc + ix;
950 : pos(1) = nyc + iy;
951 : mask.putAt( 0.0f, pos ); // mask out the inner section
952 : }
953 : }
954 : {
955 : LatticeExpr<Float> exp = (lat * mask);
956 : LatticeExprNode LEN = max( exp );
957 : myMax = LEN.getFloat();
958 : }
959 : {
960 : LatticeExpr<Float> exp = (lat * mask);
961 : LatticeExprNode LEN = min( exp );
962 : myMin = LEN.getFloat();
963 : }
964 :
965 : */
966 :
967 0 : Array<Float> arr = lat.get();
968 0 : IPosition pos( arr.shape() );
969 0 : uInt nx = lat.shape()(0);
970 0 : uInt ny = lat.shape()(1);
971 0 : uInt nxc = 0;
972 0 : uInt nyc = 0;
973 0 : Float amax = 0.0;
974 0 : for (uInt ix = 0; ix < nx; ix++) {
975 0 : for (uInt iy = 0; iy < ny; iy++) {
976 0 : if (arr(IPosition(4, ix, iy, 0, 0)) > amax) {
977 0 : nxc = ix;
978 0 : nyc = iy;
979 0 : amax = arr(IPosition(4, ix, iy, 0, 0));
980 : }
981 : }
982 : }
983 :
984 0 : uInt nxL = nxc - nCenter;
985 0 : uInt nxH = nxc + nCenter;
986 0 : uInt nyL = nyc - nCenter;
987 0 : uInt nyH = nyc + nCenter;
988 :
989 0 : for (uInt ix = 0; ix < nx; ix++) {
990 0 : for (uInt iy = 0; iy < ny; iy++) {
991 0 : if ( !(ix >= nxL && ix <= nxH && iy >= nyL && iy <= nyH) ) {
992 0 : if (arr(IPosition(4, ix, iy, 0, 0)) > myMax)
993 0 : myMax = arr(IPosition(4, ix, iy, 0, 0));
994 0 : if (arr(IPosition(4, ix, iy, 0, 0)) < myMin)
995 0 : myMin = arr(IPosition(4, ix, iy, 0, 0));
996 : }
997 : }
998 : }
999 :
1000 0 : Float absMax = max( abs(myMin), myMax );
1001 0 : return absMax;
1002 0 : };
1003 :
1004 0 : Matrix<Float>* MFCleanImageSkyModel::makeMaskMatrix(const Int& nx,
1005 : const Int& ny,
1006 : RO_LatticeIterator<Float>& maskIter,
1007 : Int& xbeg,
1008 : Int& xend,
1009 : Int& ybeg,
1010 : Int& yend) {
1011 :
1012 0 : LogIO os(LogOrigin("MFCleanImageSkyModel","makeMaskMatrix",WHERE));
1013 :
1014 0 : xbeg=0;
1015 0 : ybeg=0;
1016 0 : xend=nx-1;
1017 0 : yend=ny-1;
1018 :
1019 0 : Matrix<Float>* mask= new Matrix<Float>(maskIter.matrixCursor().shape());
1020 0 : (*mask)=maskIter.matrixCursor();
1021 : // ignore mask if none exists
1022 :
1023 0 : if(max(*mask) < 0.5) {
1024 : //If empty mask respect it for MF.... till we break WF from it
1025 0 : return mask;
1026 : }
1027 : // Now read the mask and determine the bounding box
1028 0 : xbeg=xend;
1029 0 : ybeg=yend;
1030 0 : yend=0;
1031 0 : xend=0;
1032 :
1033 0 : for (Int iy=0;iy<ny;iy++) {
1034 0 : for (Int ix=0;ix<nx;ix++) {
1035 0 : if((*mask)(ix,iy)>0.5) {
1036 0 : xbeg=min(xbeg,ix);
1037 0 : ybeg=min(ybeg,iy);
1038 0 : xend=max(xend,ix);
1039 0 : yend=max(yend,iy);
1040 : }
1041 : }
1042 : }
1043 :
1044 0 : return mask;
1045 0 : }
1046 :
1047 : } //#End casa namespace
|
