Line data Source code
1 : //# CSCleanImageSkyModel.cc: Implementation of CSCleanImageSkyModel 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 <casacore/casa/Arrays/Matrix.h>
30 : #include <synthesis/MeasurementComponents/CSCleanImageSkyModel.h>
31 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
32 : #include <casacore/images/Images/PagedImage.h>
33 : #include <casacore/casa/OS/File.h>
34 : #include <casacore/lattices/LEL/LatticeExpr.h>
35 : #include <casacore/lattices/LEL/LatticeExprNode.h>
36 : #include <casacore/lattices/Lattices/LatticeStepper.h>
37 : #include <casacore/lattices/Lattices/LatticeIterator.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 :
43 : #include <sstream>
44 :
45 : #include <casacore/casa/Logging/LogMessage.h>
46 : #include <casacore/casa/Logging/LogSink.h>
47 : #include <casacore/casa/Logging/LogIO.h>
48 :
49 : #include <msvis/MSVis/StokesVector.h>
50 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
51 : #include <synthesis/MeasurementEquations/ClarkCleanLatModel.h>
52 : #include <casacore/lattices/Lattices/SubLattice.h>
53 : #include <casacore/lattices/LRegions/LCBox.h>
54 :
55 : using namespace casacore;
56 : namespace casa {
57 :
58 0 : Int CSCleanImageSkyModel::add(ImageInterface<Float>& image,
59 : const Int maxNumXfr)
60 : {
61 0 : return CleanImageSkyModel::add(image, maxNumXfr);
62 : };
63 :
64 0 : Bool CSCleanImageSkyModel::addMask(Int image,
65 : ImageInterface<Float>& mask)
66 : {
67 0 : return CleanImageSkyModel::addMask(image, mask);
68 : };
69 :
70 0 : Bool CSCleanImageSkyModel::addResidual(Int image,
71 : ImageInterface<Float>& residual)
72 : {
73 0 : return ImageSkyModel::addResidual(image, residual);
74 : };
75 :
76 : // Clean solver
77 0 : Bool CSCleanImageSkyModel::solve(SkyEquation& se) {
78 :
79 0 : LogIO os(LogOrigin("CSCleanImageSkyModel","solve"));
80 0 : Bool converged=true;
81 0 : if(modified_p) {
82 0 : makeNewtonRaphsonStep(se, false, (numberIterations()<1)?true:False);
83 : }
84 :
85 0 : if( numberIterations() < 1)
86 0 : return true;
87 : //Make the PSFs, one per field
88 :
89 : os << LogIO::NORMAL // Loglevel PROGRESS
90 0 : << "Making approximate Point Spread Functions" << LogIO::POST;
91 0 : if(!donePSF_p)
92 0 : makeApproxPSFs(se);
93 : //
94 : // Quite a few lines of code required to pull out co-ordinate info.
95 : // from an image, one would think.
96 : //
97 0 : CoordinateSystem psfCoord=PSF(0).coordinates();
98 0 : Int dirIndex = psfCoord.findCoordinate(Coordinate::DIRECTION);
99 0 : DirectionCoordinate dc=psfCoord.directionCoordinate(dirIndex);
100 0 : Vector<Double> incr=dc.increment();
101 : //
102 : // The fitted beam params. are in arcsec. Increments returned
103 : // by the coordinate system are in radians.
104 : //
105 0 : incr *= 3600.0*180.0/M_PI;
106 0 : incr = abs(incr);
107 :
108 : // Validate PSFs for each field
109 0 : Vector<Float> psfmax(numberOfModels()); psfmax=0.0;
110 0 : Vector<Float> psfmaxouter(numberOfModels()); psfmaxouter=0.0;
111 0 : Vector<Float> psfmin(numberOfModels()); psfmin=1.0;
112 0 : Block<Vector<Float> > resmax(numberOfModels());
113 0 : Block<Vector<Float> > resmin(numberOfModels());
114 :
115 0 : Float maxSidelobe=0.0;
116 : Int model;
117 0 : os << LogIO::NORMAL1; // Loglevel INFO
118 0 : for (model=0;model<numberOfModels();model++) {
119 0 : if(isSolveable(model)) {
120 :
121 0 : IPosition blc(PSF(model).shape().nelements(), 0);
122 0 : IPosition trc(PSF(model).shape()-1);
123 0 : blc(2) = 0; trc(2) = 0;
124 0 : blc(3) = 0; trc(3) = 0;
125 :
126 0 : SubLattice<Float> subPSF;
127 0 : Int k =0;
128 0 : Int numchan= PSF(model).shape()(3);
129 : //PSF of first non zero plane
130 0 : while(psfmax(model) < 0.1 && k< numchan){
131 0 : blc(3)= k;
132 0 : trc(3)=k;
133 0 : LCBox onePlane(blc, trc, PSF(model).shape());
134 0 : subPSF=SubLattice<Float> ( PSF(model), onePlane, true);
135 : {
136 0 : LatticeExprNode node = max(subPSF);
137 0 : psfmax(model) = node.getFloat();
138 : }
139 0 : ++k;
140 : }
141 : // {
142 : // LatticeExprNode node = min(subPSF);
143 : // psfmin(model) = node.getFloat();
144 : // }
145 : // // 4 pixels: pretty arbitrary, but only look for sidelobes
146 : // // outside the inner (2n+1) * (2n+1) square
147 : // // Changed the algorithm now..so that 4 is not used
148 : // Int mainLobeSizeInPixels = (Int)(max(beam(0)[0]/incr[0],beam(0)[1]/incr[1]));
149 : // //psfmaxouter(model) = maxOuter(subPSF, 4);
150 : // psfmaxouter(model) = maxOuter(subPSF, mainLobeSizeInPixels);
151 :
152 : //
153 : // Since for CS-Clean anyway uses only a small fraction of the
154 : // inner part of the PSF matter, find the PSF outer
155 : // min/max. using only the inner quater of the PSF.
156 : //
157 0 : GaussianBeam elbeam0=beam(0)(0,0);
158 0 : Int mainLobeSizeInPixels = (Int)(max(elbeam0.getMajor("arcsec")/incr[0],elbeam0.getMinor("arcsec")/incr[1]));
159 0 : Vector<Float> psfOuterMinMax(2);
160 0 : psfOuterMinMax = outerMinMax(subPSF, mainLobeSizeInPixels);
161 0 : psfmaxouter(model)=psfOuterMinMax(1);
162 0 : psfmin(model) = psfOuterMinMax(0);
163 :
164 :
165 : os << "Model " << model+1 << ": Estimated size of the PSF mainlobe = "
166 0 : << (Int)(elbeam0.getMajor("arcsec")/incr[0]+0.5) << " X " << (Int)(elbeam0.getMinor("arcsec")/incr[1] + 0.5) << " pixels"
167 0 : << endl;
168 : os << "Model " << model+1 << ": PSF Peak, min, max sidelobe = "
169 0 : << psfmax(model) << ", " << psfmin(model) << ", " <<
170 0 : psfmaxouter(model) << endl;
171 0 : if(abs(psfmin(model))>maxSidelobe) maxSidelobe=abs(psfmin(model));
172 0 : if(psfmaxouter(model) > maxSidelobe) maxSidelobe= psfmaxouter(model );
173 : }
174 : }
175 0 : os << LogIO::POST;
176 :
177 0 : Float absmax=threshold();
178 0 : Float oldmax=absmax;
179 0 : Float cycleThreshold=0.0;
180 0 : Block< Vector<Int> > iterations(numberOfModels());
181 0 : Int maxIterations=0;
182 0 : Int oldMaxIterations=0;
183 :
184 : // Loop over major cycles
185 0 : Int cycle=0;
186 0 : Bool stop=false;
187 :
188 0 : if (displayProgress_p) {
189 0 : progress_p = new ClarkCleanProgress( pgplotter_p );
190 : } else {
191 0 : progress_p = 0;
192 : }
193 0 : Bool lastCycleWriteModel=false;
194 0 : while((absmax>=threshold())&& (maxIterations<numberIterations()) &&!stop) {
195 :
196 : os << LogIO::NORMAL << "*** Starting major cycle " << cycle + 1
197 0 : << LogIO::POST; // Loglevel PROGRESS
198 0 : cycle++;
199 :
200 : // Make the residual images. We do an incremental update
201 : // for cycles after the first one. If we have only one
202 : // model then we use convolutions to speed the processing
203 : //os << LogIO::NORMAL2 // Loglevel PROGRESS
204 : // << "Starting major cycle : making residual images for all fields"
205 : // << LogIO::POST;
206 0 : if(modified_p) {
207 0 : Bool incremental(cycle>1);
208 0 : if (incremental&&(itsSubAlgorithm == "fast")) {
209 : os << LogIO::NORMAL1 // Loglevel INFO
210 : << "Using XFR-based shortcut for residual calculation"
211 0 : << LogIO::POST;
212 0 : makeNewtonRaphsonStep(se, false);
213 : }
214 : else {
215 : os << LogIO::NORMAL1 // Loglevel INFO
216 : << "Using visibility-subtraction for residual calculation"
217 0 : << LogIO::POST;
218 0 : makeNewtonRaphsonStep(se, false);
219 : }
220 : os << LogIO::NORMAL2 // Loglevel PROGRESS
221 : << "Finished update of residuals"
222 0 : << LogIO::POST;
223 : }
224 :
225 0 : oldmax=absmax;
226 0 : absmax=maxField(resmax, resmin);
227 0 : if(cycle==1) oldmax=absmax;
228 :
229 0 : for (model=0;model<numberOfModels();model++) {
230 : os << LogIO::NORMAL // Loglevel INFO
231 : << "Model " << model+1 << ": max, min residuals = "
232 0 : << max(resmax[model]) << ", " << min(resmin[model]) << endl;
233 : }
234 0 : os << LogIO::POST;
235 :
236 : // Can we stop?
237 0 : if(absmax<threshold()) {
238 : os << LogIO::NORMAL // Loglevel INFO
239 0 : << "Reached stopping peak residual = " << absmax << LogIO::POST;
240 0 : stop=true;
241 0 : if(cycle >1)
242 0 : lastCycleWriteModel=true;
243 : }
244 : else {
245 0 : if(oldmax < absmax){
246 : //Diverging? Let's increase the cyclefactor
247 0 : cycleFactor_p=1.5*cycleFactor_p;
248 0 : oldmax=absmax;
249 : }
250 : // Calculate the threshold for this cycle. Add a safety factor
251 : // This will be fixed someday by an option for an increasing threshold
252 0 : Float fudge = cycleFactor_p * maxSidelobe;
253 0 : if (fudge > 0.8) fudge = 0.8; // painfully slow!
254 :
255 0 : cycleThreshold=max(threshold(), fudge * absmax);
256 : os << LogIO::NORMAL // Loglevel INFO
257 : << "Maximum residual = " << absmax << ", cleaning down to "
258 0 : << cycleThreshold << LogIO::POST;
259 :
260 0 : for (model=0;model<numberOfModels();model++) {
261 :
262 0 : Int nx=image(model).shape()(0);
263 0 : Int ny=image(model).shape()(1);
264 0 : Int npol=image(model).shape()(2);
265 0 : Int nchan=image(model).shape()(3);
266 :
267 0 : AlwaysAssert((npol==1)||(npol==2)||(npol==4), AipsError);
268 :
269 0 : if(cycle==1) {
270 0 : iterations[model].resize(nchan);
271 0 : iterations[model]=0;
272 : }
273 :
274 : // Initialize delta image
275 0 : deltaImage(model).set(0.0);
276 :
277 : // Only process solveable models
278 0 : if(isSolveable(model)&&psfmax(model)>0.0) {
279 :
280 0 : if((max(abs(resmax[model]))>cycleThreshold)||
281 0 : (max(abs(resmin[model]))>cycleThreshold)) {
282 :
283 : os << LogIO::NORMAL // Loglevel PROGRESS
284 0 : << "Processing model " << model+1 << LogIO::POST;
285 :
286 0 : IPosition onePlane(4, nx, ny, 1, 1);
287 :
288 0 : IPosition oneCube(4, nx, ny, npol, 1);
289 :
290 : // Loop over all channels. We only want the PSF for the first
291 : // polarization so we iterate over polarization LAST
292 :
293 : // Now clean each channel
294 0 : for (Int chan=0;chan<nchan;++chan) {
295 0 : if(nchan>1) {
296 : os << LogIO::NORMAL // Loglevel PROGRESS
297 : <<"Processing channel number "<<chan<<" of "<<nchan
298 0 : << " channels" <<LogIO::POST;
299 : }
300 0 : if((abs(resmax[model][chan])>cycleThreshold) ||
301 0 : (abs(resmin[model][chan])>cycleThreshold)) {
302 0 : LCBox psfbox(IPosition(4, 0, 0, 0, chan),
303 0 : IPosition(4, nx-1, ny-1, 0, chan),
304 0 : PSF(model).shape());
305 0 : SubLattice<Float> psf_sl (PSF(model), psfbox, false);
306 :
307 0 : LCBox imagebox(IPosition(4, 0, 0, 0, chan),
308 0 : IPosition(4, nx-1, ny-1, npol-1, chan),
309 0 : residual(model).shape());
310 :
311 :
312 0 : SubLattice<Float> residual_sl (residual(model), imagebox, true);
313 0 : SubLattice<Float> image_sl (image(model), imagebox, true);
314 0 : SubLattice<Float> deltaimage_sl (deltaImage(model), imagebox, true);
315 : // Now make a convolution equation for this
316 : // residual image and psf and then clean it
317 : {
318 0 : LatConvEquation eqn(psf_sl, residual_sl);
319 0 : ClarkCleanLatModel cleaner(deltaimage_sl);
320 0 : cleaner.setResidual(residual_sl);
321 0 : cleaner.setGain(gain());
322 0 : cleaner.setNumberIterations(numberIterations());
323 0 : cleaner.setInitialNumberIterations(iterations[model](chan));
324 0 : cleaner.setThreshold(cycleThreshold);
325 0 : cleaner.setPsfPatchSize(IPosition(2,51));
326 0 : cleaner.setMaxNumberMajorCycles(1);
327 0 : cleaner.setMaxNumberMinorIterations(100000);
328 0 : cleaner.setHistLength(1024);
329 0 : cleaner.setCycleFactor(cycleFactor_p);
330 0 : cleaner.setMaxNumPix(32*1024);
331 0 : cleaner.setChoose(false);
332 0 : if(cycleSpeedup_p >1)
333 0 : cleaner.setSpeedup(cycleSpeedup_p);
334 : else
335 0 : cleaner.setSpeedup(0.0);
336 : //Be a bit more conservative with pathologically bad PSFs
337 0 : if(maxSidelobe > 0.5)
338 0 : cleaner.setMaxNumberMinorIterations(5);
339 0 : else if(maxSidelobe > 0.35)
340 0 : cleaner.setMaxNumberMinorIterations(50);
341 :
342 : //cleaner.setSpeedup(0.0);
343 0 : if ( displayProgress_p ) {
344 0 : cleaner.setProgress( *progress_p );
345 : }
346 : os << LogIO::NORMAL // Loglevel PROGRESS
347 0 : << "Starting minor cycle of Clean" << LogIO::POST;
348 0 : SubLattice<Float> mask_sl;
349 0 : if(hasMask(model)) {
350 0 : mask_sl=SubLattice<Float> (mask(model), psfbox, true);
351 0 : cleaner.setMask(mask_sl);
352 : }
353 :
354 0 : modified_p= cleaner.singleSolve(eqn, residual_sl) || modified_p;
355 :
356 :
357 0 : if(modified_p){
358 : os << LogIO::NORMAL // Loglevel PROGRESS (See CAS-2017)
359 : << "Finished minor cycle of Clean"
360 0 : << LogIO::POST;
361 :
362 0 : if (cleaner.numberIterations()>0) {
363 : os << LogIO::NORMAL // Loglevel INFO
364 : << "Clean used " << cleaner.numberIterations()
365 : << " iterations to approach a threshold of "
366 0 : << cycleThreshold << LogIO::POST;
367 : }
368 : }
369 :
370 0 : iterations[model](chan)=cleaner.numberIterations();
371 0 : maxIterations=(iterations[model](chan)>maxIterations) ?
372 0 : iterations[model](chan) : maxIterations;
373 :
374 : os << LogIO::NORMAL2 // Loglevel PROGRESS
375 0 : << "Adding increment to existing model" << LogIO::POST;
376 0 : LatticeExpr<Float> expr=image_sl+deltaimage_sl;
377 0 : image_sl.copyData(expr);
378 : }
379 : }
380 : }// channel loop
381 0 : if(!modified_p){
382 : os << LogIO::WARN
383 : << "No clean component found below threshold of "
384 : << cycleThreshold
385 0 : << "\n in region selected to clean in model " << model << LogIO::POST;
386 :
387 : }
388 :
389 0 : if(maxIterations==0) {
390 0 : stop=true;
391 : }
392 : else{
393 0 : stop=false;
394 : }
395 : os << LogIO::NORMAL // Loglevel INFO
396 : << "Model " << model << " has "
397 0 : << LatticeExprNode(sum(image(model))).getFloat()
398 : << " Jy in clean components."
399 0 : << LogIO::POST;
400 : }
401 : else {
402 : os << LogIO::NORMAL // Loglevel INFO
403 : <<"Skipping model "<<model<<" :peak residual below threshold"
404 0 : <<LogIO::POST;
405 : }
406 : }
407 : }
408 0 : if(maxIterations != oldMaxIterations)
409 0 : oldMaxIterations=maxIterations;
410 : else {
411 : os << LogIO::NORMAL // Loglevel INFO
412 0 : << "No more cleaning occured in this major cycle - stopping now" << LogIO::POST;
413 0 : stop=true;
414 0 : converged=true;
415 : }
416 : }
417 : }
418 0 : if (progress_p) delete progress_p;
419 :
420 :
421 0 : if(modified_p || lastCycleWriteModel) {
422 : os << LogIO::NORMAL // Loglevel INFO
423 0 : << LatticeExprNode(sum(image(0))).getFloat()
424 0 : << " Jy is the sum of the clean components " << LogIO::POST;
425 : os << LogIO::NORMAL // Loglevel PROGRESS
426 0 : << "Finalizing residual images for all fields" << LogIO::POST;
427 0 : makeNewtonRaphsonStep(se, false, true);
428 0 : Float finalabsmax=maxField(resmax, resmin);
429 :
430 : os << LogIO::NORMAL // Loglevel INFO
431 0 : << "Final maximum residual = " << finalabsmax << LogIO::POST;
432 0 : converged=(finalabsmax < threshold());
433 0 : for (model=0;model<numberOfModels();model++) {
434 : os << LogIO::NORMAL // Loglevel INFO
435 : << "Model " << model+1 << ": max, min residuals = "
436 0 : << max(resmax[model]) << ", " << min(resmin[model]) << endl;
437 0 : }
438 : }
439 : else {
440 : os << LogIO::NORMAL // Loglevel INFO
441 0 : << "Residual images for all fields are up-to-date" << LogIO::POST;
442 : }
443 :
444 0 : os << LogIO::POST;
445 :
446 0 : return(converged);
447 : };
448 :
449 :
450 : // Find maximum residual
451 0 : Float CSCleanImageSkyModel::maxField(Block<Vector<Float> >& imagemax,
452 : Block<Vector<Float> >& imagemin) {
453 :
454 0 : LogIO os(LogOrigin("ImageSkyModel","maxField"));
455 :
456 0 : Float absmax=0.0;
457 :
458 : // Loop over all models
459 0 : for (Int model=0;model<numberOfModels();model++) {
460 :
461 0 : imagemax[model].resize(image(model).shape()(3));
462 0 : imagemin[model].resize(image(model).shape()(3));
463 : // Remember that the residual image can be either as specified
464 : // or created specially.
465 0 : ImageInterface<Float>* imagePtr=0;
466 0 : if(residual_p[model]) {
467 0 : imagePtr=residual_p[model];
468 : }
469 : else {
470 0 : imagePtr=(ImageInterface<Float> *)residualImage_p[model];
471 : }
472 0 : AlwaysAssert(imagePtr, AipsError);
473 0 : AlwaysAssert(imagePtr->shape().nelements()==4, AipsError);
474 0 : Int nx=imagePtr->shape()(0);
475 0 : Int ny=imagePtr->shape()(1);
476 0 : Int npol=imagePtr->shape()(2);
477 :
478 0 : AlwaysAssert((npol==1)||(npol==2)||(npol==4), AipsError);
479 :
480 : // Loop over all channels
481 0 : IPosition onePlane(4, nx, ny, 1, 1);
482 0 : LatticeStepper ls(imagePtr->shape(), onePlane, IPosition(4, 0, 1, 2, 3));
483 0 : RO_LatticeIterator<Float> imageli(*imagePtr, ls);
484 :
485 : // If we are using a mask then reset the region to be
486 : // cleaned
487 0 : Array<Float> maskArray;
488 0 : RO_LatticeIterator<Float> maskli;
489 0 : if(hasMask(model)) {
490 0 : Int mx=mask(model).shape()(0);
491 0 : Int my=mask(model).shape()(1);
492 0 : Int mpol=mask(model).shape()(2);
493 : //AlwaysAssert(mx==nx, AipsError);
494 : //AlwaysAssert(my==ny, AipsError);
495 : //AlwaysAssert(mpol==npol, AipsError);
496 0 : if((mx != nx) || (my != ny) || (mpol != npol)){
497 0 : throw(AipsError("Mask image shape is not the same as dirty image"));
498 : }
499 0 : LatticeStepper mls(mask(model).shape(), onePlane,
500 0 : IPosition(4, 0, 1, 2, 3));
501 :
502 0 : maskli=RO_LatticeIterator<Float> (mask(model), mls);
503 0 : maskli.reset();
504 0 : if (maskli.cursor().shape().nelements() > 1) maskArray=maskli.cursor();
505 : }
506 :
507 0 : Int chan=0;
508 0 : Int polpl=0;
509 : Float imax, imin;
510 0 : imax=-1E20; imagemax[model]=imax;
511 0 : imin=+1E20; imagemin[model]=imin;
512 0 : imageli.reset();
513 :
514 0 : for (imageli.reset();!imageli.atEnd();imageli++) {
515 0 : imax=-1E20;
516 0 : imin=+1E20;
517 0 : IPosition imageposmax(imageli.cursor().ndim());
518 0 : IPosition imageposmin(imageli.cursor().ndim());
519 :
520 : // If we are using a mask then multiply by it
521 0 : if (hasMask(model)) {
522 0 : Array<Float> limage=imageli.cursor();
523 : //limage*=maskArray;
524 0 : minMaxMasked(imin, imax, imageposmin, imageposmax, limage, maskArray);
525 0 : maskli++;
526 0 : if (maskli.cursor().shape().nelements() > 1) maskArray=maskli.cursor();
527 :
528 : }
529 :
530 : else {
531 0 : minMax(imin, imax, imageposmin, imageposmax, imageli.cursor());
532 : }
533 0 : if(abs(imax)>absmax) absmax=abs(imax);
534 0 : if(abs(imin)>absmax) absmax=abs(imin);
535 0 : if(imin<imagemin[model][chan]) imagemin[model][chan]=imin;
536 0 : if(imax>imagemax[model][chan]) imagemax[model][chan]=imax;
537 0 : ++polpl;
538 0 : if(polpl==npol){
539 0 : ++chan;
540 0 : polpl=0;
541 : }
542 : }
543 : }
544 0 : return absmax;
545 : };
546 :
547 :
548 0 : Vector<Float> CSCleanImageSkyModel::outerMinMax(Lattice<Float> & lat, const uInt nCenter )
549 : {
550 0 : Array<Float> arr = lat.get();
551 0 : IPosition pos( arr.shape() );
552 0 : uInt nx = lat.shape()(0);
553 0 : uInt ny = lat.shape()(1);
554 0 : uInt innerx = lat.shape()(0)/4;
555 0 : uInt innery = lat.shape()(1)/4;
556 0 : uInt nxc = 0;
557 0 : uInt nyc = 0;
558 0 : Float amax = 0.0;
559 0 : Vector<Float> amax2,minMax(2);
560 : //
561 : // First locate the location of the peak
562 : //
563 0 : for (uInt ix = 0; ix < nx; ix++)
564 0 : for (uInt iy = 0; iy < ny; iy++)
565 0 : if (arr(IPosition(4, ix, iy, 0, 0)) > amax)
566 : {
567 0 : nxc = ix;
568 0 : nyc = iy;
569 0 : amax = arr(IPosition(4, ix, iy, 0, 0));
570 : }
571 : //
572 : // Now exclude the mainlobe center on the location of the peak to
573 : // get the max. outer sidelobe.
574 : //
575 0 : Float myMax=0.0;
576 0 : Float myMin=0.0;
577 :
578 0 : uInt nxL = nxc - nCenter;
579 0 : uInt nxH = nxc + nCenter;
580 0 : uInt nyL = nyc - nCenter;
581 0 : uInt nyH = nyc + nCenter;
582 0 : uInt nx0 = nxc - innerx/2, nx1 = nxc + innerx/2;
583 0 : uInt ny0 = nyc - innery/2, ny1 = nyc + innery/2;
584 :
585 : //
586 : // Search only in the square with innerx and innery pixels on each side.
587 : //
588 0 : for (uInt ix = nx0; ix < nx1; ix++) {
589 0 : for (uInt iy = ny0; iy < ny1; iy++) {
590 0 : if ( !(ix >= nxL && ix <= nxH && iy >= nyL && iy <= nyH) ) {
591 0 : if (arr(IPosition(4, ix, iy, 0, 0)) > myMax)
592 0 : myMax = arr(IPosition(4, ix, iy, 0, 0));
593 0 : if (arr(IPosition(4, ix, iy, 0, 0)) < myMin)
594 0 : myMin = arr(IPosition(4, ix, iy, 0, 0));
595 : }
596 : }
597 : }
598 :
599 : // Float absMax = max( abs(myMin), myMax );
600 : // return absMax;
601 0 : minMax(0) = myMin;
602 0 : minMax(1) = max( abs(myMin), myMax );
603 0 : return minMax;
604 : };
605 :
606 : } //#End casa namespace
|
