Line data Source code
1 : //# ClarkCleanLatModel.cc: this defines ClarkCleanLatModel
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$
27 :
28 : #include <synthesis/MeasurementEquations/ClarkCleanLatModel.h>
29 : #include <synthesis/MeasurementEquations/ClarkCleanProgress.h>
30 : #include <casacore/casa/Arrays/Slice.h>
31 : #include <casacore/lattices/Lattices/LatticeStepper.h>
32 : #include <casacore/lattices/Lattices/LatticeIterator.h>
33 : #include <casacore/casa/Arrays/Matrix.h>
34 : #include <casacore/casa/Arrays/Vector.h>
35 : #include <casacore/casa/Arrays/Array.h>
36 : #include <casacore/casa/Arrays/ArrayError.h>
37 : #include <casacore/casa/Arrays/ArrayMath.h>
38 : #include <casacore/casa/Arrays/VectorIter.h>
39 : #include <casacore/casa/Logging/LogOrigin.h>
40 : #include <casacore/casa/Exceptions/Error.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 : #include <iostream>
43 : #include <casacore/casa/System/Choice.h>
44 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
45 : #include <synthesis/MeasurementEquations/CCList.h>
46 : #include <casacore/lattices/LEL/LatticeExpr.h>
47 : #include <casacore/lattices/LEL/LatticeExprNode.h>
48 : #include <casacore/casa/BasicSL/Constants.h>
49 :
50 : using namespace casacore;
51 : namespace casa { //# NAMESPACE CASA - BEGIN
52 :
53 : //# These are the definitions of the fortran functions
54 :
55 : #define NEED_FORTRAN_UNDERSCORES
56 :
57 : #if defined(NEED_FORTRAN_UNDERSCORES)
58 : #define getbigf getbigf_
59 : #define getbig2f getbig2f_
60 : #define getbig4f getbig4f_
61 : #define abshisf abshisf_
62 : #define abshis2f abshis2f_
63 : #define abshis4f abshis4f_
64 : #define absmaxf absmaxf_
65 : #define absmax2f absmax2f_
66 : #define absmax4f absmax4f_
67 : #define subcomf subcomf_
68 : #define subcom2f subcom2f_
69 : #define subcom4f subcom4f_
70 : #define maxabsf maxabsf_
71 : #define maxabs2f maxabs2f_
72 : #define maxabs4f maxabs4f_
73 : #define maxabmf maxabmf_
74 : #define maxabm2f maxabm2f_
75 : #define maxabm4f maxabm4f_
76 : #define abshimf abshimf_
77 : #define abshim2f abshim2f_
78 : #define abshim4f abshim4f_
79 : #define getbimf getbimf_
80 : #define getbim2f getbim2f_
81 : #define getbim4f getbim4f_
82 : #endif
83 :
84 : extern "C" {
85 : void getbigf(const Float * const pixVal, const Int * const pixPos,
86 : Int * maxPix, const Float * fluxLimit, const Float * arr,
87 : const Int * nx, const Int * ny);
88 : void getbig2f(const Float * const pixVal, const Int * const pixPos,
89 : Int * maxPix, const Float * fluxLimit, const Float * arr,
90 : const Int * nx, const Int * ny);
91 : void getbig4f(const Float * const pixVal, const Int * const pixPos,
92 : Int * maxPix, const Float * fluxLimit, const Float * arr,
93 : const Int * nx, const Int * ny);
94 : void abshisf(Int * hist, Float * minval, Float * maxval,
95 : const Int * nbins, const Float * arr, const Int * npix);
96 : void abshis2f(Int * hist, Float * minval, Float * maxval,
97 : const Int * nbins, const Float * arr, const Int * npix);
98 : void abshis4f(Int * hist, Float * minval, Float * maxval,
99 : const Int * nbins, const Float * arr, const Int * npix);
100 : void absmaxf(Float * maxelem, Float * maxval, Int * maxpos,
101 : const Float * arr, const Int * npix);
102 : void absmax2f(Float * maxelem, Float * maxval, Int * maxpos,
103 : const Float * arr, const Int * npix);
104 : void absmax4f(Float * maxelem, Float * maxval, Int * maxpos,
105 : const Float * arr, const Int * npix);
106 : void subcomf(Float * pixval, const Int * pixpos, const Int * npix,
107 : const Float * maxpix, const Int * maxpos,
108 : const Float * psf, const Int * nx, const Int * ny);
109 : void subcom2f(Float * pixval, const Int * pixpos, const Int * npix,
110 : const Float * maxpix, const Int * maxpos,
111 : const Float * psf, const Int * nx, const Int * ny);
112 : void subcom4f(Float * pixval, const Int * pixpos, const Int * npix,
113 : const Float * maxpix, const Int * maxpos,
114 : const Float * psf, const Int * nx, const Int * ny);
115 : void maxabsf(Float * maxval, const Float * arr, const Int * npix);
116 : void maxabs2f(Float * maxval, const Float * arr, const Int * npix);
117 : void maxabs4f(Float * maxval, const Float * arr, const Int * npix);
118 : void maxabmf(Float * maxval, const Float * arr, const Float * mask,
119 : const Int * npix);
120 : void maxabm2f(Float * maxval, const Float * arr, const Float * mask,
121 : const Int * npix);
122 : void maxabm4f(Float * maxval, const Float * arr, const Float * mask,
123 : const Int * npix);
124 : void abshimf(Int * hist, Float * minval, Float * maxval, const Int * nbins,
125 : const Float * arr, const Float * mask, const Int * npix);
126 : void abshim2f(Int * hist, Float * minval, Float * maxval, const Int * nbins,
127 : const Float * arr, const Float * mask, const Int * npix);
128 : void abshim4f(Int * hist, Float * minval, Float * maxval, const Int * nbins,
129 : const Float * arr, const Float * mask, const Int * npix);
130 : void getbimf(const Float * const pixVal, const Int * const pixPos,
131 : Int * maxPix, const Float * fluxLimit, const Float * arr,
132 : const Float * mask, const Int * nx, const Int * ny);
133 : void getbim2f(const Float * const pixVal, const Int * const pixPos,
134 : Int * maxPix, const Float * fluxLimit, const Float * arr,
135 : const Float * mask, const Int * nx, const Int * ny);
136 : void getbim4f(const Float * const pixVal, const Int * const pixPos,
137 : Int * maxPix, const Float * fluxLimit, const Float * arr,
138 : const Float * mask, const Int * nx, const Int * ny);
139 : };
140 :
141 : //----------------------------------------------------------------------
142 0 : ClarkCleanLatModel::ClarkCleanLatModel()
143 0 : :itsModelPtr(0),
144 0 : itsResidualPtr(0),
145 0 : itsSoftMaskPtr(0),
146 0 : itsMaxNumPix(32*1024),
147 0 : itsHistBins(1024),
148 0 : itsMaxNumberMinorIterations(10000),
149 0 : itsInitialNumberIterations(0),
150 0 : itsMaxNumberMajorCycles(-1),
151 0 : itsMaxExtPsf(0.0),
152 0 : itsPsfPatchSize(2,51,51),
153 0 : itsChoose(false),
154 0 : itsSpeedup(0.0),
155 0 : itsCycleFactor(1.5),
156 0 : itsLog(LogOrigin("ClarkCleanLatModel", "ClarkCleanLatModel()")),
157 0 : itsProgressPtr(0),
158 0 : itsJustStarting(true),
159 0 : itsWarnFlag(false),
160 0 : itsLocalResTL(false)
161 :
162 : {
163 0 : };
164 : //----------------------------------------------------------------------
165 0 : ClarkCleanLatModel::ClarkCleanLatModel(Lattice<Float> & model)
166 0 : :itsModelPtr(&model),
167 0 : itsResidualPtr(0),
168 0 : itsSoftMaskPtr(0),
169 0 : itsMaxNumPix(32*1024),
170 0 : itsHistBins(1024),
171 0 : itsMaxNumberMinorIterations(10000),
172 0 : itsInitialNumberIterations(0),
173 0 : itsMaxNumberMajorCycles(-1),
174 0 : itsMaxExtPsf(0.0),
175 0 : itsPsfPatchSize(2, 51, 51),
176 0 : itsChoose(false),
177 0 : itsSpeedup(0.0),
178 0 : itsCycleFactor(1.5),
179 0 : itsLog(LogOrigin("ClarkCleanLatModel",
180 : "ClarkCleanLatModel(const Lattice<Float> & model)")),
181 0 : itsProgressPtr(0),
182 0 : itsWarnFlag(false),
183 0 : itsLocalResTL(false)
184 : {
185 0 : AlwaysAssert(getModel().ndim() >= 2, AipsError);
186 0 : if (getModel().ndim() >= 3)
187 0 : AlwaysAssert(getModel().shape()(2) == 1 || getModel().shape()(2) == 2 ||
188 : getModel().shape()(2) == 4, AipsError);
189 0 : if (getModel().ndim() >= 4)
190 0 : for (uInt i = 3; i < getModel().ndim(); i++)
191 0 : AlwaysAssert(getModel().shape()(i) == 1, AipsError);
192 : // itsLog << LogOrigin("ClarkCleanLatModel", "ClarkCleanLatModel")
193 : // << "Model shape is:" << getModel().shape() << endl;
194 0 : };
195 : //----------------------------------------------------------------------
196 0 : ClarkCleanLatModel::ClarkCleanLatModel(Lattice<Float> & model,
197 0 : Lattice<Float> & mask)
198 0 : :itsModelPtr(&model),
199 0 : itsResidualPtr(0),
200 0 : itsSoftMaskPtr(&mask),
201 0 : itsMaxNumPix(32*1024),
202 0 : itsHistBins(1024),
203 0 : itsMaxNumberMinorIterations(10000),
204 0 : itsInitialNumberIterations(0),
205 0 : itsMaxNumberMajorCycles(-1),
206 0 : itsMaxExtPsf(0.0),
207 0 : itsPsfPatchSize(2, 51, 51),
208 0 : itsChoose(false),
209 0 : itsSpeedup(0.0),
210 0 : itsCycleFactor(1.5),
211 0 : itsLog(LogOrigin("ClarkCleanLatModel",
212 : "ClarkCleanLatModel(Lattice<Float> & model"
213 : ", Lattice<Float> & mask)")),
214 0 : itsProgressPtr(0),
215 0 : itsJustStarting(true),
216 0 : itsWarnFlag(false),
217 0 : itsLocalResTL(false)
218 :
219 : {
220 0 : AlwaysAssert(getModel().ndim() >= 2, AipsError);
221 0 : if (getModel().ndim() >= 3)
222 0 : AlwaysAssert(getModel().shape()(2) == 1 || getModel().shape()(2) == 2 ||
223 : getModel().shape()(2) == 4, AipsError);
224 0 : if (getModel().ndim() >= 4)
225 0 : for (uInt i = 3; i < getModel().ndim(); i++)
226 0 : AlwaysAssert(getModel().shape()(i) == 1, AipsError);
227 :
228 0 : AlwaysAssert(itsSoftMaskPtr->ndim() >= 2, AipsError);
229 0 : AlwaysAssert(itsSoftMaskPtr->shape()(0) == getModel().shape()(0), AipsError);
230 0 : AlwaysAssert(itsSoftMaskPtr->shape()(1) == getModel().shape()(1), AipsError);
231 0 : if (itsSoftMaskPtr->ndim() >= 3)
232 0 : for (uInt i = 2; i < itsSoftMaskPtr->ndim(); i++)
233 0 : AlwaysAssert(itsSoftMaskPtr->shape()(i) == 1, AipsError);
234 0 : };
235 0 : ClarkCleanLatModel::ClarkCleanLatModel(Lattice<Float> & model,
236 : Lattice<Float> & residual,
237 0 : Lattice<Float> & mask)
238 0 : :itsModelPtr(&model),
239 0 : itsResidualPtr(&residual),
240 0 : itsSoftMaskPtr(&mask),
241 0 : itsMaxNumPix(32*1024),
242 0 : itsHistBins(1024),
243 0 : itsMaxNumberMinorIterations(10000),
244 0 : itsInitialNumberIterations(0),
245 0 : itsMaxNumberMajorCycles(-1),
246 0 : itsMaxExtPsf(0.0),
247 0 : itsPsfPatchSize(2, 51, 51),
248 0 : itsChoose(false),
249 0 : itsSpeedup(0.0),
250 0 : itsCycleFactor(1.5),
251 0 : itsLog(LogOrigin("ClarkCleanLatModel",
252 : "ClarkCleanLatModel(Lattice<Float> & model"
253 : ", Lattice<Float> & mask)")),
254 0 : itsProgressPtr(0),
255 0 : itsJustStarting(true),
256 0 : itsWarnFlag(false),
257 0 : itsLocalResTL(false)
258 :
259 : {
260 :
261 :
262 0 : };
263 :
264 :
265 0 : ClarkCleanLatModel::~ClarkCleanLatModel() {
266 : // CAS-9268 needs this for tclean, where itsResidualPtr is not constructed outside the cleaner
267 0 : if(itsLocalResTL==True && itsResidualPtr != NULL){delete itsResidualPtr; itsResidualPtr=NULL;}
268 0 : }
269 :
270 :
271 0 : void ClarkCleanLatModel::setModel(const Lattice<Float>& model){
272 0 : AlwaysAssert(model.ndim() >= 2, AipsError);
273 0 : if (model.ndim() >= 3)
274 0 : AlwaysAssert(model.shape()(2) == 1 || model.shape()(2) == 2 ||
275 : model.shape()(2) == 4, AipsError);
276 0 : if (model.ndim() >= 4)
277 0 : for (uInt i = 3; i < model.ndim(); i++)
278 0 : AlwaysAssert(model.shape()(i) == 1, AipsError);
279 : // itsModelPtr = &model;
280 0 : throw(AipsError(
281 0 : "setModel(const Lattice<Float>& ) : only non-const version works!"));
282 : };
283 :
284 0 : void ClarkCleanLatModel::setModel(Lattice<Float>& model){
285 0 : AlwaysAssert(model.ndim() >= 2, AipsError);
286 0 : if (model.ndim() >= 3)
287 0 : AlwaysAssert(model.shape()(2) == 1 || model.shape()(2) == 2 ||
288 : model.shape()(2) == 4, AipsError);
289 0 : if (model.ndim() >= 4)
290 0 : for (uInt i = 3; i < model.ndim(); i++)
291 0 : AlwaysAssert(model.shape()(i) == 1, AipsError);
292 0 : itsModelPtr = &model;
293 0 : };
294 :
295 0 : const Lattice<Float> & ClarkCleanLatModel::getMask() const{
296 0 : return (*itsSoftMaskPtr);
297 : };
298 :
299 0 : void ClarkCleanLatModel::setMask(const Lattice<Float> & mask){
300 0 : AlwaysAssert(mask.ndim() >= 2, AipsError);
301 0 : AlwaysAssert(mask.shape()(0) == getModel().shape()(0), AipsError);
302 0 : AlwaysAssert(mask.shape()(1) == getModel().shape()(1), AipsError);
303 0 : if (mask.ndim() >= 3)
304 0 : for (uInt i = 2; i < mask.ndim(); i++)
305 0 : AlwaysAssert(mask.shape()(i) == 1, AipsError);
306 0 : itsSoftMaskPtr = &mask;
307 0 : };
308 :
309 :
310 : //----------------------------------------------------------------------
311 0 : Bool ClarkCleanLatModel::solve(LatConvEquation & eqn){
312 0 : itsLog << LogOrigin("ClarkCleanLatModel", "solve");
313 0 : AlwaysAssert(getModel().ndim() >= 2, AipsError);
314 0 : const IPosition dataShape = getModel().shape();
315 :
316 0 : Int npol = 1;
317 0 : if (getModel().ndim() >= 3)
318 0 : npol = dataShape(2);
319 0 : AlwaysAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
320 :
321 : // Determine the number of polarisations
322 : // itsLog << "Data has " << npol << " polarisations" << LogIO::POST;
323 :
324 : // compute the current residual image (using an FFT)
325 0 : if(itsResidualPtr==0){
326 0 : itsResidualPtr=new TempLattice<Float> (dataShape);
327 0 : itsLocalResTL=True;
328 0 : eqn.residual(*itsResidualPtr, *this);
329 : }
330 : // Determine the psf patch to use
331 0 : Matrix<Float> psfPatch;
332 : Float maxExtPsf;
333 :
334 0 : maxExtPsf = getPsfPatch(psfPatch, eqn);
335 : // itsLog << "PsfPatch shape is: " << psfPatch.shape()
336 : // << " and has a maximum exterior sidelobe of "
337 : // << maxExtPsf << LogIO::POST;
338 :
339 : // Declare variables needed inside the following while loop
340 : Int numPix; // the number of Active pixels
341 0 : Int maxNumPix = 0; // unclear
342 0 : uInt numIterations = itsInitialNumberIterations;
343 : // Number of Iterations done so far
344 : uInt numMinorIterations; // The number of minor iterations done/todo
345 0 : uInt numMajorCycles = 0; // The number of major cycles done
346 0 : uInt maxNumberMinorIterations = 0;// The max. number of min. iterations
347 : // ever used
348 0 : Float Fmn= 0; //1; // The "uncertainty" factor
349 : Float fluxLimit; // The fluxlimit for the current major cycle
350 0 : Float totalFlux = 0;
351 0 : Float factor=itsCycleFactor/4.5; // This factor is empirical...found to be
352 : // good without being too conservative
353 : // Note that a Matrix is used rather than say a Vector of IPositions as
354 : // it allows the inner loop (in doMinorIterations()) to be more highly
355 : // optimised (using pointers)
356 :
357 : // find its maximum value of the residual
358 0 : Float maxRes = maxResidual(*itsResidualPtr);
359 0 : Float maxResPrevious=maxRes;
360 0 : if (numIterations > 0)
361 0 : itsLog << LogIO::NORMAL << "Initial maximum residual: " << maxRes << LogIO::POST;
362 : // if flux limit or iteration limit reached then bail out.
363 0 : Bool userHalt = false;
364 0 : Int numIt=numberIterations();
365 : //Pathological PSFs
366 0 : if(maxExtPsf > 0.5)
367 0 : itsMaxNumberMinorIterations=5;
368 0 : else if(maxExtPsf > 0.35)
369 0 : itsMaxNumberMinorIterations=50;
370 :
371 0 : while ((Int(numIterations) < numIt) &&
372 0 : (maxRes > threshold()) &&
373 0 : ((itsMaxNumberMajorCycles<0)||
374 0 : (Int(numMajorCycles)<itsMaxNumberMajorCycles)) &&
375 0 : userHalt == false){
376 :
377 0 : CCList activePixels(npol, 2, 0) ; // cache of active pixel values and positions;
378 :
379 : // determine the fluxlimit for this major cycle
380 : // choose fluxlimit for residuals to be considered in minor iterations
381 : // don't consider residuals below maxRes*(value of maxPsf at outside the
382 : // patch)
383 0 : fluxLimit = maxRes * maxExtPsf * factor;
384 :
385 0 : if(factor > 1.0)
386 0 : fluxLimit = min(0.95*maxRes, fluxLimit);
387 :
388 :
389 : // itsLog << "Fluxlimit determined using the Maximum exterior Psf: "
390 : // << fluxLimit << LogIO::POST;
391 : // See if this flux limit needs to be modified because it selects too
392 : // many pixels.
393 :
394 : //*** COMMENTED below off as its giving extremely conservative limit
395 : // and making the loop very slow
396 : // minLimit = biggestResiduals(maxRes, itsMaxNumPix, fluxLimit, residual);
397 :
398 :
399 : // itsLog << "Fluxlimit determined using the maximum number active pixels: "
400 : // << minLimit << endl;
401 : //
402 : // fluxLimit = std::max(fluxLimit, minLimit);
403 :
404 :
405 : // itsLog << "Final Fluxlimit: " << fluxLimit << LogIO::POST;
406 :
407 : // Copy all the active pixels into separate areas for better memory
408 : // management and quicker indexing.
409 :
410 0 : cacheActivePixels(activePixels, *itsResidualPtr,
411 0 : std::max(fluxLimit,threshold()));
412 : // The numpix calculated here frequently differs
413 : // from the number calculated using the histogram, because of the
414 : // quantisation of the flux levels in the histogram, and the imposition
415 : // of an external fluxlevel.
416 0 : numPix = activePixels.nComp();
417 0 : if (numPix > 0) {
418 : // itsLog <<"Major cycle has "<< numPix << " active residuals, "
419 : // << "a Fluxlimit of " << std::max(fluxLimit,threshold()) << endl;
420 : // Start of minor cycles
421 :
422 0 : numMinorIterations = std::min(itsMaxNumberMinorIterations,
423 0 : numIt-numIterations);
424 0 : doMinorIterations(activePixels,
425 : psfPatch, fluxLimit, numMinorIterations,
426 : Fmn, numIterations, totalFlux);
427 0 : numIterations += numMinorIterations;
428 0 : maxNumberMinorIterations = std::max(maxNumberMinorIterations,
429 : numMinorIterations);
430 0 : maxNumPix = std::max((Int)itsMaxNumPix, numPix);
431 : // Now do a major cycle
432 0 : eqn.residual(*itsResidualPtr, *this);
433 :
434 :
435 : // find the new maximum residual
436 : // maxRes = maxResidual(*itsResidualPtr);
437 0 : maxRes = itsMaxRes;
438 0 : if(maxRes > maxResPrevious){
439 0 : if(!itsWarnFlag){
440 0 : itsLog << LogIO::WARN
441 : << "Slowing down in the minor cycle loop; "
442 : << "Could be a PSF with bad sidelobes; "
443 0 : << "Suggest using CS or Hogbom instead" << LogIO::POST;
444 : }
445 0 : factor=factor*3; //pathological PSF's go very slowly in minorcycles
446 0 : itsMaxNumberMinorIterations=10;
447 0 : itsWarnFlag=true;
448 : }
449 0 : maxResPrevious=maxRes;
450 :
451 0 : itsLog << LogIO::NORMAL
452 : << "Iteration: " << numIterations
453 : << ", Maximum residual=" << maxRes
454 0 : << " Flux limit=" << std::max(fluxLimit,threshold())
455 0 : << ", " << numPix << " Active pixels" << LogIO::POST;
456 :
457 : // Count the number of major cycles
458 0 : numMajorCycles++;
459 : }
460 : else{
461 0 : itsLog << LogIO::WARN
462 : << "Zero Pixels selected with a Fluxlimit of " << fluxLimit
463 0 : << " and a maximum Residual of " << maxRes << LogIO::POST;
464 0 : if(itsWarnFlag){
465 0 : userHalt=true;
466 0 : itsLog << LogIO::WARN
467 : << "Bailing out prior to reaching threshold as residual value is not converging "
468 0 : << LogIO::POST;
469 : }
470 : else{
471 : //lets try to increase the depth a little bit
472 0 : factor=factor*1.2;
473 : }
474 0 : itsWarnFlag=true;
475 : // userHalt = stopnow();
476 : // The above is commented off as users do not seem to find this
477 : // useful. If nobody ask for it again the function stopnow()
478 : // along with userHalt will be obliterated before the release
479 : }
480 :
481 0 : }
482 :
483 : // Is this a problem?
484 0 : if (itsProgressPtr) {
485 0 : itsProgressPtr->finalize();
486 : }
487 :
488 0 : setThreshold(maxRes);
489 0 : setNumberIterations(numIterations);
490 0 : itsMaxNumPix = maxNumPix;
491 0 : itsMaxNumberMinorIterations = maxNumberMinorIterations;
492 0 : return true;
493 0 : };
494 : //----------------------------------------------------------------------
495 : //----------------------------------------------------------------------
496 0 : Bool ClarkCleanLatModel::singleSolve(LatConvEquation & eqn, Lattice<Float>& residual){
497 0 : itsLog << LogOrigin("ClarkCleanLatModel", "singleSolve");
498 0 : AlwaysAssert(getModel().ndim() >= 2, AipsError);
499 0 : const IPosition dataShape = getModel().shape();
500 : // itsLog << "Model shape " << getModel().shape() << LogIO::POST;
501 : // itsLog << "Residual shape " << residual.shape() << LogIO::POST;
502 0 : Int npol = 1;
503 0 : if (getModel().ndim() >= 3)
504 0 : npol = dataShape(2);
505 0 : AlwaysAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
506 :
507 : // Determine the number of polarisations
508 : // itsLog << "Data has " << npol << " polarisations" << LogIO::POST;
509 :
510 : // Determine the psf patch to use
511 0 : Matrix<Float> psfPatch;
512 : Float maxExtPsf;
513 :
514 0 : maxExtPsf = getPsfPatch(psfPatch, eqn);
515 : // itsLog << "PsfPatch shape is: " << psfPatch.shape()
516 : // << " and has a maximum exterior sidelobe of "
517 : // << maxExtPsf << LogIO::POST;
518 :
519 : // Declare variables needed inside the following while loop
520 : Int numPix; // the number of Active pixels
521 0 : Int maxNumPix = 0; // unclear
522 0 : uInt numIterations = itsInitialNumberIterations;
523 : // Number of Iterations done so far
524 : uInt numMinorIterations; // The number of minor iterations done/todo
525 0 : uInt maxNumberMinorIterations = 0;// The max. number of min. iterations
526 : // ever used
527 0 : Float Fmn= 0; //1; // The "uncertainty" factor
528 : Float fluxLimit; // The fluxlimit for the current major cycle
529 0 : Float totalFlux = 0;
530 0 : Float factor=1.0/3.0; // This factor is empirical...found to be
531 : // good without being too conservative
532 : // Note that a Matrix is used rather than say a Vector of IPositions as
533 : // it allows the inner loop (in doMinorIterations()) to be more highly
534 : // optimised (using pointers)
535 :
536 : // find its maximum value of the residual
537 0 : Float maxRes = maxResidual(residual);
538 0 : itsLog << "Initial maximum residual: " << maxRes << LogIO::POST;
539 :
540 0 : CCList activePixels(npol, 2, 0) ; // cache of active pixel values and positions;
541 :
542 : // determine the fluxlimit for this major cycle
543 : // choose fluxlimit for residuals to be considered in minor iterations
544 : // don't consider residuals below maxRes*(value of maxPsf at outside the
545 : // patch)
546 0 : fluxLimit = maxRes * maxExtPsf * factor;
547 :
548 0 : if(factor > 1.0)
549 0 : fluxLimit = min(0.95*maxRes, fluxLimit);
550 :
551 :
552 : // itsLog << "Fluxlimit determined using the Maximum exterior Psf: "
553 : // << fluxLimit << LogIO::POST;
554 : // See if this flux limit needs to be modified because it selects too
555 : // many pixels.
556 :
557 : //*** COMMENTED below off as its giving extremely conservative limit
558 : // and making the loop very slow
559 : // minLimit = biggestResiduals(maxRes, itsMaxNumPix, fluxLimit, residual);
560 :
561 :
562 : // itsLog << "Fluxlimit determined using the maximum number active pixels: "
563 : // << minLimit << endl;
564 : //
565 : // fluxLimit = std::max(fluxLimit, minLimit);
566 :
567 :
568 : // itsLog << "Final Fluxlimit: " << fluxLimit << LogIO::POST;
569 :
570 : // Copy all the active pixels into separate areas for better memory
571 : // management and quicker indexing.
572 :
573 0 : cacheActivePixels(activePixels, residual, std::max(fluxLimit,threshold()));
574 : // The numpix calculated here frequently differs
575 : // from the number calculated using the histogram, because of the
576 : // quantisation of the flux levels in the histogram, and the imposition
577 : // of an external fluxlevel.
578 0 : numPix = activePixels.nComp();
579 0 : if (numPix > 0) {
580 : // itsLog <<"Major cycle has "<< numPix << " active residuals, "
581 : // << "a Fluxlimit of " << std::max(fluxLimit,threshold()) << endl;
582 : // Start of minor cycles
583 :
584 :
585 0 : numMinorIterations = std::min(itsMaxNumberMinorIterations,
586 0 : numberIterations()-numIterations);
587 0 : doMinorIterations(activePixels,
588 : psfPatch, fluxLimit, numMinorIterations,
589 : Fmn, numIterations, totalFlux);
590 0 : numIterations += numMinorIterations;
591 : // itsLog << "Clean has used " << numIterations << " Iterations" ;
592 0 : maxNumberMinorIterations = std::max(maxNumberMinorIterations,
593 : numMinorIterations);
594 0 : maxNumPix = std::max((Int)itsMaxNumPix, numPix);
595 :
596 : }
597 : else {
598 0 : itsLog << LogIO::WARN
599 : << "Zero Pixels selected with a Fluxlimit of " << fluxLimit
600 0 : << " and a maximum Residual of " << maxRes << endl;
601 0 : return false;
602 : }
603 :
604 0 : setThreshold(maxRes);
605 0 : setNumberIterations(numIterations);
606 0 : itsMaxNumPix = maxNumPix;
607 0 : itsMaxNumberMinorIterations = maxNumberMinorIterations;
608 0 : return true;
609 0 : };
610 : //----------------------------------------------------------------------
611 0 : void ClarkCleanLatModel::setResidual(Lattice<Float>& residual){
612 :
613 0 : itsResidualPtr= &residual;
614 :
615 0 : };
616 : //---------------------------------------------------------------------
617 0 : void ClarkCleanLatModel::setPsfPatchSize(const IPosition & psfPatchSize){
618 0 : if (psfPatchSize.nelements() < 2) {
619 0 : throw(AipsError("ClarkCleanLatModel::setPsfPatchSize - assumption of 2-D (or greater) not met"));
620 : } else {
621 : // only 2 dimensions used here
622 0 : itsPsfPatchSize(0)=psfPatchSize(0);
623 0 : itsPsfPatchSize(1)=psfPatchSize(1);
624 : }
625 0 : };
626 : //----------------------------------------------------------------------
627 0 : IPosition ClarkCleanLatModel::getPsfPatchSize(){
628 0 : return itsPsfPatchSize;
629 : };
630 : //----------------------------------------------------------------------
631 0 : void ClarkCleanLatModel::setHistLength(const uInt HistBins ){
632 0 : itsHistBins=HistBins;
633 0 : };
634 : //----------------------------------------------------------------------
635 0 : uInt ClarkCleanLatModel::getHistLength(){
636 0 : return itsHistBins;
637 : };
638 : //----------------------------------------------------------------------
639 0 : void ClarkCleanLatModel::setMaxNumberMinorIterations(const uInt maxNumMinorIterations){
640 0 : itsMaxNumberMinorIterations=maxNumMinorIterations;
641 0 : };
642 : //----------------------------------------------------------------------
643 0 : uInt ClarkCleanLatModel::getMaxNumberMinorIterations(){
644 0 : return itsMaxNumberMinorIterations;
645 : };
646 : //----------------------------------------------------------------------
647 0 : void ClarkCleanLatModel::setInitialNumberIterations(const uInt initialNumberIterations){
648 0 : itsInitialNumberIterations=initialNumberIterations;
649 0 : };
650 : //----------------------------------------------------------------------
651 0 : uInt ClarkCleanLatModel::getInitialNumberIterations(){
652 0 : return itsInitialNumberIterations;
653 : };
654 : //----------------------------------------------------------------------
655 0 : void ClarkCleanLatModel::setMaxNumberMajorCycles(const uInt maxNumMajorCycles){
656 0 : itsMaxNumberMajorCycles=maxNumMajorCycles;
657 0 : };
658 : //----------------------------------------------------------------------
659 0 : uInt ClarkCleanLatModel::getMaxNumberMajorCycles(){
660 0 : return itsMaxNumberMajorCycles;
661 : };
662 : //----------------------------------------------------------------------
663 0 : void ClarkCleanLatModel::setMaxNumPix(const uInt maxNumPix ){
664 0 : itsMaxNumPix=maxNumPix;
665 0 : };
666 : //----------------------------------------------------------------------
667 0 : uInt ClarkCleanLatModel::getMaxNumPix(){
668 0 : return itsMaxNumPix;
669 : };
670 : //----------------------------------------------------------------------
671 0 : void ClarkCleanLatModel::setMaxExtPsf(const Float maxExtPsf ){
672 0 : itsMaxExtPsf=maxExtPsf;
673 0 : };
674 : //----------------------------------------------------------------------
675 0 : Float ClarkCleanLatModel::getMaxExtPsf(){
676 0 : return itsMaxExtPsf;
677 : };
678 : //----------------------------------------------------------------------
679 0 : void ClarkCleanLatModel::setSpeedup(const Float speedup ){
680 0 : itsSpeedup=speedup;
681 0 : };
682 : //----------------------------------------------------------------------
683 0 : Float ClarkCleanLatModel::getSpeedup(){
684 0 : return itsSpeedup;
685 : };
686 0 : void ClarkCleanLatModel::setCycleFactor(const Float factor){
687 0 : itsCycleFactor=factor;
688 0 : }
689 :
690 : //----------------------------------------------------------------------
691 0 : void ClarkCleanLatModel::setChoose(const Bool choose ){
692 0 : itsChoose=choose;
693 0 : };
694 : //----------------------------------------------------------------------
695 0 : Bool ClarkCleanLatModel::getChoose(){
696 0 : return itsChoose;
697 : };
698 :
699 : //----------------------------------------------------------------------
700 0 : void ClarkCleanLatModel::
701 : doMinorIterations(CCList & activePixels, Matrix<Float> & psfPatch,
702 : Float fluxLimit, uInt & numberIterations, Float Fmn,
703 : const uInt totalIterations, Float& totalFlux) {
704 0 : const uInt ndim = itsModelPtr->ndim();
705 0 : DebugAssert(ndim >= 2, AipsError);
706 0 : const IPosition modelShape = itsModelPtr->shape();
707 0 : DebugAssert(modelShape(0) > 0, AipsError);
708 0 : DebugAssert(modelShape(1) > 0, AipsError);
709 0 : uInt npol = 1;
710 0 : if (ndim > 2) npol = modelShape(2);
711 0 : DebugAssert(activePixels.nPol() == npol , AipsError);
712 0 : DebugAssert(activePixels.nComp() > 0 , AipsError);
713 0 : DebugAssert( modelShape(0)*modelShape(1)*npol == uInt(modelShape.product()), AipsError);
714 0 : DebugAssert(psfPatch.nrow() > 0, AipsError);
715 0 : DebugAssert(psfPatch.ncolumn() > 0, AipsError);
716 :
717 : // itsLog << LogOrigin("ClarkCleanLatModel", "doMinorIterations");
718 0 : CCList cleanComponents(npol, 2, totalIterations);
719 :
720 : // Find the largest residual and its position.
721 0 : Block<Float> maxRes(npol);
722 : Float absRes;
723 : Float signedAbsRes;
724 : Int offRes;
725 0 : maxVect(maxRes, absRes, offRes, activePixels);
726 0 : Int * maxPosPtr = activePixels.pixelPosition(offRes);
727 0 : Block<Int> maxPos(2, maxPosPtr, false);
728 : // declare variables used inside the main loop
729 0 : Int curIter = 0;
730 0 : Float iterFluxLimit = std::max(fluxLimit, threshold());
731 0 : Float Fac = pow(fluxLimit/absRes, itsSpeedup);
732 : // itsLog << "Initial maximum residual:" << maxRes
733 : // << " (" << absRes << ") "
734 : // << " @ " << maxPos << endl;
735 :
736 : // Do the minor Iterations
737 0 : while ((curIter < Int(numberIterations)) && (absRes > iterFluxLimit)){
738 : // scale the maximum residual by the gain factor
739 0 : for (uInt i = 0; i < npol; i++) maxRes[i] *= gain();
740 0 : totalFlux += maxRes[0];
741 : // Add the new component to the clean components list
742 0 : cleanComponents.addComp(maxRes, maxPos);
743 : // Subtract the component from the current list of active pixels
744 0 : subtractComponent(activePixels, maxRes, maxPos, psfPatch);
745 : // We have now done an iteration
746 0 : curIter++;
747 : // find the next residual
748 0 : maxVect(maxRes, absRes, offRes, activePixels);
749 0 : maxPosPtr = activePixels.pixelPosition(offRes);
750 0 : maxPos.replaceStorage(2, maxPosPtr, false);
751 : // Update the uncertainty factors and fluxlimits
752 0 : Fmn += Fac/Float(totalIterations+curIter);
753 0 : iterFluxLimit = std::max(fluxLimit * Fmn, threshold());
754 :
755 0 : if (itsProgressPtr) {
756 0 : signedAbsRes = absRes * maxRes[0]/abs( maxRes[0] );
757 0 : itsProgressPtr->
758 0 : info(false, (Int)(totalIterations+curIter), (Int)numberIterations,
759 0 : signedAbsRes, IPosition(2,maxPos[0],maxPos[1]),
760 0 : totalFlux, false, itsJustStarting);
761 0 : itsJustStarting = false;
762 : }
763 : }
764 :
765 0 : itsMaxRes= absRes;
766 : // Now copy the clean components into the image.
767 0 : updateModel(cleanComponents);
768 : // Data returned to the main routine
769 0 : numberIterations = curIter;
770 0 : fluxLimit = absRes;
771 0 : }
772 : //----------------------------------------------------------------------
773 : // Find all the pixels in the residual that are greater than fluxlimit, and
774 : // store the values in the activePixels list. Increases the size of the list if
775 : // necessary, but does not decrease it. This function will weight the residual
776 : // using the mask (if set).
777 0 : void ClarkCleanLatModel::
778 : cacheActivePixels(CCList & activePixels,
779 : const Lattice<Float> & residual, Float fluxLimit) {
780 0 : const uInt ndim = residual.ndim();
781 0 : DebugAssert(ndim >= 2, AipsError);
782 0 : const IPosition residualShape = residual.shape();
783 0 : DebugAssert(residualShape == itsModelPtr->shape(), AipsError);
784 0 : if (itsSoftMaskPtr != 0) {
785 0 : DebugAssert(residualShape(0) == itsSoftMaskPtr->shape()(0), AipsError);
786 0 : DebugAssert(residualShape(1) == itsSoftMaskPtr->shape()(1), AipsError);
787 : }
788 0 : DebugAssert(residualShape(0) > 0, AipsError);
789 0 : DebugAssert(residualShape(1) > 0, AipsError);
790 0 : Int npol = 1;
791 0 : if (ndim > 2) npol = residualShape(2);
792 0 : DebugAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
793 0 : DebugAssert(npol == Int(activePixels.nPol()), AipsError);
794 0 : DebugAssert(residualShape(0)*residualShape(1)*npol == residualShape.product(), AipsError);
795 :
796 0 : IPosition cursorShape(ndim,1);
797 : Int tx, ty;
798 : {
799 0 : IPosition tileShape(residual.niceCursorShape());
800 0 : tx = cursorShape(0) = tileShape(0);
801 0 : ty = cursorShape(1) = tileShape(1);
802 0 : if (ndim > 2) cursorShape(2) = npol;
803 0 : }
804 :
805 0 : RO_LatticeIterator<Float> residualIter(residual, cursorShape);
806 0 : RO_LatticeIterator<Float> maskIter(residualIter); // There is no default ctr
807 0 : if (itsSoftMaskPtr != 0) {
808 0 : IPosition mCursorShape = cursorShape;
809 0 : mCursorShape(2) = 1;
810 0 : LatticeStepper mNav(itsSoftMaskPtr->shape(), mCursorShape, LatticeStepper::RESIZE);
811 0 : maskIter = RO_LatticeIterator<Float>(*itsSoftMaskPtr, mNav);
812 0 : maskIter.reset();
813 0 : }
814 :
815 0 : for (residualIter.reset(); !residualIter.atEnd(); residualIter++) {
816 : Bool residualCopy, maskCopy;
817 0 : const Float * residualPtr = residualIter.cursor().getStorage(residualCopy);
818 0 : const Float * maskPtr = 0;
819 0 : if (itsSoftMaskPtr != 0) {
820 0 : maskPtr = maskIter.cursor().getStorage(maskCopy);
821 : }
822 :
823 0 : Int nUsedPix = getbig(activePixels.freeFluxPtr(),
824 0 : activePixels.freePositionPtr(),
825 0 : activePixels.freeComp(), fluxLimit,
826 : residualPtr, maskPtr, npol, tx, ty);
827 :
828 0 : uInt lastLen=activePixels.nComp();
829 :
830 0 : const uInt reqLen = nUsedPix + activePixels.nComp();
831 :
832 :
833 0 : if (reqLen > activePixels.maxComp()) {
834 : // Need to resize the Pixel lists
835 0 : activePixels.resize(reqLen);
836 : // Now rescan the residual to get all the big pixels (in this iter.)
837 0 : nUsedPix = getbig(activePixels.freeFluxPtr(),
838 0 : activePixels.freePositionPtr(),
839 0 : activePixels.freeComp(), fluxLimit,
840 : residualPtr, maskPtr, npol, tx, ty);
841 :
842 0 : activePixels.nComp() = reqLen;
843 :
844 0 : AlwaysAssert(nUsedPix + activePixels.nComp() == activePixels.maxComp(),
845 : AipsError);
846 : }
847 :
848 : // Need to add the corner of the tile for the last set of components found!!!!
849 0 : Int offx=residualIter.position()(0);
850 0 : Int offy=residualIter.position()(1);
851 : // itsLog << "Adding offset " << offx << " " << offy << LogIO::POST;
852 0 : for (uInt c=lastLen;c<activePixels.nComp();c++) {
853 0 : Int* locPtr=activePixels.pixelPosition(c);
854 0 : *locPtr+=offx;
855 0 : *(locPtr+1)+=offy;
856 : }
857 :
858 0 : residualIter.cursor().freeStorage(residualPtr, residualCopy);
859 0 : if (itsSoftMaskPtr != 0) {
860 0 : maskIter.cursor().freeStorage(maskPtr, maskCopy);
861 0 : maskIter++;
862 : }
863 : }
864 : // AlwaysAssert(activePixels.nComp() > 0, AipsError);
865 0 : }
866 : //----------------------------------------------------------------------
867 : // make histogram of absolute values in array
868 0 : void ClarkCleanLatModel::absHistogram(Vector<Int> & hist,
869 : Float & minVal,
870 : Float & maxVal,
871 : const Lattice<Float> & residual) {
872 0 : const IPosition residualShape = residual.shape();
873 0 : const uInt ndim = residualShape.nelements();
874 0 : Int npol = 1;
875 0 : if (residual.ndim() > 2) npol = residualShape(2);
876 :
877 0 : IPosition cursorShape(ndim,1);
878 : {
879 0 : IPosition tileShape(residual.niceCursorShape());
880 0 : cursorShape(0) = tileShape(0);
881 0 : cursorShape(1) = tileShape(1);
882 0 : if (ndim > 2) cursorShape(2) = npol;
883 0 : }
884 0 : LatticeStepper nav(residualShape, cursorShape, LatticeStepper::RESIZE);
885 0 : RO_LatticeIterator<Float> residualIter(residual, nav);
886 0 : RO_LatticeIterator<Float> maskIter(residualIter); // There is no default ctr
887 0 : if (itsSoftMaskPtr != 0) {
888 0 : IPosition mCursorShape = cursorShape;
889 0 : mCursorShape(2) = 1;
890 0 : LatticeStepper mNav(itsSoftMaskPtr->shape(), mCursorShape, LatticeStepper::RESIZE);
891 0 : maskIter = RO_LatticeIterator<Float>(*itsSoftMaskPtr, mNav);
892 0 : maskIter.reset();
893 0 : }
894 :
895 0 : hist = 0;
896 0 : const Int nbins = hist.nelements();
897 : Bool histIsACopy;
898 0 : Int * histPtr = hist.getStorage(histIsACopy);
899 :
900 0 : Bool residualIsACopy = false;
901 0 : const Float * residualPtr = 0;
902 0 : Bool maskIsACopy = false;
903 0 : const Float * maskPtr = 0;
904 0 : for (residualIter.reset(); !residualIter.atEnd(); residualIter++) {
905 0 : residualPtr = residualIter.cursor().getStorage(residualIsACopy);
906 0 : if (itsSoftMaskPtr != 0) {
907 0 : maskPtr = maskIter.cursor().getStorage(maskIsACopy);
908 : }
909 0 : const IPosition thisCursorShape = residualIter.cursorShape();
910 0 : const Int npix = thisCursorShape(0) * thisCursorShape(1);
911 0 : switch (npol) {
912 0 : case 1:
913 0 : if (itsSoftMaskPtr == 0) {
914 0 : abshisf(histPtr, &minVal, &maxVal, &nbins, residualPtr, &npix);
915 : } else {
916 0 : abshimf(histPtr, &minVal, &maxVal, &nbins, residualPtr, maskPtr,
917 : &npix);
918 : }
919 0 : break;
920 0 : case 2:
921 0 : if (itsSoftMaskPtr == 0) {
922 0 : abshis2f(histPtr, &minVal, &maxVal, &nbins, residualPtr, &npix);
923 : } else {
924 0 : abshim2f(histPtr, &minVal, &maxVal, &nbins, residualPtr, maskPtr,
925 : &npix);
926 : }
927 0 : break;
928 0 : case 4:
929 0 : if (itsSoftMaskPtr == 0) {
930 0 : abshis4f(histPtr, &minVal, &maxVal, &nbins, residualPtr, &npix);
931 : } else {
932 0 : abshim4f(histPtr, &minVal, &maxVal, &nbins, residualPtr, maskPtr,
933 : &npix);
934 : }
935 0 : break;
936 : }
937 0 : if (itsSoftMaskPtr != 0) {
938 0 : maskIter.cursor().freeStorage(maskPtr, maskIsACopy);
939 0 : maskIter++;
940 : }
941 0 : residualIter.cursor().freeStorage(residualPtr, residualIsACopy);
942 0 : }
943 0 : hist.putStorage(histPtr, histIsACopy);
944 0 : }
945 :
946 : //----------------------------------------------------------------------
947 : // Determine the flux limit if we only select the maxNumPix biggest
948 : // residuals. Flux limit is not exact due to quantising by the histogram
949 0 : Float ClarkCleanLatModel::biggestResiduals(Float & maxRes,
950 : const uInt maxNumPix,
951 : const Float fluxLimit,
952 : const Lattice<Float> & residual) {
953 : // itsLog << LogOrigin("ClarkCleanLatModel", "biggestResiduals");
954 : // Calculate the histogram of the absolute value of the residuals
955 0 : Vector<Int> resHist(itsHistBins);
956 : // itsLog << "Created a histogram with " << resHist.nelements()
957 : // << " bins" << endl;;
958 : Float minRes;
959 0 : absHistogram(resHist, minRes, maxRes, residual);
960 : // itsLog << "Min/Max residuals are: " << minRes << " -> " << maxRes<< endl;
961 :
962 : // Deteremine how far we need to scan the histogram, before we reach the
963 : // flux cutoff imposed by the maximum exterior psf.
964 0 : Int lowbin=0;
965 0 : if (fluxLimit <= minRes)
966 0 : lowbin = 0;
967 0 : else if (fluxLimit < maxRes)
968 0 : lowbin=Int(itsHistBins*(fluxLimit-minRes)/(maxRes-minRes));
969 :
970 : // itsLog << "Select at most " << maxNumPix
971 : // << " pixels with the lowest bin being " << lowbin << endl;
972 :
973 0 : Int numPix = 0;
974 0 : Int curBin = itsHistBins - 1;
975 0 : while (curBin >= lowbin && numPix <= Int(maxNumPix)){
976 0 : numPix += resHist(curBin);
977 0 : curBin--;
978 : }
979 0 : curBin++;
980 :
981 : // Try to ensure we have maxNumPix or fewer residuals selected UNLESS
982 : // the topmost bin contains more than maxNumPix pixels. Then use all the
983 : // pixels in the topmost bin.
984 0 : if (numPix > Int(maxNumPix) && curBin != Int(itsHistBins) - 1){
985 0 : numPix -= resHist(curBin);
986 0 : curBin++;
987 : }
988 : // itsLog << "Selected " << numPix << " pixels from the top "
989 : // << itsHistBins - curBin << " bins" << LogIO::POST;
990 :
991 0 : return minRes+curBin*(maxRes-minRes)/Float(itsHistBins);
992 0 : }
993 :
994 : //----------------------------------------------------------------------
995 : // Work out the size of the Psf patch to use.
996 0 : Float ClarkCleanLatModel::getPsfPatch(Matrix<Float> & psfPatch,
997 : LatConvEquation & eqn) {
998 :
999 : // Determine the maximum possible size that should be used. Sizes greater
1000 : // than the maximum size cannot affect the cleaning and will not be used,
1001 : // even if the user requests it!
1002 0 : IPosition psfShape(eqn.psfSize());
1003 0 : uInt ndim = psfShape.nelements();
1004 0 : AlwaysAssert(ndim >= 2, AipsError);
1005 0 : IPosition modelShape(itsModelPtr->shape());
1006 0 : AlwaysAssert(modelShape.nelements() >= 2, AipsError);
1007 0 : IPosition maxShape(2, 1);
1008 0 : maxShape(0) = std::min(2*modelShape(0), psfShape(0));
1009 0 : maxShape(1) = std::min(2*modelShape(1), psfShape(1));
1010 :
1011 :
1012 : // See if the user has set a psfPatch size, and if it is less than the
1013 : // maximum size use it.
1014 0 : if (itsPsfPatchSize.nelements() != 0) {
1015 0 : IPosition psfPatchSize(2, 1, 1);
1016 0 : psfPatchSize(0) = std::min(maxShape(0), itsPsfPatchSize(0));
1017 0 : psfPatchSize(1) = std::min(maxShape(1), itsPsfPatchSize(1));
1018 0 : itsPsfPatchSize = psfPatchSize;
1019 0 : } else {
1020 0 : itsPsfPatchSize = maxShape;
1021 : }
1022 :
1023 : // Now calculate the maximum exterior psf value
1024 : // This is calculated where possible, otherwise a user supplied value is
1025 : // used.
1026 :
1027 : // Check if Psf is big enough to do a proper calculation
1028 : // This is the full psf, not the patch
1029 0 : Lattice<Float> *psf_p = 0;
1030 0 : Float maxExtPsf(0);
1031 0 : IPosition beyond(psfShape.nelements(), 0);
1032 0 : beyond(0) = itsPsfPatchSize(0)/2;
1033 0 : beyond(1) = itsPsfPatchSize(1)/2;
1034 0 : if (max((2*modelShape-psfShape).asVector()) <= 0){
1035 0 : if ( (itsPsfPatchSize(0) == (2*modelShape(0))) &&
1036 0 : (itsPsfPatchSize(1) == (2*modelShape(1))) ) {
1037 0 : maxExtPsf = Float(0); // Here the PsfPatch is used is big enough so
1038 : // that exterior sidelobes are irrelevant
1039 : }
1040 : else { // Calculate the exterior sidelobes
1041 0 : psf_p = eqn.evaluate(psfShape/2, Float(1), psfShape);
1042 0 : maxExtPsf = absMaxBeyondDist(beyond, psfShape/2, *psf_p);
1043 : }
1044 : }
1045 : else { // psf is not big enough so try and estimate something sensible
1046 :
1047 0 : if ((itsPsfPatchSize(0) == psfShape(0)) &&
1048 0 : (itsPsfPatchSize(1) == psfShape(1)) )
1049 : {
1050 0 : maxExtPsf = itsMaxExtPsf; // must use the user supplied value as it is
1051 : // impossible to estimate anything
1052 : } else { // try and estimate the ext. Psf and use the maximum of this
1053 : // value and the user supplied value
1054 0 : psf_p = eqn.evaluate(psfShape/2, Float(1), psfShape);
1055 0 : Float tempMax = absMaxBeyondDist(beyond, psfShape/2, *psf_p);
1056 0 : maxExtPsf = std::max(tempMax, itsMaxExtPsf);
1057 :
1058 : }
1059 : }
1060 0 : if(psf_p) delete psf_p; psf_p=0;
1061 : // set the max external psf Value to what is actually used.
1062 : // So the user can find out.
1063 0 : itsMaxExtPsf = maxExtPsf;
1064 :
1065 : // Now get a psf of the required size
1066 :
1067 :
1068 : {
1069 0 : Array<Float> a_psfPatch;
1070 0 : IPosition a_PsfPatchSize(psfShape.nelements(), 1);
1071 0 : a_PsfPatchSize(0) = itsPsfPatchSize(0);
1072 0 : a_PsfPatchSize(1) = itsPsfPatchSize(1);
1073 0 : AlwaysAssert(a_PsfPatchSize.product() == itsPsfPatchSize.product(), AipsError);
1074 0 : eqn.evaluate(a_psfPatch, beyond, Float(1), a_PsfPatchSize);
1075 0 : psfPatch = a_psfPatch.reform(itsPsfPatchSize);
1076 0 : }
1077 0 : return maxExtPsf;
1078 0 : }
1079 : //----------------------------------------------------------------------
1080 : // The maximum residual is the absolute maximum.
1081 0 : Float ClarkCleanLatModel::maxResidual(const Lattice<Float> & residual) {
1082 0 : const IPosition residualShape = residual.shape();
1083 0 : const uInt ndim = residualShape.nelements();
1084 0 : Int npol = 1;
1085 0 : if (ndim > 2) npol = residualShape(2);
1086 :
1087 0 : IPosition cursorShape(ndim, 1);
1088 : {
1089 0 : IPosition tileShape(residual.niceCursorShape());
1090 0 : cursorShape(0) = tileShape(0);
1091 0 : cursorShape(1) = tileShape(1);
1092 0 : if (ndim > 2) cursorShape(2) = npol;
1093 0 : }
1094 0 : LatticeStepper nav(residualShape, cursorShape, LatticeStepper::RESIZE);
1095 0 : RO_LatticeIterator<Float> residualIter(residual, nav);
1096 0 : RO_LatticeIterator<Float> maskIter(residualIter); // There is no default ctr
1097 0 : if (itsSoftMaskPtr != 0) {
1098 0 : IPosition mCursorShape = cursorShape;
1099 0 : mCursorShape(2) = 1;
1100 0 : LatticeStepper mNav(itsSoftMaskPtr->shape(), mCursorShape, LatticeStepper::RESIZE);
1101 0 : maskIter = RO_LatticeIterator<Float>(*itsSoftMaskPtr, mNav);
1102 0 : maskIter.reset();
1103 0 : }
1104 :
1105 0 : Float maxVal = 0, iterMaxVal = 0;
1106 0 : Bool residualIsACopy = false;
1107 0 : const Float * residualPtr = 0;
1108 0 : Bool maskIsACopy = false;
1109 0 : const Float * maskPtr = 0;
1110 0 : for (residualIter.reset(); !residualIter.atEnd(); residualIter++) {
1111 0 : residualPtr = residualIter.cursor().getStorage(residualIsACopy);
1112 0 : if (itsSoftMaskPtr != 0) {
1113 0 : maskPtr = maskIter.cursor().getStorage(maskIsACopy);
1114 : }
1115 0 : const IPosition thisCursorShape = residualIter.cursorShape();
1116 0 : const Int npix = thisCursorShape(0) * thisCursorShape(1);
1117 0 : switch (npol) {
1118 0 : case 1:
1119 0 : if (itsSoftMaskPtr == 0) {
1120 0 : maxabsf(&iterMaxVal, residualPtr, &npix);
1121 : } else {
1122 0 : maxabmf(&iterMaxVal, residualPtr, maskPtr, &npix);
1123 : }
1124 0 : break;
1125 0 : case 2:
1126 0 : if (itsSoftMaskPtr == 0) {
1127 0 : maxabs2f(&iterMaxVal, residualPtr, &npix);
1128 : } else {
1129 0 : maxabm2f(&iterMaxVal, residualPtr, maskPtr, &npix);
1130 : }
1131 0 : break;
1132 0 : case 4:
1133 0 : if (itsSoftMaskPtr == 0) {
1134 0 : maxabs4f(&iterMaxVal, residualPtr, &npix);
1135 : } else {
1136 0 : maxabm4f(&iterMaxVal, residualPtr, maskPtr, &npix);
1137 : }
1138 : }
1139 0 : if (iterMaxVal > maxVal) maxVal = iterMaxVal;
1140 0 : if (itsSoftMaskPtr == 0) {
1141 0 : maskIter.cursor().freeStorage(maskPtr, maskIsACopy);
1142 0 : maskIter++;
1143 : }
1144 0 : residualIter.cursor().freeStorage(residualPtr, residualIsACopy);
1145 0 : }
1146 0 : return maxVal;
1147 0 : }
1148 : //----------------------------------------------------------------------
1149 0 : void ClarkCleanLatModel::
1150 : maxVect(Block<Float> & maxVal, Float & absVal, Int & offset,
1151 : const CCList & activePixels) {
1152 0 : const Int npol = activePixels.nPol();
1153 0 : const Int numComp = activePixels.nComp();
1154 0 : DebugAssert(numComp > 0 , AipsError);
1155 0 : DebugAssert((maxVal.nelements() == 1 || maxVal.nelements() == 2
1156 : || maxVal.nelements() == 4) , AipsError);
1157 :
1158 0 : const Float * dataPtr = activePixels.fluxPtr();
1159 0 : Float * maxPtr = maxVal.storage();
1160 0 : switch (npol) {
1161 0 : case 1:
1162 0 : absmaxf(maxPtr, &absVal, &offset, dataPtr, &numComp);
1163 0 : break;
1164 0 : case 2:
1165 0 : absmax2f(maxPtr, &absVal, &offset, dataPtr, &numComp);
1166 0 : break;
1167 0 : case 4:
1168 0 : absmax4f(maxPtr, &absVal, &offset, dataPtr, &numComp);
1169 : }
1170 0 : }
1171 :
1172 : //----------------------------------------------------------------------
1173 0 : void ClarkCleanLatModel::
1174 : subtractComponent(CCList & activePixels, const Block<Float> & maxVal,
1175 : const Block<Int> & maxPos, const Matrix<Float> & psfPatch) {
1176 0 : const Int nx = psfPatch.nrow();
1177 0 : const Int ny = psfPatch.ncolumn();
1178 :
1179 0 : const Int numComp = activePixels.nComp();
1180 0 : Float * pixFluxPtr = activePixels.fluxPtr();
1181 0 : const Int * pixPosPtr = activePixels.positionPtr();
1182 0 : const Float * maxValPtr = maxVal.storage();
1183 0 : const Int * maxPosPtr = maxPos.storage();
1184 : Bool psfIsACopy;
1185 0 : const Float * psfPtr = psfPatch.getStorage(psfIsACopy);
1186 0 : switch (activePixels.nPol()){
1187 0 : case 1:
1188 0 : subcomf(pixFluxPtr, pixPosPtr, &numComp, maxValPtr, maxPosPtr,
1189 : psfPtr, &nx, &ny);
1190 0 : break;
1191 0 : case 2:
1192 0 : subcom2f(pixFluxPtr, pixPosPtr, &numComp, maxValPtr, maxPosPtr,
1193 : psfPtr, &nx, &ny);
1194 0 : break;
1195 0 : case 4:
1196 0 : subcom4f(pixFluxPtr, pixPosPtr, &numComp, maxValPtr, maxPosPtr,
1197 : psfPtr, &nx, &ny);
1198 0 : break;
1199 : }
1200 0 : psfPatch.freeStorage(psfPtr, psfIsACopy);
1201 0 : }
1202 :
1203 : //----------------------------------------------------------------------
1204 : // For an Array make a vector which gives the peak beyond distance n, p(n):
1205 : // p(0)= central value, p(n)=max value outside hypercube with side 2n-1
1206 : // Distance is measured from the point centre in the array
1207 0 : Float ClarkCleanLatModel::absMaxBeyondDist(const IPosition & maxDist,
1208 : const IPosition & centre,
1209 : const Lattice<Float> & psf) {
1210 0 : const IPosition psfShape = psf.shape();
1211 0 : const uInt ndim = psfShape.nelements();
1212 0 : const Int nx = psfShape(0);
1213 0 : const Int ny = psfShape(1);
1214 0 : AlwaysAssert(psfShape.product() == Int(nx*ny), AipsError);
1215 0 : IPosition cursorShape(ndim, 1);
1216 0 : cursorShape(0) = nx;
1217 0 : LatticeStepper nav(psfShape, cursorShape);
1218 0 : RO_LatticeIterator<Float> iter(psf, nav);
1219 0 : Int miny = centre(1) - maxDist(1);
1220 0 : Int maxy = centre(1) + maxDist(1);
1221 0 : Int nleft = centre(0) - maxDist(0);
1222 0 : Int nright = nx - (centre(0) + maxDist(0));
1223 0 : Int rightoffset = centre(0) + maxDist(0);
1224 0 : Float psfMax = 0, rowMax = 0;
1225 0 : Bool psfIsACopy = false;
1226 0 : const Float * psfPtr = 0;
1227 0 : const Float * endPtr = 0;
1228 0 : for (iter.reset(); !iter.atEnd(); iter++) {
1229 0 : psfPtr = iter.cursor().getStorage(psfIsACopy);
1230 0 : const IPosition cursorPos = iter.position();
1231 0 : if (cursorPos(1) < miny || cursorPos(1) > maxy) { // do the whole row
1232 0 : maxabsf(&rowMax, psfPtr, &nx);
1233 : } else { // just the ends of the row
1234 0 : maxabsf(&rowMax, psfPtr, &nleft);
1235 0 : if (rowMax > psfMax) psfMax = rowMax;
1236 0 : endPtr = psfPtr + rightoffset;
1237 0 : maxabsf(&rowMax, endPtr, &nright);
1238 0 : endPtr = 0;
1239 : }
1240 0 : iter.cursor().freeStorage(psfPtr, psfIsACopy);
1241 0 : if (rowMax > psfMax) psfMax = rowMax;
1242 0 : }
1243 0 : return psfMax;
1244 0 : }
1245 :
1246 0 : Bool ClarkCleanLatModel::stopnow() {
1247 0 : if (itsChoose == true) {
1248 0 : Vector<String> choices(2);
1249 0 : choices(0) = "Continue";
1250 0 : choices(1) = "Stop Now";
1251 0 : choices(2) = "Don't ask again";
1252 : String choice = Choice::choice("Do you want to continue or stop?",
1253 0 : choices);
1254 0 : if (choice == choices(1)) {
1255 0 : itsLog << "Clark clean stopped at user request" << LogIO::POST;
1256 0 : return true;
1257 : }
1258 0 : if (choice == choices(2)) {
1259 0 : setChoose(false);
1260 0 : itsLog << "Continuing: won't ask again" << LogIO::POST;
1261 : }
1262 0 : }
1263 0 : return false;
1264 : }
1265 :
1266 0 : Int ClarkCleanLatModel::
1267 : getbig(Float const * pixValPtr, Int const * pixPosPtr, const Int nPix,
1268 : const Float fluxLimit,
1269 : const Float * const residualPtr, const Float * const maskPtr,
1270 : const uInt npol, const Int nx, const Int ny)
1271 : {
1272 0 : Int nUsedPix = nPix;
1273 0 : switch (npol) {
1274 0 : case 1:
1275 0 : if (maskPtr == 0) {
1276 0 : getbigf(pixValPtr, pixPosPtr, &nUsedPix, &fluxLimit, residualPtr,
1277 : &nx, &ny);
1278 : } else {
1279 0 : getbimf(pixValPtr, pixPosPtr, &nUsedPix, &fluxLimit,
1280 : residualPtr, maskPtr, &nx, &ny);
1281 : }
1282 0 : break;
1283 0 : case 2:
1284 0 : if (maskPtr == 0) {
1285 0 : getbig2f(pixValPtr, pixPosPtr, &nUsedPix, &fluxLimit, residualPtr,
1286 : &nx, &ny);
1287 : } else {
1288 0 : getbim2f(pixValPtr, pixPosPtr, &nUsedPix, &fluxLimit,
1289 : residualPtr, maskPtr, &nx, &ny);
1290 : }
1291 0 : break;
1292 0 : case 4:
1293 0 : if (maskPtr == 0) {
1294 0 : getbig4f(pixValPtr, pixPosPtr, &nUsedPix, &fluxLimit, residualPtr,
1295 : &nx, &ny);
1296 : } else {
1297 0 : getbim4f(pixValPtr, pixPosPtr, &nUsedPix, &fluxLimit,
1298 : residualPtr, maskPtr, &nx, &ny);
1299 : }
1300 0 : break;
1301 : }
1302 0 : return nUsedPix;
1303 : }
1304 :
1305 0 : void ClarkCleanLatModel::
1306 : updateModel(CCList & cleanComponents) {
1307 0 : const IPosition modelShape(itsModelPtr->shape());
1308 0 : const uInt ndim = modelShape.nelements();
1309 0 : IPosition cursorShape(ndim,1);
1310 0 : if (ndim > 2) cursorShape(2) = modelShape(2);
1311 0 : IPosition cs = itsModelPtr->niceCursorShape();
1312 0 : cleanComponents.tiledSort(cs);
1313 : const Int * compPositionPtr;
1314 0 : IPosition compPosition(ndim,0);
1315 0 : Array<Float> compFlux(cursorShape), modelFlux(cursorShape);
1316 0 : for(uInt c = 0; c < cleanComponents.nComp(); c++) {
1317 0 : compPositionPtr = cleanComponents.pixelPosition(c);
1318 0 : compPosition(0) = *compPositionPtr;
1319 0 : compPosition(1) = *(compPositionPtr+1);
1320 0 : compFlux.takeStorage(cursorShape, cleanComponents.pixelFlux(c), SHARE);
1321 0 : itsModelPtr->getSlice(modelFlux, compPosition, cursorShape);
1322 : // itsLog << "Model " << modelFlux << " @ " << compPosition;
1323 0 : modelFlux += compFlux;
1324 : // itsLog << " -> " << modelFlux;
1325 : // itsLog << " Subtracting:" << compFlux << " @ " << compPosition;
1326 0 : itsModelPtr->putSlice(modelFlux, compPosition);
1327 : }
1328 0 : }
1329 : // Local Variables:
1330 : // compile-command: "gmake OPTLIB=1 ClarkCleanLatModel"
1331 : // End:
1332 :
1333 : } //# NAMESPACE CASA - END
1334 :
|
