Line data Source code
1 : //# ClarkCleanModel.cc: this defines ClarkCleanModel
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,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/ClarkCleanModel.h>
29 : #include <synthesis/MeasurementEquations/ClarkCleanProgress.h>
30 : #include <casacore/casa/Arrays/IPosition.h>
31 : #include <casacore/casa/Arrays/Slice.h>
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/Vector.h>
34 : #include <casacore/casa/Arrays/Array.h>
35 : #include <casacore/casa/Arrays/ArrayError.h>
36 : #include <casacore/casa/Arrays/ArrayMath.h>
37 : #include <casacore/casa/Arrays/ArrayLogical.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 <casacore/casa/iostream.h>
43 : #include <casacore/casa/System/Choice.h>
44 :
45 : using namespace casacore;
46 : namespace casa { //# NAMESPACE CASA - BEGIN
47 :
48 : //# These are the definitions of the fortran functions
49 :
50 : #define NEED_FORTRAN_UNDERSCORES
51 :
52 : #if defined(NEED_FORTRAN_UNDERSCORES)
53 : #define getbigf getbigf_
54 : #define getbig2f getbig2f_
55 : #define getbig4f getbig4f_
56 : #define abshisf abshisf_
57 : #define abshis2f abshis2f_
58 : #define abshis4f abshis4f_
59 : #define absmaxf absmaxf_
60 : #define absmax2f absmax2f_
61 : #define absmax4f absmax4f_
62 : #define subcomf subcomf_
63 : #define subcom2f subcom2f_
64 : #define subcom4f subcom4f_
65 : #define maxabsf maxabsf_
66 : #define maxabs2f maxabs2f_
67 : #define maxabs4f maxabs4f_
68 : #define maxabmf maxabmf_
69 : #define maxabm2f maxabm2f_
70 : #define maxabm4f maxabm4f_
71 : #define abshimf abshimf_
72 : #define abshim2f abshim2f_
73 : #define abshim4f abshim4f_
74 : #define getbimf getbimf_
75 : #define getbim2f getbim2f_
76 : #define getbim4f getbim4f_
77 : #endif
78 :
79 : extern "C" {
80 : void getbigf(Float * pixVal, Int * pixPos, Int * maxPix,
81 : const Float * fluxLimit, const Float * arr,
82 : const Int * nx, const Int * ny);
83 : void getbig2f(Float * pixVal, Int * pixPos, Int * maxPix,
84 : const Float * fluxLimit, const Float * arr,
85 : const Int * nx, const Int * ny);
86 : void getbig4f(Float * pixVal, Int * pixPos, Int * maxPix,
87 : const Float * fluxLimit, const Float * arr,
88 : const Int * nx, const Int * ny);
89 : void abshisf(Int * hist, Float * minval, Float * maxval,
90 : const Int * nbins, const Float * arr, const Int * npix);
91 : void abshis2f(Int * hist, Float * minval, Float * maxval,
92 : const Int * nbins, const Float * arr, const Int * npix);
93 : void abshis4f(Int * hist, Float * minval, Float * maxval,
94 : const Int * nbins, const Float * arr, const Int * npix);
95 : void absmaxf(Float * maxelem, Float * maxval, Int * maxpos,
96 : const Float * arr, const Int * npix);
97 : void absmax2f(Float * maxelem, Float * maxval, Int * maxpos,
98 : const Float * arr, const Int * npix);
99 : void absmax4f(Float * maxelem, Float * maxval, Int * maxpos,
100 : const Float * arr, const Int * npix);
101 : void subcomf(Float * pixval, const Int * pixpos, const Int * npix,
102 : const Float * maxpix, const Int * maxpos,
103 : const Float * psf, const Int * nx, const Int * ny);
104 : void subcom2f(Float * pixval, const Int * pixpos, const Int * npix,
105 : const Float * maxpix, const Int * maxpos,
106 : const Float * psf, const Int * nx, const Int * ny);
107 : void subcom4f(Float * pixval, const Int * pixpos, const Int * npix,
108 : const Float * maxpix, const Int * maxpos,
109 : const Float * psf, const Int * nx, const Int * ny);
110 : void maxabsf(Float * maxval, const Float * arr, const Int * npix);
111 : void maxabs2f(Float * maxval, const Float * arr, const Int * npix);
112 : void maxabs4f(Float * maxval, const Float * arr, const Int * npix);
113 : void maxabmf(Float * maxval, const Float * arr, const Float * mask,
114 : const Int * npix);
115 : void maxabm2f(Float * maxval, const Float * arr, const Float * mask,
116 : const Int * npix);
117 : void maxabm4f(Float * maxval, const Float * arr, const Float * mask,
118 : const Int * npix);
119 : void abshimf(Int * hist, Float * minval, Float * maxval, const Int * nbins,
120 : const Float * arr, const Float * mask, const Int * npix);
121 : void abshim2f(Int * hist, Float * minval, Float * maxval, const Int * nbins,
122 : const Float * arr, const Float * mask, const Int * npix);
123 : void abshim4f(Int * hist, Float * minval, Float * maxval, const Int * nbins,
124 : const Float * arr, const Float * mask, const Int * npix);
125 : void getbimf(Float * pixVal, Int * pixPos, Int * maxPix,
126 : const Float * fluxLimit, const Float * arr, const Float * mask,
127 : const Int * nx, const Int * ny);
128 : void getbim2f(Float * pixVal, Int * pixPos, Int * maxPix,
129 : const Float * fluxLimit, const Float * arr, const Float * mask,
130 : const Int * nx, const Int * ny);
131 : void getbim4f(Float * pixVal, Int * pixPos, Int * maxPix,
132 : const Float * fluxLimit, const Float * arr, const Float * mask,
133 : const Int * nx, const Int * ny);
134 : };
135 :
136 :
137 : //----------------------------------------------------------------------
138 0 : ClarkCleanModel::ClarkCleanModel()
139 : :ArrayModel<Float>(),
140 0 : theHistBins(1024),
141 0 : theMaxExtPsf(0.0),
142 0 : theMaxNumberMinorIterations(10000),
143 0 : theInitialNumberIterations(0),
144 0 : theMaxNumberMajorCycles(-1),
145 0 : theMaxNumPix(32*1024),
146 0 : thePsfPatchSize(2,51,51),
147 0 : theSpeedup(0.0),
148 0 : theCycleSpeedup(-1.0),
149 0 : theChoose(false),
150 0 : theMask(),
151 0 : theLog(LogOrigin("ClarkCleanModel", "ClarkCleanModel()")),
152 0 : theIterCounter(0),
153 0 : itsProgressPtr(0),
154 0 : itsJustStarting(true)
155 : {
156 0 : };
157 : //----------------------------------------------------------------------
158 0 : ClarkCleanModel::ClarkCleanModel(Array<Float> & model)
159 : :ArrayModel<Float>(model),
160 0 : theHistBins(1024),
161 0 : theMaxExtPsf(0.0),
162 0 : theMaxNumberMinorIterations(10000),
163 0 : theInitialNumberIterations(0),
164 0 : theMaxNumberMajorCycles(-1),
165 0 : theMaxNumPix(32*1024),
166 0 : thePsfPatchSize(2, 51, 51),
167 0 : theSpeedup(0.0),
168 0 : theCycleSpeedup(-1.0),
169 0 : theChoose(false),
170 0 : theMask(),
171 0 : theLog(LogOrigin("ClarkCleanModel",
172 : "ClarkCleanModel(const Array<Float> & model)")),
173 0 : theIterCounter(0),
174 0 : itsProgressPtr(0),
175 0 : itsJustStarting(true)
176 : {
177 0 : AlwaysAssert(theModel.ndim() >= 2, AipsError);
178 0 : if (theModel.ndim() >= 3)
179 0 : AlwaysAssert(theModel.shape()(2) == 1 || theModel.shape()(2) == 2 ||
180 : theModel.shape()(2) == 4, AipsError);
181 0 : if (theModel.ndim() >= 4)
182 0 : for (uInt i = 3; i < theModel.ndim(); i++)
183 0 : AlwaysAssert(theModel.shape()(i) == 1, AipsError);
184 : // theLog << LogOrigin("ClarkCleanModel", "ClarkCleanModel")
185 : // << "Model shape is:" << theModel.shape() << endl;
186 0 : };
187 : //----------------------------------------------------------------------
188 0 : ClarkCleanModel::ClarkCleanModel(Array<Float> & model,
189 0 : Array<Float> & mask)
190 : :ArrayModel<Float>(model),
191 0 : theHistBins(1024),
192 0 : theMaxExtPsf(0.0),
193 0 : theMaxNumberMinorIterations(10000),
194 0 : theInitialNumberIterations(0),
195 0 : theMaxNumberMajorCycles(-1),
196 0 : theMaxNumPix(32*1024),
197 0 : thePsfPatchSize(2, 51, 51),
198 0 : theSpeedup(0.0),
199 0 : theCycleSpeedup(-1.0),
200 0 : theChoose(false),
201 0 : theMask(mask),
202 0 : theLog(LogOrigin("ClarkCleanModel",
203 : "ClarkCleanModel(Array<Float> & model"
204 : ", Array<Float> & mask)")),
205 0 : theIterCounter(0),
206 0 : itsProgressPtr(0),
207 0 : itsJustStarting(true)
208 : {
209 0 : AlwaysAssert(theModel.ndim() >= 2, AipsError);
210 0 : if (theModel.ndim() >= 3)
211 0 : AlwaysAssert(theModel.shape()(2) == 1 || theModel.shape()(2) == 2 ||
212 : theModel.shape()(2) == 4, AipsError);
213 0 : if (theModel.ndim() >= 4)
214 0 : for (uInt i = 3; i < theModel.ndim(); i++)
215 0 : AlwaysAssert(theModel.shape()(i) == 1, AipsError);
216 :
217 0 : AlwaysAssert(theMask.ndim() >= 2, AipsError);
218 0 : AlwaysAssert(theMask.shape()(0) == theModel.shape()(0), AipsError);
219 0 : AlwaysAssert(theMask.shape()(1) == theModel.shape()(1), AipsError);
220 0 : if (theMask.ndim() >= 3)
221 0 : for (uInt i = 2; i < theMask.ndim(); i++)
222 0 : AlwaysAssert(theMask.shape()(i) == 1, AipsError);
223 0 : };
224 :
225 0 : void ClarkCleanModel::getModel(Array<Float>& model) const{
226 0 : model = theModel;
227 0 : };
228 0 : void ClarkCleanModel::setModel(const Array<Float>& model){
229 0 : AlwaysAssert(model.ndim() >= 2, AipsError);
230 0 : if (model.ndim() >= 3)
231 0 : AlwaysAssert(model.shape()(2) == 1 || model.shape()(2) == 2 ||
232 : model.shape()(2) == 4, AipsError);
233 0 : if (model.ndim() >= 4)
234 0 : for (uInt i = 3; i < model.ndim(); i++)
235 0 : AlwaysAssert(model.shape()(i) == 1, AipsError);
236 0 : theModel = model;
237 0 : };
238 0 : void ClarkCleanModel::setModel(Array<Float> & model){
239 0 : AlwaysAssert(model.ndim() >= 2, AipsError);
240 0 : if (model.ndim() >= 3)
241 0 : AlwaysAssert(model.shape()(2) == 1 || model.shape()(2) == 2 ||
242 : model.shape()(2) == 4, AipsError);
243 0 : if (model.ndim() >= 4)
244 0 : for (uInt i = 3; i < model.ndim(); i++)
245 0 : AlwaysAssert(model.shape()(i) == 1, AipsError);
246 0 : theModel.reference(model);
247 0 : };
248 :
249 0 : void ClarkCleanModel::getMask(Array<Float>& mask) const{
250 0 : mask = theMask;
251 0 : };
252 0 : void ClarkCleanModel::setMask(const Array<Float>& mask){
253 0 : AlwaysAssert(mask.ndim() >= 2, AipsError);
254 0 : AlwaysAssert(mask.shape()(0) == theModel.shape()(0), AipsError);
255 0 : AlwaysAssert(mask.shape()(1) == theModel.shape()(1), AipsError);
256 0 : if (mask.ndim() >= 3)
257 0 : for (uInt i = 2; i < mask.ndim(); i++)
258 0 : AlwaysAssert(mask.shape()(i) == 1, AipsError);
259 0 : theMask = mask;
260 0 : };
261 0 : void ClarkCleanModel::setMask(Array<Float> & mask){
262 0 : AlwaysAssert(mask.ndim() >= 2, AipsError);
263 0 : AlwaysAssert(mask.shape()(0) == theModel.shape()(0), AipsError);
264 0 : AlwaysAssert(mask.shape()(1) == theModel.shape()(1), AipsError);
265 0 : if (mask.ndim() >= 3)
266 0 : for (uInt i = 2; i < mask.ndim(); i++)
267 0 : AlwaysAssert(mask.shape()(i) == 1, AipsError);
268 0 : theMask.reference(mask);
269 0 : };
270 :
271 :
272 : //----------------------------------------------------------------------
273 0 : Bool ClarkCleanModel::solve(ConvolutionEquation & eqn){
274 0 : theLog << LogOrigin("ClarkCleanModel", "solve");
275 0 : AlwaysAssert(theModel.ndim() >= 2, AipsError);
276 0 : const IPosition dataShape = theModel.shape();
277 0 : Int npol = 1;
278 0 : if (theModel.ndim() >= 3)
279 0 : npol = dataShape(2);
280 0 : AlwaysAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
281 :
282 : // Determine the number of polarisations
283 : // theLog << "Data has " << npol << " polarisations" << LogIO::POST;
284 :
285 : // compute the current residual image (using an FFT)
286 0 : Array<Float> residual;
287 0 : eqn.residual(residual, *this);
288 :
289 : // Determine the psf patch to use
290 0 : Matrix<Float> psfPatch;
291 : Float maxExtPsf;
292 0 : maxExtPsf = getPsfPatch(psfPatch, eqn);
293 : // theLog << "PsfPatch shape is: " << psfPatch.shape()
294 : // << " and has a maximum exterior sidelobe of "
295 : // << maxExtPsf << LogIO::POST;
296 :
297 : // Declare variables needed inside the following while loop
298 : Float minLimit; // the min flux limit when using the
299 : // maximum number of active pixels
300 : Int numPix; // the number of Active pixels
301 0 : Int maxNumPix = 0; // The max. number of active pixels ever used
302 0 : uInt numIterations = theInitialNumberIterations;
303 : // Number of Iterations done so far
304 : uInt numMinorIterations; // The number of minor iterations done/todo
305 0 : uInt numMajorCycles = 0; // The number of major cycles done
306 0 : uInt maxNumberMinorIterations = 0;// The max. number of min. iterations
307 : // ever used
308 0 : Matrix<Float> pixelValue; // cache of "active" pixel values
309 0 : Matrix<Int> pixelPos; // cache of "active" pixel positions
310 0 : Float Fmn=1; // The "uncertainty" factor
311 : Float fluxLimit; // The fluxlimit for the current major cycle
312 0 : Float totalFlux = 0;
313 :
314 : // Note that a Matrix is used rather than say a Vector of IPositions as
315 : // it allows the inner loop (in doMinorIterations()) to be more highly
316 : // optimised (using pointers)
317 :
318 : // find its maximum value of the residual
319 0 : Float maxRes = maxResidual(residual);
320 :
321 0 : theLog << "Initial maximum residual: " << maxRes << LogIO::POST;
322 : // if flux limit or iteration limit reached then bail out.
323 0 : Bool userHalt = false;
324 0 : while ((Int(numIterations) < numberIterations()) &&
325 0 : (maxRes > threshold()) &&
326 0 : ((theMaxNumberMajorCycles<0)||(numMajorCycles<(uInt)theMaxNumberMajorCycles)) &&
327 0 : userHalt == false){
328 :
329 : // determine the fluxlimit for this major cycle
330 : // choose fluxlimit for residuals to be considered in minor iterations
331 : // don't consider residuals below maxRes*(value of maxPsf at outside the
332 : // patch)
333 0 : fluxLimit = maxRes * maxExtPsf;
334 : // theLog << "Fluxlimit determined using the Maximum exterior Psf: "
335 : // << fluxLimit << LogIO::POST;
336 : // See if this flux limit needs to be modified because it selects too
337 : // many pixels.
338 0 : minLimit = biggestResiduals(maxRes, theMaxNumPix, fluxLimit, residual);
339 : // theLog << "Fluxlimit determined using the maximum number active pixels: "
340 : // << minLimit << endl;
341 0 : fluxLimit = max(fluxLimit, minLimit);
342 0 : fluxLimit /=3.0; //This factor was found empirically as
343 : // as it it was too conservative in the minor loop
344 : // theLog << "Final Fluxlimit: " << fluxLimit << LogIO::POST;
345 :
346 : // Copy all the active pixels into separate areas for better memory
347 : // management and quicker indexing.
348 0 : numPix = cacheActivePixels(pixelValue, pixelPos, residual,
349 0 : max(fluxLimit,threshold()));
350 : // The numpix calculated here frequently differs
351 : // from the number calculated using the histogram, because of the
352 : // quantisation of the flux levels in the histogram, and the imposition
353 : // of an external fluxlevel.
354 0 : if (numPix > 0) {
355 : // theLog <<"Major cycle has "<< numPix << " active residuals, "
356 : // << "a Fluxlimit of " << max(fluxLimit,threshold()) << endl;
357 : // Start of minor cycles
358 0 : numMinorIterations = min(theMaxNumberMinorIterations,
359 0 : numberIterations()-numIterations);
360 0 : doMinorIterations(theModel, pixelValue, pixelPos, numPix,
361 : psfPatch, fluxLimit, numMinorIterations,
362 : Fmn, numIterations, totalFlux);
363 0 : numIterations += numMinorIterations;
364 : // theLog << "Clean has used " << numIterations << " Iterations" ;
365 0 : maxNumberMinorIterations = max(maxNumberMinorIterations,
366 : numMinorIterations);
367 0 : maxNumPix = max(maxNumPix, numPix);
368 : // Now do a major cycle
369 0 : eqn.residual(residual, *this);
370 :
371 : // find the new maximum residual
372 0 : maxRes = maxResidual(residual);
373 0 : theLog << "Iteration: " << numIterations
374 0 : << ", Maximum residual=" << maxRes << LogIO::POST;
375 : // << " Flux limit=" << max(fluxLimit,threshold())
376 : // << ", " << numPix << " Active pixels" << LogIO::POST;
377 :
378 : // theLog << " to get to a maximum residual of " << maxRes << LogIO::POST;
379 :
380 : // Count the number of major cycles
381 0 : numMajorCycles++;
382 : }
383 : else{
384 0 : theLog << LogIO::WARN
385 : << "Zero Pixels selected with a Fluxlimit of " << fluxLimit
386 0 : << " and a maximum Residual of " << maxRes << LogIO::POST;
387 0 : userHalt=true;
388 : }
389 0 : userHalt = userHalt || stopnow();
390 : }
391 0 : setThreshold(maxRes);
392 0 : setNumberIterations(numIterations);
393 0 : theMaxNumPix = maxNumPix;
394 0 : theMaxNumberMinorIterations = maxNumberMinorIterations;
395 0 : return true;
396 0 : };
397 :
398 : //----------------------------------------------------------------------
399 0 : Bool ClarkCleanModel::singleSolve(ConvolutionEquation & eqn,
400 : Array<Float>& residual){
401 0 : theLog << LogOrigin("ClarkCleanModel", "singleSolve");
402 0 : AlwaysAssert(theModel.ndim() >= 2, AipsError);
403 0 : const IPosition dataShape = theModel.shape();
404 0 : Int npol = 1;
405 0 : if (theModel.ndim() >= 3)
406 0 : npol = dataShape(2);
407 0 : AlwaysAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
408 :
409 : // Determine the number of polarisations
410 : // theLog << "Data has " << npol << " polarisations" << LogIO::POST;
411 :
412 : // Determine the psf patch to use
413 0 : Matrix<Float> psfPatch;
414 : Float maxExtPsf;
415 0 : maxExtPsf = getPsfPatch(psfPatch, eqn);
416 : // theLog << "PsfPatch shape is: " << psfPatch.shape()
417 : // << " and has a maximum exterior sidelobe of "
418 : // << maxExtPsf << LogIO::POST;
419 :
420 : // Declare variables needed inside the following while loop
421 : Float minLimit; // the min flux limit when using the
422 : // maximum number of active pixels
423 : Int numPix; // the number of Active pixels
424 0 : Int maxNumPix = 0; // The max. number of active pixels ever used
425 0 : uInt numIterations = theInitialNumberIterations;
426 : // Number of Iterations done so far
427 : uInt numMinorIterations; // The number of minor iterations done/todo
428 0 : uInt maxNumberMinorIterations = 0;// The max. number of min. iterations
429 : // ever used
430 0 : Matrix<Float> pixelValue; // cache of "active" pixel values
431 0 : Matrix<Int> pixelPos; // cache of "active" pixel positions
432 0 : Float Fmn=1; // The "uncertainty" factor
433 : Float fluxLimit; // The fluxlimit for the current major cycle
434 0 : Float totalFlux = 0;
435 :
436 : // Note that a Matrix is used rather than say a Vector of IPositions as
437 : // it allows the inner loop (in doMinorIterations()) to be more highly
438 : // optimised (using pointers)
439 :
440 : // find its maximum value of the residual
441 0 : Float maxRes = maxResidual(residual);
442 :
443 : // move this to later TT
444 : //theLog << "Initial maximum residual: " << maxRes << LogIO::POST;
445 : //
446 : // if flux limit or iteration limit reached then bail out.
447 :
448 : // determine the fluxlimit for this major cycle
449 : // choose fluxlimit for residuals to be considered in minor iterations
450 : // don't consider residuals below maxRes*(value of maxPsf at outside the
451 : // patch)
452 0 : fluxLimit = maxRes * maxExtPsf;
453 : // theLog << "Fluxlimit determined using the Maximum exterior Psf: "
454 : // << fluxLimit << LogIO::POST;
455 : // See if this flux limit needs to be modified because it selects too
456 : // many pixels.
457 0 : minLimit = biggestResiduals(maxRes, theMaxNumPix, fluxLimit, residual);
458 : // theLog << "Fluxlimit determined using the maximum number active pixels: "
459 : // << minLimit << endl;
460 0 : fluxLimit = max(fluxLimit, minLimit);
461 0 : fluxLimit /= 8.0; //emepirically found that fluxlimit was too
462 : // too conservative
463 : // theLog << "Final Fluxlimit: " << fluxLimit << LogIO::POST;
464 :
465 : // Copy all the active pixels into separate areas for better memory
466 : // management and quicker indexing.
467 0 : numPix = cacheActivePixels(pixelValue, pixelPos, residual,
468 0 : max(fluxLimit,threshold()));
469 : // The numpix calculated here frequently differs
470 : // from the number calculated using the histogram, because of the
471 : // quantisation of the flux levels in the histogram, and the imposition
472 : // of an external fluxlevel.
473 0 : if (numPix > 0) {
474 : // theLog <<"Major cycle has "<< numPix << " active residuals, "
475 : // << "a Fluxlimit of " << max(fluxLimit,threshold()) << endl;
476 : // Start of minor cycles
477 0 : theLog << "Initial maximum residual: " << maxRes << LogIO::POST;
478 0 : numMinorIterations = min(theMaxNumberMinorIterations,
479 0 : numberIterations()-numIterations);
480 0 : doMinorIterations(theModel, pixelValue, pixelPos, numPix,
481 : psfPatch, fluxLimit, numMinorIterations,
482 : Fmn, numIterations, totalFlux);
483 0 : numIterations += numMinorIterations;
484 : // theLog << "Clean has used " << numIterations << " Iterations" ;
485 0 : maxNumberMinorIterations = max(maxNumberMinorIterations,
486 : numMinorIterations);
487 0 : maxNumPix = max(maxNumPix, numPix);
488 : }
489 : /*
490 : else
491 : theLog << LogIO::WARN
492 : << "Zero Pixels selected with a Fluxlimit of " << fluxLimit
493 : << " and a maximum Residual of " << maxRes << LogIO::POST;
494 : */
495 0 : setNumberIterations(numIterations);
496 0 : theMaxNumPix = maxNumPix;
497 0 : theMaxNumberMinorIterations = maxNumberMinorIterations;
498 :
499 0 : return true;
500 0 : };
501 : //----------------------------------------------------------------------
502 0 : void ClarkCleanModel::setPsfPatchSize(const IPosition & psfPatchSize){
503 0 : thePsfPatchSize=psfPatchSize;
504 0 : };
505 : //----------------------------------------------------------------------
506 0 : IPosition ClarkCleanModel::getPsfPatchSize(){
507 0 : return thePsfPatchSize;
508 : };
509 : //----------------------------------------------------------------------
510 0 : void ClarkCleanModel::setHistLength(const uInt HistBins ){
511 0 : theHistBins=HistBins;
512 0 : };
513 : //----------------------------------------------------------------------
514 0 : uInt ClarkCleanModel::getHistLength(){
515 0 : return theHistBins;
516 : };
517 : //----------------------------------------------------------------------
518 0 : void ClarkCleanModel::setMaxNumberMinorIterations(const uInt maxNumMinorIterations){
519 0 : theMaxNumberMinorIterations=maxNumMinorIterations;
520 0 : };
521 : //----------------------------------------------------------------------
522 0 : uInt ClarkCleanModel::getMaxNumberMinorIterations(){
523 0 : return theMaxNumberMinorIterations;
524 : };
525 : //----------------------------------------------------------------------
526 0 : void ClarkCleanModel::setInitialNumberIterations(const uInt initialNumberIterations){
527 0 : theInitialNumberIterations=initialNumberIterations;
528 0 : };
529 : //----------------------------------------------------------------------
530 0 : uInt ClarkCleanModel::getInitialNumberIterations(){
531 0 : return theInitialNumberIterations;
532 : };
533 : //----------------------------------------------------------------------
534 0 : void ClarkCleanModel::setMaxNumberMajorCycles(const uInt maxNumMajorCycles){
535 0 : theMaxNumberMajorCycles=maxNumMajorCycles;
536 0 : };
537 : //----------------------------------------------------------------------
538 0 : uInt ClarkCleanModel::getMaxNumberMajorCycles(){
539 0 : return theMaxNumberMajorCycles;
540 : };
541 : //----------------------------------------------------------------------
542 0 : void ClarkCleanModel::setMaxNumPix(const uInt maxNumPix ){
543 0 : theMaxNumPix=maxNumPix;
544 0 : };
545 : //----------------------------------------------------------------------
546 0 : uInt ClarkCleanModel::getMaxNumPix(){
547 0 : return theMaxNumPix;
548 : };
549 : //----------------------------------------------------------------------
550 0 : void ClarkCleanModel::setMaxExtPsf(const Float maxExtPsf ){
551 0 : theMaxExtPsf=maxExtPsf;
552 0 : };
553 : //----------------------------------------------------------------------
554 0 : Float ClarkCleanModel::getMaxExtPsf(){
555 0 : return theMaxExtPsf;
556 : };
557 : //----------------------------------------------------------------------
558 0 : void ClarkCleanModel::setSpeedup(const Float speedup ){
559 0 : theSpeedup=speedup;
560 0 : };
561 : //----------------------------------------------------------------------
562 0 : Float ClarkCleanModel::getSpeedup(){
563 0 : return theSpeedup;
564 : };
565 : //----------------------------------------------------------------------
566 0 : void ClarkCleanModel::setCycleSpeedup(const Float speedup ){
567 0 : theCycleSpeedup=speedup;
568 0 : };
569 : //----------------------------------------------------------------------
570 0 : Float ClarkCleanModel::getCycleSpeedup(){
571 0 : return theCycleSpeedup;
572 : };
573 : //----------------------------------------------------------------------
574 0 : void ClarkCleanModel::setChoose(const Bool choose ){
575 0 : theChoose=choose;
576 0 : };
577 : //----------------------------------------------------------------------
578 0 : Bool ClarkCleanModel::getChoose(){
579 0 : return theChoose;
580 : };
581 : //----------------------------------------------------------------------
582 0 : void ClarkCleanModel::doMinorIterations(Array<Float> & model,
583 : Matrix<Float> & pixVal,
584 : const Matrix<Int> & pixPos,
585 : const Int numPix,
586 : Matrix<Float> & psfPatch,
587 : Float fluxLimit,
588 : uInt & numberIterations,
589 : Float Fmn,
590 : const uInt totalIterations,
591 : Float &totalFlux){
592 0 : DebugAssert(model.ndim() >= 2, AipsError);
593 0 : DebugAssert(model.shape()(0) > 0, AipsError);
594 0 : DebugAssert(model.shape()(1) > 0, AipsError);
595 0 : Int npol = 1;
596 0 : if (model.ndim() >= 3)
597 0 : npol = model.shape()(2);
598 0 : DebugAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
599 0 : DebugAssert(model.shape()(3) == 1, AipsError);
600 0 : DebugAssert(Int(pixVal.nrow()) == npol, AipsError);
601 0 : DebugAssert(numPix <= Int(pixVal.ncolumn()), AipsError);
602 0 : DebugAssert(0 < numPix , AipsError);
603 0 : DebugAssert(pixPos.nrow() == 2, AipsError);
604 0 : DebugAssert(pixPos.ncolumn() == pixVal.ncolumn(), AipsError);
605 0 : DebugAssert(psfPatch.nrow() > 0, AipsError);
606 0 : DebugAssert(psfPatch.ncolumn() > 0, AipsError);
607 :
608 : // theLog << LogOrigin("ClarkCleanModel", "doMinorIterations");
609 : // Find the largest residual and its position.
610 0 : Vector<Float> maxRes(npol);
611 0 : Vector<Int> maxPos(2);
612 : Float absRes;
613 : Float signedAbsRes;
614 : Int offRes;
615 0 : maxVect(maxRes, absRes, offRes, pixVal, numPix);
616 0 : maxPos = pixPos.column(offRes);
617 : // declare variables used inside the main loop
618 0 : Int curIter = 0;
619 0 : Float iterFluxLimit = max(fluxLimit, threshold());
620 0 : Float Fac = pow(fluxLimit/absRes, theSpeedup);
621 0 : IPosition position(model.ndim(), 0);
622 : // theLog << "Initial maximum residual:" << maxRes
623 : // << " (" << absRes << ") "
624 : // << " @ " << maxPos << endl;
625 :
626 : // Do the minor Iterations
627 0 : while ((curIter < Int(numberIterations)) && (absRes > iterFluxLimit)){
628 0 : iterFluxLimit = max(fluxLimit, threshold()); // threshold() changes now!
629 0 : maxRes *= gain();
630 0 : totalFlux += maxRes(0);
631 : // Add the new component to the current model
632 0 : position(0) = maxPos(0);
633 0 : position(1) = maxPos(1);
634 0 : if (model.ndim() >= 3)
635 0 : for (Int p = 0; p < npol; p++){
636 0 : position(2) = p;
637 : // theLog << "Model " << model(position) << " @ " << position;
638 0 : model(position) += maxRes(p);
639 : // theLog << " -> " << model(position);
640 : }
641 : else {
642 : // theLog << "Model " << model(position) << " @ " << position;
643 0 : model(position) += maxRes(0);
644 : // theLog << " -> " << model(position);
645 : }
646 :
647 : // theLog << " Subtracting:" << maxRes
648 : // << " @ " << position;
649 : // Subtract the component from the current list of active pixels
650 0 : subtractComponent(pixVal, pixPos, numPix,
651 : maxRes, maxPos, psfPatch);
652 : // We have now done an iteration
653 0 : curIter++;
654 0 : theIterCounter++;
655 : // find the next residual
656 0 : maxVect(maxRes, absRes, offRes, pixVal, numPix);
657 0 : maxPos = pixPos.column(offRes);
658 : // theLog << " After Iteration: " << curIter
659 : // << " the Maximum residual is:" << maxRes
660 : // << " (" << absRes << ") "
661 : // << " @ " << maxPos << LogIO::POST;
662 : // Update the uncertainty factors and fluxlimits
663 0 : Fmn += Fac/Float(totalIterations+curIter);
664 0 : iterFluxLimit = max(fluxLimit * Fmn, threshold());
665 :
666 0 : if (itsProgressPtr) {
667 : try {
668 : // if this does not throw an exception, we are in business
669 0 : itsProgressPtr->hasPGPlotter();
670 0 : signedAbsRes = absRes * maxRes(0)/abs( maxRes(0) );
671 0 : itsProgressPtr->
672 0 : info(false, (Int)(totalIterations+curIter), (Int)numberIterations,
673 0 : signedAbsRes, IPosition(2,maxPos(0),maxPos(1)),
674 0 : totalFlux, false, itsJustStarting );
675 0 : itsJustStarting = false;
676 0 : } catch (AipsError x) {
677 : // if it throw an exception, do nothing
678 0 : }
679 : }
680 : }
681 : // Data returned to the main routine
682 0 : numberIterations = curIter;
683 0 : fluxLimit = absRes;
684 0 : itsMaxRes = maxRes(0);
685 0 : };
686 : //----------------------------------------------------------------------
687 0 : Int ClarkCleanModel::
688 : cacheActivePixels(Matrix<Float> & pixVal, Matrix<Int> & pixPos,
689 : const Array<Float> & data, const Float fluxLimit){
690 0 : DebugAssert(data.ndim() >= 2, AipsError);
691 0 : const IPosition dataShape = data.shape();
692 0 : const Int nx = dataShape(0);
693 0 : const Int ny = dataShape(1);
694 0 : Int npol = 1;
695 0 : if (data.ndim() >= 3)
696 0 : npol = dataShape(2);
697 0 : DebugAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
698 0 : DebugAssert(nx > 0, AipsError);
699 0 : DebugAssert(ny > 0, AipsError);
700 0 : DebugAssert(pixVal.ncolumn() == pixPos.ncolumn(), AipsError);
701 :
702 0 : Int nBigPix = pixVal.ncolumn();
703 : Bool dataCopy, valCopy, posCopy;
704 0 : const Float * dataPtr = data.getStorage(dataCopy);
705 0 : Float * valPtr = pixVal.getStorage(valCopy);
706 0 : Int * posPtr = pixPos.getStorage(posCopy);
707 :
708 0 : if (theMask.nelements() == 0) {
709 0 : switch (npol){
710 0 : case 1:
711 0 : getbigf(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, &nx, &ny);
712 0 : break;
713 0 : case 2:
714 0 : getbig2f(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, &nx, &ny);
715 0 : break;
716 0 : case 4:
717 0 : getbig4f(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, &nx, &ny);
718 : }
719 0 : if (nBigPix > 0){ // I could be more efficient about this
720 0 : nBigPix += pixVal.ncolumn();
721 0 : pixVal.putStorage(valPtr, valCopy); pixPos.putStorage(posPtr, posCopy);
722 0 : pixVal.resize(npol, nBigPix);
723 0 : pixPos.resize(2, nBigPix);
724 0 : valPtr = pixVal.getStorage(valCopy); posPtr = pixPos.getStorage(posCopy);
725 0 : switch (npol){
726 0 : case 1:
727 0 : getbigf(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, &nx, &ny);
728 0 : break;
729 0 : case 2:
730 0 : getbig2f(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, &nx, &ny);
731 0 : break;
732 0 : case 4:
733 0 : getbig4f(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, &nx, &ny);
734 : }
735 0 : AlwaysAssert(nBigPix == 0, AipsError);
736 : }
737 : }
738 : else {
739 : Bool maskCopy;
740 0 : const Float * maskPtr = theMask.getStorage(maskCopy);
741 0 : switch (npol){
742 0 : case 1:
743 0 : getbimf(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, maskPtr, &nx, &ny);
744 0 : break;
745 0 : case 2:
746 0 : getbim2f(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, maskPtr, &nx, &ny);
747 0 : break;
748 0 : case 4:
749 0 : getbim4f(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, maskPtr, &nx, &ny);
750 : }
751 0 : if (nBigPix > 0){ // I could be more effecient about this
752 0 : nBigPix += pixVal.ncolumn();
753 0 : pixVal.putStorage(valPtr, valCopy); pixPos.putStorage(posPtr, posCopy);
754 0 : pixVal.resize(npol, nBigPix);
755 0 : pixPos.resize(2, nBigPix);
756 0 : valPtr = pixVal.getStorage(valCopy); posPtr = pixPos.getStorage(posCopy);
757 0 : switch (npol){
758 0 : case 1:
759 0 : getbimf(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, maskPtr, &nx, &ny);
760 0 : break;
761 0 : case 2:
762 0 : getbim2f(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, maskPtr, &nx, &ny);
763 0 : break;
764 0 : case 4:
765 0 : getbim4f(valPtr, posPtr, &nBigPix, &fluxLimit, dataPtr, maskPtr, &nx, &ny);
766 : }
767 0 : AlwaysAssert(nBigPix == 0, AipsError);
768 : }
769 0 : theMask.freeStorage(maskPtr, dataCopy);
770 : }
771 0 : pixVal.putStorage(valPtr, valCopy);
772 0 : pixPos.putStorage(posPtr, posCopy);
773 0 : data.freeStorage(dataPtr, dataCopy);
774 0 : DebugAssert(nBigPix <= 0 && (nBigPix + Int(pixVal.ncolumn())) >= Int(0), AipsError);
775 0 : return pixVal.ncolumn() + nBigPix;
776 0 : };
777 : //----------------------------------------------------------------------
778 : // make histogram of absolute values in array
779 0 : void ClarkCleanModel::absHistogram(Vector<Int> & hist,
780 : Float & minVal,
781 : Float & maxVal,
782 : const Array<Float> & data) {
783 0 : DebugAssert(data.ndim() >= 2, AipsError);
784 0 : const IPosition dataShape = data.shape();
785 0 : const Int nx = dataShape(0);
786 0 : const Int ny = dataShape(1);
787 0 : Int npol = 1;
788 0 : if (data.ndim() >= 3)
789 0 : npol = dataShape(2);
790 0 : const Int npix = nx*ny;
791 0 : DebugAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
792 0 : DebugAssert(nx > 0, AipsError);
793 0 : DebugAssert(ny > 0, AipsError);
794 : Bool dataCopy, histCopy;
795 0 : const Int nbins = hist.nelements();
796 0 : const Float * dataPtr = data.getStorage(dataCopy);
797 0 : Int * histPtr = hist.getStorage(histCopy);
798 0 : hist = 0;
799 0 : if (theMask.nelements() == 0)
800 0 : switch (npol){
801 0 : case 1:
802 0 : abshisf(histPtr, &minVal, &maxVal, &nbins, dataPtr, &npix);
803 0 : break;
804 0 : case 2:
805 0 : abshis2f(histPtr, &minVal, &maxVal, &nbins, dataPtr, &npix);
806 0 : break;
807 0 : case 4:
808 0 : abshis4f(histPtr, &minVal, &maxVal, &nbins, dataPtr, &npix);
809 : }
810 : else{
811 : Bool maskCopy;
812 0 : const Float * maskPtr = theMask.getStorage(maskCopy);
813 0 : switch (npol){
814 0 : case 1:
815 0 : abshimf(histPtr, &minVal, &maxVal, &nbins, dataPtr, maskPtr, &npix);
816 0 : break;
817 0 : case 2:
818 0 : abshim2f(histPtr, &minVal, &maxVal, &nbins, dataPtr, maskPtr, &npix);
819 0 : break;
820 0 : case 4:
821 0 : abshim4f(histPtr, &minVal, &maxVal, &nbins, dataPtr, maskPtr, &npix);
822 : }
823 0 : theMask.freeStorage(maskPtr, dataCopy);
824 : }
825 0 : data.freeStorage(dataPtr, dataCopy);
826 0 : hist.putStorage(histPtr, histCopy);
827 :
828 0 : };
829 : //----------------------------------------------------------------------
830 : // Determine the flux limit if we only select the maxNumPix biggest
831 : // residuals. Flux limit is not exact due to quantising by the histogram
832 0 : Float ClarkCleanModel::biggestResiduals(Float & maxRes,
833 : const uInt maxNumPix,
834 : const Float fluxLimit,
835 : const Array<Float> & residual) {
836 : // theLog << LogOrigin("ClarkCleanModel", "biggestResiduals");
837 : // Calculate the histogram of the absolute value of the residuals
838 0 : Vector<Int> resHist(theHistBins);
839 : // theLog << "Created a histogram with " << resHist.nelements()
840 : // << " bins" << endl;;
841 : Float minRes;
842 0 : absHistogram(resHist, minRes, maxRes, residual);
843 : // theLog << "Min/Max residuals are: " << minRes << " -> " << maxRes<< endl;
844 :
845 : // Deteremine how far we need to scan the histogram, before we reach the
846 : // flux cutoff imposed by the maximum exterior psf.
847 : Int lowbin;
848 0 : if (fluxLimit <= minRes)
849 0 : lowbin = 0;
850 0 : else if (fluxLimit >= maxRes)
851 0 : lowbin = theHistBins - 1;
852 : else
853 0 : lowbin=Int(theHistBins*(fluxLimit-minRes)/(maxRes-minRes));
854 :
855 : // theLog << "Select at most " << maxNumPix
856 : // << " pixels with the lowest bin being " << lowbin << endl;
857 :
858 0 : Int numPix = 0;
859 0 : Int curBin = theHistBins - 1;
860 0 : while (curBin >= lowbin && numPix <= Int(maxNumPix)){
861 0 : numPix += resHist(curBin);
862 0 : curBin--;
863 : }
864 0 : curBin++;
865 :
866 : // Try to ensure we have maxNumPix or fewer residuals selected UNLESS
867 : // the topmost bin contains more than maxNumPix pixels. Then use all the
868 : // pixels in the topmost bin.
869 0 : if (numPix > Int(maxNumPix) && curBin != Int(theHistBins - 1)){
870 0 : numPix -= resHist(curBin);
871 0 : curBin++;
872 : }
873 : // theLog << "Selected " << numPix << " pixels from the top "
874 : // << theHistBins - curBin << " bins" << LogIO::POST;
875 :
876 0 : return minRes+curBin*(maxRes-minRes)/Float(theHistBins);
877 0 : }
878 : //----------------------------------------------------------------------
879 : // Work out the size of the Psf patch to use.
880 0 : Float ClarkCleanModel::getPsfPatch(Array<Float>& psfPatch,
881 : ConvolutionEquation& eqn) {
882 :
883 : // Determine the maximum possible size that should be used. Sizes greater
884 : // than the maximum size cannot affect the cleaning and will not be used,
885 : // even if the user requests it!
886 0 : IPosition psfSize(eqn.psfSize());
887 0 : uInt ndim = psfSize.nelements();
888 0 : IPosition modelSize(theModel.shape().getFirst(ndim));
889 0 : IPosition maxSize(min(2*modelSize.asVector(),
890 0 : psfSize.asVector()));
891 : // See if the user has set a psfPatch size, and if it is less than the
892 : // maximum size use it.
893 0 : IPosition psfPatchSize;
894 0 : if (thePsfPatchSize.nelements() != 0) {
895 0 : psfPatchSize = casacore::min(maxSize.asVector(),
896 0 : thePsfPatchSize.asVector());
897 : }
898 : else {
899 0 : psfPatchSize = maxSize;
900 : }
901 : // set the psf Patch size to what is actually used. So the user can find out.
902 0 : thePsfPatchSize = psfPatchSize;
903 :
904 : // Now calculate the maximum exterior psf value
905 :
906 : // This is calculated where possible, otherwise a user supplied value is
907 : // used.
908 :
909 : // Check if Psf is big enough to do a proper calculation
910 0 : Array<Float> psf;
911 0 : Float maxExtPsf(0);
912 0 : if (max((2*modelSize-psfSize).asVector()) <= 0){
913 0 : if (psfPatchSize.isEqual(2*modelSize)) {
914 0 : maxExtPsf = Float(0); // Here the PsfPatch is used is big enough so
915 : // that exterior sidelobes are irrelevant
916 : }
917 : else { // Calculate the exterior sidelobes
918 0 : eqn.evaluate(psf, psfSize/2, Float(1), psfSize);
919 0 : maxExtPsf = absMaxBeyondDist(psfPatchSize/2, psfSize/2, psf);
920 : }
921 : }
922 : else { // psf is not big enough so try and estimate something sensible
923 0 : if (psfPatchSize.isEqual(psfSize)) {
924 0 : maxExtPsf = theMaxExtPsf; // must use the user supplied value as it is
925 : // impossible to estimate anything
926 : }
927 : else { // try and estimate the ext. Psf and use the maximum of this
928 : // value and the user supplied value
929 0 : eqn.evaluate(psf, psfSize/2, Float(1), psfSize);
930 0 : maxExtPsf = max(absMaxBeyondDist(psfPatchSize/2, psfSize/2, psf),
931 : theMaxExtPsf);
932 : }
933 : }
934 0 : eqn.flushPsf(); // Tell the convolution equation to release the cached psf
935 : // set the max external psf Value to what is actually used.
936 : // So the user can find out.
937 0 : theMaxExtPsf = maxExtPsf;
938 : // Now get a psf of the required size
939 0 : eqn.evaluate(psfPatch, psfPatchSize/2, Float(1), psfPatchSize);
940 0 : return maxExtPsf;
941 0 : };
942 : //----------------------------------------------------------------------
943 : // The maximum residual is the absolute maximum.
944 0 : Float ClarkCleanModel::maxResidual(const Array<Float> & residual) {
945 0 : DebugAssert(residual.ndim() >= 2, AipsError);
946 0 : const IPosition dataShape = residual.shape();
947 0 : const Int nx = dataShape(0);
948 0 : const Int ny = dataShape(1);
949 0 : Int npol = 1;
950 0 : if (residual.ndim() >= 3)
951 0 : npol = dataShape(2);
952 0 : const Int npix = nx*ny;
953 0 : DebugAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
954 0 : DebugAssert(nx > 0, AipsError);
955 0 : DebugAssert(ny > 0, AipsError);
956 :
957 : Float maxVal;
958 : Bool dataCopy;
959 0 : const Float * dataPtr = residual.getStorage(dataCopy);
960 0 : if (theMask.nelements() == 0)
961 0 : switch (npol){
962 0 : case 1:
963 0 : maxabsf(&maxVal, dataPtr, &npix);
964 0 : break;
965 0 : case 2:
966 0 : maxabs2f(&maxVal, dataPtr, &npix);
967 0 : break;
968 0 : case 4:
969 0 : maxabs4f(&maxVal, dataPtr, &npix);
970 : }
971 : else {
972 : Bool maskCopy;
973 0 : const Float * maskPtr = theMask.getStorage(maskCopy);
974 0 : switch (npol){
975 0 : case 1:
976 0 : maxabmf(&maxVal, dataPtr, maskPtr, &npix);
977 0 : break;
978 0 : case 2:
979 0 : maxabm2f(&maxVal, dataPtr, maskPtr, &npix);
980 0 : break;
981 0 : case 4:
982 0 : maxabm4f(&maxVal, dataPtr, maskPtr, &npix);
983 : }
984 0 : theMask.freeStorage(maskPtr, dataCopy);
985 : }
986 0 : residual.freeStorage(dataPtr, dataCopy);
987 0 : return maxVal;
988 0 : };
989 : //----------------------------------------------------------------------
990 0 : void ClarkCleanModel::maxVect(Vector<Float> & maxVal,
991 : Float & absVal,
992 : Int & offset,
993 : const Matrix<Float> & pixVal,
994 : const Int numPix) {
995 0 : const Int npol = pixVal.nrow();
996 0 : DebugAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
997 0 : DebugAssert(numPix <= Int(pixVal.ncolumn()), AipsError);
998 0 : DebugAssert(0 < numPix , AipsError);
999 :
1000 : Bool dataCopy, maxCopy;
1001 0 : const Float * dataPtr = pixVal.getStorage(dataCopy);
1002 0 : Float * maxPtr = maxVal.getStorage(maxCopy);
1003 0 : switch (npol){
1004 0 : case 1:
1005 0 : absmaxf(maxPtr, &absVal, &offset, dataPtr, &numPix);
1006 0 : break;
1007 0 : case 2:
1008 0 : absmax2f(maxPtr, &absVal, &offset, dataPtr, &numPix);
1009 0 : break;
1010 0 : case 4:
1011 0 : absmax4f(maxPtr, &absVal, &offset, dataPtr, &numPix);
1012 : }
1013 0 : pixVal.freeStorage(dataPtr, dataCopy);
1014 0 : maxVal.putStorage(maxPtr, dataCopy);
1015 0 : };
1016 : //----------------------------------------------------------------------
1017 0 : void ClarkCleanModel::subtractComponent(Matrix<Float> & pixVal,
1018 : const Matrix<Int> & pixPos,
1019 : const Int numPix,
1020 : const Vector<Float> & maxVal,
1021 : const Vector<Int> & maxPos,
1022 : const Matrix<Float> & psf){
1023 0 : const Int npol = pixVal.nrow();
1024 0 : DebugAssert(npol == 1 || npol == 2 || npol == 4, AipsError);
1025 0 : DebugAssert(numPix <= Int(pixVal.ncolumn()), AipsError);
1026 0 : DebugAssert(0 < numPix , AipsError);
1027 0 : DebugAssert(pixPos.nrow() == 2, AipsError);
1028 0 : DebugAssert(pixPos.ncolumn() == pixVal.ncolumn(), AipsError);
1029 0 : DebugAssert(Int(maxVal.nelements()) == npol, AipsError);
1030 0 : DebugAssert(maxPos.nelements() == 2, AipsError);
1031 0 : const Int nx = psf.nrow();
1032 0 : const Int ny = psf.ncolumn();
1033 0 : DebugAssert(nx > 0, AipsError);
1034 0 : DebugAssert(ny > 0, AipsError);
1035 :
1036 : Bool pixValCopy, pixPosCopy, maxValCopy, maxPosCopy, psfCopy;
1037 0 : Float * pixValPtr = pixVal.getStorage(pixValCopy);
1038 0 : const Int * pixPosPtr = pixPos.getStorage(pixPosCopy);
1039 0 : const Float * maxValPtr = maxVal.getStorage(maxValCopy);
1040 0 : const Int * maxPosPtr = maxPos.getStorage(maxPosCopy);
1041 0 : const Float * psfPtr = psf.getStorage(psfCopy);
1042 0 : switch (npol){
1043 0 : case 1:
1044 0 : subcomf(pixValPtr, pixPosPtr, &numPix, maxValPtr, maxPosPtr,
1045 : psfPtr, &nx, &ny);
1046 0 : break;
1047 0 : case 2:
1048 0 : subcom2f(pixValPtr, pixPosPtr, &numPix, maxValPtr, maxPosPtr,
1049 : psfPtr, &nx, &ny);
1050 0 : break;
1051 0 : case 4:
1052 0 : subcom4f(pixValPtr, pixPosPtr, &numPix, maxValPtr, maxPosPtr,
1053 : psfPtr, &nx, &ny);
1054 : }
1055 0 : psf.freeStorage(psfPtr, psfCopy);
1056 0 : maxPos.freeStorage(maxPosPtr, maxPosCopy);
1057 0 : maxVal.freeStorage(maxValPtr, maxValCopy);
1058 0 : pixPos.freeStorage(pixPosPtr, pixPosCopy);
1059 0 : pixVal.putStorage(pixValPtr, pixValCopy);
1060 0 : };
1061 : //----------------------------------------------------------------------
1062 : // For an Array make a vector which gives the peak beyond distance n, p(n):
1063 : // p(0)= central value, p(n)=max value outside hypercube with side 2n-1
1064 : // Distance is measured from the point centre in the array
1065 0 : Float ClarkCleanModel::absMaxBeyondDist(const IPosition &maxDist,
1066 : const IPosition ¢re,
1067 : const Array<Float> &array){
1068 0 : if (maxDist.nelements() != array.ndim()) {
1069 0 : throw(ArrayError("Vector<Float> absMaxBeyondDist("
1070 : "const IPosition &maxDist,const IPosition ¢re,"
1071 : "const Array<Float> &array) - "
1072 0 : "maxDist dimensionality inconsistent with array"));
1073 : }
1074 0 : if (array.nelements() == 0) {
1075 0 : throw(ArrayError("Vector<Float> absMaxBeyondDist("
1076 : "const IPosition &maxDist,const IPosition ¢re,"
1077 : "const Array<Float> &array) - "
1078 0 : "Array has no elements"));
1079 : }
1080 : {
1081 0 : Vector<Int> tmp1 = (centre-maxDist).asVector();
1082 0 : Vector<Int> tmp2 = (centre+maxDist-array.endPosition()).asVector();
1083 0 : if (min(tmp1)<0 || max(tmp2)>0) {
1084 0 : throw(ArrayError("Vector<Float> absMaxBeyondDist("
1085 : "const IPosition &maxDist,const IPosition ¢re,"
1086 : "const Array<Float> &array) - "
1087 0 : "maxDist too large for Array"));
1088 : }
1089 0 : }
1090 : // Initialize
1091 0 : ReadOnlyVectorIterator<Float> ai(array);
1092 0 : Float maxVal(0);
1093 0 : uInt n = ai.vector().nelements();
1094 0 : uInt start = centre(0) - maxDist(0);
1095 0 : uInt end = centre(0) + maxDist(0) + 1;
1096 0 : IPosition vecPos(array.ndim());
1097 0 : IPosition vecDist(array.ndim());
1098 : uInt i;
1099 :
1100 : // loop though array accumulating maxima for each distance
1101 0 : while (! ai.pastEnd()) {
1102 : // find the distance of the current vector to the midpoint of the array
1103 0 : vecPos=ai.pos(); vecPos(0)=centre(0);
1104 0 : vecDist=abs((vecPos-centre).asVector());
1105 : // skip if current vector too far from midpoint in any dimension
1106 0 : if (max((vecDist-maxDist).asVector()) > 0)
1107 0 : for (i=0; i<n; i++) {
1108 0 : maxVal=max(maxVal, Float(abs(ai.vector()(i))));
1109 : }
1110 : else {
1111 0 : for (i=0; i<start; i++) {
1112 0 : maxVal=max(maxVal, Float(abs(ai.vector()(i))));
1113 : }
1114 0 : for (i=end; i<n; i++) {
1115 0 : maxVal=max(maxVal, Float(abs(ai.vector()(i))));
1116 : }
1117 : }
1118 0 : ai.next();
1119 : }
1120 0 : return maxVal;
1121 0 : };
1122 :
1123 0 : Bool ClarkCleanModel::stopnow() {
1124 0 : if(theChoose) {
1125 0 : Vector<String> choices(2);
1126 0 : choices(0)="Continue";
1127 0 : choices(1)="Stop Now";
1128 0 : choices(2)="Don't ask again";
1129 : String choice = Choice::choice("Do you want to continue or stop?",
1130 0 : choices);
1131 0 : if (choice==choices(0)) {
1132 0 : return false;
1133 : }
1134 0 : else if (choice==choices(2)) {
1135 0 : setChoose(false);
1136 0 : theLog << "Continuing: won't ask again" << LogIO::POST;
1137 0 : return false;
1138 : }
1139 : else {
1140 0 : theLog << "Clark clean stopped at user request" << LogIO::POST;
1141 0 : return true;
1142 : }
1143 0 : }
1144 : else {
1145 0 : return false;
1146 : }
1147 : }
1148 :
1149 0 : Float ClarkCleanModel::threshold()
1150 : {
1151 0 : Float thresh = Iterate::threshold();
1152 0 : if (theCycleSpeedup > 0.0) {
1153 0 : thresh = thresh * pow(2.0, ((Double)(theIterCounter)/theCycleSpeedup) );
1154 : }
1155 0 : return thresh;
1156 : };
1157 :
1158 : } //# NAMESPACE CASA - END
1159 :
|
